25.2 PLAIN Monte Carlo

The plain Monte Carlo algorithm samples points randomly from the integration region to estimate the integral and its error. Using this algorithm the estimate of the integral E(f; N) for N randomly distributed points x_i is given by,

E(f; N) = =  V <f> = (V / N) \sum_i^N f(x_i)

where V is the volume of the integration region. The error on this estimate \sigma(E;N) is calculated from the estimated variance of the mean,

\sigma^2 (E; N) = (V^2 / N^2) \sum_i^N (f(x_i) -  <f>)^2.

For large N this variance decreases asymptotically as \Var(f)/N, where \Var(f) is the true variance of the function over the integration region. The error estimate itself should decrease as \sigma(f)/\sqrt{N}. The familiar law of errors decreasing as 1/\sqrt{N} applies—to reduce the error by a factor of 10 requires a 100-fold increase in the number of sample points.

The functions described in this section are declared in the header file gsl_monte_plain.h.

Function: gsl_monte_plain_state * gsl_monte_plain_alloc (size_t dim)

This function allocates and initializes a workspace for Monte Carlo integration in dim dimensions.

Function: int gsl_monte_plain_init (gsl_monte_plain_state* s)

This function initializes a previously allocated integration state. This allows an existing workspace to be reused for different integrations.

Function: int gsl_monte_plain_integrate (gsl_monte_function * f, const double xl[], const double xu[], size_t dim, size_t calls, gsl_rng * r, gsl_monte_plain_state * s, double * result, double * abserr)

This routines uses the plain Monte Carlo algorithm to integrate the function f over the dim-dimensional hypercubic region defined by the lower and upper limits in the arrays xl and xu, each of size dim. The integration uses a fixed number of function calls calls, and obtains random sampling points using the random number generator r. A previously allocated workspace s must be supplied. The result of the integration is returned in result, with an estimated absolute error abserr.

Function: void gsl_monte_plain_free (gsl_monte_plain_state * s)

This function frees the memory associated with the integrator state s.