Monte Carlo 方法与 MCMC 简介

Monte Carlo 方法也即随机模拟方法的别称，它的基本原理是：当求解随机事件方式的概率或随机变量的数学期望时，通过设计某种实验，得出某个特定事件发生的频率，使用这个频率来近似表示这一事件发生的概率，从而得到问题的数值解。 可以看出，Monte Carlo 方法包含三个核心问题：构造概率过程、从已知概率分布中抽样、建立估计量。

大数定律

首先我们讨论下概率与频率的关系问题。独立地重复一个成功概率为 $p$ 的 Bernoulli 试验，用 $\xi_{n}$ 表示前 $n$ 次试验成功的次数，那么 $\frac{\xi_{n}}{n}$ 是前 $n$ 次试验成功的频率，也是一个随机变量，这个频率在某种意义下收敛于成功概率 $p$，即大数定律。

Bernoulli 定理: 对任何 $\epsilon>0$，有

$$\lim_{n}P\left( \left| \frac{\xi_{n}}{n}-p \right| > \epsilon \right) = 0$$
大数定律说明当样本很大时，不确定性就消失了。Monte Carlo 方法是大数定律的直接结果。

静态 Monte Carlo 方法

如果我们要算一个单位方块 $\Omega$ 内任意区域 $D$ 的面积，我们重复地做随机试验：在 $\Omega$ 上任取一点，用 $\xi_{n}$ 表示第 $n$ 次取点落在 $D$ 中这个时间的指标，那么 $D$ 的面积 $|D|=p|\Omega|$，其中 $p=P(\xi_{n}=1)$，这里 $p$ 是未知数，用大数定律，我们可以用频率来估计它，这样就可以得到 $D$ 的面积近似值。这就是一种简单的静态 Monte Carlo 方法的案例。

静态 Monte Carlo 方法通过构造独立同分布的随机数来计算积分，有频率法和期望法两种。

频率法

举个例子说明：计算任意定义域为 $[0,1]$ 的函数 $g(x)$ 在区间 $[0,1]$ 上的积分 $\int_{0}^{1}g(x)\mathrm{d}x$。

假设随机变量 $X$ 和 $Y$ 服从 $[0,1]$ 上的均匀分布，且相互独立，则二位均匀分布 $(X,Y)$ 的联合概率密度为：

$$f(x,y) = \left\{ \begin{aligned} & 1,\quad 0<x<1,0<y<1 \\ & 0,\quad\text{otherwise} \end{aligned} \right.$$ 现用 $B$ 表示事件 $\{w:Y\leq g(x)\}$，也即我们向矩形区域 $[0,1]\times[0,1]$ 随机投点，其中点落在以 $[0,1]$ 为底，以函数 $g(x)$ 为曲边的曲边梯形内。那么事件 $B$ 发生的概率为：
$$P(B)=\iint_{Y\leq g(x)}f(x,y)\mathrm{d}x\mathrm{d}y=\int_{0}^{1}\left[ \int_{0}^{g(x)}1\mathrm{d}y \right]\mathrm{d}x = \int_{0}^{1}g(x)\mathrm{d}x$$ 可以看出，积分的计算转化为求事件 $B$ 发生的概率。由大数定律可知，可重复试验中事件 $B$ 发生的频率近似表示事件 $B$ 发生的概率。

总结一下，具体做法如下：

1. 产生服从 $[0,1]$ 上均匀分布的随机数；
   • 模拟试验，考察 $n$ 次投点试验，记录时间 $B$ 发生的概率，用来近似表示其发生的概率，即得到积分值 $\int_{0}^{1}g(x)\mathrm{d}x$。

   • 期望法

     我们依然用上一小节的例子来说明。

     假设随机变量 $X$ 服从 $[0,1]$ 上的均匀分布，那么 $Y=g(X)$ 的期望为：
     

     $$E[Y]=E[g(X)]=\int_{0}^{1}g(x)\mathrm{d}x$$ 也就是说，积分的计算转化为计算 $g(X)$ 的数学期望值。由大数定律可知，若 $X_{n}$ 是独立同分布的随机变量序列，则 $\frac{1}{n}\sum_{i=1}^{n}X_{i}$ 依概率收敛到 $E[X_{i}]$，也就是说可用 $g(x)$ 的观察值的均值估计 $g(X)$ 的期望值。具体做法如下：
     1. 产生服从 $[0,1]$ 上均匀分布的随机数 $x_{i}$；
     2. 对每个 $x_{i}$ 计算 $g(x_{i})$ ，即可得到积分 $\int_{0}^{1}g(x)\mathrm{d}x$ 的估计值 $\frac{1}{n}\sum_{i=1}^{n}X_{i}$。

     动态 Monte Carlo 方法（MCMC）

     维数非常高的情况下，由于计算量太大，使用静态 Monte Carlo 方法处理速度太慢。动态 Monte Carlo 即 Markov Chain Monte Carlo 方法 （简称 MCMC）主要用于对维度非常高的随机向量取样。

     MCMC 方法首先建立一个 Markov 链，使得其极限分布是平稳分布。从目标分布中产生随机样本，就是从达到平稳状态的 Markov 链中产生样本路径。一个好的 Markov 链应满足从任意位置出发都能快速达到平稳分布这一性质。MCMC的理论依据是几个极限定理，下面简要介绍。

     遍历的 Markov 链（即不可约、正常返、非周期）的极限分布是平稳分布且是唯一平稳分布。考虑一个状态空间为 $S=\{ 1,2,\dots,N\}$ 的Markov链，记其转移概率为
     

     $$P = \begin{pmatrix} p_{11} & p_{12} &\dots & p_{1N} \\ p_{21} & p_{22} & \dots & p_{2N} \\ \vdots & \vdots & \ddots & \vdots \\ p_{N1} & p_{N2} & \dots & p_{NN} \end{pmatrix}$$ 对于遍历Markov链，极限分布为
     $$\pi_{j}=\lim_{n\to\infty}p_{ij}^{n}$$ 其中 $\pi_{j}$ 表示长时间运行后状态 $j$ 出现的时间比例，称为平稳概率。令 $\pi=(\pi_{1},\pi_{2},\dots, \pi_{N})$ ，那么 $\pi_{j}$ 是满足下列线性方程组的唯一解：
     $$\begin{split} \pi &=\pi P \\ \sum_{j=1}^{N}\pi_{j} &=1 \end{split}$$ 第一个方程表示 Markov 链处于状态 $j$ 的时间所占的比例等于 Markov 链从状态 $i$ 转移到状态 $j$ 所占的比例对 $i$ 求和；第二个方程表示 Markov 链处于状态 $j$ 的时间所占的比例对 $j$ 求和为 $1$。这个定理说明当 Markov 链运行足够长时间后的分布（极限分布）和初始分布无关。同时，从任意时刻开始以相反方向考察系统的变化情况，仍是一个转移概率为 $P$ 的Markov 链。其次，Markov链还有如下性质：假设 $\{X_{n},n=0,1,2,\dots\}$ 为一平稳分布为 $\pi$ 的遍历的Markov链，则 $X_{n}$ 依分布收敛到分布为 $\pi$ 的随机变量 $X$ ，且对任意函数 $g$ 当 $E_{\pi}[g(X)]$ 存在且 $n\to\infty$ 时，有
     $$\bar{g}_{n}=\frac{1}{n}\sum_{i=1}^{n}g(X_{i}) \to E_{\pi}[g(X)]$$ 换句话说， Markov链的实值函数的遍历均值几乎处处收敛到极限分布下的均值。若一个Markov链是一致几何遍历（转移速度以速度 $\lambda^{t}(0<\lambda<1)$ 收敛）的，并且 $f$ 相对于平稳分布 $\pi$ 是平方可积，则有
     $$\sqrt{n}\frac{\bar{f}_{n}-E_{\pi}[f(x)]}{\Gamma} \to N(0,1), n\to\infty$$ 这说明 Markov链的遍历均值做合适的变化后依分布收敛到标准正态分布。

     综上，我们可建立一个以 $\pi$ 为平稳分布的Markov链，则在运行此链足够长时间后，该Markov链会达到平稳状态。此时Markov链的值就相当于从分布 $\pi$ 中抽取的样本。

     Metropolis-Hastings 算法

     MCMC方法的重点在于构造合适的Markov链，Metropolis-Hastings 算法（简称 M-H算法）就是构造一个给定概率分布作为极限分布的 Markov 链的方法。

     M-H算法的原理此处暂时省略，这里只给出关键步骤。

     1. 构造合适的建议分布 $g(\cdot | X_{t})$，并产生服从该分布的 $X_{0}$。
     2. 从 $g(\cdot | X_{t})$ 中产生 $Y$，从 $U(0,1)$ 中产生 $U$，若 $$U\leq \frac{f(Y)g(X_{t}|Y)}{f(X_{t})g(Y|X_{t})}$$ 则接受 $Y$ 并令 $X_{t+1} = Y$，否则令 $X_{t+1} = X_{t}$。
     3. 重复步骤 2 中的过程直至Markov链达到平稳状态。

     Gibbs 抽样

     Gibbs抽样是M-H算法的一个特例，它将高维问题转化为一维问题。Gibbs抽样方法的重要特点是分别逐一对每个分类进行抽取。在对每个分量进行抽取时，是对其它所有分量的条件分布进行抽样的。

     这里也暂时仅介绍其关键步骤，以二维分布 $(X_{1},X_{2})$ 为例：

     1. 令 $(x_{1},x_{2})=X(t-1)$；
     2. 从 $f(x_{1}|x_{2})$ 中产生候选点 $X_{1}^{*}(t)$，更新 $x_{1}=X_{1}^{*}(t)$；
     3. 从 $f(x_{2}|x_{1})$ 中产生候选点 $X_{2}^{*}(t)$，更新 $x_{2}=X_{2}^{*}(t)$；
     4. 令 $X(t) = (X_{1}^{*}(t), X_{2}^{*}(t))$ 。

     贝叶斯 MCMC 估计

     应用贝叶斯方法分析问题的时候，一个主要的困难是得到的后验分布需要进行高维积分函数计算，通常高维积分的计算非常困难且耗时。MCMC方法是一种避免直接进行高纬积分计算的方法。

     假设 $p(y|\theta)$ 为抽样概率密度函数，其中 $\theta$ 是待估计未知参数向量。设未知参数的先验概率密度函数为 $\pi(\theta)$，则后验概率密度为
     

     $$\pi(\theta | y) = \frac{p(y|\theta)\pi(\theta)}{\int_{\Theta}p(y|\theta)\pi(\theta)\mathrm{d}\theta}$$ 实际问题中上述后验密度通常是比较复杂的未知形式。这些困难可以使用MCMC方法来解决。具体做法此处暂时不做详细介绍，敬请期待。