蒙特卡洛（Monte Carlo）方法

蒙特卡洛方法可以近似计算某个概率值，计算结果随着实验次数增加而愈加精确，本文记录相关内容。

简介

1. 蒙特卡洛方法Monte Carlo 可以通过采用随机投点法来求解不规则图形的面积。 求解结果并不是一个精确值，而是一个近似值。当投点的数量越来越大时，该近似值也越接近真实值。
2. 蒙特卡洛方法也可以用于根据概率分布来随机采样的任务。

布丰投针

布丰投针问题是1777年法国科学家布丰提出的一种计算圆周率的方法：随机投针法。

执行步骤

• 首先取一张白纸，在上面绘制许多条间距为d 的平行线。
• 取一根长度为l , l \lt d的针，随机地向纸上投掷n次，观测针与直线相交的次数，记做 m。
• 计算针与直线相交的概率 p=\frac{m}{n}。因此有：

相交概率证明

由于向纸上投针是完全随机的, 因此用二维随机变量 (X, Y) 来确定针在纸上的具体位置。其中:

• X 表示针的中点到平行线的距离，它是 [0, d / 2] 之间的均匀分布。
• Y 表示针与平行线的夹角， 它是 \left[0, \frac{\pi}{2}\right] 之间的均匀分布。

• 当 X<\frac{l}{2} \sin Y
• 由于 X, Y 相互独立， 因此有概率密度函数：

• 因此，针与直线相交的概率为：

• 根据 \frac{2 l}{\pi d}=\frac{m}{n} 即可得证。

蒙特卡洛积分

对于函数 f(x) , 其在区间 [a, b] 上的积分 \int_{a}^{b} f(x) d x 可以采用两种方法来求解: 投点法、期望法。

投点法

• 对函数 f(x) , 对其求积分等价于求它的曲线下方的面积。
• 定义一个常数 M，使得 M>\max _{a \leq x \leq b} f(x) [a, b] 上的面积就是矩形面积 M(b-a) .
• 随机向矩形框中随机的、均匀的投点，设落在函数 f(x) 下方的点为绿色，落在 f(x) 和M之间的点为红色。
• 则有：落在 f(x) 下方的点的概率等于 f(x) 的面积比上矩形框的面积 。
• 从 [a, b] 之间的均匀分布中采样 x_{0} , 从 [0, M] 之见的均匀分布中采样 y_{0}, \quad\left(x_{0}, y_{0}\right) 构成一个 随机点。
• 若 y_{0} \leq f\left(x_{0}\right) , 则说明该随机点在函数 f(x) 下方，染成绿色。
• 若 f\left(x_{0}\right)<y_{0} \leq M f(x) 上方，染成红色。
• 假设绿色点有n_1个，红色点有n_2个，总的点数为n_1 + n_2 ，因此有：

期望法

• 假设需要求解积分 I=\int_{a}^{b} f(x) d x ，则任意选择一个概率密度函数 p(x) ，其中 p(x) 满足条件:

• 令:

• 则有:

• I刚好是一个期望：设随机变量 X 服从分布 X \sim p(x) , 则 I=\mathbb{E}_{X \sim p}\left[f^{*}(X)\right]
• 则期望法求积分的步骤是： 任选一个满足条件的概率分布 p(x) 。 根据 p(x) , 生成一组服从分布 p(x) 的随机数 x_{1}, x_{2}, \cdots, x_{N} 。 计算均值 \bar{I}=\frac{1}{N} \sum_{i=1}^{N} f^{*}\left(x_{i}\right) , 并将 \bar{I} 作为 I 的近似。
• 在期望法求积分中, 如果 a, b 均为有限值, 则 p(x) 可以取均匀分布的概率密度函数：

	此时 $ f^{*}(x)=(b-a) f(x), \quad \bar{I}=\frac{b-a}{N} \sum_{i=1}^{N} f\left(x_{i}\right) $ 。

• 其物理意义为: \frac{\sum_{i}^{N} f\left(x_{i}\right)}{N} 为在区间 [a, b] 上函数的平均高度, 乘以区间宽度 b-a 就是平均面积。

蒙特卡洛采样

采样问题的主要任务是：根据概率分布 p(x) , 生成一组服从分布 p(x) 的随机数 x_{1}, x_{2}, \cdots .

均匀分布模拟$p(x)$采样

• 如果 p(x) 就是均匀分布，则均匀分布的采样非常简单。
• 如果 p(x) 是非均匀分布，则可以通过均匀分布的采样来实现。

1. 首先根据均匀分布 U(0,1) 随机生成一个样本 z_{i} 。
2. 设 \tilde{P}(x) 为概率分布 p(x) 的累计分布函数:

1. 令 z_{i}=\tilde{P}\left(x_{i}\right) , 计算得到 x_{i}=\tilde{P}^{-1}\left(z_{i}\right) , 其中 \tilde{P}^{-1} 为反函数, 则 x_{i} 为对 p(x) 的采样。

1. 通过均匀分布的采样的原理：假设随机变量 Z, X 满足 Z=\tilde{P}(X) , 则 X 的概率分布为:

• 因为 Z 是 [0,1] 上面的均匀分布，因此 p_{Z}(z)=1 ; \tilde{P}(x) 为概率分布 p(x) 的累计分布函数, 因此 \frac{d}{d x} \tilde{P}(x)=p_{X}(x) 。因此上式刚好等于 p(x) , 即： x_{i} 服从概率分布 p(x) 。

1. 这其中有两个关键计算：

• 根据 p(x) , 计算累计分布函数 \tilde{P}(x)=\int_{-\infty}^{\pi} p(z) d z 。
• 根据 z=\tilde{P}(x) 计算反函数 x=\tilde{P}^{-1}(z) 。

如果累计分布函数无法计算，或者反函数难以求解，则该方法无法进行。

接受-拒绝采样

对于复杂的概率分布p(x) ，难以通过均匀分布来实现采样。此时可以使用接受-拒绝采样 策略。

• 首先选定一个容易采样的概率分布 q(x) ，选择一个常数 k , 使得在定义域的所有位置都满足 p(x) \leq k \times q(x) 。
• 然后根据概率分布 q(x) 随机生成一个样本 x_{i} 。
• 计算 \alpha_{i}=\frac{p\left(x_{i}\right)}{k q\left(x_{i}\right)} , 以概率 \alpha_{i} 接受该样本。

方法一

根据均匀分布 U(0,1) 随机生成一个点 u_{i} ， 如果 u_{i} \leq \alpha_{i} ，则接受该样本；否则拒绝该样本。

方法二

根据均匀分布 U\left(0, k q\left(x_{i}\right)\right) 生成一个随机点, 如果该点落在灰色区间((p(x_{i}), k q(x_{i})]) 则拒绝该样本；如果该点落在白色区间 \left(\left[0, p\left(x_{i}\right)\right]\right) 则接受该样本。

不足

接受-拒绝采样 在高维的情况下会出现两个问题：

• 合适的q 分布比较难以找到。
• 难以确定一个合理的k值。

这两个问题会导致拒绝率很高，无效计算太多。

参考资料

• http://www.huaxiaozhuan.com/数学基础/chapters/4_monte_carlo.html