% Monte Carlo computation of the D-dimensional integral:
% I(f) = int_0^1 sin(pi*(x1+x2+...+xD)/4)^2 dx1 dx2... dxD
% using the Sample-Mean method.
% K. Ming Leung, 02/18/03

clear all; hold off;
repeat = 25;     % Repeat calculation for a fixed N 25 times & averaged
D = 5;              % Dimension 5
exactI = (1-(2*sqrt(2)/pi)^D*cos(D*pi/4))/2;     % Exact value of I

nPts = 100*10.^(0:3);         % Number of points used in simulation
for k = 1:4
% Get sample points and sample function values
s = sin((pi/4)*sum(rand(repeat,nPts(k),D),3));
s = s.*s;
% Sample mean method
Mean = sum(s,2)/nPts(k);
Integral = Mean;
ActualError = abs(Integral-exactI);
AvgError(k) = sum(ActualError)/length(ActualError);
end

t = -log(nPts);
y = log(AvgError);
plot(t,y,'*'); hold on;
xlabel('-log(N)');
ylabel('log(Error)');
title('Sample-Mean Monte Carlo Method for a 5-Dimensional Integral');

% Linear fit of the data
A = [ones(4,1) t'];
AT = A';
AA = AT*A;
b = AT*y';
x = AA\b;
disp('The slope is: '); disp(x(2));
plot(t,x(1)+x(2)*t,'r-');