function [value] = n_d_Monte_Carlo_integral (n, m, n_cases, f)
% Accompanying program for the text
%
% Classical and Modern Numerical Analysis:
% Theory, Methods and Practice
% by Azmy S. Ackleh, Edward J. Allen,
% R. Baker Kearfott, and Padmanabhan Seshaiyer
%
% (Taylor and Francis / CRC Press, 2009)
%
% [value] = n_d_Monte_Carlo_integral (n, m, n_cases, f) returns a guess
% for the integral of a function f of n variables over the unit n-cube
% [0,1]^n, using m sample points. The mean of these tries is returned in
% value and the estimated standard deviation is returned in deviation.
% This routine also prints the estimated value and estimated
% for each try standard deviation.
%
% The string "f" is the name of an m-file containing the function to be
% integrated. For example, if sum_of_squares.m contains
%
% function [y] = sum_of_squares(x)
% y = x'*x;
%
% then issuing
%
% value = n_d_Monte_Carlo_integral (10, 5000, 5, 'sum_of_squares')
%
% will return an approximation to 10/3 in value by guessing at a
% 10-dimensional integral, using 5000 sample points in each of 5 tries,
% resulting in a standard deviation of deviation.
%
% The state of the random nuber generator is reset at each call of this
% function, so the results are reproducible.
rand('state',1);
value = 0;
disp(' try no. m est. mean est. dev. ')
for icase=1:n_cases
s1=0.0;
s2=0.0;
for i=1:m
x=rand(n,1);
fval = feval(f,x);
s1=s1+fval/m;
s2=s2+fval^2/m;
end
estmean=s1;
value = value + estmean;
eststdev=sqrt(s2-s1*s1)/sqrt(m);
disp((sprintf(' %5.0f %10.0f %10.6f %10.6f', icase, m, estmean, eststdev)))
end
value = value / n_cases;