Name

intprb_crudeld — Estimates a multidimensional integral using a Low Discrepancy sequence.

Calling 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 )
   
   

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.

ldseq :

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.

callf :

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

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)])

integr :

a 1 x 1 matrix of doubles, the approximate value of the integral, the mean of the function.

Description

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.

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
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 )

   

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