Name

intprb_montecarlo — Estimates a multidimensional integral using iterated Monte Carlo.

Calling Sequence

   integr = intprb_montecarlo ( func , n )
   integr = intprb_montecarlo ( func , n , evalf )
   integr = intprb_montecarlo ( func , n , evalf , bounds )
   integr = intprb_montecarlo ( func , n , evalf , bounds , tol )
   integr = intprb_montecarlo ( func , n , evalf , bounds , tol , randgen )
   integr = intprb_montecarlo ( func , n , evalf , bounds , tol , randgen , verbose )
   [ integr , accur , varf , evalf ] = intprb_montecarlo ( ... )
   
   

Parameters

func :

a function or a list, the name of the function to be evaluated.

n:

a 1 x 1 matrix of floating point integers, the spatial dimension.

bounds :

a n x 2 matrix of doubles, where bounds(:,1) = lower bound and bounds(:,2) = upper bound of integration (default = [zeros(n,1) ones(n,1)])

evalf :

a 2 x 1 matrix of floating point integers, evalf(1) = the maximum number of function evaluations and evalf(2) = the minimum number of function evaluations (default maxevalf = 10^4, minevalf = 10^2)

tol :

a 2 x 1 matrix of floating point integers, tol(1) = relative tolerance on the accuracy, tol(2) = absolute tolerance on the accuracy (default tolrel = 10^-2, tolabs = 10^-2)

randgen :

a function, the random number generator. (default = grand)

verbose :

a 1 x 1 matrix of boolean, set to %t to enable verbose messages. (default = %f)

integr :

a 1 x 1 matrix of doubles, the approximate value of the integral, the mean of the function.

accur :

a 1 x 1 matrix of doubles, the estimated error on the integral.

varf :

a 1 x 1 matrix of doubles, the approximate value of the variance.

evalf :

a 1 x 1 matrix of floating point integers, the actual number of calls to the function.

iter :

a 1 x 1 matrix of floating point integers, the number of loops in the algorithm

Description

The algorithm uses a crude Monte-Carlo method to estimate the integral. It tries to minimize the number of calls to the function, while maintaining good performance using vectorization. The performance is achieved by calling the function and giving it never less than minevalf experiments to perform by call. Minimization of the number of calls is achieved by trying to get the expected accuracy with the lowest possible amount of calls to the function, never exceeding maxevalf total function calls.

If an optional input argument is equal to the empty matrix, it is set to its default value.

The function should have header y = func( m , n , x ) where m is a floating point integer representing the number of experiments to perform and x is a m x n matrix of doubles. The uniform random number generator must have header y = randgen ( m , n ) where m is a 1 x 1 matrix of floating point integers, the number of experiments to be performed and y is a m x n matrix of doubles in the interval [0,1].

On output, the number of calls to the function, evalf, is between minevalf and maxevalf. In order to get good performance, the number of calls to the function is generally not lower than minevalf.

The error is estimated as accur = 2 * stdy / sqrt(evalf), where stdy is the standard deviation associated with the outputs of the function. This error ensures that the probability that the actual integral is in the interval [integr-accur,integr+accur] is greater than 95%. This is because the Monte-Carlo estimate follows a Normal (Laplace-Gauss) law, with variance variance(f)/sqrt(evalf). The probability 95% comes as the approximate probability that a normal variable is in the range [mu-2*sigma,mu+2*sigma].

If the randgen input argument is the empty matrix, the grand function is used to produce the uniform random numbers.

The algorithm stops if the condition accur < tolrel * integr + tolabs is met. The default value of tolrel and tolabs is 10^-2. This corresponds to the order of magnitude of the expected accuracy for a function with mean zero and variance one, if the number of simulations is equal to 10^4 and varf = 2 (in this case, we have 2 * 1 / sqrt(10^2) = 2.10^-2).

Examples

// A basic example
function r = sumfunction (m,n,x)
// sum of n variables
// Expectation = 0, Variance = 1
v = n/12
e = n/2
r = sqrt(1/v) * (sum(x,"c") - e)
endfunction
// Test with default settings
grand("setsd",123456);
dim_num = 10;
[ integr , accur , varf , evalf , iter ] = intprb_montecarlo ( sumfunction , dim_num )

// Test with customized maximum number of evaluations
grand("setsd",123456);
dim_num = 10;
maxevalf = 1.e4;
minevalf = 1.e3;
evalf = [maxevalf minevalf];
[ integr , accur , varf , evalf , iter ] = intprb_montecarlo ( sumfunction , dim_num , evalf )

// Test verbose
grand("setsd",123456);
dim_num = 10;
bounds = [];
evalf = [];
tol = [];
randgen = [];
integr = intprb_montecarlo ( sumfunction , dim_num , evalf , bounds , tol , randgen , %t )

// Test tolrel, tolabs
grand("setsd",123456);
dim_num = 10;
bounds = [];
evalf = [1.e5 1.e3];
tol = [1.e-3 1.e-2];
[ integr , accur , varf , evalf , iter ] = intprb_montecarlo ( sumfunction , dim_num , evalf , ..
bounds , tol )

// Customize the bounds
// integral_10^30 sqrt(x)  = 88.46266043324404
// E(f) = integral_10^30 sqrt(x)  * 1/20 = 4.423133021662201
// V(f) = integral_10^30 (sqrt(x) - 4.423133021662201 )^2 * 1/20 = 0.4358942726814049
function r = sqrtfunction (m,n,x)
r = sqrt(x)
endfunction
grand("setsd",123456);
dim_num = 1;
evalf = [];
bounds = [10 30];
[ integr , accur ] = intprb_montecarlo ( sqrtfunction , dim_num , evalf , bounds )

// Customize the random number generator
function r = shrand ( m , n )
global shrageseed
r = zeros(m,n)
for k = 1 : m
for i = 1 : n
[ value , shrageseed ] = intprb_shragerandom ( shrageseed );
r(k,i) = value
end
end
endfunction
global shrageseed;
shrageseed = 123456;
dim_num = 3;
evalf = [1000 1000];
bounds = [];
tol = [];
[ integr , accur ] = intprb_montecarlo ( sumfunction , dim_num , evalf , bounds , tol , shrand )

// Customize the random number generator :
// configure a random number generator which has optionnal arguments.
function u = mygen ( m , n , seed , gentype )
grand ( "setgen" , gentype )
grand ( "setsd" , seed )
u = grand ( m , n , "def" )
endfunction
backgen = grand("getgen");
dim_num = 3;
evalf = [1000 1000];
bounds = [];
tol = [];
randgen = list(mygen,123456,"urand");
[ integr , accur ] = intprb_montecarlo ( sumfunction , dim_num , evalf , bounds , tol , randgen )
// Configure back the default generator.
grand("setgen",backgen);

// Customize the function :
// use a function which has additionnal parameters.
function r = sumfunction2 (m,n,x,e,v)
r = sqrt(1/v) * (sum(x,"c") - e)
endfunction
dim_num = 3;
func = list(sumfunction2,dim_num/2,dim_num/12);
[ integr , accur , varf , evalf , iter ] = intprb_montecarlo ( func , dim_num )

   

Authors

Michael Baudin - 2010 - DIGITEO

Bibliography

"Methods of Numerical Integration", Second Edition, Philip Davis, Philip Rabinowitz, Dover, 2007,

GSL/src/monte/plain.c, "Gnu Scientific Library"

genz/software/fort77/mvtexoh.f/RCRUDE, "Crude Monte-Carlo Algorithm with simple antithetic variates and weighted results on restart", Alan Genz