References：

马尔可夫链蒙特卡罗算法（MCMC）

走进贝叶斯统计（四）—— 蒙特卡洛方法

蒙特卡罗方法详解

【概述】

虽然共轭先验能够极大地简化贝叶斯统计推断过程中的计算难度，但它不一定能够很好地描绘对目标参数 $\theta$ 的经验判断

然而，能够准确描绘对目标参数 $\theta$ 的经验判断的非共轭先验分布 $p(\theta)$ 又往往容易产生极其复杂、难以计算性质的后验分布

这时候，蒙特卡洛方法（Monte Carlo Method）的重要性便得以凸显，其是一种从概率模型的随机抽样中进行近似数值计算的方法，它帮助贝叶斯统计摆脱了计算上的束缚，能够随心所欲地选取先验分布

【蒙特卡罗法】

假设有一个半径 $r=1$ 的圆和一个边长为 $1$ 的正方形，圆的面积为 $\pi r^2=\pi$，那么正方形内部的相切圆的面积为整个圆的 $\frac{1}{4}$，即 $\frac{\pi}{4}$，正方形的面积为 $1$

下面不断的向正方形中随机打点，这些点有一定的概率落在相切圆中，一定的概率落在相切圆外，显然落在相切圆中的概率是：

$P(落在相切圆中的概率) = \frac{S_{圆}}{S_{正方}} = \frac{\pi}{4}$

于是，即可得到圆周率的计算公式：

$\pi = 4 * P(落在相切圆中的概率) = 4 * \frac{红色点数}{总点数}$

如下图所示，随着随机点数的增加，$\pi$ 的近似度也在不断增大

上述的这个例子实际上道出了蒙特卡洛法最基本的思想：从目标分布中采样大量的样本，并用这些样本不断去逼近这个分布

也就是说，蒙特卡罗法的应用场景是：假设概率分布已知，那么通过抽样获得概率分布的随机样本，然后通过得到的随机样本对概率分布的特征进行分析

这种做法的理论支撑来源于格里文科定理与大数定律：

格里文科定理：当样本量趋于无穷时，经验分布函数收敛于真实的总体分布函数
大数定律：样本均值在样本量趋于无穷时依分布收敛于总体期望

因此，在这两个定理的支撑下，当采样的样本越多时，对分布的模拟就越加准确

如下图所示，当采样的样本量越大，便越能够逼近目标的分布，因此就可以通过样本来得到目标分布的期望、方差等相关信息

但关键问题是：假设已知目标分布 $X\sim p(x)$ 的具体形式，那么该如何对目标分布进行采样？

对于那些通过计算得到的形式较为罕见的后验分布 $p(x)$，一个直观而可行的方法是利用该分布的累积分布函数 $F(X)$ 与均匀分布 $U(0,1)$ 的采样函数进行采样

对于随机变量 $X$，其分布函数为 $F(X)$，那么对于任意 $u\in [0,1]$，有：

$\begin{align*} P(U\leq u) &= P(F(X)\leq u) \\ &= P(X\leq F^{-1}(u)) \\ &= F(F^{-1}(u)) \\ &= u \end{align*}$

即，$F(X)\sim U(0,1)$

也就是说，通过该分布的累积分布函数的反函数 $F^{-1}$，与从均匀分布 $U(0,1)$ 中采样得到的样本 $u$，可以得到目标分布 $p(x)$ 的样本，即：

$x = F^{-1}(u)$

但在大多数情况下，连目标分布的累积分布函数 $F(X)$ 的形式都难以写出，就更别说求出累积分布函数 $F(X)$ 以及它的反函数 $F^{-1}$ 了

因此，找到一个能够对任意分布都有效的采样方法非常重要，所以蒙特卡罗法的核心是随机抽样（Random Sampling）

一般的蒙特卡罗法有直接采样法、接受-拒绝采样法、重要性采样法等，这里仅介绍基于接受-拒绝采样法的蒙特卡罗法

【接受-拒绝采样法】

基本思路

接受-拒绝采样法（Accept-reject Sampling Method）是一种随机抽样方法，其能够对任意分布进行有效的采样

若已知目标分布 $X\sim p(x)$ 的具体形式，现在想要得到该概率分布的随机样本，以对这个概率分布进行分析

假设 $p(x)$ 无法直接进行抽样，那么需要找一个可以直接抽样的分布，其被称为提议分布（Proposal Distribution）

设 $q(x)$ 为提议分布的概率密度函数，且满足 $c\cdot q(x)\geq p(x),c>0$，即：

$\frac{p(x)}{q(x)}\leq c,\quad c>0$

即常数 $c$ 一定大于等于 $\frac{p(x)}{q(x)}$ 的上界

之后，对 $q(x)$ 进行抽样，假设得到的结果是 $x^{*}$，按照 $\frac{p(x^{*})}{c\cdot q(x^{*})}$ 的比例，随机决定是否接受 $x^{*}$

如图所示，直观上落在 $p(x)$ 的范围内就接受，落在 $p(x)$ 的范围外就拒绝，也就是说，接受-拒绝采样法实际上是按照 $p(x)$ 的涵盖面积占 $c\cdot q(x)$ 的涵盖面积的比例进行抽样

在接受-拒绝采样法中，需要注意的有以下几点：

理论上可能出现 $c=\infty$ 的情况，但在常规的应用场景下，这种情况很难发生
对于提议分布 $q(x)$ 的选取可以是任意的，一般选择最简单的均匀分布即可
提议分布 $q(x)$ 的选取能够影响采样的效率，若 $p(x)$ 的涵盖体积占 $c\cdot q(x)$ 的涵盖体积的比例很低，会导致拒绝的比例很高，采样效率很低
在高维空间进行抽样时，即使目标分布 $p(x)$ 和提议分布 $q(x)$ 很接近，但两者涵盖的高维体积差异可能会很大

算法流程

基于接受-拒绝采样法的蒙特卡罗法算法流程如下

输入：抽样的目标分布的概率密度函数 $p(x)$

输出：满足目标分布 $p(x)$ 的随机样本 $x_1,x_2,\cdots,x_n$

参数：样本数 $n$

算法流程：

Step 1：选择概率密度为 $q(x)$ 的概率分布作为提议分布，使其对任一 $x$ 满足

$\frac{p(x)}{q(x)}\leq c,\quad c>0$

Step 2：在提议分布 $q(x)$ 中随机抽样，得到样本 $x^{*}$

Step 3：从分布 $U(0,1)$ 中采样一个随机数 $u$

Step 4：若有 $u\leq \frac{p(x^{*})}{c\cdot q(x^{*})}$，则将 $x^{*}$ 视为从目标分布 $p(x)$ 中采样得到的样本，否则视为采样失败，回到 Step 2，重新开始采样

Step 5：当得到 $n$ 个随机样本时，算法结束

证明

从基于接受-拒绝采样法的蒙特卡罗法的算法流程来看，存在一个极为重要的点，即：为什么从任意选取的提议分布 $q(x)$ 中采样得到的 $x^{*}$，一旦满足 $u\leq \frac{f(x^{*})}{c \cdot g(x^{*})}$ 时，即将 $x^{*}$ 视为从目标分布 $p(x)$ 中采取的样本？

下面开始证明从上述的接受-拒绝采样法生成的样本 $x^{*}$ 服从目标分布 $p(x)$

若要证明从上述采样方法生成的样本 $x^{*}$ 服从目标分布 $p(x)$，就等价于在已知 $x^{*}\sim q(x)$ 的情况下，证明：

$F(x^{*}) = P(X\leq x^{*}| u\leq \frac{p(x^{*})}{c\cdot q(x^{*})})$

其中，$F(X)$ 为 $p(x)$ 的累积分布函数，$u$ 是从均匀分布 $U(0,1)$ 中采样一个随机数

根据贝叶斯公式 $P(A|B)=\frac{P(B|A)P(A)}{P(B)}$，可以将 $P(X\leq x^{*}| u\leq \frac{p(x^{*})}{c\cdot q(x^{*})})$ 转化为：

$P(X\leq x^{*}| u\leq \frac{p(x^{*})}{c\cdot q(x^{*})}) = \frac{P( u\leq \frac{p(x^{*})}{c\cdot q(x^{*})} | X\leq x^{*}) P(X\leq x^{*})}{P(u\leq \frac{p(x^{*})}{c\cdot q(x^{*})})}$

下面，只需分析等式右边的三个函数，以证明其组合等于 $F(x^{*})$ 即可

1）对于 $P(u\leq \frac{p(x^{*})}{c\cdot q(x^{*})})$，根据全概率公式，有：

$P(u\leq \frac{p(x^{*})}{c\cdot q(x^{*})}) = \int P(u\leq \frac{p(x^{*})}{c\cdot q(x^{*})}|X=x^{*})P(X=x^{*}) dx$

由于 $x^{*}$ 是从提议分布 $q(x^{*})$ 中抽样得到的，故 $P(X=x^{*})=q(x^{*})$

同时，在抽样得到一个具体数值 $X=x^{*}$ 后，采样成功的条件概率即为：

$P(u\leq \frac{p(x^{*})}{c\cdot q(x^{*})}|X=x^{*}) = \frac{p(x^{*})}{c\cdot q(x^{*})}$

故而可以得到：

$\begin{align*} P(u\leq \frac{p(x^{*})}{c\cdot q(x^{*})}) &= \int P(u\leq \frac{p(x^{*})}{c\cdot q(x^{*})}|X=x^{*})P(X=x^{*})dx \\ &= \int \frac{p(x^{*})}{c\cdot q(x^{*})} q(x^{*}) dx \\ &= \int \frac{p(x^{*})}{c} dx \\ &= \frac{1}{c} \end{align*} \tag{1}$

2）对于 $P(X\leq x^{*})$，已知 $x^{*} \sim q(x^{*})$，故有：

$P(X\leq x^{*}) = G(x^{*}) \tag{2}$

其中，$G(Y)$ 为 $q(x)$ 的累积分布函数

3）对于 $P( u\leq \frac{p(x^{*})}{c\cdot q(x^{*})} | X\leq x^{*}) $，根据条件概率公式，有：

$\begin{align*} P( u\leq \frac{p(x^{*})}{c\cdot q(x^{*})} | X\leq x^{*}) &= \frac{P(u\leq \frac{p(x^{*})}{c\cdot q(x^{*})},X\leq x^{*})}{P(X\leq x^{*})} \\ &= \frac{P(u\leq \frac{p(x^{*})}{c\cdot q(x^{*})},X\leq x^{*})}{G(x^{*})} \\ &= \int_{-\infty}^{x^{*}} \frac{P(u\leq \frac{p(x^{*})}{c\cdot q(x^{*})}|X=w\leq x^{*})}{G(x^{*})} g(w) dw \\ &= \frac{1}{G(x^{*})} \int_{-\infty}^{x^{*}} \frac{p(w)}{c\cdot q(w)} g(w) dw \\ &= \frac{1}{c\cdot G(x^{*})}F(x^{*}) \end{align*} \tag{3}$

结合式 $(1),(2),(3)$，即可得到：

$\begin{align*} P(X\leq x^{*}| u\leq \frac{p(x^{*})}{c\cdot q(x^{*})}) &= \frac{P( u\leq \frac{q(x^{*})}{c\cdot q(x^{*})} | X\leq x^{*}) P(X\leq x^{*})}{P(u\leq \frac{p(y_i)}{c\cdot q(x^{*})})} \\ &= \frac{G(x^{*}) \frac{1}{c\cdot G(x^{*})} F(x^{*})}{\frac{1}{c} } \\ &= F(x^{*}) \end{align*}$

即：

$P(X\leq x^{*}| u\leq \frac{p(x^{*})}{c\cdot q(x^{*})}) = F(x^{*})$

即通过拒绝采样法生成的样本 $x^{*}$，服从目标分布 $p(x)$

【应用】

数学期望估计

通过蒙特卡罗法采样得到的样本，可用于进行数学期望估计（Estimation of Mathematical Expectation）

设随机变量 $X\sim p(x)$，$f(x)$ 是定义在 $X\in \mathcal{X}$ 上的函数，现要求函数 $f(x)$ 关于概率密度函数 $p(x)$ 的数学期望 $E_{p(x)}[f(x)]$

针对该问题，使用蒙特卡罗法独立地抽取服从分布 $p(x)$ 的 $n$ 个样本 $x_1,x_2,\cdots,x_n$，然后计算函数 $f(x)$ 的样本均值

$\hat{f_n} = \frac{1}{n}\sum_{i=1}^n f(x_i)$

作为数学期望 $E_{p(x)}[f(x)]$ 的近似值

根据大数定理可知，当样本容量增大时，样本均值以概率 $1$ 收敛于数学期望，即：

$\hat{f_n} \rightarrow E_{p(x)}[f(x)],\quad n\rightarrow \infty$

由此，即可得到数学期望的近似计算方法

$E_{p(x)}[f(x)] \approx \frac{1}{n}\sum_{i=1}^n f(x_i)$

定积分计算

蒙特卡罗法也可用于定积分的近似计算，称为蒙特卡罗积分（Monte Carlo Integration）

假设有一个函数 $h(x)$，现要计算该函数的积分

$\int_{\mathcal{X}} h(x) dx$

若能将 $h(x)$ 分解为一个函数 $f(x)$ 和一个概率密度函数 $p(x)$ 的乘积，那么有：

$\int_{\mathcal{X}} h(x) dx = \int_{\mathcal{X}} f(x)p(x)dx =E_{p(x)}[f(x)]$

即函数 $h(x)$ 的积分可表示为函数 $f(x)$ 关于概率密度函数 $p(x)$ 的数学期望

实际上，给定一个概率密度函数 $p(x)$，只需取 $f(x)=\frac{h(x)}{p(x)}$ 即可得到上述公式，也就是说，任何一个函数的积分都可以表示为某个函数的数学期望的形式

而函数的数学期望又可以通过函数的样本均值估计，于是，就可以利用样本均值来近似计算积分

$\int_{\mathcal{X}} h(x) dx =E_{p(x)}[f(x)] \approx \frac{1}{n}\sum_{i=1}^n f(x_i)$

例如：用蒙特卡罗积分求 $\int_{0}^1 e^{-\frac{x^2}{2}}dx$

假设随机变量 $X$ 在 $(0,1)$ 区间内遵循均匀分布，则有：

$\begin{align*} f(x) &= e^{-\frac{x^2}{2}} \\ p(x) &= 1(0<x<1) \end{align*}$

使用蒙特卡罗法，在 $(0,1)$ 区间内按均匀分布抽取 $10$ 个样本 $x_1,x_2,\cdots,x_{10}$，然后计算样本的函数均值

$\hat{f_{10}} = \frac{1}{10}\sum_{i=1}^{10} e^{-\frac{x_i^2}{x}} =0.832$

作为积分的近似，随机样本数越大，计算就越精确