蒙特卡罗积分

    本章描述了多维蒙特卡罗积分的例程。其中包括传统的蒙特卡罗方法和自适应算法，如VEGAS和MISER，这些算法使用重要抽样和分层抽样技术。每个算法计算一个多维定积分的估计值，

     在一个超立方区域((xl, xu)，(yl, yu)，…)上使用固定数量的函数调用。例程还提供了结果误差的统计估计。这个误差估计应该作为一个指导，而不是作为一个严格的误差界-随机抽样的区域，可能不会发现函数所有重要特性，导致误差的低估。

    这些函数定义在每个例程的独立头文件中，分别是gsl_monte_plain.h、gsl_monte_master .h和gsl_monte_vegas.h。

27.1 接口

所有的蒙特卡罗积分程序都使用相同的通用接口形式。有一个分配器来为控制变量和工作空间分配内存，有一个例程来初始化这些控制变量，积分器本身，以及一个函数来释放空间。

    每个积分函数都需要提供一个随机数生成器，并返回积分及其标准差的估计值。结果的准确性由用户指定的函数调用次数决定。如果需要已知的精度水平，可以通过多次调用积分器并对单个结果进行平均，直到获得所需的精度来实现。

    蒙特卡罗例程中使用的随机抽样点总是严格地选择在积分区域内，从而自动避免了端点奇点。

要积分的函数有自己的数据类型，定义在头文件gsl_mon.h中。

gsl_monte_function

    本数据类型定义了一个带有参数的通用函数用于蒙特卡罗积分。

double (* f) (double * x, size_t dim, void * params)	本函数将返回参数x和参数params的值f(x, params)，其中x是一个大小为dim的数组，给出了函数将要计算的点的坐标。
size_t dim	x的维数。
void * params	一个指向函数参数的指针。

这是一个二维二次函数的例子，

     其中a=3，b=2，c=1。下面的代码定义了一个gsl_monte_function F，可以传递给一个积分器:

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; 

函数f(x)可以通过下面的宏定义赋值，

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

27.2 普通蒙特卡洛

普通蒙特卡罗算法从积分区域随机采样点来估计积分及其误差。使用该算法，可以得出N个随机分布点xi的积分E(f;N)的估计值，

     其中V是积分区域的体积。此估计的误差σ(E;N)由均值的估计方差计算，

     对于大N，这个方差随着Var(f)/N渐近减小，其中Var(f)是函数在积分区域上的真实方差。误差估计本身应减小为σ(f)/N。我们熟悉的误差随1/N的增加而减少的规律——要将误差减少10倍，就需要增加100倍的样本点。

本节中描述的函数声明在头文件gsl_monte_plain.h中。

gsl_monte_plain_state

    这是普通蒙特卡洛积分的工作空间。

gsl_monte_plain_state * gsl_monte_plain_alloc(size_t dim)

    本函数分配并初始化一个工作空间，用于在dim维数下进行蒙特卡罗积分。

int gsl_monte_plain_init(gsl_monte_plain_state* s)

    本函数初始化先前分配的积分状态。这允许现有的工作区用于不同的积分。

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)

本函数使用普通蒙特卡罗算法,在由数组xl和xu的上下限定义的dim维超立方区域上对函数f进行积分，每个数组的大小为dim。该积分使用固定数量的函数调用calls，并使用随机数生成器r获得随机采样点。必须提供先前分配的工作空间s。积分的结果以result的形式返回，并带有一个估计的绝对误差abserr。

void gsl_monte_plain_free(gsl_monte_plain_state * s)

    本函数释放与积分器状态s相关的内存。

27.3 MISER

    Press和Farrar的MISER算法是基于递归分层抽样的。该方法的目的是通过将积分点集中在方差最大的区域来减小总体积分误差。