intprb_crudeld — 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 )
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.
a 1 x 1 matrix of floating point integers, the spatial dimension.
a 1 x 1 matrix of strings, the name of the sequence (Default = "sobolf"). The list of available sequences can be computed by lowdisc_methods(). The ldseq variable might also be a function, or a list.
a 1 x 1 matrix of floating point integers, the number of calls to the function (default = 1.e3)
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)])
a 1 x 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 "lowdisc" module must be installed to use this function. To install this module, run
atomsInstall("lowdisc")
then restart Scilab.
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.
If an optional input argument is equal to the empty matrix, it is set to its default value.
The low discrepancy generator may also be a function, with header u = ldseq ( 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 low discrepancy generator requires additionnal arguments to be evaluated. In this case, we can use the following feature. The function ldseq can also be a list, with 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 = 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 n = 10; // The exact answer is 0. integr = intprb_crudeld ( sumfunction , n ) // Customize the Low Discrepancy sequence dim_num = 10; ldseq = "haltonf"; // The exact answer is 0. integr = intprb_crudeld ( sumfunction , 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) - 4.423133021662201 )^2 * 1/20 = 0.4358942726814049 function r = sqrtfunction (m,n,x) r = sqrt(x) endfunction n = 1; callf = []; ldseq = []; bounds = [10 30]; // The exact answer is 88.46266043324404. integr = intprb_crudeld ( sqrtfunction , 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 = haltonsequence ( 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 ( sumfunction , dim_num , haltonsequence , callf ) // Customize ldseq, which has additionnal parameters function next = haltonsequence2 ( 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(haltonsequence2,skip,leap,primematrix); callf = 100; integr = intprb_crudeld ( sumfunction , dim_num , ldseq , callf )