
function demo_QMC_1D

NV=2.^(1:18);
f=@(x) exp(-x.^2); I=sqrt(pi)*erf(1)/2;
trials=1000;
AE_mat=[];

QMC_type=5; 

for k=1:trials
    AE=[];
    for n=NV
        P=QMC_pointset(QMC_type,n);
        fP=feval(f,P);
        In=sum(fP)/n;
        AE(end+1,1)=abs(I-In);
    end
    AE_mat(:,k)=AE;
end

AE_mat_aver=abs(sum(AE_mat,2)/trials);

fprintf('\n \n \t \t * AVERAGE ABSOLUTE ERRORS \n \n');
for k=1:length(NV)
    fprintf('\n \t n: %7.0f ae: %1.5e',NV(k),AE_mat_aver(k));
end

fprintf('\n \n \t \t \t * AVERAGE RATIOS \n \n');
rat=AE_mat_aver(1:end-1)./AE_mat_aver(2:end);
for k=1:length(NV)-1
    fprintf('\n \t n1: %7.0f n2: %7.0f rat: %1.5e',NV(k),NV(k+1),rat(k));
end
fprintf('\n \n');






function P=QMC_pointset(QMC_type,n)

switch QMC_type

    case 0
        P=rand(n,1);

    case 1

        Pset = sobolset(1);
        P=Pset(1:n,:);

    case 2
        Pset = haltonset(1);
        P=Pset(1:n,:);

    case 3
        Pset = sobolset(1);
        Pset = scramble(Pset,'RR2');
        P=Pset(1:n,:);

    case 4
        Pset = haltonset(1);
        Pset = scramble(Pset,'RR2');
        P=Pset(1:n,:);

    case 5
        P = Hammersley(n,1); P=P';

end

