clear all; close all;

% Initializations
s = 1000;
true = 2^s * cos(2*pi+s/2)*sin(1/2)^s;
fun = @(x) cos(2*pi+sum(x));
maxiter = 20;
sums = [];
means = [];

% Compute the Monte Carlo error
for ii = 0:maxiter
    disp(['iteration: ',num2str(ii)]);
    nincr = max(1,2^(ii-1)); % reuse the previous 2^(ii-1) points
    ntot = 2^ii; % total number of points
    sample = rand(s,nincr);
    sums = [sums;sum(fun(sample))];
    means = [means;sum(sums)/ntot];
end
errors = [2.^(0:maxiter)',abs(means-true)];

% Least squares fit
lsq = [ones(length(errors),1),log(errors(:,1))]\log(errors(:,2));
lsq(1) = exp(lsq(1));

% Visualize the results
figure;
loglog(errors(:,1),errors(:,2),'r.','MarkerSize',18), hold on
x = errors(:,1);
loglog(x,lsq(1)*x.^(lsq(2)),'b','LineWidth',2), hold off;
set(gca,'FontSize',14,'FontWeight','bold')
legend('data',['slope: ',num2str(lsq(2))])
title(['Monte Carlo (s = ',num2str(s),')'])
xlabel('n')
ylabel('absolute error')