Estimates a multidimensional integral using iterated Monte Carlo.
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 ( ... )
a function or a list, 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. See below for details.
a 1-by-1 matrix of floating point integers, the spatial dimension.
a n-by-2 matrix of doubles, where bounds(:,1) = lower bound and bounds(:,2) = upper bound of integration (default = [zeros(n,1) ones(n,1)])
a 2-by-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)
a 2-by-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)
a function or a list, the random number generator (default = grand). See below for details.
a 1-by-1 matrix of boolean, set to %t to enable verbose messages. (default = %f)
a 1-by-1 matrix of doubles, the approximate value of the integral, the mean of the function.
a 1-by-1 matrix of doubles, the estimated error on the integral.
a 1-by-1 matrix of doubles, the approximate value of the variance.
a 1-by-1 matrix of floating point integers, the actual number of calls to the function.
a 1-by-1 matrix of floating point integers, the number of loops in the algorithm
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.
It might happen that the function requires additionnal arguments to be evaluated. In this case, we can use the following feature. The argument func can also be a list, where the first item is a function f. The function f must then have the 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 f.
The uniform random number generator must have header
y = randgen ( m , n )
where m is a 1-by-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].
It might happen that the random number generator requires additionnal arguments to be evaluated. In this case, we can use the following feature. The argument randgen can also be a list, where the first element rg is a function. This function rg must have the 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.
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).
// A basic example function r=sumfun(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(sumfun,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(sumfun,dim_num,evalf) // Test verbose grand("setsd",123456); dim_num = 10; bounds = []; evalf = []; tol = []; randgen = []; integr = intprb_montecarlo ( sumfun,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(sumfun,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) - E(f) )^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(sumfun,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(sumfun,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=sumfun2(m, n, x, e, v) r = sqrt(1/v) * (sum(x,"c") - e) endfunction dim_num = 3; func = list(sumfun2,dim_num/2,dim_num/12); [integr,accur,varf,evalf,iter]=intprb_montecarlo(func,dim_num ) | ![]() | ![]() |
"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