蒙特卡洛法求积分

问题一：我们如何用蒙特卡洛方法求积分？问题二：如何近似求一个随机变量的数学期望？问题三：估计的误差是多少？问题四：如何从理论上对蒙特卡洛估计做分析？结论

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('whitegrid')

问题一：我们如何用蒙特卡洛方法求积分？

你眼中的蒙特卡洛方法求积分，可能是这样子的：

Image Name 最最经典的例子就是求 的近似值了，生成若干个均匀的点，然后统计在圆内的点的个数的比例，这个比例就是 的近似了！ 但是这种方法的计算量非常大，而且随着维数的增长需要的计算量也会成倍上升，收敛速度也并不快。

但是我要讲的蒙特卡洛求积跟这个些许不一样。

假如我想求 ，我们来用概率的语言表达一下它。 设随机变量 ，即 上的均匀分布， 具有密度函数 。 那么就有： ，这个公式是下面推导中非常重要的一环。 事实上，借助这个公式，我们将求积分转化为求某个随机变量的数学期望！

接下来我们就可以用统计学的手段来处理了。

问题二：如何近似求一个随机变量的数学期望？

通常情况下， 都是一个比较复杂的函数。我们想要近似求期望，只能用统计学的手段。

设随机变量 ，一个常用的办法是，如果我们找到 个随机变量 的样本 那么 就是 一个好的近似！

容易知道，上式中的 服从 上的均匀分布。

所以我们的做法可以总结如下：

1. 生成 个 上均匀分布的随机数 。
2. 计算函数值
3. 积分 具有近似值

接下来看一个例子来验证我们的结果：

, 用蒙特卡洛方法近似计算 ，我们知道真实值是 ，介于2.6到2.7之间 取

N = 1000
x = 2 * np.random.random(size=N)  # [0, 2] 上的均匀分布
ans = 2 / N * np.sum(x ** 2)
print(ans)

2.6431046604313226

目前看来效果还算不错，但是问题就这样结束了么？

问题三：估计的误差是多少？

凡估计必有误差

每一次采样都可以得到一个估计值，我们多次采样，得到多个估计值，画出多个估计值的分布图，从图上就可以近似看出估计的误差了。

record = []  # 记录多次采样的估计值
N = 1000  # 每次采样取1000个点
for _ in range(1000):
    x = 2 * np.random.random(size=N)
    ans = 2 / N * np.sum(x ** 2)
    record.append(ans)
plt.title('N=1000')
sns.distplot(record); plt.show()

惊讶的观察到：

• 虽然我们真实值是2.67左右，但是估计出2.5、2.9这种偏离程度大的值还是有可能的。
• 多次采样的结果分布像是正态分布。（这是巧合吗？）

把单次采样点增加到2000个，直觉告诉我们，误差会减小！

record = []  # 记录多次采样的估计值
N = 2000  # 每次采样取2000个点
for _ in range(1000):
    x = 2 * np.random.random(size=N)
    ans = 2 / N * np.sum(x ** 2)
    record.append(ans)
plt.title('N=2000')
sns.distplot(record); plt.show()

明显看到分布更加紧凑了一些。

问题四：如何从理论上对蒙特卡洛估计做分析？

在统计学上，我们评价一个估计好不好的标准有哪些呢？ 统计量有三大基本性质：

• 无偏性、有效性、相合性（一致性）

无偏性表示这个估计有没有 bias；有效性指这个估计的方差够不够小；相合性或者说一致性，说的是当样本容量非常大的时候，估计值是否趋近于真实值。

接下来我们从理论上来讨论这三点： 设 是 的一个估计量，则：

1. 是无偏的，这点很好理解，对于样本 , 。
2. 是相合估计量，根据大数定律，估计量 几乎绝对收敛与 。
3. 有效性没法说明，但是我们可以知道 的方差为 （样本之间独立） 通过中心极限定理，还可以知道，当样本容量很大的时候， 渐进服从以 为均值， 为方差的正态分布。

最后，我想展示一下，本文所述的转化为估计随机变量期望的蒙特卡洛方法 与 传统的往正方形内投点计算落在圆内的点个数来估计 值的方法的不同。

Image Name

plt.figure(figsize=(12, 4))
record = []
for _ in range(1000):
    x = np.random.random(size=(2, 200))
    ans = 4 * np.mean(np.linalg.norm(x, axis=0) < 1)  # 是否落在圆内
    record.append(ans)
plt.subplot(121)
sns.distplot(record, kde=False)

record = []
for _ in range(1000):
    x = np.random.random(size=200)
    ans = 4 * np.mean(np.sqrt(1 - x ** 2))  # y = sqrt(1 - x ^ 2)
    record.append(ans)
plt.subplot(122)
sns.distplot(record, kde=False); plt.show()

同样是取了2000个点（做200次计算），统计1000次结果。左图为传统方法，右图为本文所述转化为求期望的方法。

明显右边的效果更好！

结论

本文简要介绍了蒙特卡洛方法求积分的思路，以及相应的理论推导。蒙特卡洛求积分的本质是利用随机模拟估计一个随机变量的期望。理解好蒙特卡洛求积的思想有助于进一步学习MCMC方法。

进一步还可以思考：

• 如何用蒙特卡洛估计重积分？这种方法会随着维数的增大而出现计算困难吗？（即：是多项式级别的更困难还是指数级别的更困难）
• 说到用随机模拟估计一个随机变量的期望，如果我们能够获取到一定数量的某随机变量的样本，但是这个随机变量的分布非常畸形（考虑中国10亿劳动人口有6亿月收入不到1千，福布斯中国排行榜上前100的富豪能对中国人均收入产生巨大的影响），如何改进我们模拟的方法？