Next: , Up: Monte Carlo Integration

23.1 Interface

All of the Monte Carlo integration routines use the same general form of interface. There is an allocator to allocate memory for control variables and workspace, a routine to initialize those control variables, the integrator itself, and a function to free the space when done.

Each integration function requires a random number generator to be supplied, and returns an estimate of the integral and its standard deviation. The accuracy of the result is determined by the number of function calls specified by the user. If a known level of accuracy is required this can be achieved by calling the integrator several times and averaging the individual results until the desired accuracy is obtained.

Random sample points used within the Monte Carlo routines are always chosen strictly within the integration region, so that endpoint singularities are automatically avoided.

The function to be integrated has its own datatype, defined in the header file gsl_monte.h.

— Data Type: gsl_monte_function

This data type defines a general function with parameters for Monte Carlo integration.

double (* f) (double * x, size_t dim, void * params)
this function should return the value f(x,params) for the argument x and parameters params, where x is an array of size dim giving the coordinates of the point where the function is to be evaluated.
size_t dim
the number of dimensions for x.
void * params
a pointer to the parameters of the function.

Here is an example for a quadratic function in two dimensions, 

     f(x,y) = a x^2 + b x y + c y^2

with a = 3, b = 2, c = 1. The following code defines a gsl_monte_function F which you could pass to an integrator: 

     struct my_f_params { double a; double b; double c; };
     
     double
     my_f (double x[], size_t dim, void * p) {
        struct my_f_params * fp = (struct my_f_params *)p;
     
        if (dim != 2)
           {
             fprintf (stderr, "error: dim != 2");
             abort ();
           }
     
        return  fp->a * x[0] * x[0]
                  + fp->b * x[0] * x[1]
                    + fp->c * x[1] * x[1];
     }
     
     gsl_monte_function F;
     struct my_f_params params = { 3.0, 2.0, 1.0 };
     
     F.f = &my_f;
     F.dim = 2;
     F.params = &params;

The function f(x) can be evaluated using the following macro, 

     #define GSL_MONTE_FN_EVAL(F,x)
         (*((F)->f))(x,(F)->dim,(F)->params)