如何做Gibbs采样(how to do gibbs-sampling)

原文地址：《如何做Gibbs采样(how to do gibbs-sampling)》

随机模拟

　　随机模拟（或者统计模拟）方法最早有数学家乌拉姆提出，又称做蒙特卡洛方法。蒙特卡洛是一个著名的赌场，赌博总是和统计有着密切的关系，所以这个命名风趣而贴切，被广为接受。
　　随机模拟的一个重要问题就是给定一个概率分布\(p(x)\)，如何生成它的样本。均匀分布\(Uniform(0,1)\)的样本可以通过线性同余发生器生成的伪随机数来模拟。常见的概率分布，无论是离散的还是连续的分布，都可以基于\(Uniform(0,1)\)的样本经过变换生成。
　　然而，当概率分布\(p(x)\)的形式很复杂或者\(p(x)\)是个高维分布的时候，样本的生成就会变得困难。此时就需要一些更加复杂的随机模拟方法来生成样本。接下来将会介绍下MCMC(Markov Chain Monte Carlo)和Gibbs Sampling算法。

马尔科夫链及其平稳分布

马尔科夫链的数学定义如下：
\[P(X_{t+1}=x| X_t, X_{t-1},...) = P(X_{t+1}| X_t)\]
定理 如果一个非周期马尔科夫链具有转移概率矩阵\(P\)，且它的任何两个状态是连通的，那么

1. $\lim_{n\rightarrow \infty}{P^n}= \left[ \begin{array}{ccccc} \pi(1) & \pi(2) &...&\pi(j)&...\\ \pi(1) & \pi(2) &...&\pi(j)&...\\ ...&...&...&...&...\\ \pi(1) & \pi(2) &...&\pi(j)&...\\ ...&...&...&...&...\\ \end{array} \right]. $

2. \(\pi(j)=\sum_{i=0}^{\infty}{\pi(i)P_{ij}}\).
3. \(\pi\)是方程\(\pi P=\pi\)的唯一非负解，其中\(\pi =[\pi(1), \pi(2), ..., \pi(j), ...], \sum_{i=0}^{\infty}{\pi_i} = 1.\)
   \(\pi\)称为马尔科夫链的平稳分布。

收敛到平稳分布

从初始概率分布\(\pi_0\)出发，我们在马尔科夫链上做状态转移，记\(X_i\)的概率分布为\(\pi_i\)，则有

\begin{array}{cc} X_0\sim \pi_0(x)& \\ X_i\sim \pi_i(x),& \pi_i(x) = \pi_{i-1}(x)P=\pi_0(x)P^i \end{array}

由马尔科夫链收敛的定理，概率分布\(\pi_i(x)\)将收敛到平稳分布\(\pi(x)\)。假设到第\(n\)步的时候马尔科夫链收敛，则有

\begin{array}{l} X_0 \sim \pi_0(x)\\ X_1\sim \pi_1(x)\\ ...\\ X_n\sim \pi_n(x) = \pi(x)\\ X_{n+1}\sim \pi(x)\\ X_{n+2}\sim \pi(x)\\ ... \end{array}

所以\(X_n,X_{n+1},X_{n+2},...\sim \pi(x)\)都是属于同分布的随机变量，当然他们并不是独立的。如果从\(x_0\)开始沿着马尔科夫链按照概率转移矩阵跳转，那么我们将会得到一个转移序列\(\{x_0,x_1,...,x_n,x_{n+1},...\}\)，由于马尔科夫链具有收敛性，所以\(\{x_{n},x_{n+1},...\}\)都是平稳分布\(\pi(x)\)的样本。

MCMC

定理(细致平稳条件) 如果非周期马尔科夫链的转移矩阵P和分布\(\pi(x)\)满足
\[\pi(i)P_{ij}=\pi(j)P_{ji}, ~for~all~ i,j\]
则\(\pi(x)\)是马尔科夫链的平稳分布，上式被称为细致平稳条件。
假设我们已经有了一个转移矩阵为\(Q\)的马尔科夫链(\(q(i\rightarrow j)\) 表示从状态\(i\)转移到状态\(j\)的概率)，\(p(x)\)表示当前的概率分布\(\pi\)，显然通常情况下：
\[p(i)q(i\rightarrow j) \ne p(j)q(j\rightarrow i)\]
也就是细致平稳条件不成立，所以\(p(x)\)不太可能是马尔科夫链的平稳分布。为了能够对\(p(x)\)进行抽样，我们必须对马尔科夫链进行改造，使得细致平稳条件成立。比如，引入一个接受率\(\alpha(i,j)\)使得
\[p(i)q(i\rightarrow j)\alpha(i,j) = p(j)q(j\rightarrow i)\alpha(j,i)\]
最简单的按照对称性，我们可以取如下\(\alpha(i,j),\alpha(j,i)\)使得上式成立
\[\alpha(i,j)=p(j)q(j\rightarrow i), ~ \alpha(j,i)=p(i)q(i\rightarrow j)\]

我们可以理解这个改造过程为在原来的马尔科夫链上，从状态\(i\)以概率\(q(i\rightarrow j)\)转移到状态\(j\)时，以概率\(\alpha(i,j)\)接受这个转移。所以新的马尔科夫链\(Q^\prime\)的转移概率为\(q(i\rightarrow j)\alpha(i,j)\)。

上面的改造还有一个小小的问题：马尔科夫链的接受率\(\alpha(i,j)\)在转移过程中可能会偏小，如果这样子采样过程中马尔科夫链容易原地踏步，拒绝大量跳转，致使收敛速度太慢。为了提高采样的接受率，我们可以对接受率进行细微的修改：\[\alpha(i,j)=\min \left \\{ \frac{p(j)q(j\rightarrow i) }{p(i)q(i\rightarrow j )}, 1 \right \\}.\]

Metropolis-Hastings采样算法

1. 初始化马尔科夫链初始状态\(X_0=x_0\)
2. 对t=0,1,2...,循环以下过程进行采样

• 第t时刻马尔科夫链状态\(X_t=x_t\)，采样\(y\sim q(x|x_t)\)
• 从均匀分布中采样\(u\sim Uniform(0,1)\)
• 如果\(u\lt \alpha(x_t,y)\)则接受转移\(x_t \rightarrow y\)，即\(X_{t+1} = y\)
• 否则不接受转移，即\(X_{t+1}=x_t\)

Gibbs Sampling

上述的Metropolis-Hasting算法也可以被应用与高维的场景，但是由于接受率的存在，算法的效率仍然较低。Gibbs Sampling，找到了一个接受率为1的转移矩阵。考察二维的情形，假设概率分布为\(p(x,y)\)。给定\(x\)坐标相同的两个点\(A(x_1,y_1),B(x_1,y_2)\)不难发现
\[ \begin{array}{c} 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)\\\\ p(x_1,y_1)p(y_2|x_1) = p(x_1,y_2)p(y_1|x_1)\\\\ \end{array} \]
也就是\(p(A)p(y_2|x_1)=p(B)p(y_1|x_1)\)。也就是说给定\(x=x_1\)这条平行于\(y\)轴的直线，使用分布\(p(y|x_1\)作为这条直线上任意两点的转移概率，该分布满足细致平稳条件。同样，对于点\(A(x_1,y_1),C(x_2,y_1)\)我们也可以得到\(p(A)p(x_2|y1)=p(C)p(x_1|y_1)\)。

构造如下转移矩阵
\[ \begin{array}{l} Q(A\rightarrow B) = p(y_B|x_1)\\\\ Q(A\rightarrow C) = p(x_C|y_1)\\\\ Q(A\rightarrow D) = 0 \end{array} \]
显然上面的转移矩阵\(Q\)满足细致平稳条件。于是这个二维的马尔科夫链将会收敛到平稳分布\(p(x,y)\)。这个算法被称为Gibbs Sampling，是由德国物理学家Gibbs首先提出的。
二维的算法很容易扩展到\(n\)维空间。对于\(n\)维空间的概率分布\(p(x_1,x_2,...,x_n)\)的转移矩阵如下：

1. 对于当前状态\((x_1,x_2,...,x_n)\)，马尔科夫链只能沿着坐标轴转移。比如，但状态沿着\(x_i\)转移时，其马尔科夫链转移概率为\(p(x_i|x_1,...,x_{i-1},x_{i+1},...,x_n)\)。
2. 其它无法沿着坐标轴进行跳转的转移概率均为0。

n维Gibbs Sampling算法

1. 随机选取一个初始状态\((x_1,x_2,...,x_n)\)
2. 对t=0,1,2... 循环采样：

• \(x_1^{(t+1)}\sim p(x_1|x_2^{(t)},x_3^{(t)},...,x_n^{(t)})\)
• \(x_2^{(t+1)}\sim p(x_2|x_1^{(t+1)},x_3^{(t)},...,x_n^{(t)})\)
• ...
• \(x_j^{(t+1)}\sim p(x_j|x_1^{(t+1)},x_2^{(t+1)},...,x_{j-1}^{(t+1)},x_{j+1}^{(t)},...,x_n^{(t)})\)
• ...
• \(x_n^{(t+1)}\sim p(x_n|x_1^{(t+1)},x_2^{(t+1)},...,x_{n-1}^{(t+1)})\)

转载于:https://www.cnblogs.com/forzhongyou/p/4573082.html