随着计算机技术的发展，随机模拟成为了人类科学发展的重要工具。从1946年最早的简单版本的Monte-Carlo算法诞生开始，随机模拟方法就开始吸引人们的注意，并在20世纪产生了丰富且深刻的应用结果。被评为20世纪的十大算法之一。随着普通Monte-Carlo算法的不断发展，特别是1970年提出的Metropolis算法促进了这一算法的极大繁荣。1980年代，两个重要的MCMC算法出现在计算机视觉和人工智能领域的研究当中。直至1990年，MCMC算法开始在统计学中有重要影响。现在MCMC算法已经广泛应用在统计学，经济学，计算机科学等各个学科当中. MCMC算法经常被用来解决在高维空间中的积分和优化问题。这两类问题也在现代的机器学习，物理学，统计学，经济学，决策分析等领域扮演着基础性的角色. 本文将对MCMC算法的动机、原理做一些基本介绍，主要参考的是C.Andrieu,N.De Freitas, A. Doucet, M.I.Jordan的一个综述性的介绍文章 An introduction to MCMC for machine learning.
MCMC算法的动机和基本原理
Monte-Carlo原理
Monte-Carlo方法的基本思想是：为了求解一个问题，建立一定的概率模型或随机过程模型，使得其参量（如事件的概率，随即变量的期望等）等于所求问题的解，然后对模型或过程进行重复多次的随机抽样试验，并对结果进行统计分析，计算参量的值并得到所要求解问题的近似解.
一个简单的例子是利用MC方法计算积分的值.例如积分\[I=\iint_Df(x)dx\]，一个非常基本的MC做法是：生成D内的均匀随机数\(x_1,x_2,\cdots,x_n\),设\(D\)的面积为\(S_D\)，则 \[ E(f(X))=\iint_Df(x)\frac{1}{S_D}dx \] 由大数定律知\[ E(f(X))\approx \frac{1}{n}\sum_{k=1}^nf(x_k) \] 从而可以得到近似解 \[ \iint_D f(x)dx \approx \frac{1}{n}\sum_{i=1}^nf(x_i)S_D \] 一般认为MC方法的重要好处是估计值的精确程度与\(x\)的维度无关.但是当用MC算法处理更多问题时就会发现，单纯的事实上，生成特定分布的样本也是MC算法中的一个重要的问题.
重要度抽样与舍选法
在前面计算积分的例子中，实际上也可以不必采用均匀分布.事实上对任意一个在D上有支撑的连续型随机变量的密度函数\(p(x)\)有 \[ E[\frac{1}{N}\sum_{i=1}^N\frac{f(X_i)}{p(X_i)}]=\frac{1}{N}\sum_{i=1}^NE[\frac{f(X_i)}{p(X_i)}] =\frac{1}{N}\sum_{i=1}^N \frac{f(x)}{p(x)}p(x)dx =\int_a^b f(x)dx \] 这种抽样方法的好处是，当\(p\)取的比较合适时，可以减小随机数的方差，提高计算效率.然而当维数较高时，就很难找到合适的\(p\)
另一方面，为了解决生成特定分布随机数的问题，有一种间接的方法舍选法.其基本思想是给定一个与目标函数具有相同支撑集的候选密度函数\(f_c(y)\)，并计算\(M=\sup_y \frac{f(y)}{f_c(y)}\).若\(M<\infty\),则进行舍选法的核心迭代步骤：
- 生成独立的两个随机变量\(U\sim U(0,1),V\sim f_c(v)\)
- 如果\(U<\frac{1}{M}f(V)/f_c(V),\text{令}Y=V\)即为生成的一个符合要求的伪随机数；否则返回步骤\((1)\)
这个算法的有效性是因为可以证明 \[ P(Y\leq y)=P(V\leq y | U<\frac{1}{M}\frac{f(V)}{f_c(V)})=\frac{\int_{-\infty}^y \int_0^{\frac{1}{M}\frac{f(V)}{f_c(v)}}duf_c(v)dv}{\int_{-\infty}^{\infty} \int_0^{\frac{1}{M}\frac{f(V)}{f_c(V)}}duf_c(v)dv}=\int_{-\infty}^y f(v)dv \] 故\(Y\)的分布即服从所要求的分布函数.
下面考虑算法的迭代次数.计算\(N=\text{生成一个符合要求的}Y\text{需要的}(U,V)\text{对的数量}.\)注意到\(U,V\)相互独立，从而 \[ P(U<\frac{1}{M}\frac{f(V)}{f_c(V)})=\int_{-\infty}^{\infty} \int_0^{\frac{1}{M}\frac{f(V)}{f_c(V)}}duf_c(v)dv=\frac{1}{M} \] 从而\(N\sim G(\frac{1}{M})\)的几何分布.故\(EN=M\)因此调整候选密度的重要意义在于可以减少产生达到符合要求的伪随机数的迭代次数，从而提高算法效率.
但是对于维数较高的随机模拟问题，舍选法和重要度抽样一样，都会面临难以选择较为合适的候选密度的问题.因此在高维空间中生成服从给定密度的随机样本成为制约MC方法发展的一个关键问题.
MCMC算法的基本原理
在以上静态的模拟方法难以奏效的时候，基于马尔可夫链的MCMC采样方法开始展露头角.其关键思想在于：通过平稳的马尔可夫链进行转移计算，等效于从指定的密度函数\(P(x)\)的分布进行采样.
首先回顾马尔可夫链的定义.在最简单的有限状态空间中，设随机过程\(X_n\)的概率分布满足马尔可夫性，即 \[ P(X_n|X_{n-1},\cdots,X_1)=P(X_n|X_{n-1}) \] 则称该过程是一个马尔可夫过程.在\(X_n\)为离散时间时，称为马尔可夫链.描述一个马尔可夫链的常用方法是用它的转移概率矩阵.一步转移概率矩阵的元素\(p_{ij}^{(n)}=P(X_{n+1}=j|X_{n}=i)\),特别地当这个条件分布的值与\(n\)无关时，称该链为齐次的马尔可夫链.
而MCMC方法利用最多的马尔可夫链性质是它的平稳分布.对于一个有限状态的马尔可夫链，若存在一个概率分布\(\mathbf{\pi}=(\pi_1,\pi_2,\cdots,\pi_k)\),满足\(\mathbf{\pi}=\mathbf{\pi}P\)，其中\(P\)是单步概率转移矩阵，则称\(\mathbf{\pi}\)为该马尔可夫链的平稳分布.例如，假设有一马尔可夫链的单步转移概率矩阵为 \[ T=\begin{pmatrix} 0&1&0 \\ 0&0.1&0.9 \\ 0.6& 0.4 & 0 \end{pmatrix} \] 这个马尔可夫链存在一个极限分布\(\mathbf{\pi}=(0.2,0.4,0.4)\)事实上，由马尔可夫链相关知识知道，有限状态的不可约非周期马尔可夫链一定存在平稳分布
从任意的初始概率分布\(\pi_0\)出发，在马尔可夫链上做状态转移，记\(X_i\)的概率分布为\(\pi_i\)，则\(X_0\sim \pi_0(x),X_i\sim \pi_i(x),\pi_i(x)=\pi_{i-1}(x)P=\pi_0(x)P^n\) 因此概率分布\(\pi_i(x)\)将收敛到平稳分布\(\mathbf{x}\)
MCMC取样器一般选择目标分布函数座位马尔可夫链的平稳分布.在设计抽样器时，出了要满足平稳条件，产生平稳分布以外，也需要使得收敛速度能够尽可能的快.
常见的MCMC算法
Metropolis-Hastings 算法
根据马尔可夫链平稳分布的性质，对于一个给定的概率分布\(p(x)\),如果能构造一个转移矩阵为\(P\)的马尔可夫链，使得其平稳分布恰好是\(p(x)\),那么从任何状态\(x_0\)出发沿马尔可夫链进行转移，当它在第\(n\)步收敛后，即可得到一个样本\(x_n,x_{n+1},\cdots\)这就是1953年提出的Metropolis算法.
为了便于构造出能够产生所要求的分布\(p(x)\)的马尔可夫链的概率转移矩阵，需要满足所谓细致平稳条件条件:
若非周期马尔可夫链的转移矩阵\(P\)和分布\(\pi(x)\)满足 \[ \pi(i)P_{ij}=\pi(j)P_{ji} \] 则\(\pi(x)\)是马尔可夫链的平稳分布，上式被称为细致平稳条件.
证明：若细致平稳条件满足，则 \[ \sum_i \pi(i)P_{ij}=\sum_i \pi(j)P_{ji}=\pi(j)\sum_i P_{ji}=\pi(j) \] 此即\(\pi P = \pi\),也就说明了\(\pi\)是平稳分布. \end{proof} 一般情况下，细致平稳条件不成立，假设一个具有转移矩阵\(Q\)的普通马尔可夫链，考虑引入\(\alpha(i,j)\)满足 \[ p(i)q(i,j)\alpha(i,j)=p(j)q(j,i)\alpha(j,i) \quad(\star) \] 并取 \[ \alpha(i,j)=p(j)q(j,i),\quad\alpha(j,i)=p(i)q(i,j) \] 则\((\star)\)式成立，这样令\(Q‘(i,j)=q(i,j)\alpha(i,j)\),就把\(Q\)转化成了\(Q'\)，并且\(p(x)\)恰好就是其平稳分布.把引入的\(\alpha(i,j)\)称为接受率.这个过程理解为，在原来马式链上，从状态\(i\)以\(q(i,j)\)的概率转移到状态\(j\)的时候，以\(\alpha(i,j)\)的概率接受这个转移.整理以上过程，得到以下计算采样概率分布的算法.
- 初始化马式链初始状态\(X_0＝x_0\)
- 对\(t=0,1,2,\cdots\)，循环以下过程进行采样：
- 第\(t\)个时刻马式链的状态为\(X_t=x_t\),采样\(y\sim q(x|x_t)\)
- 从均匀分布采样\(u\sim U(0,1)\)
- 若\(u<\alpha(x_t,y)=p(y)q(x_t|y)\)，则接受转移\(x_t\rightarrow y\)，即\(X_{t+1}=y\)，否则\(X_{t+1}=x_t\)
上述过程中的\(p(x),q(x|y)\)都是离散情形.但是若它们连续，\(p(x),q(x|y)\)表示相应的密度函数和条件密度函数.
上述过程中仍有一个问题在于，接受率\(\alpha(i,j)\)有可能会严重偏小，使得收敛速度过慢.因此考虑将\(\alpha(i,j),\alpha(j,i)\)同比例放大，其中一个达到1，因此可以采用\(\alpha(i,j)=\min\{\frac{p(j)q(j,i)}{p(i)q(i,j)},1\}\)代替原来的接受率，即可得到常见的Metropolis-Hastings算法.
Gibbs抽样
对于高维的抽样过程，Gibbs取样器是一种常用的抽样方法.它实际上可以看作是MH算法的一个特例.例如，在\(n=2\)的情形下，我们希望构造马式链使得二维空间上的点\((x,y)\)能够具有平稳分布(目标分布函数\(p(x,y)\)).考虑两点\(A(x_1,y_1),B(x_2,y_2)\),则 \[ \begin{align*} p(x_1,y_1)p(y_2|x_1)&=p(x_1)p(y_1|x_1)p(y_2|x_1) \\ p(x_1,y_2)p(y_1|x_1)&=p(x_1)p(y_2|x_1)p(y_1|x_1) \\ \end{align*} \] \[ p(x_1,y_1)p(y_2|x_1)=p(x_1,y_2)p(y_1|x_1) \] 同理有 \[ p(x_1,y_1)p(x_2|y_1)=p(x_2,y_1)p(x_1|y_1) \] 在二维情形，可以构造如下的转移概率矩阵\(Q\) \[ \begin{align*} Q(B|A)=p(y_B|x_1), &\quad x_A=x_B=x_1 \\ Q(C|A)=p(x_C|y_1),&\quad y_A=y_C=y_1 \\ Q(D|A)=0,&\quad \text{其它} \end{align*} \] 这给出了二维情形下Gibbs Sampling的基本思想：轮换坐标轴循环采样.将这一思想推广到高维.,其中转移矩阵的定义满足：
- 如果当前状态为\((x_1,\cdots,x_n)\),在马式链转移的过程中，只能进行沿着坐标轴的转移，沿着\(x_i\)坐标轴转移时，转移概率由条件概率\(p(x_i|x_1,\cdots,x_{i-1},x_{i+1},\cdots,x_n)\)定义
- 其它无法沿着一根坐标轴进行的转移，转移概率为0.
由此可整理出Gibbs Sampling算法的流程如下：
- 初始化\(\{x_i:i=1,2,\cdots,n\}\)
对\(t=0,1,2,\cdots\)，循环采样：
- \(x_1^{(t+1)}\sim p(x_1|x_2^{(t)},x_3^{(t)},\cdots,x_n^{(t)})\)
- \(x_2^{(t+1)}\sim p(x_2|x_1^{(t+1)},x_3^{(t)},\cdots,x_n^{(t)})\)
- \(\cdots\)
- \(x_j^{(t+1)}\sim p(x_j|x_1^{(t+1)},\cdots,x_{j-1}^{(t+1)},x_{j+1}^{(t)},\cdots,x_n^{(t)})\)
- \(\cdots\)
- \(x_n^{(t+1)}\sim p(x_n|x_1^{(t+1)},x_2^{(t+1)},\cdots,x_{n-1}^{(t+1)})\)
需要指出的是Gibbs Sampling给出的样本并非是独立的.
总结
MCMC方法是20世纪计算机模拟技术的一个重大突破，它的诞生在人类许多科技项目上都有重大意义.其实MCMC方法利用的还是马尔可夫链的平稳分布的相关性质，从这一点我们也可以窥见随机过程对于现代科技发展的重要作用.当然MCMC算法本身的收敛效率等也是有赖于进一步的研究的.