蒙特卡罗积分中的准随机序列与自适应方法

### 蒙特卡罗积分中的准随机序列与自适应方法 #### 1. 准随机序列概述 在蒙特卡罗积分中，在n维空间均匀随机选择N个点会导致误差项以$1/\sqrt{N}$的速度减小。然而，这并非不可避免，例如在笛卡尔网格上选择样本点，其分数误差至少以$N^{-1}$的速度减小。但网格存在问题，需要提前确定其精细程度，且难以根据收敛或终止准则进行采样。 于是，人们提出了准随机序列的概念，它是n元组序列，比不相关的随机点更均匀地填充n维空间。实际上，准随机序列并非真正的随机，而是经过精心设计的次随机序列，其样本点在某种精确意义上“最大限度地相互避开”。 #### 2. 哈尔顿序列 哈尔顿序列是一个概念上简单的准随机序列示例。在一维情况下，序列中的第j个数$H_j$通过以下步骤获得： 1. 将j写成以素数b为基数的数。例如，当$b = 3$，$j = 17$时，写成122。 2. 反转数字并在序列前加上小数点（以b为基数的小数点）。在上述例子中，得到$0.221_{(3)}$，结果即为$H_j$。 要在n维空间中获得n元组序列，每个分量使用具有不同素数基数b的哈尔顿序列，通常使用前n个素数。 哈尔顿序列的工作原理是，每当j的位数增加一位时，j的反转数字分数会变得更精细。这一过程相当于在一系列越来越精细的笛卡尔网格上填充所有点，并且在每个网格上以一种最大程度分散的顺序进行。 #### 3. 索博尔序列 除了哈尔顿序列，还有其他生成准随机序列的方法，如Faure、Sobol’、Niederreiter等人提出的方法。这里重点介绍Sobol’序列的Antonov - Saleev变体。 ##### 3.1 基本原理 Sobol’序列直接将0到1之间的数字生成长度为w位的二进制分数，基于一组w个特殊的二进制分数$V_i$（$i = 1, 2, ..., w$），称为方向数。在Sobol的原始方法中，第j个数$X_j$通过对满足“j的第i位非零”这一条件的$V_i$进行按位异或（XOR）运算得到。随着j的增加，不同的$V_i$在不同时间尺度上在$X_j$中进出。$V_1$交替出现和消失的速度最快，而$V_k$每$2^{k - 1}$步才从出现变为消失（或反之）。 Antonov和Saleev的贡献是指出，可以使用j的格雷码$G(j)$的位来选择方向数，而不是使用整数j的位。由于$G(j)$和$G(j + 1)$仅在j的二进制表示中最右边的零位位置不同，因此第$j + 1$个Sobol’ - Antonov - Saleev数可以通过将第j个数与单个$V_i$进行异或运算得到，其中i是j中最右边零位的位置，这使得序列的计算非常高效。 ##### 3.2 方向数的生成 每个不同的Sobol’序列（或n维序列的分量）基于模2整数上的不同本原多项式。设P是一个q次的本原多项式： $P = x^q + a_1x^{q - 1} + a_2x^{q - 2} + ... + a_{q - 1} + 1$ 定义整数序列$M_i$的q项递推关系为： $M_i = 2a_1M_{i - 1} \oplus 2^2a_2M_{i - 2} \oplus ... \oplus 2^{q - 1}M_{i - q + 1}a_{q - 1} \oplus (2^qM_{i - q} \oplus M_{i - q})$ 其中按位异或用$\oplus$表示。递推的起始值$M_1, ..., M_q$可以分别是小于$2, ..., 2^q$的任意奇数。然后，方向数$V_i$由下式给出： $V_i = M_i/2^i$，$i = 1, ..., w$ 以下是模2本原多项式（次数$q \leq 10$）的表格： | 次数 | 本原多项式 | | ---- | ---- | | 1 | 0 (即$x + 1$) | | 2 | 1 (即$x^2 + x + 1$) | | 3 | 1, 2 (即$x^3 + x + 1$和$x^3 + x^2 + 1$) | | 4 | 1, 4 (即$x^4 + x + 1$和$x^4 + x^3 + 1$) | | 5 | 2, 4, 7, 11, 13, 14 | | 6 | 1, 13, 16, 19, 22, 25 | | 7 | 1, 4, 7, 8, 14, 19, 21, 28, 31, 32, 37, 41, 42, 50, 55, 56, 59, 62 | | 8 | 14, 21, 22, 38, 47, 49, 50, 52, 56, 67, 70, 84, 97, 103, 115, 122 | | 9 | 8, 13, 16, 22, 25, 44, 47, 52, 55, 59, 62, 67, 74, 81, 82, 87, 91, 94, 103, 104, 109, 122, 124, 137, 138, 143, 145, 152, 157, 167, 173, 176, 181, 182, 185, 191, 194, 199, 218, 220, 227, 229, 230, 234, 236, 241, 244, 253 | | 10 | 4, 13, 19, 22, 50, 55, 64, 69, 98, 107, 115, 121, 127, 134, 140, 145, 152, 158, 161, 171, 181, 194, 199, 203, 208, 227, 242, 251, 253, 265, 266, 274, 283, 289, 295, 301, 316, 319, 324, 346, 352, 361, 367, 382, 395, 398, 400, 412, 419, 422, 426, 428, 433, 446, 454, 457, 472, 493, 505, 508 | ##### 3.3 索博尔序列的实现 以下是实现Sobol’序列的代码： ```c #include "nrutil.h" #define MAXBIT 30 #define MAXDIM 6 void sobseq(int *n, float x[]) { int j,k,l; unsigned long i,im,ipp; static float fac; static unsigned long in,ix[MAXDIM+1],*iu[MAXBIT+1]; static unsigned long mdeg[MAXDIM+1]={0,1,2,3,3,4,4}; static unsigned long ip[MAXDIM+1]={0,0,1,1,2,1,4}; static unsigned long iv[MAXDIM*MAXBIT+1]={ 0,1,1,1,1,1,1,3,1,3,3,1,1,5,7,7,3,3,5,15,11,5,15,13,9}; if (*n < 0) { for (k=1;k<=MAXDIM;k++) ix[k]=0; in=0; if (iv[1] != 1) return; fac=1.0/(1L << MAXBIT); for (j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM) iu[j] = &iv[k]; for (k=1;k<=MAXDIM;k++) { for (j=1;j<=mdeg[k];j++) iu[j][k] <<= (MAXBIT-j); for (j=mdeg[k]+1;j<=MAXBIT;j++) { ipp=ip[k]; i=iu[j-mdeg[k]][k]; i ^= (i >> mdeg[k]); for (l=mdeg[k]-1;l>=1;l--) { if (ipp & 1) i ^= iu[j-l][k]; ipp >>= 1; } iu[j][k]=i; } } } else { im=in++; for (j=1;j<=MAXBIT;j++) { if (!(im & 1)) break; im >>= 1; } if (j > MAXBIT) nrerror("MAXBIT too small in sobseq"); im=(j-1)*MAXDIM; for (k=1;k<=IMIN(*n,MAXDIM);k++) { ix[k] ^= iv[im+k]; x[k]=ix[k]*fac; } } } ``` 当n为负数时，该函数内部初始化一组MAXBIT个方向数，用于MAXDIM个不同的Sobol’序列。当n为正数（但$\leq$MAXDIM）时，返回这些序列中n个序列的下一个值。 ##### 3.4 索博尔序列的性能 对于n维空间中光滑函数的蒙特卡罗积分，Sobol’序列的分数误差随样本数量N的增加以$(\ln N)^n/N$的速度减小，几乎与$1/N$一样快。例如，对三维空间中圆环内非零函数的积分实验表明，Sobol’序列在几千个样本时就能达到1%的精度，而伪随机采样则需要近100,000个样本。 然而，当被积函数在采样区域内有硬边界（不连续边界）时，Sobol’序列的优势会减弱。例如，圆环内为1，圆环外为0的函数，Sobol’序列的分数误差约为$N^{-2/3}$，但与独立随机点相比，在中等精度（如1%）下仍有显著优势，在更高精度下优势更大。 #### 4. 拉丁超立方体采样 拉丁超立方体采样是一种在N维空间中极其稀疏地采样M个点的技术。例如，要同时测试汽车的4个不同设计参数的耐撞性，但只有3辆可消耗的汽车。 其思想是将每个设计参数（维度）划分为M个段，使整个空间划分为$M^N$个单元格。然后通过以下算法选择M个单元格包含样本点： 1. 随机选择$M^N$个单元格中的一个作为第一个点。 2. 消除与该点在任何参数上相同的所有单元格（即划掉同一行、列等中的所有单元格），留下$(M - 1)^N$个候选单元格。 3. 随机选择其中一个，消除新的行和列，继续此过程，直到只剩下一个单元格，该单元格包含最终的样本点。 这种构造的结果是每个设计参数的每个子范围都得到了测试。如果系统的响应由某个设计参数主导，该参数可以通过这种采样技术找到。但如果不同设计参数之间存在重要的相互作用，拉丁超立方体采样则没有特别的优势，使用时需谨慎。 以下是上述内容的mermaid流程图： ```mermaid graph TD; A[蒙特卡罗积分] --> B[均匀随机选点]; A --> C[准随机序列]; C --> D[哈尔顿序列]; C --> E[索博尔序列]; E --> F[基本原理]; E --> G[方向数生成]; E --> H[实现代码]; E --> I[性能评估]; C --> J[拉丁超立方体采样]; ``` #### 5. 自适应和递归蒙特卡罗方法概述 接下来介绍更高级的蒙特卡罗积分技术，包括重要性采样、分层采样、混合策略以及具体的自适应蒙特卡罗算法VEGAS和递归分层采样算法MISER。这些技术都属于方差缩减的范畴，但各有特点。 #### 6. 重要性采样 重要性采样在之前的方程中已有隐含体现。假设被积函数$f$可以写成几乎恒定的函数$h$与另一个正函数$g$的乘积，即$\int f dV = \int (f/g) g dV = \int h g dV$。 可以通过采样$h$来计算$f$的积分，但不是以均匀概率密度$dV$采样，而是以非均匀密度$g dV$采样。更一般地，对于在体积$V$内以概率密度$p$选择的点$x_i$，满足$\int p dV = 1$，积分的估计公式为： $I \equiv \int f dV = \int \frac{f}{p} p dV \approx \left\langle\frac{f}{p}\right\rangle \pm \sqrt{\frac{\left\langle f^2/p^2\right\rangle - \left\langle f/p\right\rangle^2}{N}}$ 其中角括号表示对$N$个点的算术平均值。 最佳采样密度$p$的选择是使$h = f/p$尽可能接近常数。通过泛函变分可以得到最优的$p$与$|f|$成正比，即$p = \frac{|f|}{\int |f| dV}$。当$f$在积分区域内只有一个符号时，方差可降为零。但实际中，精确知道$p$需要先知道$\int |f| dV$，这相当于已经知道要计算的积分。 如果函数$f$在大部分体积$V$中取已知的常数值，添加一个常数使该值为零是个好主意。重要性采样的实际精度取决于可实现的$p$与理想$p$的接近程度。 #### 7. 分层采样 分层采样的思想与重要性采样不同。设$\langle\langle f\rangle\rangle$表示函数$f$在体积$V$上的真实平均值（即积分除以$V$），$\langle f\rangle$表示最简单的（均匀采样）蒙特卡罗估计值： $\la