统计计算四|蒙特卡罗方法（Monte Carlo Method）

系列文章目录

统计计算一|非线性方程的求解
统计计算二|EM算法（Expectation-Maximization Algorithm，期望最大化算法）
统计计算三|Cases for EM

---

一、基本概念

蒙特卡洛方法：为了解决某确定性问题，把它变成一个概率模型的求解问题，然后产生符合模型的大量随机数，对产生的随机数进行分析从而求解问题的方法，又称为随机模拟方法。

• 随机数：设 $X$ 是具有分布函数 $F(x)$ 的随机变量，从分布 $F(x)$ 中随机抽样得到的序列 { x i , i = 1 , 2 , . . . } \{x_i, i = 1, 2, ...\} {xi​,i=1,2,...} 称为该分布的随机数序列，$x_i$ 称为分布 $F(x)$ 的随机数。

（一）估算$\pi$

向正方形$D=\{(x,y):x\in [0,1],y\in [0,1]\}$内随机等可能投点，落入四分之一圆$C=\{(x,y):x^2+y^2\le 1,x>0,y>0\}$的概率为面积之比$p=\frac{\pi }{4}$。如果独立重复地投了$n$个点，落入$C$中的点的个数为$\xi$，则有：
$$\frac{\xi }{n}\approx \frac{\pi }{4},\pi \approx \hat{\pi }=\frac{4\xi }{n}$$

由于$\xi$服从$Binomial(n,\frac{\pi }{4})$分布，有：
$$Var(\hat{\pi })=\frac{\pi (4-\pi )}{16n}$$
由中心极限定理，$\hat{\pi }$近似服从$N(\pi ,\frac{\pi (4-\pi )}{16n})$分布，所以随机模拟误差的幅度大约在$\pm 2\sqrt{\frac{\pi (4-\pi )}{16n}}$（随机模拟误差95%以上落入此区间）

（二）求积分

将积分转化为期望来计算：对于$Q={\int }_a^{b}h(x)dx$，取$U\sim U(a,b)$，有：
$$Q=(b-a){\int }_a^{b}h(u)\frac{1}{b-a}du=(b-a)E[h(U)]$$
若取 { U i , i = 1 , . . . , N } \{U_i,i=1,...,N\} {Ui​,i=1,...,N}独立同$U(a,b)$分布，并设$Y_i=h(U_i),i=1,2,...,N$是$iid$随机变量列，由强大数律：
$$\bar{Y}=\frac{1}{N}\sum _{i=1}^{N}h(U_i)\to Eh(U)=\frac{Q}{b-a},a.s.(N\to \infty )$$
于是有：
$$\hat{Q}=(b-a)\bar{Y}=\frac{b-a}{N}\sum _{i=1}^{N}h(U_i)$$

由中心极限定理：
$$\sqrt{N}(\hat{Q}-Q)\stackrel{d}{\to }N(0,(b-a{)}^{2}Var(h(U)))$$
$$Var[h(U)]={\int }_a^b[h(u)-Eh(U){]}^2\frac{1}{b-a}du$$
$Var[(h(U)]$可以用模拟样本 { Y i = h ( U i ) } \{Y_i=h(U_i)\} {Yi​=h(Ui​)}估计为：
$$Var(h(U))\approx \frac{1}{N}\sum _{i=1}^N(Y_i-\bar{Y}{)}^2$$

（三）使用步骤

蒙特卡洛方法的理论基础是大数定律。样本数量越多，则随机数的平均值就越接近期望，也就是要计算的真实值。

• 将实际问题转化为求期望，并定义要采样的随机变量
• 计算机模拟采样过程，处理产生的随机数得到期望

二、变换采样

（一）基本概念

如果随机变量 $\eta$ 不容易抽样，但是存在另一个容易抽样的随机变量 $\xi$ 和随机变量 $\eta$ 间具有一一对应关系，即 $\eta =h(\xi )$ 或 $\xi =h^{-1}(\eta )$，同分布。那么可以先产生随机变量 $\xi$，再由函数关系 $h(\cdot )$ 得到随机变量 $\eta$，这种产生随机数的方法称为变换抽样法。

（二）相关定理

设随机变量$\xi$具有概率密度函数$f(x)$，另有一函数$h(\cdot )$严格单调，其反函数记为$h^{-1}(\cdot )$且导函数存在，则$\eta =h(\xi )$是随机变量$\xi$的函数，其概率密度函数为：
p ( z ) = f ( h − 1 ( z ) ) ⋅ ∣ { h − 1 ( z ) } ′ ∣ p(z)=f(h^{-1}(z))·|\{h^{-1}(z)\}'| p(z)=f(h−1(z))⋅∣{h−1(z)}′∣

证明：

（三）示例

• 用变换抽样法产生分布为 $N(\mu ,\sigma 2)$ 的随机数。
  

• 用变换抽样法产生分布为 $Gamma(1/2,3)$ 的随机数。
  

三、逆变换采样

（一）相关定理

假设$X$为一个连续随机变量，其累计分布函数为$F_X$，此时可证明随机变量$Y=F_X(X)$服从区间$[0,1]$上的均匀分。逆变换采样就是将上述过程反过来进行。

设连续型随机变量$\eta$的分布函数$F(x)$是连续且严格单调上升的分布函数，其反函数存在且记为$F^{-1}(x)$。则有：

• 随机变量$F(\eta )$服从$(0,1)$上的均匀分布，即$F(\eta )\sim U(0,1)$
• 对于随机变量$U\sim U(0,1)$，$F^{-1}(U)$的分布函数为$F(x)$

证明：

（二）逆变换法

逆变换法：当随机变量$\eta$的分布函数$F(x)$的反函数存在，且容易计算时，可通过产生均匀分布的随机数来产生$\eta$的随机数序列 { η i , i = 1 , 2 , . . . } \{\eta_i,i=1,2,...\} {ηi​,i=1,2,...}。这种产生非均匀分布随机数的方法称为逆变换法或反函数法。

（三）采样步骤

• 产生$U(0,1)$的随机数序列 { u i , i = 1 , 2 , . . . } \{u_i,i=1,2,...\} {ui​,i=1,2,...}
• $\eta$的随机数序列为：
  $${\eta }_i=F^{-1}(u_i),i=1,2,...$$

（四）示例

• 产生概率密度函数为 f(x) 的随机数，其中：
  $$f(x)=\{\begin{array}{ll}\frac{x}{{\sigma }^2}{e}^{-\frac{x^2}{2{\sigma }^2}},& \text{x>0}\\ 0,& \text{z≤0}\end{array}$$
  
• 产生分布函数为$F(x)$ 的随机数 $\eta$，其中
  $$F(x)=\frac{x^2+x}{2},0\le x\le 1$$
  

四、接受拒绝采样

（一）基本概念

拒绝抽样是基于以下观察而提出的：要在一维中抽样一个随机变量，可以对二维笛卡尔图进行均匀随机抽样，并将样本保留在其密度函数图形下的区域中。

想象将一个随机变量的密度函数绘制在一个大矩形板上，并向其投掷飞镖。假设这些飞镖在整个板上均匀分布。现在移除所有落在曲线下方以外区域的飞镖。剩下的飞镖将在曲线下方的区域内均匀分布，并且这些飞镖的 x 坐标将按照随机变量的密度分布。这是因为在曲线最高的地方，也就是概率密度最大的地方，飞镖着陆的空间最多。

拒绝抽样的一般形式假设板子的形状不一定是矩形，而是根据某个提议分布的密度来确定（该分布不一定归一化为 1）。通常情况下将其视为某个已知的分布的倍数。提议分布中的每个点至少与想要抽样的分布一样高，以便前者完全包围后者。（否则，想要抽样的曲线区域中的某些部分可能永远无法到达。）

（二）工作原理

拒绝抽样的工作原理：

• 从提议分布中在 x 轴上抽样一个点。
• 在该 x 位置上画一条竖直线，直到提议分布的概率密度函数的 y值。
• 在这条线上从 0 到提议分布的概率密度函数的 y 值之间均匀抽样。如果抽样值大于该竖直线上所需分布的密度函数值，则拒绝该 x 值并返回第 1 步；否则，该 x 值就是所需分布的一个样本。

拒绝抽样算法可以用于从任何曲线下方进行抽样，无论函数是否积分为 1。事实上，通过常数缩放函数对抽样的 x 位置没有影响。因此，该算法可以用于从归一化常数未知的分布中进行抽样。

（三）采样步骤

提案分布$g$：为了从密度为$f$的分布$X$中获取样本，利用了容易采样的密度函数为$g$的分布$Y$，$g$就是提案分布

设$M$为似然比$f(x)/g(x)$的上界，即一个常数满足$1\le M<\infty$。也就是说$M$必须满足$f(x)\le Mg(x)$对任意$x$都成立，因此$Y$分布的支撑要包含$X$的支撑

• 从分布$Y$获取样本$y\sim g$，并从$Unif(0,1)$（单位区间上的均匀分布）获取样本$u$
• 检查是否$u<f(y)/Mg(y)$.($M\ge 1$)
  • 成立则接受$y$作为从$f$中抽取的样本
  • 不成立则拒绝$y$的值并重新获取样本

保留样本不大于值 y 的概率为：

其中$P[U\le \frac{f(Y)}{Mg(Y)}]=\frac{1}{M}$为接受率，接受率越大，采样效率就越高。接受拒绝算法平均需要 M 次迭代才能获得样本，并且M 越小越好，即包络线越贴近目标分布越好，可设
$$M=\underset{x}{\max }\frac{f(x)}{g(x)}$$

（四）示例

试用接受拒绝采样产生服从均值为 0，方差为 1 的半正态分布的随机数 $\eta$，该分布的概率密度函数为
$$p(z)=\{\begin{array}{ll}\sqrt{2/\pi }{e}^{-\frac{z^2}{2}},& \text{z≥0}\\ 0,& \text{z<0}\end{array}$$

五、重要性采样

（一）基本概念

接受拒绝采样完美的解决了累积分布函数不可求时的采样问题。但是接受拒绝采样非常依赖于提议分布的选择，如果提议分布选择的不好，可能采样时间很长却获得很少满足分布的粒子。而重要性采样就解决了这一问题。

重要性采样是使用蒙特卡洛方法估算积分 (期望) 时，提高对积分计算重要区域的抽样，从而达到减少方差的目的。
$$\mu =E_{X\sim f}(h(X))=\int h(x)f(x)dx=\int h(x)\frac{f(x)}{g(x)}g(x)dx$$

或若 $\int f(x)\ne 1$, 也就是只知道分布成比例于某个函数，差一个归一化常数，则
$$\mu =E_{X\propto f}(h(X))=\int h(x)\frac{f(x)}{\int f(x)dx}dx=\frac{\int h(x)\frac{f(x)}{g(x)}g(x)dx}{\int \frac{f(x)}{g(x)}g(x)dx}$$

（二）蒙特卡洛估计

上式建议用来估计 $Eh(X)$ 的一种 Monte Carlo 方法：

• 从$g$中抽取独立同分布的样本$X_1,...,X_n$，并采用估计：
  $${\hat{\mu }}_{IS}^{\ast }=\frac{1}{n}\sum _{i}h(x_i)\frac{f(x_i)}{g(x_i)}\to E_{X\sim g}(h(x)\frac{f(x)}{g(x)})$$

• 可以写成：（$w^{\ast }(X_i)=f(X_i)/g(X_i)$是未标准化权重，称为重要性比率）
  $${\hat{\mu }}_{IS}^{\ast }=\frac{1}{n}\sum _{i}h(X_i)w^{\ast }(X_i)$$

• 若是差一个比例常数的 f，则（$w^{\ast }(X_i)=w^{\ast }(X_i)/{\sum }_{i=1}^n{w}^{\ast }(X_i)$是标准化权重）
  $${\hat{\mu }}_{IS}=\frac{1}{n}\sum _{i}h(X_i)w(X_i)$$

• 估计值的方差是
  V a r ( μ ^ ) = 1 n V a r h ( X ) f ( X ) g ( X ) = 1 n ∫ ( h ( x ) f ( x ) g ( x ) − μ 2 ) 2 g ( x ) d x = 1 n { ∫ ( h 2 ( x ) f 2 ( x ) g ( x ) d x − μ 2 } \begin{aligned} Var(\hat{\mu})=&\frac{1}{n}Var\frac{h(X)f(X)}{g(X)} \\ =&\frac{1}{n}\int (\frac{h(x)f(x)}{g(x)}-\mu^2)^2g(x)dx\\ =&\frac{1}{n}\{ \int(\frac{h^2(x)f^2(x)}{g(x)}dx-\mu^2\}\\ \end{aligned} Var(μ^​)===​n1​Varg(X)h(X)f(X)​n1​∫(g(x)h(x)f(x)​−μ2)2g(x)dxn1​{∫(g(x)h2(x)f2(x)​dx−μ2}​
  如果$h^2(x)f^2(x)/g^2(x)={\mu }^2$，就有$Var(\hat{\mu })=0$，即方差达到最小。此时：
  $$g(x)=\frac{|h(x)|f(x)}{\mu }=\frac{|h(x)|f(x)}{\int |h(x)|f(x)dx}$$
  密度函数 $g(x)$ 的最佳选择就是和被积函数 $|h(x)|f(x)$ 具有相同的形状。对积分值贡献越大的区域，希望以较大的概率抽取到随机数。
  在实际中， $\mu$ 是未知量，因此无法选取 $g(x)$，使得 $Var(\hat{\mu })=0$。通常情况下，我们会选取一个形状接近 $|h(x)|f(x)$ 的函数作为 $g(x)$。
  $$\mu =\int h(x)f(x)dx=\int h(x)\frac{f(x)}{g(x)}g(x)dx$$
  

（三）示例

例：用重要抽样法估计$\mu ={\int }_0^{1}xf(x)dx={\int }_0^1{e}^{x}dx$的估计值

参考：
MCMC入门（一）蒙特卡罗方法与拒绝-接受采样