CWYAlpha

Just another WordPress.com site

Thought this was cool: LDA-math-MCMC 和 Gibbs Sampling(2)

leave a comment »


3 LDA-math-MCMC 和 Gibbs Sampling(2)
3.2 Markov Chain Monte Carlo

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

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

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

定理:[细致平稳条件] 如果非周期马氏链的转移矩阵$P$和分布$\pi(x)$ 满足
\begin{equation}
\pi(i)P_{ij} = \pi(j)P_{ji} \quad\quad \text{for all} \quad i,j
\end{equation}
则 $\pi(x)$ 是马氏链的平稳分布,上式被称为细致平稳条件(detailed balance condition)。

其实这个定理是显而易见的,因为细致平稳条件的物理含义就是对于任何两个状态$i,j$, 从 $i$ 转移出去到$j$ 而丢失的概率质量,恰好会被从 $j$ 转移回$i$ 的概率质量补充回来,所以状态$i$上的概率质量$\pi(i)$是稳定的,从而$\pi(x)$是马氏链的平稳分布。数学上的证明也很简单,由细致平稳条件可得
\begin{align*}
& \sum_{i=1}^\infty \pi(i)P_{ij} = \sum_{i=1}^\infty \pi(j)P_{ji}
= \sum_{i=1}^\infty \pi(j)P_{ji} = \pi(j) \\
& \Rightarrow \pi P = \pi
\end{align*}
由于$\pi$ 是方程 $\pi P = \pi$的解,所以$\pi$是平稳分布。

假设我们已经有一个转移矩阵为$Q$马氏链($q(i,j)$表示从状态 $i$转移到状态$j$的概率,也可以写为 $q(j|i)$或者$q(i\rightarrow j)$), 显然,通常情况下 $$ p(i) q(i,j) \neq p(j) q(j,i) $$ 也就是细致平稳条件不成立,所以 $p(x)$ 不太可能是这个马氏链的平稳分布。我们可否对马氏链做一个改造,使得细致平稳条件成立呢?譬如,我们引入一个 $\alpha(i,j)$, 我们希望
\begin{equation}
\label{choose-alpha}
p(i) q(i,j)\alpha(i,j) = p(j) q(j,i)\alpha(j,i)  \quad (*)
\end{equation}
取什么样的 $\alpha(i,j)$ 以上等式能成立呢?最简单的,按照对称性,我们可以取
$$ \alpha(i,j)= p(j) q(j,i), \quad \alpha(j,i) = p(i) q(i,j)$$
于是(*)式就成立了。所以有
\begin{equation}
\label{detailed-balance}
p(i) \underbrace{q(i,j)\alpha(i,j)}_{Q’(i,j)}
= p(j) \underbrace{q(j,i)\alpha(j,i)}_{Q’(j,i)}  \quad (**)
\end{equation}
于是我们把原来具有转移矩阵$Q$的一个很普通的马氏链,改造为了具有转移矩阵$Q’$的马氏链,而 $Q’$恰好满足细致平稳条件,由此马氏链$Q’$的平稳分布就是$p(x)$!

在改造 $Q$ 的过程中引入的 $\alpha(i,j)$称为接受率,物理意义可以理解为在原来的马氏链上,从状态 $i$ 以$q(i,j)$ 的概率转跳转到状态$j$ 的时候,我们以$\alpha(i,j)$的概率接受这个转移,于是得到新的马氏链$Q’$的转移概率为$q(i,j)\alpha(i,j)$。

mcmc-transition马氏链转移和接受概率

假设我们已经有一个转移矩阵Q(对应元素为$q(i,j)$), 把以上的过程整理一下,我们就得到了如下的用于采样概率分布$p(x)$的算法。

mcmc-algo-1

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

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

假设 $\alpha(i,j)=0.1, \alpha(j,i)=0.2$, 此时满足细致平稳条件,于是
$$ p(i)q(i,j)\times 0.1 = p(j)q(j,i) \times 0.2 $$
上式两边扩大5倍,我们改写为
$$ p(i)q(i,j) \times 0.5 = p(j)q(j,i) \times 1 $$
看,我们提高了接受率,而细致平稳条件并没有打破!这启发我们可以把细致平稳条件(**) 式中的$\alpha(i,j),\alpha(j,i)$ 同比例放大,使得两数中最大的一个放大到1,这样我们就提高了采样中的跳转接受率。所以我们可以取
$$ \alpha(i,j) = \min\left\{\frac{p(j)q(j,i)}{p(i)q(i,j)},1\right\} $$
于是,经过对上述MCMC 采样算法中接受率的微小改造,我们就得到了如下教科书中最常见的 Metropolis-Hastings 算法。

mcmc-algo-2

对于分布 $p(x)$,我们构造转移矩阵 $Q’$ 使其满足细致平稳条件
$$ p(x) Q’(x\rightarrow y) = p(y) Q’(y\rightarrow x) $$
此处 $x$ 并不要求是一维的,对于高维空间的 $p(\mathbf{x})$,如果满足细致平稳条件
$$ p(\mathbf{x}) Q’(\mathbf{x}\rightarrow \mathbf{y}) = p(\mathbf{y}) Q’(\mathbf{y}\rightarrow \mathbf{x}) $$
那么以上的 Metropolis-Hastings 算法一样有效。

3.2 Gibbs Sampling

对于高维的情形,由于接受率 $\alpha$的存在(通常 $\alpha < 1$), 以上 Metropolis-Hastings 算法的效率不够高。能否找到一个转移矩阵Q使得接受率 $\alpha=1$ 呢?我们先看看二维的情形,假设有一个概率分布 $p(x,y)$, 考察$x$坐标相同的两个点$A(x_1,y_1), B(x_1,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*}
所以得到
\begin{equation}
\label{gibbs-detailed-balance}
p(x_1,y_1)p(y_2|x_1) = p(x_1,y_2)p(y_1|x_1)  \quad (***)
\end{equation}

$$ p(A)p(y_2|x_1) = p(B)p(y_1|x_1) $$
基于以上等式,我们发现,在 $x=x_1$ 这条平行于 $y$轴的直线上,如果使用条件分布 $p(y|x_1)$做为任何两个点之间的转移概率,那么任何两个点之间的转移满足细致平稳条件。同样的,如果我们在 $y=y_1$ 这条直线上任意取两个点 $A(x_1,y_1), C(x_2,y_1)$,也有如下等式
$$ p(A)p(x_2|y_1) = p(C)p(x_1|y_1). $$

gibbs-transition平面上马氏链转移矩阵的构造

于是我们可以如下构造平面上任意两点之间的转移概率矩阵Q
\begin{align*}
Q(A\rightarrow B) & = p(y_B|x_1) & \text{如果} \quad x_A=x_B=x_1 & \\
Q(A\rightarrow C) & = p(x_C|y_1) & \text{如果} \quad y_A=y_C=y_1 & \\
Q(A\rightarrow D) & = 0 & \text{其它} &
\end{align*}

有了如上的转移矩阵 Q, 我们很容易验证对平面上任意两点 $X,Y$, 满足细致平稳条件
$$ p(X)Q(X\rightarrow Y) = p(Y) Q(Y\rightarrow X) $$
于是这个二维空间上的马氏链将收敛到平稳分布 $p(x,y)$。而这个算法就称为 Gibbs Sampling 算法,由物理学家 Gibbs 首先给出的。

gibbs-algo-1

two-stage-gibbs二维Gibbs Sampling 算法中的马氏链转移

以上采样过程中,如图所示,马氏链的转移只是轮换的沿着坐标轴 $x$轴和$y$轴做转移,于是得到样本 $(x_0,y_0), (x_0,y_1), (x_1,y_1), (x_1,y_2),(x_2,y_2), \cdots $ 马氏链收敛后,最终得到的样本就是 $p(x,y)$ 的样本,而收敛之前的阶段称为 burn-in period。额外说明一下,我们看到教科书上的 Gibbs Sampling 算法大都是坐标轴轮换采样的,但是这其实是不强制要求的。最一般的情形可以是,在$t$时刻,可以在$x$轴和$y$轴之间随机的选一个坐标轴,然后按条件概率做转移,马氏链也是一样收敛的。轮换两个坐标轴只是一种方便的形式。

以上的过程我们很容易推广到高维的情形,对于(***) 式,如果$x_1$ 变为多维情形$\mathbf{x_1}$,可以看出推导过程不变,所以细致平稳条件同样是成立的
\begin{equation}
\label{gibbs-detailed-balance-n-dimen}
p(\mathbf{x_1},y_1)p(y_2|\mathbf{x_1}) = p(\mathbf{x_1},y_2)p(y_1|\mathbf{x_1})
\end{equation}
此时转移矩阵 Q 由条件分布 $p(y|\mathbf{x_1})$ 定义。上式只是说明了一根坐标轴的情形,和二维情形类似,很容易验证对所有坐标轴都有类似的结论。所以$n$维空间中对于概率分布 $p(x_1,x_2,\cdots, x_n)$ 可以如下定义转移矩阵

  1. 如果当前状态为$(x_1,x_2,\cdots, x_n)$,马氏链转移的过程中,只能沿着坐标轴做转移。沿着 $x_i$ 这根坐标轴做转移的时候,转移概率由条件概率 $p(x_i|x_1, \cdots, x_{i-1}, x_{i+1}, \cdots, x_n)$ 定义;
  2.  其它无法沿着单根坐标轴进行的跳转,转移概率都设置为 0。

于是我们可以把Gibbs Smapling 算法从采样二维的 $p(x,y)$ 推广到采样$n$ 维的 $p(x_1,x_2,\cdots, x_n)$

gibbs-algo-2

以上算法收敛后,得到的就是概率分布$p(x_1,x_2,\cdots, x_n)$的样本,当然这些样本并不独立,但是我们此处要求的是采样得到的样本符合给定的概率分布,并不要求独立。同样的,在以上算法中,坐标轴轮换采样不是必须的,可以在坐标轴轮换中引入随机性,这时候转移矩阵 $Q$ 中任何两个点的转移概率中就会包含坐标轴选择的概率,而在通常的 Gibbs Sampling 算法中,坐标轴轮换是一个确定性的过程,也就是在给定时刻$t$,在一根固定的坐标轴上转移的概率是1。

相关文章:

  1. LDA-math-MCMC 和 Gibbs Sampling(1)
  2. LDA-math-认识Beta/Dirichlet分布(1)
  3. LDA-math-神奇的Gamma函数(3)
  4. LDA-math-认识Beta/Dirichlet分布(3)
  5. LDA-math-认识Beta/Dirichlet分布(2)
  6. LDA-math-神奇的Gamma函数(1)
  7. 概率语言模型及其变形系列-LDA及Gibbs Sampling
  8. LDA-math-神奇的Gamma函数(2)
  9. 正态分布的前世今生(四)
  10. 正态分布的前世今生(二)


from 我爱自然语言处理: http://www.52nlp.cn/lda-math-mcmc-%e5%92%8c-gibbs-sampling2?utm_source=feedburner&utm_medium=feed&utm_campaign=Feed%3A+52nlp+%28%E6%88%91%E7%88%B1%E8%87%AA%E7%84%B6%E8%AF%AD%E8%A8%80%E5%A4%84%E7%90%86%29

Written by cwyalpha

一月 13, 2013 在 8:23 上午

发表在 Uncategorized

发表评论

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / 更改 )

Twitter picture

You are commenting using your Twitter account. Log Out / 更改 )

Facebook photo

You are commenting using your Facebook account. Log Out / 更改 )

Google+ photo

You are commenting using your Google+ account. Log Out / 更改 )

Connecting to %s

%d 博主赞过: