intprb_montecarlo — 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 name of the function to be evaluated.
a 1 x 1 matrix of floating point integers, the spatial dimension.
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)])
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)
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)
a function, the random number generator. (default = grand)
a 1 x 1 matrix of boolean, set to %t to enable verbose messages. (default = %f)
a 1 x 1 matrix of doubles, the approximate value of the integral, the mean of the function.
a 1 x 1 matrix of doubles, the estimated error on the integral.
a 1 x 1 matrix of doubles, the approximate value of the variance.
a 1 x 1 matrix of floating point integers, the actual number of calls to the function.
a 1 x 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. 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).
// 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 )