http://blog.csdn.net/pipisorry/article/details/51373090

吉布斯采样算法详解

为什么要用吉布斯采样

通俗解释一下什么是sampling。

sampling就是以一定的概率分布，看发生什么事件。举一个例子。甲只能E：吃饭、学习、打球，时间T：上午、下午、晚上，天气W：晴朗、刮风、下雨。现在要一个sample，这个sample可以是：打球+下午+晴朗。。。

问题是我们不知道p(E,T,W)，或者说，不知道三件事的联合分布。当然，如果知道的话，就没有必要用gibbs sampling了。但是，我们知道三件事的conditional distribution。也就是说，p(E|T,W),p(T|E,W),p(W|E,T)。现在要做的就是通过这三个已知的条件分布，再用gibbs sampling的方法，得到joint distribution。
具体方法。首先随便初始化一个组合,i.e. 学习+晚上+刮风，然后依条件概率改变其中的一个变量。具体说，假设我们知道晚上+刮风，我们给E生成一个变量，比如，学习-》吃饭。我们再依条件概率改下一个变量，根据学习+刮风，把晚上变成上午。类似地，把刮风变成刮风（当然可以变成相同的变量）。这样学习+晚上+刮风-》吃饭+上午+刮风。同样的方法，得到一个序列，每个单元包含三个变量，也就是一个马尔可夫链。然后跳过初始的一定数量的单元（比如100个），然后隔一定的数量取一个单元（比如隔20个取1个）。这样sample到的单元，是逼近联合分布的。

马氏链收敛定理

马氏链定理： 如果一个非周期马氏链具有转移概率矩阵,且它的任何两个状态是连通的，那么 存在且与无关，记, 我们有

1. 是方程 的唯一非负解

其中,   称为马氏链的平稳分布。

所有的 MCMC(Markov Chain Monte Carlo) 方法都是以这个定理作为理论基础的。 定理的证明相对复杂。

定理内容的一些解释说明

1. 该定理中马氏链的状态不要求有限，可以是有无穷多个的；
2. 定理中的“非周期“这个概念不解释，因为我们遇到的绝大多数马氏链都是非周期的；
3. 两个状态     是连通并非指 可以直接一步转移到   (),而是指   可以通过有限的步转移到达()。马氏链的任何两个状态是连通的含义是指存在一个, 使得矩阵 中的任何一个元素的数值都大于零。
4. 我们用 表示在马氏链上跳转第步后所处的状态，如果 存在，很容易证明以上定理的第二个结论。由于
   

   上式两边取极限就得到

[马尔科夫模型]

Monte Carlo随机采样方法

产生独立样本
基本方法:概率积分变换(第一部分已讲)
接受—拒绝(舍选)采样
重要性采样
产生相关样本:Markov Chain Monte Carlo
Metropolis-Hastings算法
Gibbs Sampler

这里主要讲马尔科夫链蒙特卡罗( Markov chain Monte Carlo, MCMC )，它使得我们可以从一大类概率分布中进行采样，并且可以很好地应对样本空间维度的增长。

皮皮blog

Metropolis-Hastings采样算法

马氏链的收敛定理非常重要，所有的 MCMC(Markov Chain Monte Carlo) 方法都是以这个定理作为理论基础的。

Idea

对于给定的概率分布p(x),我们希望能有便捷的方式生成它对应的样本。由于马氏链能收敛到平稳分布， 于是一个很的漂亮想法是：如果我们能构造一个转移矩阵为P的马氏链，使得该马氏链的平稳分布恰好是p(x), 那么我们从任何一个初始状态x0出发沿着马氏链转移, 得到一个转移序列 x0,x1,x2,⋯xn,xn+1⋯,， 如果马氏链在第n步已经收敛了，于是我们就得到了 π(x) 的样本xn,xn+1⋯。

这个绝妙的想法在1953年被 Metropolis想到了，为了研究粒子系统的平稳性质， Metropolis 考虑了物理学中常见的波尔兹曼分布的采样问题，首次提出了基于马氏链的蒙特卡罗方法，即Metropolis算法，并在最早的计算机上编程实现。Metropolis 算法是首个普适的采样方法，并启发了一系列 MCMC方法，所以人们把它视为随机模拟技术腾飞的起点。 Metropolis的这篇论文被收录在《统计学中的重大突破》中， Metropolis算法也被遴选为二十世纪的十个最重要的算法之一。

我们接下来介绍的MCMC 算法是 Metropolis 算法的一个改进变种，即常用的 Metropolis-Hastings 算法。由上一节的例子和定理我们看到了，马氏链的收敛性质主要由转移矩阵P 决定, 所以基于马氏链做采样的关键问题是如何构造转移矩阵P,使得平稳分布恰好是我们要的分布p(x)。如何能做到这一点呢？

细致平稳条件

定理：[细致平稳条件] 如果非周期马氏链的转移矩阵和分布     满足

则 是马氏链的平稳分布，上式被称为细致平稳条件(detailed balance condition)。

其实这个定理是显而易见的，因为细致平稳条件的物理含义就是对于任何两个状态, 从 转移出去到 而丢失的概率质量，恰好会被从 转移回 的概率质量补充回来，所以状态上的概率质量       是稳定的，从而     是马氏链的平稳分布。数学上的证明也很简单，由细致平稳条件可得

由于 是方程     的解，所以是平稳分布。

Note: 细致平稳条件是达到稳态的充分条件，并不是必要条件（也就是说要想达到稳态，必须要满足细致平衡条件，但也不是说满足细致平衡条件时的点就是稳态，而是说细致平衡条件是达到稳态的必经之路）。如在[马尔科夫模型]示例中0.28*0.286=0.08008，0.15*0.489=0.07335不相等，并不符合细致平稳条件。

构建满足条件的马氏链

假设我们已经有一个转移矩阵为马氏链(表示从状态转移到状态的概率，也可以写为或者), 显然，通常情况下

也就是细致平稳条件不成立，所以 不太可能是这个马氏链的平稳分布。我们可否对马氏链做一个改造，使得细致平稳条件成立呢？譬如，我们引入一个 , 我们希望

取什么样的 以上等式能成立呢？最简单的，按照对称性，我们可以取

于是(*)式就成立了。所以有

于是我们把原来具有转移矩阵的一个很普通的马氏链，改造为了具有转移矩阵的马氏链，而恰好满足细致平稳条件，由此马氏链的平稳分布就是！

在改造 的过程中引入的 称为接受率，物理意义可以理解为在原来的马氏链上，从状态 以     的概率转跳转到状态 的时候，我们以       的概率接受这个转移，于是得到新的马氏链的转移概率为。

Note：

1. 也就是说之前具有转移矩阵Q的马氏链可能收敛于另一个分布，所以我们要改造这个转移矩阵，使其转移矩阵Q’能够使新的马氏链收敛于p(x)。

2. 当按照上面介绍的构造方法把Q–>Q’后，就不能保证Q’是一个转移矩阵了，即Q’的每一行加和为1。这时应该在当 j != i 的时候概率Q'(i, j) 就如上处理， 当j = i 的时候， Q'(i, i) 应该设置Q'(i, i) = 1- 其它概率之和，归一化概率转移矩阵。

马氏链转移和接受概率：

MCMC采样算法

采样概率分布p(x)的算法，假设我们已经有一个转移矩阵Q(对应元素为q(i,j))

Note: 一开始的采样还没有收敛，并不是平稳分布（p(x)）的样本，只有采样了多次（如20次）可能收敛了，其样本才算是平稳分布的样本。

Metropolis-Hastings算法

{提高接受率alpha}

上述过程中      说的都是离散的情形，事实上即便这两个分布是连续的，以上算法仍然是有效，于是就得到更一般的连续概率分布  的采样算法，而   就是任意一个连续二元概率分布对应的条件分布。

以上的 MCMC 采样算法已经能很漂亮的工作了，不过它有一个小的问题：马氏链在转移的过程中的接受率 可能偏小，这样采样过程中马氏链容易原地踏步，拒绝大量的跳转，这使得马氏链遍历所有的状态空间要花费太长的时间，收敛到平稳分布 的速度太慢。有没有办法提升一些接受率呢?

假设       , 此时满足细致平稳条件，于是

上式两边扩大5倍，我们改写为

看，我们提高了接受率，而细致平稳条件并没有打破！这启发我们可以把细致平稳条件(**) 式中的               同比例放大，使得两数中最大的一个放大到1，这样我们就提高了采样中的跳转接受率。所以我们可以取

于是，经过对上述MCMC 采样算法中接受率的微小改造，我们就得到了如下教科书中最常见的 Metropolis-Hastings 算法。

对于分布 ,我们构造转移矩阵 使其满足细致平稳条件

此处 并不要求是一维的，对于高维空间的 ，如果满足细致平稳条件

那么以上的 Metropolis-Hastings 算法一样有效。

皮皮blog

Gibbs Sampling算法

{提高高维随机变量采样接受率alpha}

MCMC：Markov链通过转移概率矩阵可以收敛到稳定的概率分布。这意味着MCMC可以借助Markov链的平稳分布特性模拟高维概率分布；当Markov链经过burn-in阶段，消除初始参数的影响，到达平稳状态后，每一次状态转移都可以生成待模拟分布的一个样本。

而Gibbs抽样是MCMC的一个特例，它交替的固定某一维度，然后通过其他维度的值来抽样该维度的值，注意，gibbs采样只对z是高维（2维以上）（Gibbs sampling is applicable in situations where Z has at least two dimensions）情况有效。

基本算法如下：

1. 选择一个维度，可以随机选择；
2. 根据分布 抽样。

吉布斯采样满足细致平稳条件的转移矩阵的构造

对于高维的情形，由于接受率 的存在(通常), 以上 Metropolis-Hastings 算法的效率不够高。能否找到一个转移矩阵Q使得接受率 呢？我们先看看二维的情形，假设有一个概率分布 , 考察坐标相同的两个点，我们发现

所以得到

即

基于以上等式，我们发现，在 这条平行于 轴的直线上，如果使用条件分布  作为任何两个点之间的转移概率，那么任何两个点之间的转移满足细致平稳条件。

同样的，如果我们在 这条直线上任意取两个点 ,也有如下等式

于是我们可以如下构造平面上任意两点之间的转移概率矩阵Q

有了如上的转移矩阵 Q, 我们很容易验证对平面上任意两点 , 满足细致平稳条件

于是这个二维空间上的马氏链将收敛到平稳分布    。而这个算法就称为 Gibbs Sampling 算法,是 Stuart Geman 和Donald Geman 这两兄弟于1984年提出来的，之所以叫做Gibbs Sampling 是因为他们研究了Gibbs random field, 这个算法在现代贝叶斯分析中占据重要位置。

二维吉布斯采样算法

Note: 采样算法中右边的概率我们是知道的，例如你要采样的是二维高斯分布，那么固定xt后就是二维高斯分布固定xt后的一维高斯分布，且每次采样的坐标不同，这样这个一维高斯分布概率密度函数也就不一样了。

Gibbs Sampling 算法中的马氏链转移

Note：2.1步是从 (x0,y0)转移到(x0,y1)，满足Q(A→B)=p(yB|x1)的细致平稳条件，所以会收敛到平稳分布；同样2.2步是从(x0,y1)转移到(x1,y1)，也会收敛到平稳分布。也就是整个第2步是从 (x0,y0)转移到(x1,y1)，满足细致平稳条件，在循环多次后会收敛于平稳分布，采样得到的(xn,yn),(xn+1,yn+1)...就是平稳分布的样本（注意，并不是(xn,yn),(xn,yn+1),(xn+1,yn+1)...lz理解的，或许可以？）。

以上采样过程中，如图所示，马氏链的转移只是轮换的沿着坐标轴 轴和轴做转移，于是得到样本                 马氏链收敛后，最终得到的样本就是   的样本，而收敛之前的阶段称为 burn-in period。也就是说马氏链跳转过程就是采样的过程(采样的意思就是说在xt下，我们先计算出p(y|xt)的概率分布，并依这个概率分布采样某个yt+1)， 马氏链任何一个时刻 i 到达的状态 x_i 都是一个样本。 只是要等到 i 足够大( i > K) ， 马氏链收敛到平稳分布后, 那么 x_K, x_{K+1}, …. 这些样本就都是平稳分布的样本。lz：收敛到的平稳分布pi(x)是我们看不到的，我们看到的只是第t个循环的采样结果，在某个循环后，到达平稳分布pi(x)，在这之后的采样结果(xn,yn),(xn+1,yn+1)...统计一下计算不同样本的概率，就组成了我们可以看到的近似平稳分布（平稳分布的采样分布）。

坐标轴轮换采样

一般地，Gibbs Sampling 算法大都是坐标轴轮换采样的，但是这其实是不强制要求的。最一般的情形可以是，在  时刻，可以在轴和 轴之间随机的选一个坐标轴，然后按条件概率做转移，马氏链也是一样收敛的。轮换两个坐标轴只是一种方便的形式。

n维吉布斯采样算法

以上的过程我们很容易推广到高维的情形，对于(***) 式，如果变为多维情形，可以看出推导过程不变，所以细致平稳条件同样是成立的

此时转移矩阵 Q 由条件分布 定义。上式只是说明了一根坐标轴的情形，和二维情形类似，很容易验证对所有坐标轴都有类似的结论。

所以维空间中对于概率分布  可以如下定义转移矩阵

1. 如果当前状态为，马氏链转移的过程中，只能沿着坐标轴做转移。沿着 这根坐标轴做转移的时候，转移概率由条件概率     定义；
2. 其它无法沿着单根坐标轴进行的跳转，转移概率都设置为 0。

于是我们可以把Gibbs Sampling 算法从采样二维的    推广到采样 维的

   

Note:

1 Algorithm 8中的第2步的第6小步，x2(t)应为x2(t+1)。

2 要完成Gibbs抽样，需要知道如下条件概率：

Gibbs sampling makes it possible to obtain samples from probability distributions without having to explicitly calculate the values for their marginalizing integrals, e.g. computing expected values, by defining a conceptually straightforward approximation.

以上算法收敛后，每次完整的循环完成后得到的(xn+1,yn+1)...就是概率分布p(x1,x2,⋯,xn)的一个样本。

同样的，在以上算法中，坐标轴轮换采样不是必须的，可以在坐标轴轮换中引入随机性，这时候转移矩阵 Q 中任何两个点的转移概率中就会包含坐标轴选择的概率，而在通常的 Gibbs Sampling 算法中，坐标轴轮换是一个确定性的过程，也就是在给定时刻t，在一根固定的坐标轴上转移的概率是1。

吉布斯采样的相关性和独立性

马尔科夫链中的连续的样本是高度相关的（they were produced by a process that conditions on the previous point in the chain to generate the next point. This is referred to as autocorrelation (sometimes serial autocorrelation), i.e. correlation between successive values in the data.），这些样本并不独立，但是我们此处要求的是采样得到的样本符合给定的概率分布，并不要求独立。

吉布斯采样的接受率计算

吉布斯采样可以看做是M-H算法的一个特例，即接受率的情况。证明如下：
考虑一个M-H采样的步骤，它涉及到变量，剩余变量保持不变。同时，从z到的转移概率为。我们注意到，因为再采样步骤中，向量的各个元素都不变。又，因此，确定M-H算法中的接受概率为
=1

吉布斯采样收敛的判断

无论metropolis-hasting算法还是gibbs算法，都需要一个burn in过程，只有在达到平衡状态时候得到的样本才能是平衡状态时候的目标分布的样本，因此，在burn in过程中产生的样本都需要被舍弃。如何判断一个过程是否达到了平衡状态还没有一个成熟的方法来解决，目前常见的方法是看是否状态已经平稳（例如画一个图，如果在较长的过程中，变化已经不大，说明很有可能已经平衡）当然这个方法并不能肯定一个状态是否平衡，你可以举出反例，但是却是实际中没有办法的办法。

关于链的收敛有这样一些检验方法
（1）图形方法 这是简单直观的方法。我们可以利用这样一些图形：
          （a）迹图（trace plot）：将所产生的样本对迭代次数作图，生成马氏链的一条样本路径。如果当t足够大时，路径表现出稳定性没有明显的周期和趋势，就可以认为是收敛了。
           (b)自相关图（Autocorrelation plot）：如果产生的样本序列自相关程度很高，用迹图检验的效果会比较差。一般自相关随迭代步长的增加而减小，如果没有表现出这种现象，说明链的收敛性有问题。
          （c）遍历均值图（ergodic mean plot）：MCMC的理论基础是马尔科夫链的遍历定理。因此可以用累积均值对迭代步骤作图，观察遍历均值是否收敛。
（2）蒙特卡洛误差
（3）Gelman-Rubin方法

关于gibbs sampling的收敛，可以采用R^统计量。同时，可以多开几个chain进行模拟。判断收敛的时候，应该掐头去尾计算R^。

在工程实践中我们更多的靠经验和对数据的观察来指定 Gibbs Sampling 中的 burn-in 的迭代需要多少次。

Robert and Casella ( 1999 )总结了马尔科夫链蒙特卡罗算法的收敛性检测。

Jordan Boyd-Graber (personal communication) also recommends looking at Neal’s [15] discussion of likelihood as a metric of convergence.

[15] R. M. Neal. Probabilistic inference using markov chain monte carlo methods. Technical Report CRGTR-93-1, University of Toronto, 1993. http://www.cs.toronto.edu/∼radford/ftp/review.pdf.

burn-in阶段可能不需要

当然可以一条链跑到黑，但是一条链跑到黑只能是写学术 paper 的做法，在工程上可能还是要考虑很实际的速度和效率的问题，做 LDA 的时候我们就得考虑每秒钟能处理多少个请求,这时候不得不设置 burn-in。

Charles Geyer:为什么burn-in不是必要的Burn-In is Unnecessary

[Charles Geyer:为什么一条链走到黑是正确的One Long Run in MCMC]

Jason Eisner (personal communication) argues, with some support from the literature, that burn-in, lag, and multiple chains are in fact unnecessary and it is perfectly correct to do a single long sampling run and keep all samples. See [4, 5], MacKay ([13], end of section 29.9, page 381) and Koller and Friedman ([10], end of section 12.3.5.2, page 522).

[4] C. Geyer. Burn-in is unnecessary, 2009. http://www.stat.umn.edu/∼charlie/mcmc/burn.html, Downloaded October 18, 2009.
[5] C. Geyer. One long run, 2009. http://www.stat.umn.edu/∼charlie/mcmc/one.html.

[10] D. Koller and N. Friedman. Probabilistic Graphical Models: Principles and Techniques. MIT Press, 2009.
[13] D. Mackay. Information Theory, Inference, and Learning Algorithms. Cambridge University Press, 2003.

得到近似的独立样本

为了得到近似独立的样本可以对序列进行采样。

In order to avoid autocorrelation problems (so that the chain “mixes well”), many implementations of Gibbs sampling average only every L th value, where L is referred to as thelag.

the choice of L seems to be more a matter of art than science: people seem to look at plots of the autocorrelation for different values of L and use a value for which the autocorrelation drops off quickly. The autocorrelation for variable Z i with lag L is simply the correlation between the sequence Z i(t) and the sequence Z i(t−L) . Which correlation function is used seems to vary.

[Philip Resnik : GIBBS SAMPLING FOR THE UNINITIATED]

皮皮blog

吉布斯采样高级主题

隐变量Gibbs抽样公式

如果模型包含隐变量，通常需要知道后验概率分布  （如LDA中要推断的目标分布p(z|w)），但是这个分布因涉及很多离散随机变量等原因很难求解。

我们期望Gibbs抽样器可以通过Markov链利用全部的条件分布p(zi|z¬i,w)  来模拟p(z|w)  。所以包含隐变量的Gibbs抽样器公式如下：

Note：lz提示一下，这里分母实际只是分子对zi的一个积分而已，将变量zi积分掉，就得到p(z-i, x)，所以重点在联合分布p(z,w)公式的推导上，一般先推出联合分布公式再积分就可以使用上面的隐变量gibbs采样公式了。

多链Multiple chains

As is the case for many other stochastic algorithms (e.g. expectation maximization as used in the forward-backward algorithm for HMMs), people often try to avoid sensitivity to the starting point chosen at initialization time by doing multiple runs from different starting points. For Gibbs sampling and other Markov Chain Monte Carlo methods, these are referred to as “multiple chains”.

超参数的采样Hyperparameter sampling

Rather than simply picking hyperparameters, it is possible, and in fact often critical, to assign their values via sampling (Boyd-Graber, personal communication). See, e.g., Wallach et al. [20] and Escobar and West [3].

[3] M. D. Escobar and M. West. Bayesian density estimation and inference using mixtures. Journal of the American Statistical Association, 90:577–588, June 1995.

[20] H. Wallach, C. Sutton, and A. McCallum. Bayesian modeling of dependency trees using hierarchical Pitman-Yor priors. In ICML Workshop on Prior Knowledge for Text and Language Processing, 2008.

其它蒙特卡罗方法

Adaptive Rejection Sampling for Gibbs Sampling

吉布斯采样的实现

github Search · Gibbs Sampling

from: http://blog.csdn.net/pipisorry/article/details/51373090

ref: 随机采样方法整理与讲解（MCMC、Gibbs Sampling等）*

各种参数估计方法+Gibbs 采样精准细致的论述：Gregor Heinrich.Parameter estimation for text analysis*

MCMC案例学习

LDA-math-MCMC 和 Gibbs Sampling*

PRML

从随机过程到马尔科夫链蒙特卡洛方法

An Introduction to MCMC for Machine Learning，2003

Introduction to Monte Carlo Methods

吉布斯采样

吉布斯采样分布式实现

吉布斯采样与lda