Neil Zhu，简书ID Not_GOD，University AI 创始人 & Chief Scientist，致力于推进世界人工智能化进程。制定并实施 UAI 中长期增长战略和目标，带领团队快速成长为人工智能领域最专业的力量。
作为行业领导者，他和UAI一起在2014年创建了TASA（中国最早的人工智能社团）, DL Center（深度学习知识中心全球价值网络），AI growth（行业智库培训）等，为中国的人工智能人才建设输送了大量的血液和养分。此外，他还参与或者举办过各类国际性的人工智能峰会和活动，产生了巨大的影响力，书写了60万字的人工智能精品技术内容，生产翻译了全球第一本深度学习入门书《神经网络与深度学习》，生产的内容被大量的专业垂直公众号和媒体转载与连载。曾经受邀为国内顶尖大学制定人工智能学习规划和教授人工智能前沿课程，均受学生和老师好评。
原文
贝叶斯方法的广泛实践应用的一个主要限制就是在获取后验分布时进行的高维度函数积分运算。这是一个相当困难的问题，但是也有了几种方式已经提出来解决了这个问题（Smith 1991，Evans and Swartz 1995，Tanner 1996）。本文聚焦Markov Chain Monte Carlo方法，该方法尝试从复杂的目标分布中直接进行取样。之所以称为MCMC，是因为这个方法生成了一个Markov Chain（新生成的样本只和前面一个样本有关系，也被称作无记忆性）。
在1990年由Gelfand和Smith提出的Gibbs sampler是MCMC方法的特例，该方法可以用于广泛的一类贝叶斯问题，于是导致了贝叶斯方法的迅速发展，这种强烈的兴趣推动了某个时代的到来。MCMC方法源自Metropolis算法（Metropolis and Ulam 1949，Metropolis et al. 1953）这是由物理学家做出的一个计算复杂积分的尝试得到的，他们将积分表示成某个分布的期望，通过从分布中取样来预测这个期望。Gibbs sampler诞生于图片处理领域。略带讽刺意味的是，MCMC方法最初对于统计领域没有太大的影响。MCMC方法在Tanner（1996）和Draper（2000）的工作中才得到深刻和详尽的研究。后面会给出额外的参考文献。
Monte Carlo Integration
最早的Monte Carlo思想就是由物理学家发展出来的使用随机数生成来计算积分的一套方法。假设我们要算一个复杂的积分
\int_a^b h(x) dx
如果我们可以将h(x)分解成定义在区间(a, b)一个函数f(x)和一个概率密度函数p(x)，注意到
\int_a^b h(x) dx = \int_a^b f(x)p(x) dx = E_{p(x)} [f(x)]
所以这个积分可以被表示成密度函数p(x)的函数f(x)的期望。因此如果我们从密度函数p(x)中取样出大量随机变量x_1,...,x_n，那么就有：
\int_a^b h(x) dx = E_{p(x)} [f(x)] ~ \frac{1}{n} \sum_{i=1}^{n} f(x_i)
这就是传说中的Monte Carlo Integration。
Monte Carlo积分可以用来近似后验分布（或者边缘后验分布）——这对于贝叶斯分析相当重要。考虑一下积分I(y) = \int f(y|x) p(x) dx，这个我们可以使用下面的方式近似：
\hat{I}(y) = \frac{1}{n} \sum_{i=1}^{n} f(y|x_i)
其中x_i就是从概率分布p(x)中抽样出来的。估计的Monte Carlo标准误差就是：
SE^2[\hat{I}(y)] = \frac{1}{n} (\frac{1}{n-1} \sum_{i=1}^{n} (f(y|x_i)-\hat(y))^2)
重要性取样（importance sampling）
假设概率密度p(x)粗略地近似密度函数q(x)（我们感兴趣的），那么
\int f(x) q(x) dx = \int f(x) (\frac{q(x)}{p(x)}) p(x) dx = E_{p(x)} [f(x) (\frac{q(x)}{p(x)})]
这是重要性取样的基础，
\int f(x) q(x) dx ~ \frac{1} \sum_{i=1}^n f(x_i) (\frac{q(x_i)}{p(x_i)})
其中x_i是从p(x)分布中取样出来的。例如，如果我们对一个关于y的边缘密度感兴趣，J(y) = \int f(y|x) q(x) dx，我们使用下面的方式近似它，
J(n) ~ \frac{1}{n} \sum_{i=1}^n f(y|x_i) (\frac{q(x_i)}{p(x_i)})
其中x_i是从近似密度函数p中取样出来的。
另一个对重要性取样的刻画方式是使用
\int f(x) q(x) dx ~ \hat{I} = \sum_{i=1}^n w_i f(x_i) / \sum_{i=1}^n w_i，其中w_i = \frac{q(x_i)}{p(x_i)}，x_i从p(x)中取样。与之相关的还有一个Monte Carlo方差：
Var(\hat(I)) = \sum_{i=1}^n w_i (f(x_i) - \hat{I})^2 / \sum_{i=1}^n w_i
Markov Chain介绍
略
The Metropolis-Hasting algorithm
应用Monte Carlo积分的一个问题是如何从复杂的概率分布p(x)中取样。对解决这个方法的各种尝试便是MCMC方法的开始。尤其是，由数学物理学家通过随机取样的方式对复杂函数进行积分（Metropolis and Ulam 1949，Metropolis et al. 1953，Hasting 1970）的方法最终成为了现在的Metropolis-Hastings算法。有关此方法更加详细参见Chib和Greenberg（1995）。
假设我们的目标是从某个分布p(\theta)中采样使得p(\theta) = f(\theta) / K，其中正规化常量K可能未知，并且很难计算。Metropolis算法从这个分布产生样本的过程如下：
- 从任何满足条件
f(\theta_0)>0的初始值\theta_0开始 - 使用当前的
theta值，从某个jumping分布q(\theta_1,\theta_2)中取样出候选样本点\theta^*，q是给定\theta_1前一个值获得\theta_2的值的概率。这个分布同样也被称为proposal或者candidate-generating分布。对于jump密度的唯一限制就是对称性（即q(\theta_1, \theta_2)=q(\theta_2,\theta_1)） - 给定候选样本点
\theta^*，计算密度函数在候选点\theta^*和当前点\theta_{t-1}的比值，
\alpha = \frac{p(\theta^*)}{p(\theta_{t-1})} = \frac{\theta^*}{f(\theta_{t-1})}，注意此处K已被消去了。 - 如果jump提高了密度（
\alpha > 1），接受这个候选样本点（设置\theta_t=\theta^*）并返回到step 2。如果jump降低了密度（\alpha < 1），则以概率\alpha接受候选点，否则拒绝他并返回step 2。
总结一下Metropolis采样：首先计算
\alpha = min(\frac{f(\theta^*)}{f(\theta){t-1})}, 1)
接着以概率\alpha接受候选样本点。这样产生一个Markov chain(\theta_0, \theta_1, ..., \theta_k, ...)，因为从\theta_t至\theta_{t+1}的转移概率只依赖于\theta_t。在一个充分的burn-id period后（k步以后），这个chain达到自身的稳定分布，接着就可以从p(x)中采样出n个样本了(n_{k+1},..., n_{k+n})。
Hastings在1970年推广了Metropolis算法，使用一个任意的转移概率函数q(\theta_1, \theta_2)=Pr(\theta_1 \rightarrow \theta_2)，并设置候选样本点的接受概率为
\alpha = min(\frac{f(\theta^*)q(\theta^*,\theta_{t-1})}{f(\theta_{t-1})q(\theta_{t-1},\theta^*)}, 1)
这就是Metropolis-Hastings算法。
