【Practical】蒙特卡罗法及其应用

蒙特卡罗法.

• Monte Carlo Method，也称之为统计模拟方法Statistical Simulation Method，是通过从概率模型的随机抽样进行近似数值计算的方法。根据所使用概率模型的不同可以细化出多种不同的蒙特卡罗法，最简洁也最为人熟知的 —— 蒙特卡罗法计算圆周率：在单位正方形内依据二维均匀分布对点进行随机抽样，当抽样次数足够多时，最终落在单位圆内的样本点数量与全部样本点数量就近似等于圆面积与正方形面积之比：$$a\approx \frac{\pi \cdot (1/2{)}^2}{1^2}=\frac{\pi }{4}$$于是可以反推出圆周率 $\pi$ 的近似表示。
• 蒙特卡罗法主要的应用有：随机抽样、数学期望估计以及定积分计算。

随机抽样.

• 随机抽样是蒙特卡洛法的核心，假定我们已知概率分布的定义，那么通过随机抽样获得随机样本后，可以通过随机样本的样本特征得到对于概率分布总体特征的估计，例如以样本均值估计总体期望。除了最为暴力简洁的直接抽样法外，当概率密度函数过于复杂 —— 多变量、变量之间不独立时，可以使用接受拒绝抽样法。
• 接受拒绝抽样法大多应用在原始概率分布 $p(x)$ 复杂，难以进行直接抽样时，我们选取另一个易于直接抽样的分布 $q(x)$，并且 $q(x)$ 满足：$$c\cdot q(x)\ge p(x),c>0$$这里的分布 $q(x)$ 被称为建议分布，按照 $q(x)$ 进行直接抽样的结果是 $x^{\ast }$，再对均匀分布 $U(0,1)$ 进行直接抽样得到 $u$，我们将 $u$ 与 $\frac{p(x^{\ast })}{c\cdot q(x^{\ast })}$ 进行比较，当且仅当 $u\le \frac{p(x^{\ast })}{c\cdot q(x^{\ast })}$ 时接受样本 $x^{\ast }$，否则拒绝它。
• 直观上来看，接受拒绝采样法在样本落到 $p(x)$ 范围内时就接受，而落到 $p(x)$ 范围外时就拒绝，其实质是按照 $p(x)$ 的涵盖面积占 $c\cdot q(x)$ 的涵盖面积之比来进行抽样的。
  
• 【证明】我们记按照 $q(x)$ 进行采用的样本被接受为事件 $A$，将原分布记为 $\pi (x)$ 以示区分，那么需要证明 $p(x|A)=\pi (x).$ 基于贝叶斯公式得到：$$p(x|A)=\frac{p(x,A)}{{\int }_{X}p(x,A)dx}$$其中联合概率分布 $p(x,A)$ 可以写为：$$p(A|x)\cdot p(x)$$并且其中的先验分布 $p(x)$ 就是我们用于直接采样的建议分布 $q(x)$，而样本 $x$ 是根据独立的均匀分布 $U(0,1)$ 来决定的，因此：$$p(A|x)={\int }_{U}p(A|x,u)p(u|x)du={\int }_{U}p(A|x,u)p(u)du$$根据接受拒绝采样法的接受原则我们知道，当且仅当 $u\le \frac{\pi (x)}{c\cdot q(x)}$ 时 $p(A|x,u)=1$，其它情况均为 $0$，因此：$$p(A|x)={\int }_0^{\frac{\pi (x)}{c\cdot q(x)}}p(u)du=\frac{\pi (x)}{c\cdot q(x)}$$上述积分是均匀分布 $U(0,1)$ 的积分。于是我们得到：$$p(x,A)=p(A|x)\cdot q(x)=\frac{\pi (x)}{c}$$因此：$$p(A|x)=\frac{\pi (x)/c}{{\int }_X(\pi (x)/c)dx}=\frac{\pi (x)}{{\int }_X\pi (x)dx}=\pi (x)$$该式说明我们基于接受拒绝采样法进行采样的分布 $p(A|x)$ 在概率意义上等于原始的采样概率分布 $\pi (x).$

数学期望估计.

• 【问题描述】假设有随机变量 $x$，取值 $x\in \mathcal{X}$，其概率密度函数为 $p(x)$，$f(x)$ 是定义在 $\mathcal{X}$ 上的函数，目标是求函数 $f(x)$ 关于密度函数 $p(x)$ 的数学期望 $E_{p(x)}[f(x)].$
• 处理该问题的常见做法就是，按照概率分布独立地随机抽取 $n$ 个样本 { x i ∣ i = 1 , 2 , ⋯   , n } \{x_i|i=1,2,\cdots,n\} {xi​∣i=1,2,⋯,n}，之后计算关于函数 $f(x)$ 的样本均值 $\hat{f}(n)$：$${\hat{f}}_n=\frac{1}{n}\sum _{i=1}^{n}f(x_i)$$将其作为数学期望 $E_{p(x)}[f(x)]$ 的近似值。依据弱大数定律，也称为辛钦大数定律，当样本容量趋于无穷时，样本均值的极限等于总体数学期望：$$\underset{n\to \infty }{\lim }{\hat{f}}_n=E_{p(x)}[f(x)]$$
• 上述过程是极为自然的处理策略，这种估计总体数学期望的技巧在下面使用蒙特卡罗方法进行数值定积分时也得到应用。

定积分.

• 数值积分方法可以对定积分进行近似计算，使用蒙特卡罗方法进行的定积分计算也被称作蒙特卡罗积分，对于函数 $h(x)$，计算它在区域 $\mathcal{X}$ 内的定积分：$${\int }_{\mathcal{X}}h(x)dx$$
• 如果能够将 $h(x)$ 进行分解：$$h(x)=f(x)\cdot p(x)$$其中 $p(x)$ 构成一个概率密度函数，那么上述定积分可以表示为：$${\int }_{\mathcal{X}}h(x)dx={\int }_{\mathcal{X}}f(x)\cdot p(x)dx=E_{p(x)}[f(x)]$$即我们将定积分问题转换成了 $f(x)$ 关于概率分布 $p(x)$ 的数学期望问题。实际上，我们选定一个概率分布函数 $p(x)$，而后取 $f(x)=\frac{h(x)}{p(x)}$ 即可满足分解要求。
• 参照上部分进行数学期望估计的方法，我们按照分布 $p(x)$ 进行随机采样，当 $n$ 足够大时，$f(x)$ 的样本均值能够很好的估计总体的数学期望，即：$${\int }_{\mathcal{X}}h(x)dx=E_{p(x)}[f(x)]\approx {\hat{f}}_n=\frac{1}{n}\sum _{i=1}^{n}f(x_i)$$

---

• 【例】计算定积分：$${\int }_0^1\exp (-x^2/2)dx$$
• 【解】在区间 $(0,1)$ 内取 $p(x)$ 为均匀分布，即：$$p(x)=1,0<x<1$$从而 $f(x)=\exp (-\frac{x^2}{2}).$
• 而后按照蒙特卡罗积分法，在 $(0,1)$ 区间内随机抽取 $n$ 个样本点，他们的样本均值即为我们要求的定积分结果近似值，并且辛钦大数定律指出随着 $n$ 趋向于无穷，这一近似程度总体上是越来越好的，不保证局部不会有震荡。
• 下面是对上述例子使用蒙特卡罗积分法的Python代码：

# -*- coding: utf-8 -*-
"""
Spyder Editor

This is a temporary script file.
"""

# In[]
import numpy as np
import matplotlib.pyplot as plt

# In[]
def func(x):
    return np.exp(-x**2/2)

# In[]
n = np.linspace(10,1000,1000).astype(int)

res = np.zeros(len(n))

# In[]
for i in range(len(n)):
    x = np.random.rand(n[i])
    x = func(x)
    
    res[i] = np.mean(x)
    
# In[]
plt.plot(n,res)

• 其运行结果如下，过程中存在震荡，但总体收敛到一个稳定值附近。
  

---