12 重要抽样法

\(I = \int_C h(\boldsymbol x) \,d \boldsymbol x\)中积分区域\(C\)可能是任意形状的, 也可能无界，\(h(\boldsymbol x)\)在\(C\)内各处的取值大小差异可能很大， 使得直接用平均值法估计\(I\)时， 很多样本点处于\(|h(\boldsymbol x)|\)接近于零的地方， 造成浪费，另外使得\(\hat I_2\)的渐近方差(见式(11.23))中的 \(\text{Var}(h(\boldsymbol U))\)很大(\(\boldsymbol U \sim \text{U}(C)\))。 为此，考虑用非均匀抽样: \(|h(\boldsymbol x)|\)大的地方密集投点， \(|h(\boldsymbol x)|\)小的地方稀疏投点。 这样可以有效利用样本。 设\(g(\boldsymbol x), \boldsymbol x \in C\)是一个密度， 其形状与\(|h(\boldsymbol x)|\)相近， 且当\(g(\boldsymbol x)=0\)时\(h(\boldsymbol x)=0\)， 当\(\| \boldsymbol x \| \to \infty\)时 \(h(\boldsymbol x) = o(g(\boldsymbol x))\)。 称\(g(\boldsymbol x)\)为试投密度或重要抽样密度。

设\(\boldsymbol X_i\) iid \(\sim g(\boldsymbol x)\), \(i=1,2,\ldots, N\)。 令 \[\begin{aligned} \eta_i = \frac{h(\boldsymbol X_i)}{g(\boldsymbol X_i)}, \ i=1,2,\dots,N, \end{aligned}\] 则 \[\begin{aligned} E \eta_1 = \int_C \frac{h(\boldsymbol x)}{g(\boldsymbol x)} g(\boldsymbol x) \,d \boldsymbol x = \int_C h(\boldsymbol x) \,d \boldsymbol x = I, \end{aligned}\] 因此可以用\(\{ \eta_i, i=1,2,\dots, N \}\)的样本平均值来估计\(I\)，即 \[\begin{align} \hat I_3 = \bar\eta = \frac{1}{N} \sum_{i=1}^N \frac{h(\boldsymbol X_i)}{g(\boldsymbol X_i)}. \tag{12.1} \end{align}\] \(\hat I_3\)是\(I\)的无偏估计和强相合估计。 \(\hat I_3\)的渐近方差为 \[\begin{align} \text{Var}(\hat I_3) = \text{Var}\left( \frac{h(\boldsymbol X)}{g(\boldsymbol X)} \right) \frac{1}{N} = O(\frac{1}{N}), \tag{12.2} \end{align}\] 当\(g(\boldsymbol x)\)和\(|h(\boldsymbol x)|\)形状很接近时\(\frac{|h(\boldsymbol x)|}{g(\boldsymbol x)}\)近似为常数， 方差\(\text{Var}(\hat I_3)\)很小, 这时\(\hat I_3\)的随机误差可以很小。 用(12.1)式估计\(I\)与舍选法II有类似的想法， 这种方法叫做重要抽样法(importance sampling)， 是随机模拟的重要方法。

为什么当\(g(\boldsymbol x)\)和\(|h(\boldsymbol x)|\)形状很接近时\(\text{Var}(\hat I_3)\)很小？ 事实上， \[\begin{align} \text{Var}\left( \frac{h(\boldsymbol X)}{g(\boldsymbol X)} \right) =& E \left( \frac{h^2(\boldsymbol X)}{g^2(\boldsymbol X)} \right) - \left( E \frac{h(\boldsymbol X)}{g(\boldsymbol X)} \right)^2 = E \left( \frac{h^2(\boldsymbol X)}{g^2(\boldsymbol X)} \right) - \left( \int_C h(\boldsymbol x) \,d\boldsymbol x \right)^2 \nonumber \\ \geq& \left( E \frac{|h(\boldsymbol X)|}{g(\boldsymbol X)} \right)^2 - I^2 \quad \text{(Jensen不等式)} \tag{12.3} \\ =& \left( \int_C |h(\boldsymbol x)| \,d\boldsymbol x \right)^2 - I^2 , \end{align}\] 当且仅当\(\frac{|h(\boldsymbol X)|}{g(\boldsymbol X)}\)为常数时(12.3)的不等式中的等号成立， 即\(g(\boldsymbol x) = \frac{1}{\int_C |h(\boldsymbol t)| \,d\boldsymbol t} |h(\boldsymbol x)|, \boldsymbol x \in C\) 时\(\text{Var}(\hat I_3)\)达到最小。

上述求积分的问题也可以表述为求随机变量函数的期望问题。 设\(\boldsymbol Y \sim f(\boldsymbol y)\)， 要求\(\boldsymbol Y\)的函数\(h(\boldsymbol Y)\) 的期望值\(E h(\boldsymbol Y) = \int h(\boldsymbol y) f(\boldsymbol y) \,d\boldsymbol y\), 可以抽取\(\boldsymbol Y\)的随机数\(\boldsymbol Y_i, i=1,2,\dots,N\), 然后用平均值法 \(\frac{1}{N} \sum_{i=1}^N h(\boldsymbol Y_i)\)估计\(E h(\boldsymbol Y)\)。 但是，如果很难生成\(\boldsymbol Y\)的随机数， 或者\(\boldsymbol Y\)的分布集中于\(h(\boldsymbol y)\)接近于零的位置以至于积分效率很低， 可以找一个试投密度\(g(\boldsymbol x)\)， 从\(g(\boldsymbol x)\)产生随机数\(\boldsymbol X_i, i=1,2,\dots,N\), 用如下重要抽样法 \[\begin{align} \hat I_{3.1} =& \frac{1}{N} \sum_{i=1}^N h(\boldsymbol X_i) \frac{f(\boldsymbol X_i)}{g(\boldsymbol X_i)} \tag{12.4} \end{align}\] 估计\(E h(\boldsymbol Y)\), 易见\(E \hat I_{3.1} = E h(\boldsymbol Y)\)。 称\(W_i = \frac{f(\boldsymbol X_i)}{g(\boldsymbol X_i)}\)为重要性权重, \(E h(\boldsymbol Y)\)可以用重要性权重估计为 \[\begin{align} \hat I_{3.1} =& \frac{1}{N} \sum_{i=1}^N W_i h(\boldsymbol X_i) . \tag{12.5} \end{align}\]

选取试投密度\(g(\boldsymbol x)\)时， 要求当\(h(\boldsymbol x) f(\boldsymbol x) \neq 0\)时\(g(\boldsymbol x) \neq 0\), 当\(\boldsymbol x\)趋于无穷时 \(h(\boldsymbol x) f(\boldsymbol x) = o(g(\boldsymbol x))\)（\(g(\boldsymbol x)\)相对厚尾）, 一般还要求\(g(\boldsymbol x)\)的形状与\(|h(\boldsymbol x)| f(\boldsymbol x)\)形状接近。 如果\(g(\boldsymbol x)\)的相对厚尾性难以确定， 可以使用如下保险的试投密度 \[\begin{align} \tilde g(\boldsymbol x) = \rho g(\boldsymbol x) + (1-\rho) r(\boldsymbol x), \tag{12.6} \end{align}\] 其中\(\rho\)接近于1， \(r(\boldsymbol x)\)是柯西分布或Pareto分布这样的重尾分布。 要生成\(N\)个\(\tilde g(\boldsymbol x)\)的随机数， 只要生成\(N \rho\)个\(g(\boldsymbol x)\)的随机数和\(N(1-\rho)\)个\(r(\boldsymbol x)\)的随机数。

选取不适当的试投密度会把绝大多数样本点投到了对计算积分不重要的位置， 使得样本点中只有极少数点是真正有作用的。 在多维问题中合适的试投密度尤其难找， 经常需要反复试验。

例12.1 用MC积分法计算 \(I = \int_0^1 e^x \,dx = e-1 \approx 1.718\)。 对被积函数\(h(x)=e^x\)做泰勒展开得 \[\begin{aligned} e^x = 1 + x + \frac{x^2}{2!} + \cdots \end{aligned}\]

取 \[\begin{aligned} g(x)=c(1+x) = \frac{2}{3}(1+x) \end{aligned}\] 要产生\(g(x)\)的随机数可以用逆变换法， 密度\(g(x)\)的分布函数\(G(x)\)的反函数为 \[\begin{aligned} G^{-1}(y) = \sqrt{1+3y}-1, \ 0