Name

intprb_crudemc — Estimates a multidimensional integral using Monte Carlo.

Calling Sequence

   integr = intprb_crudemc ( func , n )
   integr = intprb_crudemc ( func , n , callf )
   integr = intprb_crudemc ( func , n , callf , bounds )
   integr = intprb_crudemc ( func , n , callf , bounds , randgen )
   [ integr , accur , varf ] = intprb_crudemc ( ... )
   
   

Parameters

func :

the function to be evaluated. If func is a list, the first element is expected to be a function and the remaining elements of the list are input arguments of the function, which are appended at the end.

n:

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

callf :

a 1 x 1 matrix of floating point integers, the number of calls to the function (default = 1.e4)

bounds :

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

randgen :

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

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 of f.

Description

The algorithm uses a crude Monte-Carlo method to estimate the integral.

The function should have header y = func( m , n , x ) where n is a floating point representing the spatial dimension, m is a floating point integer representing the number of experiments to perform and x is a m x n matrix of doubles.

It might happen that the function requires additionnal arguments to be evaluated. In this case, we can use the following feature. The function func can also be a list, with header y = f ( m , n , x , a1 , a2 , ... ). In this case the func variable should hold the list (f,a1,a2,...) and the input arguments a1, a2, will be automatically be appended at the end of the calling sequence of rg.

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

It might happen that the random number generator requires additionnal arguments to be evaluated. In this case, we can use the following feature. The function randgen can also be a list, with header u = rg ( m , n , a1 , a2 , ... ). In this case the randgen variable should hold the list (rg,a1,a2,...) and the input arguments a1, a2, will be automatically be appended at the end of the calling sequence of rg.

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

The error is estimated as accur = volume * varf / sqrt(callf), where varf is the standard deviation associated with the outputs of the function. The Monte-Carlo estimate follows a Normal (Laplace-Gauss) law, with standard deviation accur. The value of accur ensures that the probability that the actual integral is in the interval [integr-accur,integr+accur] is greater than 68%. The probability is 95% that the exact integral is in the range [I-2*accur,I+2*accur].

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

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
grand("setsd",123456);
n = 10;
// Exact integral is 0.
[ integr , accur , varf ] = intprb_crudemc ( sumfunction , n )

// 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);
n = 1;
callf = [];
lowb = 10;
uppb = 30;
bounds = [lowb uppb];
// Exact integral is 88.46266043324404.
[ integr , accur , varf ] = intprb_crudemc ( sqrtfunction , n , callf , 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;
n = 3;
callf = 100;
bounds = [];
[ integr , accur , varf ] = intprb_crudemc ( sumfunction , n , callf , bounds , 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
n = 3;
callf = 1000;
bounds = [];
randgen = list(mygen,123456,"urand");
[ integr , accur , varf ] = intprb_crudemc ( sumfunction , n , callf , bounds , randgen )

// 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
n = 3;
func = list(sumfunction2,n/2,n/12);
[ integr , accur , varf ] = intprb_crudemc ( sumfunction , n )

   

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