蒙特卡洛(Monte Carlo)方法

本文最后更新于:2022年8月5日 晚上

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

简介

  1. 蒙特卡洛方法Monte Carlo 可以通过采用随机投点法来求解不规则图形的面积。

    求解结果并不是一个精确值,而是一个近似值。当投点的数量越来越大时,该近似值也越接近真实值。

  2. 蒙特卡洛方法也可以用于根据概率分布来随机采样的任务。

布丰投针

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

执行步骤

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

$$
\pi = 2 \frac { n \times l } { m \times d }
$$

相交概率证明

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

  • $ X $ 表示针的中点到平行线的距离,它是 $ [0, d / 2] $ 之间的均匀分布。

  • $ Y $ 表示针与平行线的夹角, 它是 $ \left[0, \frac{\pi}{2}\right] $ 之间的均匀分布。

  • 当 $ X<\frac{l}{2} \sin Y $ 时,针与直线相交。
  • 由于 $ X, Y $ 相互独立, 因此有概率密度函数:
$$ p(X=x, Y=y)=\left\{\begin{array}{ll} \frac{4}{\pi d}, & 0 \leq x \leq d / 2,0 \leq y \leq \pi / 2 \\ 0, & \text { else } \end{array}\right. $$
  • 因此,针与直线相交的概率为:
$$ \begin{array}{c} P \left\{X<\frac{l}{2} \sin Y\right\} &= & \iint_{x<\frac{l}{2} \sin y} p(x, y) d x d y\\ &=& \int _ { x = 0 } ^ { x = \frac{1}{2} \sin y } \int_{y=0}^{ y = \pi / 2} \frac {4}{ \pi d} d x d y\\ &=& \frac{ 2l }{ \pi d } \end{array} $$
  • 根据 $ \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$ ,因此有:

$$
\int_{a}^{b} f(x) d x=\frac{n_{1}}{n_{1}+n_{2}} \times M(b-a)
$$

期望法

  • 假设需要求解积分$ I=\int_{a}^{b} f(x) d x $,则任意选择一个概率密度函数 $ p(x) $,其中$ p(x) $满足条件:
$$ \begin{array}{c} \int_{a}^{b} p(x) d x=1 \\ \text { if } f(x) \neq 0 \text { then } p(x) \neq 0, \quad a \leq x \leq b \end{array} $$
  • 令:
$$ f^{*}(x)=\left\{\begin{array}{ll} \frac{f(x)}{p(x)}, & p(x) \neq 0 \\ 0, & p(x)=0 \end{array}\right. $$
  • 则有:

$$
I=\int_{a}^{b} f(x) d x=\int_{a}^{b} f^{*}(x) p(x) d 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) $ 可以取均匀分布的概率密度函数:
$$ p(x)=\left\{\begin{array}{ll}\frac{1}{b-a}, & a \leq x \leq b \\ 0, & \text { else }\end{array}\right. $$
	此时 $ 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) $ 的累计分布函数:

$$
\tilde{P}(x)=\int_{-\infty}^{x} p(z) d z
$$

  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 $ 的概率分布为:

$$
p_{Z}(z) \frac{d}{d x} \tilde{P}(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$值。

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

参考资料


蒙特卡洛(Monte Carlo)方法
https://www.zywvvd.com/notes/study/math/monte-carlo/monte-carlo/
作者
Yiwei Zhang
发布于
2021年7月26日
许可协议