Estimates a multidimensional integral using a Low Discrepancy sequence.
integr = intprb_crudeld ( func , n ) integr = intprb_crudeld ( func , n , ldseq ) integr = intprb_crudeld ( func , n , ldseq , callf ) integr = intprb_crudeld ( func , n , ldseq , callf , bounds )
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 1-by-1 matrix of strings, a function or a list, the name of the sequence (Default = "sobol"). The list of available sequences can be computed by lowdisc_methods(). The ldseq variable might also be a function, or a list. See below for details.
a 1-by-1 matrix of floating point integers, the number of calls to the function (default = 1.e3)
a n-by-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)])
a 1-by-1 matrix of doubles, the approximate value of the integral, the mean of the function.
The algorithm uses a crude Low Discrepancy sequence to estimate the integral.
The function func
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-by-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 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 rg.
If an optional input argument is equal to the empty matrix, it is set to its default value.
The low discrepancy generator ldseq may also be a function, with header
u = ldseq ( m , n )
where m is a 1-by-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 low discrepancy generator requires additionnal arguments to be evaluated. In this case, we can use the following feature. The argument ldseq can also be a list, where the first item is a function rg. The function rg must have the header :
u = rg ( m , n , a1 , a2 , ... )
In this case the ldseq 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.
// 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 n = 10; // The exact answer is 0. integr = intprb_crudeld ( sumfun , n ) // Customize the Low Discrepancy sequence dim_num = 10; ldseq = "haltonf"; // The exact answer is 0. integr = intprb_crudeld ( sumfun , dim_num , ldseq ) // 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=sqrtfun(m, n, x) r = sqrt(x) endfunction n = 1; callf = []; ldseq = []; bounds = [10 30]; // The exact answer is 88.46266043324404. integr=intprb_crudeld(sqrtfun,n,callf,ldseq,bounds) // Customize the low Discrepancy sequence. // Uses our own generator (here, a basic random number generator). // Caution : this is a wrong use of the function, which expects // a Low Discrepancy sequence. // Customize ldseq with a home-made Halton sequence // restricted to dimension 10. function result=vdc(i, basis) current = i ib = 1.0 / basis result = 0.0 while (current>0) digit = modulo ( current , basis ) current = int ( current / basis ) result = result + digit * ib ib = ib / basis end endfunction function next=haltonseq(m, n) next = zeros ( m , n ) primematrix = [2,3,5,7,11,13,17,19,23,29] for k = 1 : m for i = 1 : n basis = primematrix ( i ) next(k,i) = vdc ( k , basis ) end end endfunction dim_num = 3; callf = 100; integr=intprb_crudeld(sumfun,dim_num,haltonseq,callf) // Customize ldseq, which has additionnal parameters function next=haltonseq2(m, n, skip, leap, prmat) next = zeros ( m , n ) k = skip + 1 index = 0 for index = 1 : m for i = 1 : n basis = prmat ( i ) next(index,i) = vdc ( k , basis ) end k = k + leap end endfunction dim_num = 3; primematrix = [2,3,5,7,11,13,17,19,23,29]; skip = 10; // Kocis and Whiten recommend leap = 409. leap = 409; ldseq = list(haltonseq2,skip,leap,primematrix); callf = 100; integr=intprb_crudeld(sumfun,dim_num,ldseq,callf) // 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 n = 3; func = list(sumfun2,n/2,n/12); integr = intprb_crudeld ( func , n ) | ![]() | ![]() |
"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