问题一:我们如何用蒙特卡洛方法求积分?问题二:如何近似求一个随机变量的数学期望?问题三:估计的误差是多少?问题四:如何从理论上对蒙特卡洛估计做分析?结论
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('whitegrid')
问题一:我们如何用蒙特卡洛方法求积分?
你眼中的蒙特卡洛方法求积分,可能是这样子的:
![Image Name](https://i-blog.csdnimg.cn/blog_migrate/0550e8db43f34e1187c8ee45b0559861.png)
最最经典的例子就是求 的近似值了,生成若干个均匀的点,然后统计在圆内的点的个数的比例,这个比例就是 的近似了!
但是这种方法的计算量非常大,而且随着维数的增长需要的计算量也会成倍上升,收敛速度也并不快。
但是我要讲的蒙特卡洛求积跟这个些许不一样。
假如我想求
,我们来用概率的语言表达一下它。
设随机变量
,即
上的均匀分布,
具有密度函数
。
那么就有:
,这个公式是下面推导中非常重要的一环。
事实上,借助这个公式,我们将求积分转化为求某个随机变量的数学期望!
接下来我们就可以用统计学的手段来处理了。
问题二:如何近似求一个随机变量的数学期望?
通常情况下, 都是一个比较复杂的函数。我们想要近似求期望,只能用统计学的手段。
设随机变量
,一个常用的办法是,如果我们找到
个随机变量
的样本
那么
就是
一个好的近似!
容易知道,上式中的 服从 上的均匀分布。
所以我们的做法可以总结如下:
生成 个 上均匀分布的随机数 。
计算函数值
积分 具有近似值
接下来看一个例子来验证我们的结果:
, 用蒙特卡洛方法近似计算
,我们知道真实值是
,介于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;有效性指这个估计的方差够不够小;相合性或者说一致性,说的是当样本容量非常大的时候,估计值是否趋近于真实值。
接下来我们从理论上来讨论这三点:
设
是
的一个估计量,则:
是无偏的,这点很好理解,对于样本 , 。
是相合估计量,根据大数定律,估计量 几乎绝对收敛与 。
有效性没法说明,但是我们可以知道 的方差为 (样本之间独立)
通过中心极限定理,还可以知道,当样本容量很大的时候, 渐进服从以 为均值, 为方差的正态分布。
最后,我想展示一下,本文所述的转化为估计随机变量期望的蒙特卡洛方法 与 传统的往正方形内投点计算落在圆内的点个数来估计 值的方法的不同。
![Image Name](https://i-blog.csdnimg.cn/blog_migrate/66914809d61e7fde4cb739ad6fdcf0d7.png)
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的富豪能对中国人均收入产生巨大的影响),如何改进我们模拟的方法?