
function demo_montecarlo_2D

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

for k=1:trials
    AE=[];
    for n=NV
        P=rand(n,2);
        fP=feval(f,P(:,1),P(:,2));
        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);

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

AE_mat_aver(1:end-1)./AE_mat_aver(2:end)


% semilogy(NV,AE_mat_aver);
% hold on;
% semilogy(NV,1./sqrt(NV),'r');
% hold off
% 


fprintf('\n \n');