从零开始学LDA（Latent Dirichlet allocation）

去年初学LDA，看完了Rickjin老师的《LDA数学八卦》，觉得不是很懂，查阅了很多资料之后才对LDA有了更深入地认识。一直想加入自己的理解后更简单的讲述这个模型，今日补上。

注：文章的行文思路与大多数公式参考了Rickjin《LDA数学八卦》。有兴趣更深入了解的同学可以在看完本篇后继续阅读。

LDA(Latent Dirichlet allocation)是主题模型的一种，最早是由Blei D M^[1]等人在2003年提出的，常被用于文档的主题识别。时至今日，在大大小小的会议上仍能看到不少对LDA（或其衍生、或类似的主题模型）提出改进的文章。其自然的思路，优雅的求解和拓展性强的模型结构，实在是非常适合作为对无监督学习一个入门模型来学习和改进。

本文主要介绍LDA背后的相关数学知识和模型的构建与求解思路，力求只有少量统计知识的初学者也能看懂。文中略去了大量的公式推导细节，如果希望深入的了解，请根据需要查找相关的资料。

¶从扔硬币说起

介绍模型之前，我们先来简单回顾一些知识点。本章主要讲解二项分布，贝叶斯理论，beta分布，multi分布和Dirichlet分布。如果你熟悉这些知识点，可以直接跳过这一章。

¶二项分布

二项分布（英语：Binomial distribution）是n个独立的是/非试验中成功的次数的离散概率分布，其中每次试验的成功概率为p。这样的单次成功/失败试验又称为伯努利试验。实际上，当n = 1时，二项分布就是伯努利分布。二项分布是显著性差异的二项试验的基础。

摘自二项分布维基百科

扔一枚硬币，一共扔N次，其中k次都是正面朝上的概率是多少？这个概率服从二项分布(Bernoulli distribution) $$P(k|N,p)=C^k_N\cdot p^k\cdot (1-p)^{N-k}$$ 其中p是硬币扔一次正面朝上的概率。真实场景中，我们面对的问题往往是，我想知道这个公式中的p是多少。一个我们熟悉且自然的思路是，把这个硬币扔N次，得到的结果中k次朝上的那么可以估计

$$p=\frac{k}{N}$$

如果你学过概率论，那么你就知道这是对二项分布中参数p的一个极大似然估计。

¶贝叶斯估计

当我们像上面去估计p值的时候，其实已经默认了一个假设，就是p是一个定值，它在每次实验中都是相同的。p是一个客观存的量，只是我们不知道它是多少。这样的想法很容易理解，但有时候会带来一些问题。当我们扔了10次硬币，10次都朝上的时候（尽管概率很小，但这一现象也会发生），用上面的方法估计，p=1。这一结论和我们的常识不符，当全部出现正面朝上这种小概率事件时，我么可能会倾向认为p值比0.5大，但没有大到等于1。

这时候，有些人就想，或许p并不是一个定值，而是服从某一种概率分布，对p的一个较好的估计是这个概率分布的期望，那么即使发生了每次都正面朝上的情况，只要这个概率分布不是始终等于1，就不会出现违背常识的情况了。

怎么理解？这就好比说，你扔出的硬币在落下的一瞬间之前，在你无法感知的空间中散落成无数的硬币，这些硬币的密度不尽相同，也就是落下时的p是不同的，但这些p是的有约束的，都是服从某一概率分布。在降落的一瞬间，有且仅有一枚硬币投射到你感知到的硬币上，也就是这一次抛硬币对应的p。由于无数硬币的p值并不是都等于1，所以即使出现扔硬币多次全部正面朝上，也可以理解成这些硬币中p值较大的硬币比较多，而不会是所有硬币的p都等于1。

如果你觉得上面这略显中二的说法也不太好理解的话😑，我们就直接看知识点好了。根据伯努利试验中p值是否固定，划分出了两个学派——频率学派和贝叶斯学派。这两个学派的思路差异主要体现在：

频率派把需要推断的参数θ看做是固定的未知常数，即概率虽然是未知的，但是确定的一个值，同时，样本X是随机的，所以频率派重点研究样本空间，大部分的概率计算都是针对样本X的分布；

贝叶斯派的观点则截然相反，他们认为参数是随机变量，而样本X 是固定的，由于样本是固定的，所以他们重点研究的是参数的分布。

在贝叶斯学派看来，所有的参数都不是固定值，而是服从一个概率分布，而对这个概率分布的估计有一个经典的模式：

先验分布𝜋(𝜃)+样本信息Χ⇒后验分布𝜋(𝜃|𝑥)

上述思考模式意味着，新观察到的样本信息将修正人们以前对事物的认知。换言之，在得到新的样本信息之前，人们对的认知是先验分布$𝜋(𝜃)$，在得到新的样本信息Χ后，人们对的认知为$𝜋(𝜃|𝑥)$。而这里求解$𝜋(𝜃|𝑥)$最重要依据就是贝叶斯公式（条件概率公式） $$P(𝜃,𝑋)=P(𝜃│𝑋)∙P(𝑋)=P(𝑋|𝜃)∙P(𝜃)$$

举个例子来说明两个学派分析的思路差异好了。加入我们扔一枚硬币，连续扔10次，有6次出现了正面朝上，那么扔一次正面朝上的概率是多少？

按照频率学派的观点， $$p=\frac{6}{10}=0.6$$

按照贝叶斯学派，首先p有一个先验分布，也就是你对p的一个经验估计，为了计算简单和表达的方便，咱们这里取一个离散分布,令$𝜃=p=0.4,0.5,0.6$，且

$P(𝜃=0.4)=0.1$

$P(𝜃=0.5)=0.8$

$P(𝜃=0.6)=0.1$

先验分布的期望 $$E(𝜃)=0.4\cdot0.1+0.5\cdot0.8+0.6\cdot0.1=0.5$$

根据样本信息，10次中有6次朝上，由贝叶斯公式变形可以得到，在该样本条件下，$𝜃=6$的概率为

$$𝑝(𝜃=0.6│𝑋=6)=\frac{𝑝(𝑋=6|𝜃=0.6)\cdot𝑝(𝜃=0.6)}{𝑝(𝑋=6)} $$

其中：

$𝜃$就是抛一次硬币是正面的概率，X是正面的次数，在本例中就是6；
$𝑝(𝜃=0.6│𝑋=6)$就是在6次正面的情况下，$𝜃=0.6$的概率，即后验概率；
$𝑝(𝑋=6|𝜃=0.6)$是𝜃=0.6时可以抛出6次正面的概率，即样本信息，这里求法按二项分布来求；
$𝑝(𝜃=0.6)$就是𝜃的先验分布，此例子中$𝑝(𝜃=0.6)=0.1$；
$𝑝(𝑋)$称为Normalizer，或者叫做marginal probability，是一个定值，$𝜃$分布如果是离散的，就是$𝑝(𝑋|𝜃)\cdot 𝑝(𝜃)$在所有𝜃情况下的和；如果$𝜃$是连续的，$𝑝(𝑋)=\int𝑝(𝑋|𝜃)\cdot 𝑝(𝜃)𝑑𝜃$。

我们就可以分别求出在样本信息为10次实验6次正面的情况下，后验分布

$P(𝜃=0.4|X=6)=0.056$

$P(𝜃=0.5|X=6)=0.819$

$P(𝜃=0.6|X=6)=0.125$

后验分布的期望 $$E(𝜃|X=6)=0.4\cdot0.056+0.5\cdot0.819+0.6\cdot0.125=0.5069$$

可以看到，后验分布的概率变得比0.5大了，而且并没有大很多。由于这里先验分布取的很特殊，哪怕出现无穷多次都是正面朝上的情况，后验分布的期望也只会逼近0.6。

¶Beta分布与共轭

上面的例子中，我们选取了一个离散分布作为先验分布，可以看出这个离散分布限定了p只能取三个值，得到的后验分布期望所在的区间也十分有限，选取这个分布并不太合适。数学家们给出了一个很漂亮的分布——Beta分布，来作为二项分布的先验分布。为什么选取Beta分布呢？主要是因为两点：

Beta分布形态多变，可以满足多种情况下的分布。

Beta分布做先验分布，样本信息服从二项分布时，两种分布符合良好的特性可以简化计算。

为了更好的理解Beat分布与二项分布的关系，我们先来看一下Beta分布的表达式。 $$p(\theta)=Beta(\theta|\alpha,\beta)=\frac{\theta^{\alpha-1}\cdot (1-\theta)^{\beta-1}}{B(\alpha,\beta)},$$ $$B(\alpha,\beta)=\int^1_0 t^{\alpha-1} (1-t)^{\beta-1}dt$$

公式中$\alpha$和$\beta$是Beta分布的两个超参数，通过调整$\alpha$和$\beta$的大小可以调整Beta分布的形态，如下图所示。

由图像我们不难看出，改变$\alpha$和$\beta$可以得到丰富的Beta分布的形式，得到不同的期望与分布的均匀程度，可以满足很多的分析场景。这是我们提到的选取Beta分布作为先验分布的第一点原因。

那么，它与二项分布的关系呢？我么不妨来求一下先验分布是Beta分布，样本服从二项分布的情况下，后验分布是什么样子的。如果先验分布是Beta分布

$$p(\theta)=Beta(\theta|\alpha,\beta)=\frac{\theta^{\alpha-1}\cdot (1-\theta)^{\beta-1}}{B(\alpha,\beta)}$$

样本服从二项分布

$$P(X|\theta)=C^k_N\cdot \theta^k\cdot (1-\theta)^{N-k}$$

那么根据贝叶斯公式，后验分布满足

$$\begin{align} p(\theta|X) & =\frac{p(X|\theta)\cdot p(\theta)}{p(X)}\\\\ & =\frac{p(X|\theta)\cdot p(\theta)}{\int p(X|\theta)\cdot p(\theta)d\theta}\\\\ & =\frac{(C^k_N\cdot \theta^k\cdot (1-\theta)^{N-k})\cdot(\frac{\theta^{\alpha-1}\cdot (1-\theta)^{\beta-1}}{B(\alpha,\beta)})}{\int (C^k_N\cdot \theta^k\cdot (1-\theta)^{N-k})\cdot (\frac{\theta^{\alpha-1}\cdot (1-\theta)^{\beta-1}}{B(\alpha,\beta)})d\theta}\\\\ & =\frac{\frac{C^k_N}{B(\alpha,\beta)}\cdot \theta^{\alpha+k-1}\cdot (1-\theta)^{\beta+N-k-1}}{\frac{C^k_N}{B(\alpha,\beta)}\cdot \int \theta^{\alpha+k-1}\cdot (1-\theta)^{\beta+N-k-1}d\theta}\\\\ & =\frac {\theta^{\alpha+k-1}\cdot (1-\theta)^{\beta+N-k-1}}{B(\theta|\alpha+k,\beta+(N-k))}\\\\ &= Beta(\theta|\alpha+k,\beta+(N-k)) \end{align}$$

得到的仍然是一个Beta分布！也就是说，在先验分布是Beta分布的情况下，如果样本服从二项分布，那么后验分布也是Beta分布。在贝叶斯估计中，将这种根据样本求得的后验分布与先验分布的分布形式一致时，样本服从的分布与先验分布共称为共轭分布。将上面的式子我们简化的记为

$$Beta(\theta|\alpha,\beta)+BinomCount(k,N-k)=Beta(\theta|\alpha+k,\beta+N-k)$$

共轭是一个非常好的特性。在之前的先验分布是离散分布的例子中，我们需要把每一步都求解出来才能知道后验分布的期望是多少，但是如果先验分布是Beta分布，我们只需要会求Beta分布的期望，就能很快的得到后验分布的期望了。

事实上，$Beta(\theta|\alpha,\beta)$的期望非常好求（这里就不推导公式了）

$$E(Beta(\theta|\alpha,\beta))=\frac{\alpha}{\alpha+\beta}$$

等下，这个结果是不是有点眼熟？如果我们另$\alpha=k,\beta=N-k$，在扔硬币的试验中，Beta的期望，也就是对p值的估计就变成了

$$E(Beta(\theta|k,N-k))=\frac{k}{N}$$

和频率学派对p值得估计统一起来了！从这里我们可以看出超参数$\alpha,\beta$的物理意义实际上就是扔硬币实验中正面朝上负面朝上的次数。根据我们的经验，我们认为p值应该是0.5，那么我们可以设先验Beta分布中$\alpha=\beta=1$，十次实验全部正面朝上，后验Beta分布的期望，也就是对p值的估计$E(\theta)=11/12$，一个接近1但不等于1的数。甚至，如果你对p值等于0.5非常有自信，你可以先验Beta分布中$\alpha=\beta=100$，那么十次实验全部正面朝上的话，$E(\theta)=101/102$，对p值的估计仍在0.5左右。记住这一点，对理解共轭的特性和Beta分布的含义有较大的帮助。

Beta分布是分布的分布，它与二项分布共轭。在Beta分布中，超参数$\alpha,\beta$物理含义是伯努利试验中两种结果发生的次数，它们比值决定了Beta分布的期望，他们的大小决定了样本对后验分布的影响程度。

¶多项分布与Dirichlet共轭

好的，下面我们来把二项分布推广吧。想象硬币就是一个“二面体”，它的每一面朝上的事件都会发生，而且发生的概率并不相同，对二项分布

$$P(k|N,p)=C^k_N\cdot p^k\cdot (1-p)^{N-k}$$ 令$x_1=k$，$x_2=N-k$，$p_1=p$，$p_2=1-p$，$\vec{x}=(x_1.x_2)$，$\vec{p}=(p_1,p_2)$那么上式可以改写为

$$p(\vec{x}|N,\vec{p})=\frac{N!}{x_1!\cdot x_2!}\cdot p^{x_1}_1\cdot p^{x_2}_2$$

这里式子的形式已经很清楚了。如果我们抛一个k面的凸多面体，一共抛N次，每个面朝上的次数为$\vec{x}$，每个面朝上的概率为$\vec{p}$,那么对应的概率就是多项分布

$$p(\vec{x}|N,\vec{p})=\frac{N!}{\prod_i^ kx_i}\cdot \coprod_{i}^ {k}p ^{x_i}_i$$

与多项分布共轭的分布就是Dirchlet分布

$$p(\vec{x})=Dir(\vec{x}|\vec{\alpha})=\frac{1}{\Delta (\vec{\alpha})}\coprod_{i}^ {k}x_i ^{\alpha_i-1}$$ $$\Delta (\vec{\alpha})=\int\coprod_{i} ^ {k}x ^{\alpha-1}_id \vec {x}$$

Dirchlet分布与多项分布的共轭也满足我们之前说的漂亮的性质

$$Dir(\vec{x}|\vec{\alpha})+MultCount(\vec{m})=Dir(\vec{x}|\vec{\alpha}+\vec{m})$$

而Dirchlet分布的期望等于

$$E(\vec{p})=(\frac{\alpha_1}{\sum_{k}^{i}\alpha_i},\frac{\alpha_2}{\sum_{k} ^ {i}\alpha_i}, \ldots ,\frac{\alpha_N}{\sum_{k} ^{i}\alpha_i})$$

多项分布与Dirihlet分布是整个LDA模型构建的基础。了解它们的特性对于理解和改进LDA模型会有很大的帮助。至于LDA模型到底是什么，我们下节继续。

¶Latent Dirichlet Allocation

OK，正餐来了。那么什么是LDA呢？让我们从Unigram Model开始说起。

¶Unigram Model

一篇文档是由许多词语构成的，事实上，抛开语法结构不谈，如果一篇文章是由一些没有顺序的词语构成，我们仍能从一定程度理解文档在探讨的主题。就好像那句有趣的话一样

研表究明，汉字序顺并不定一影阅响读

那么我们不妨将文档看做是一些词的集合，那么最简单的生成文档的方式就是，先选定这篇文档的词语个数N，然后选N个词出来。这一过程就好像是我们有一个巨大的多面体，它有字典里的词那么多面，每一面对应一个词。当我们要创作一篇有N个词的文档时，我们就把这个多面体扔N次，然后文档就生成了。

是的，上面这个简单的文档生成模型就叫做Unigram Model。它所对应的一篇文档的生成概率就是一个多项分布，根据词语的数量和总词数就可以估计每个词出现的概率。

Unigram model中文档生成的过程如下

Unigram model文档生成过程

假设：只有一个多面体，这个多面体有N面，每个面代表一个词；每个面的朝上的概率不尽相同。

过程：

1.选定文档的词语数量M；

2.抛掷多面体，记录朝上的面对应的词；

3.重复步骤共M次，直到所有词都生成。

由生成过程我们知道，当一篇文档每一个词确定下来之后得到$\vec{w}=(w_1,\cdots,w_V)$，则该文档的生成概率为

$$p(\vec{w})=p(w_1)\cdots p(w_V)$$

如果语料库中有M篇文档，那么语料库$W=(\vec{w}_1)\cdots \vec{w}_M$生成的概率为

$$p(W)=p(\vec{w}_1)\cdots p(\vec{w}_M)$$

假如语料库中的无重复的词语数量是N，每一个词在所有文档中出现的次数为$n_i$，那么$\vec{n}=(n_1,\cdots,n_N)$正好是一个多项分布，整个语料库生成的概率为 $$p(W)=p(\vec{n})=Mult(\vec{n}|\vec{p},N)=\coprod _{i=1} ^{N} p_i ^{n_i}$$

这里，对每一个参数$p_k$的一个估计值就是

$$\hat{p_i}=\frac{n_i}{N}$$

当然，别忘记了贝叶斯估计，如果给这个多面体每个面朝上的概率一个先验分布$Dir(\vec{p}| \vec{\alpha})$，那么就可以轻松地根据样本信息 $\vec{n}$ 来估计后验分布$Dir(\vec{p}| \vec{\alpha}+\vec{n})$了。它的概率图模型如下图所示。

这里稍微提一下概率图模型。我没系统学过，面试的时候老是被问，但是总也画不标准。后来看了这一篇，大概了解了一些基本规则。

这种记图的方式叫做plate notation，我们会用到符号包括

单圆圈表示隐变量

双圆圈表示观察到的变量

箭头表示采样

相同参数采样得到的独立同分布的变量放在一个方框里，并在方框中标注变量数量

这里，我们不难根据之前Normalizer的定义，得出先验分布是超参数为$\vec{alpha}$的Dirichlet分布的Unigram Model，文档生成概率如下（这个公式比较重要，我们后面会用到，记为公式3.1。）

$$\begin{align} p(\vec{W}|\vec{\alpha})& =\int p(\vec{W}|\vec{p}) \cdot p(\vec{p}|\vec{\alpha})d\vec{p}\\\\ &=\int \prod_{i=1}^{N}p_i^{n_i}\cdot Dir(\vec{p}|\vec{\alpha}))d\vec{p} \\\\ &=\int \prod_{i=1}^{N}p_i^{n_i} \cdot \frac{1}{\Delta (\vec{\alpha})} \prod_{i=1}^{N}p_i^{\alpha_i-1}d\vec{p} \\\\ &=\frac{1}{\Delta (\vec{\alpha})}\int \prod_{i=1}^{N}p_i^{n_i+\alpha_i-1}\vec{p}\\\\ &=\frac{\Delta(\vec{n}+\vec{\alpha})}{\Delta (\vec{\alpha})} \end{align}$$

¶主题模型

Unigram Model够简单，但不太符合我们日常的写作习惯。写作的人在写文档时，往往会先挑选一些主题。比如我在写这篇博客的时候，我想写一些关于机器学习，关于统计，关于工程实现方面的东西，那么我的这个多面体恐怕“煎饼果子”，“热干面”这样的词出现的概率就得很小，可当我写一些美食相关的博客时，这样一个多面体或许又显得有些不够用。

主题模型很好地解决了这一问题，相比Unigram Model，它添加了主题层，让文档生成的方式更加自然合理。一篇文章往往由多个主题构成，一个主题可以由多个与之相关的词语来描述。简单来说，当我们现在要创作一篇文当时，我们并不是只有一个大的多面体了，而是有多个“主题多面体”。在计算机主题相关的多面体上，内存、硬盘、复杂度、编程这些词出现的概率会更大；而在音乐相关的主题上。莫扎特、钢琴、二分音符、C大调这样的词出现的概率会更大。

当然，这时我们生成一篇文档，就需要两种多面体。第一种多面体有一个，它有K个面，每个面代表一个主题；第二种多面体有K个，每一个有W革面，每一个代表一个词。当我们生成一篇有N个词文档时，对每一个词，我们先扔第一种多面体，选定一个主题，然后扔对应的第二种多面体，得到最终的词。

这样直观的想法最初由Hoffmm在199年给出的PLSA（Probabilistic Latent Semantic Analysis）模型中进行了明确的数学化。PLSA在对文档主题建模时，除了词袋模型的大前提外，主要有三点假设：

一篇文章可以由多个主题构成

每一个主题由一组词的概率分布来表示

一篇文章中的每一个具体的词都来自于一个固定的主题

对于PLSA模型而言一篇文档生成的方式如下

PLSA文档生成过程

假设：有两类多面体分别是

doc-topic多面体，每一个doc-topic多面体有K个面，K代表主题个数，每个面对应一个主题编号。

topic-word多面体，一共有K个，每个topic-word多面体有V个面，每个面对应一个词。

过程：

1.选定文档的词语数量N；

2.抛掷doc-topic多面体，记录朝上的面对应的主题k；

3.抛掷k对应的topic-word多面体，记录朝上的面对应的词v

3.重复步骤2,3, N次，直到所有词都生成。

每篇文档生成的概率和整个语料库生成的概率就可以由多个多项分布的联合分布求出来，具体公式这里略去，我们直接开始讲LDA。

¶LDA文本主题建模

好了，终于到LDA了。

细心的你一定会想，unigram model中多面体的概率能用贝叶斯估计，那PLSA是否也能用贝叶斯估计能。当然！如果把PLSA中doc-topic多面体和topic-word多面体每一面朝上的概率加上Drichlet先验分布，就得了LDA模型！简单来说，可以理解为LDA就是把PLSA中的概率估计改造成了一个贝叶斯估计，选用的先验分布是Drichlet分布。（下面这个图我个人觉得把$Dir(\vec{\beta)}$）指向右边的topic-word多面体会更好一些。)

每一篇文章的主题分布服从Dirichlet-Multi共轭结构，这样的结构一共有M个；每一个主题的词分布也符合Dirichlet-Multi共轭结构，这样的结构一共有K个。LDA建模的本质么就是这个M+K个共轭结构和他们之间的关系，LDA模型的求解也就是求这M+K个共轭结构构成的联合分布的期望。

如果觉得这段话不好理解，也没关系，下面我们开始列公式，我们可以直接从公式中看出LDA模型的本质。

为了更好的讲解公式，统一用法，减少歧义，我们现在把即将用到的公式中的每一个变量和参数列在下面，供读者查阅。

K:主题个数

V：语料库中所有的、无重复的词语个数

M: 语料库中文档数量

$\vec{\alpha}$：doc-topic多面体的先验分布的超参数，维数为K。

$\vec{\beta}$：topic-word多面体的先验分布的超参数，维数为V。

$\vec{\theta}_m$：第m篇文档的主题分布，即这篇文档的doc-topic多面体投掷时出现每一个topic的概率，维度为K。

$\vec{\varphi}_k$：第k个主题的词分布，即这个主题的topic-word多面体投掷时出现每一个word的概率，维度为V

$N_m$：语料库中第m篇文档的词数量

$\vec{z}$：整个语料库的每一个词的主题构成的向量

$\vec{z}_m$：第m篇文档的每一个词的主题构成的向量，维度为N_m

$z_{m,n}$：第m篇文档的第n个词的主题

$w_{m,n}$：第m篇文档的第n个词的具体词语

$\vec{n}_m=(n_m ^1,\cdots,n_m ^K)$：$n_m ^K$表示在语料库中，第m篇文档中属于主题k的词语的个数。

$\vec{w}_k=(v_k ^1,\cdots,v_k ^V)$：第k个主题下每个词的个数。$v_k ^t$表示在语料库中，第k个主题包含的某一个具体的词语V的个数。

$\vec{w}=(\vec{w}_1,\cdots,\vec{w}_k)$：所有主题下的每个词的个数。

先上概率图。

由概率图和文档生成过程我们可以将LDA中每一个词生成的过程分解为两个大的物理过程，分别是

1.生成每一个词对应的主题，$\vec{\alpha} \rightarrow \vec{\theta_m} \rightarrow z_{m,n}$。在这个过程中，首先根据超参数$\vec{\alpha}$采样得到第m篇文档的主题分布$\vec{\theta_m}$，相当于的到了一个doc-topic骰子；然后投掷这个骰子，得到第n个词对应的topic编号$z_{m,n}$。这里M篇文档对应M个Dirichlet-Multi共轭结构。

2.根据主题生成具体的词，$\vec{\beta} \rightarrow \vec{\varphi_k} \rightarrow w_{m,n}$。在这个过程中，根据超参数$\vec{\beta}$采样得到第k个主题的词分布$\vec{\varphi_k}$，相当于的到了一个编号是$k=z_{m,n}$的topic-doc骰子；然后投掷这个骰子，得到第n个词对应的词$w_{m,n}$。这里K个主题对应K个Dirichlet-Multi共轭结构。

需要注意的是，在这两个过程中，每一个Dirichlet-Multi共轭结构都是独立的。这也就意味着，不仅可以一次就把$w_{m,n}$具体的词语生成出来，也可以先生成所有的$z_{m,n}$,再生成所有的$w_{m,n}$，这两个过程是完全等价的。根据这一特性，我们可以把LDA的文档生成过程定义为如下形式。

LDA文档生成过程

假设：有两类多面体分别是

doc-topic多面体，一共有M个，每个doc-topic多面体有K个面，每个面对应一个主题编号。这个多面体每个面朝上的概率分布服从超参数为$\vec{\alpha}$的Dirichlet分布。

topic-word多面体，一共有K个，每个topic-word多面体有V个面，每个面对应一个词。这个多面体每个面朝上的概率分布服从超参数为$\vec{\beta}$的Dirichlet分布。

过程：

选定每一篇文档的词语数量$N_m$；

对第m篇文档，根据超参数$\vec{\alpha}$采样的到对应的doc-topic多面体。

对第m篇文档的第n个词，投掷doc-topic多面体，得到该词的主题$z_{m,n}$。

重复步骤3共$N_m$次，得到第m篇文档的每一个词的主题编号$\vec{z_m}$。

重复步骤2-4共M次，得到每一篇文档的的每一个词的主题编号。（过程1结束）

对第k个主题，根据超参数$\vec{\beta}$采样的到对应的topic-word多面体。

对第m篇文档的第n个词，根据主题编号选择对应的topic-word多面体，投掷该多面体，得到对应的词$w_{m,n}$。

重复步骤7，直到所有词都得到具体的词。（过程2结束）

从这一过程就可以很清晰的看出，LDA模型求解的过程实际上就是求这M+K个Dirichlet-Multi共轭结构下文档生成概率的联合分布了。那么Dirichlet-Multi共轭结构如何如何求解呢？

在过程1中，超参数为$\vec{\alpha}$的第m篇篇文档下每一个词的主题编号生成概率，我们可以由公式3.1得到

$$p(\vec{z}_m|\vec{\alpha}) =\frac{\Delta(\vec{n}_m+\vec{\alpha})}{\Delta (\vec{\alpha})}$$

同时，我们得到隐含参数$\vec{\theta}_m$的后验分布恰好是

$$Dir(\vec{\theta}_m|\vec{n}_m+\vec{\alpha})$$

整个语料库下的所有词的主题编号的生成概率就是每一篇文档的联合分布

$$p(\vec{z}|\vec{\alpha}) =\prod_{m} ^{M} \frac{\Delta(\vec{n}_m+\vec{\alpha})}{\Delta (\vec{\alpha})}$$

同样的道理，我们可以得到过程二中，每个主题下的词生成概率为 $$p(\vec{w}_k|\vec{\beta}) =\frac{\Delta(\vec{v}_k+\vec{\beta})}{\Delta (\vec{\beta})}$$

同时，我们得到隐含参数$\vec{\varphi}_k$的后验分布恰好是

$$Dir(\vec{\varphi}_k|\vec{v}_k+\vec{\beta})$$

对于整个语料库的所有词的生成过程来说，每一个词具体的主题编号也是最终决定词语的先决条件，但是无论选择哪一个主题编号，并不会改变着K个主题下的词分布(这一点很重要，可以将过程1和过程2分开来看)，所以有 $$p(\vec{w}|\vec{z},\vec{\beta})=p(\vec{w}|\vec{\beta})=\prod_{k} ^{V} \frac {\Delta(\vec{v}_k+\vec{\beta})}{\Delta (\vec{\beta})}$$

结合上面的式子，我们最终可以得到，两个过程加在一起，LDA模型的语料库的生成概率为

$$p(\vec{w},\vec{z}|\vec{\alpha},\vec{\beta})=p(\vec{w},\vec{z}|\vec{\beta})\cdotp(\vec{z}|\vec{\alpha})=\prod_{k} ^{V} \frac {\Delta (\vec{v} _k + \vec{\beta}) } { \Delta ( \vec{\beta}) } \cdot \prod _{m} ^{M} \frac{\Delta(\vec{n}_m+\vec{\alpha})}{\Delta (\vec{\alpha})}$$

到这里，整个LDA模型的文档生成过程就介绍完了。LDA模型在实际应用中，就是要估计隐含变量$\vec{\theta}_m$和$\vec{\varphi}_k$的值，也就是要求出doc-topic多面体和topic-word多面体的每个面朝上的概率大小，由此我们可以去分析文档的主题分布和主题的词分布。那么，这两个概率大小究竟该怎么求呢，我们继续往下探讨。

¶LDA模型的求解

OK，模型的基本构成我们了解完了，接下来就要看模型的求解了。对于LDA模型来说，最基本的目标有两个：

Training。根据一个已有的语料库，我们需要这个语料库用LDA方式生成过程中的参数$\vec{\theta}_m$和$\vec{\varphi}_k$，由此可以把生成这个语料库的LDA模型完全确定下来。
Inference。对于一篇不在这个语料库中的文档，我们能够根据已经确定的LDA模型来计算这篇文档的主题分布。

当然，要先进行Training，然后才能进行Inference。Training的过程实际上就是对于上一节中我们得到的$p(\vec{w},\vec{z}|\vec{\alpha},\vec{\beta})$分布求期望，由于$p(\vec{w})$是我们观测到的已知变量，实际上就是要求分布$p(\vec{z}|\vec{w}$的期望，也就是每一个词对应的最可能的主题编号。

如何求这个期望是一个难题。早期的PLSA和LDA的求解都是基于最大期望（Expectation Maximization）算法做的，要对$p(\vec{w},\vec{z}|\vec{\alpha},\vec{\beta})$直接积分求解，这一过程的时间复杂度和空间复杂度都比较高，好处是可以得到精确解。在工程上和较新的一些关于主题模型的研究中，则往往采用Gibbs Sampling的方式求解。这种方式得到的是近似解，但是求解速度快，工程实现也比较简单。下面，我们就讲解如何通过Gibbs Sampling的方式求解LDA模型。

¶采样

Sampling，就是采样。理解Gibbs Sampling的基础，是要明白采样的目的和原理。所谓采样，是指获得符合某一概率分布的一些列样本，根据大数定律，当这些样本足够多时，样本的期望会逐渐趋近于概率分布的期望。也就是说，采样的目的就是通过样本期望去估计分布的期望。

举个例子，有一个每个面大小不相同的骰子，现在想知道每个面朝上的概率大小是多少，怎么做？最先想到的最简单的方法，就是采样。不断投掷这个骰子，记录下每次投掷朝上的面，然后计算这些面出现的次数占总次数比，就可以大致的估计每个面朝上的概率。

对LDA求解也是一样，直接求解期望很麻烦，那么我们不妨通过采样得到一系列符合联合分布的样本，计算样本的期望就可以估计分布的期望。如何得到符合联合分布的样本呢，这就要引入另一组数学知识——马氏链和蒙特卡洛方法。

¶马氏链

马氏链（Markov Chain）的数学定义很简单， $$P(X_{t+1}=x|X_t,X_{t-1},\cdots)=P(X_{t+1}=x|X_t\cdots)$$ 即在一个状态转移链中，如果下一时刻状态转移的概率只和当前状态有关，那么这个转移链就是一个马氏链。

用一个例子来说明马氏链和马氏链产生的过程。（这个例子我实在《LDA数学八卦》里看到的，不过网上好像用烂了，最早应该也不是这本书中提出的）社会学家经常把人按其经济状况分成3类：下层(lower-class)、中层(middle-class)、上层(upper-class)，我们用1,2,3 分别代表这三个阶层。社会学家们发现决定一个人的收入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别，那么他的孩子属于下层收入的概率是 0.65, 属于中层收入的概率是 0.28, 属于上层收入的概率是 0.07。事实上，从父代到子代，收入阶层的变化的转移概率如下。

我们用矩阵的形式表示，状态转移矩阵记为

$$P=\begin{bmatrix} 0.65 & 0.28 & 0.07\\ 0.15 & 0.67 & 0.18\\ 0.12 & 0.36 & 0.52 \end{bmatrix}$$

我们假设当前这一代下、中、上层人口的比例是$\pi_0=[0.21,0.68,0.11]$，那么下一代的比例分布式$\pi_0=\pi_1\cdot P$，以此类推，可以计算前n代人中下、中、上层人口比例为

从第7代开始，这个比例就不变了。这是偶然吗？换一个初始概率$\pi_0=[0.75, 0.15, 0.1]]$再试一次，前n代人中下、中、上层人口比例为

第9代人的时候, 分布又收敛了。奇特的是，两次给定不同的初始概率分布，最终都收敛到概率分布$\pi=[0.286,0.489,0.225]$，也就是说收敛的行为和初始概率分布 $\pi_0$ 无关。这说明这个收敛行为主要是由概率转移矩阵P决定的。我们计算一下这个矩阵的n次方，可知

$$ P^{20} =P^{21}=\cdots =P100=\cdots =\begin{bmatrix} 0.286 & 0.489 & 0.225\\ 0.286 & 0.489 & 0.225\\ 0.286 & 0.489 & 0.225 \end{bmatrix} $$

当n足够大的时候，P^n的每一行都会收敛于分布$\pi=[0.286,0.489,0.225]$。这并不是这个过程特有的，这里我们引入一个非常重要的定理。

马氏链收敛定理：如果一个非周期马氏链具有转移概率矩阵P,且它的任何两个状态是连通的，那么$\lim_{n\rightarrow \infty }P_{ij}$ 存在且与i无关，记 $\lim_{n\rightarrow \infty }P_{ij}=\pi(j)$, 我们有

1.$\lim_{n\rightarrow \infty }P^n =\begin{bmatrix} \pi(1) & \pi(2) & \cdots &\pi(j) &\cdots \\ \pi(1) & \pi(2) & \cdots &\pi(j) &\cdots \\ \cdots & \cdots & \cdots &\cdots &\cdots \\ \pi(1) & \pi(2) & \cdots &\pi(j) &\cdots \\ \cdots & \cdots & \cdots &\cdots &\cdots \\ \end{bmatrix} $

2.$\pi(j)=\sum_{i=0}^{\infty }\pi(i)P_{ij}$

3.$\pi$是方程$\pi P=\pi$的唯一非负解

其中 $$\pi=[\pi(1),\pi(2),\cdots,\pi(j),\cdots], \sum_{i=0}^{\infty} \pi_i =1$$ $$P^n_{ij}指状态矩阵n次方的第i行第j列$$

公式的证明以及更详细的说明这里不说了，在LDA中会用到的马氏链都是非周期的。关于马氏链的收敛定律，对理解LDA来说，最重要的是记住两点：

平稳分布可以是离散的，也可以是连续的
对于一个转移矩阵为P的马氏链，在足够长（达到收敛）的步数之后，对象处在不同状态的概率是一定的，与时间无关。

什么意思？其实我最早看马氏链的时候，觉得马氏链收敛定律的思想是好理解的，但是一直不是很明白这和采样有什么关系。我们还是用之前的社会阶层的例子来说明这个问题。

现在我们不观察整个社会在不同阶层的分布了，而是集中到一个人和他的后代上。初始的分布$\pi(0)$对应的就是最初的这一个人落在不同阶层的概率。让这个人按这个概率落在一个阶层上，比如阶层1，那么他的后代就会按照状态转移矩阵的第一列对应的概率落在不同的阶层上。按照第一种情况，在第7代时整个社会的分布收敛到了平稳分布，也就是说第七代后代落在三个阶层的概率与第八代后代落在三个阶层的概率是一样的，与后面所有后代落在三个阶层的概率也是一样。那么如果把这个人第七代后代之后很多代的子孙所处的阶层记录下来，我们就得到了一组服从平稳分布的样本。

根据马氏链收敛定律，不管初始分布是怎样的，在足够长（达到收敛）的步数之后，对象处在不同状态的概率是一定的，与时间无关。那么我们最初让观察对象随机的处在一个状态，让对象按照这个马氏链进行转移，达到足够长的步数，记录之后时刻对象的状态，就得到了服从平稳分布的一组样本。

这一点意义重大！因为如果我们能构造一个转移矩阵P，使得马氏链的平稳分布正好符合我们预期的分布，那么从任意的初始状态出发，在足够长的步数之后，我们就能得到符合预期分布的样本，从而估计样本的期望。而MCMC算法，就正好运用了这一思想。

¶MCMC

Markov Chain Monte Carlo，从名字就能看出来，MCMC算法就是马氏链+蒙特卡洛。所谓蒙特卡洛，是对一类具有相同特性的随机算法的统称。这类方法的特点是，可以在随机采样上计算得到近似结果，随着采样的增多，得到的结果是正确结果的概率逐渐加大，但在（放弃随机采样，而采用类似全采样这样的确定性方法）获得真正的结果之前，无法知道目前得到的结果是不是真正的结果（摘自知乎@陈义的回答）

举个例子，一个数组有一亿个数，我想知道他们的均值，就随机抽取其中的一些，把这些样本的均值当做总体的均值。我抽的越多，这个值越准，但是除非我把所有数遍历一边，不然我还不知道真正的结果。

我们用MCMC方法来求解Beta分布期望的例子来说明MCMC算法的思路。这里蒙特卡洛就体现在，我们会去采样符合这个Beta分布的样本，样本越多，这个期望符合正确解的概率就越大，但是无论样本多大，只要是有限个，就不能认定是真实值。

根据之前的思路，现在的目标是要构造一个转移矩阵P，使得平稳分布是Beta分布。如何构造呢，我们要引入下面一个定理：

马氏链收敛的细致平稳条件

如果非周期马氏链的转移矩阵P和分布$\pi(x)$满足 $$\pi(i)P_{ij}=\pi(j)P_{ji}$$ 则 $\pi(x)$是马氏链的平稳分布，上式被称为细致平稳条件(detailed balance condition)。

证明略去。公式很好理解，说的就是两个状态之间互相转移消耗的代价是相等的。总之这里重要的是我们要构造一个转移矩阵，让上面的公式成立。

假设我们已经有一个转移矩阵为Q的马氏链，其中$q(i,j)$表示从状态i转移到状态j的概率，也就是矩阵的第i行第j列，也可以写成$q(i|j)$或者$q(i \rightarrow j)$，显然，通常情况下细致平稳条件不成立，

$$p(i)q(i,j)\neq p(j)q(j,i)$$

也就是说此时的p(x)不太可能是这个马氏链的平稳分布，为了满足这个平稳条件，我们不妨对马氏链进行简单的改造，引入一个参数$\alpha(i,j)$，使得

$$p(i)q(i,j)\alpha(i,j)=p(j)q(j,i)\alpha(j,i)$$

这样平稳细致条件就成立了诶，我们就可以做采样了。那么$\alpha(i,j)$又该怎么选呢？观察上式，不难发现，直接按对称性，可以取

$$\alpha(i,j)=p(j)q(j,i)\quad \alpha(j,i)=p(i)q(i,j)$$

此时有

$$p(i)\underbrace{q(i,j)\alpha(i,j)} _{ {q}’ (i,j)} =p(j)\underbrace{q(j,i)\alpha(j,i)} _{ {q}’ (j,i)}$$

这样转移矩阵为${Q}’$的马氏链满足平稳细致条件，平稳分布就是p(x)啦，就可以做采样啦~

卧槽这是在干嘛？怎么这么随便就添加一个东西,然后这个东西其实还是原来的式子啊，怎么就能采样了😑😑😑

我刚开始也这么想，但其实这正是这个算法的厉害之处。这里需要理解和注意的是，最开始的转移矩阵Q可以是任意构造的（只要满足任意两状态联通，保证马氏链可以收敛）！也就是说，实际上在整个过程中，我们真正重要的就是我们要采样的分布p(x)，其他的都是构造出来了。在知道p(x)的情况下，可以用这个方法对任意的分布做采样！牛逼之极！

这里引入的$\alpha(i,j)$成为接受率，物理意义可以理解为在原来的马氏链上，从状态 i 以q(i,j) 的概率转跳转到状态j 的时候，我们以α(i,j)的概率接受这个转移，于是得到新的马氏链Q′的转移概率为q(i,j)，α(i,j)。一条马氏链可以看做是一个状态量在不断转移状态的过程，这里接受率可以理解为转移的过程改变成，每次转移时，先根据转移概率确定要转移到那个状态，然后再根据接受概率决定是否要发生转移。

把整个采样的过程整理出来就是

MCMC采样算法

已知分布p(x)，要获得满足p(x)分布的样本

下面过程中，$X_i$表示第i时刻的状态值

1.随机初始化状态量。取$X _ 0$为p(x)分布中x可取到的一个随机值。

2.构造转移矩阵Q。

3.对之后的时刻，循环以下过程采样

第t时刻状态量的值为$X _ t$,采样$y\sim q(x|X_t)$

从均匀分布采样u~Uniform[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$

这里要多做一些解释，帮助大家理解。

1.这个方法既可以用于离散的情况，也可以用于连续的情况。

2.注意平稳细致条件保证的是样本会服从特定分布，而不是马氏链收敛的判断依据。根据马氏链收敛定理，样本点最终一定会收敛，但具体在哪一步收敛，很难判断。通常的做法是，取一个较大的步长n，从第n+1步才开始记录样本的状态。这里的n就叫做Burnin time。

3.如何构造转移矩阵Q？只要保证任意两个状态都是联通的就行。这里联通不是指任意两状态都可以一步到达，只要没有孤立的状态就行。最简单的，我们可以令所有的$q(i,j)$相等,也就是任意两个状态之间的转移概率相同。

4.采样$y\sim q(x|X_t)$是如何实现的？这里实际上就是当已知t时刻状态$X_t$，根据$X_t$转移到其他状态的概率做一个采样。比如何p(x)是之前例子中的$\pi=[0.286,0.489,0.225]$，转移矩阵Q满足所有$q(i,j)$相等，那么就是按转移到状态任意状态的概率都是1/3。在计算机中我们要做这样一个采样，通常会这样做：按照概率累加划定区间，比如三个状态概率都是0.33，那么第一个状态对应的区间是[0,0.33]，第二个状态对应的区间是[0.33,0.66]，第三个状态对应的区间是[0.66,1]。从均匀分布采样u~Uniform[0,1]，u落在那个区间，就取相应区间对应的状态值作为采样的值。

5.为什么要引入接受率的概念？其实本质上，循环采样时，我们直接按$y\sim q(x|X_t)’=q(x|X_t)\alpha(x|X_t)$做采样，就可以得到符合p(x)分布的样本点了。但是$y\sim q(x|X_t)’$采样其实很难做，因为$\alpha(x|X_t)$包含$p(x)$，本质上还是要对p(x)做采样。p(x)只有在x确定的情况下比较好求值，而不好直接采样，不然咱用这方法干嘛是吧？这里，先从q采样再算rejection, 实际上等价于用 rejection sampling 从${q}’(x|X_t)$中直接采样了。

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

观察我们构造的平稳细致等式$p(i)q(i,j)\alpha(i,j)=p(j)q(j,i)\alpha(j,i)$，如果我们等比例的扩大$\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 算法。M-H算法和上述MCMC算法的思路都是一样的，差异就只有接受率而已。

M-H采样算法

已知分布p(x)，要获得满足p(x)分布的样本

下面过程中，$X_i$表示第i时刻的状态值

1.随机初始化状态量。取$X _ 0$为p(x)分布中x可取到的一个随机值。

2.构造转移矩阵Q。

3.对之后的时刻，循环以下过程采样

第t时刻状态量的值为$X _ t$,采样$y\sim q(x|X_t)$

从均匀分布采样u~Uniform[0,1]

如果$u<alpha(x_t,y)=min \left\{ \frac{p(y)q(y,x_t)}{p(x_t)q(x_t,y)},1 \right\}$,则接受转移$X_t \rightarrow y$,即$X_{t+1}=y$

否则，不接受转移，$X_{t+1}=X_t$

OK,我们来看如何用M-H算法对Beta分布做采样。这里稳态分布为Beta分布 $$p(\theta)=Beta(\theta|\alpha,\beta)=\frac{\theta^{\alpha-1}\cdot (1-\theta)^{\beta-1}}{B(\alpha,\beta)}$$，取转移矩阵q为各状态间的转移概率相同，则接受率

$$\begin{align} alpha(i,j) &=min \left\{ \frac{p(j)q(j,i)}{p(i)q(i,j)},1 \right\}\\\\ &=min \left\{ \frac{p(j)}{p(i)},1 \right\}\\\\ &=min \left\{ \frac{j^{\alpha-1}(1-j)^{\beta-1}}{i^{\alpha-1}(1-i)^{\beta-1}},1 \right\}\\\\ \end{align}$$

得到接受率的表达式后，我们就可以按上面提到的步骤开始做采样，具体步骤如下。

1.随机初始化状态量。在Beta分布的取值范围[0,1]之间随机取一个值x_0,具体方法可以是采样x_0~Uniform[0,1]。

2.对之后的时刻，循环以下过程采样

第t时刻状态量的值为$X _ t$,采样$y\sim q(x|X_t)$。由于这里所有的状态间的转移概率都是相同的，实际上等价于在所有取值上做随机一个值，具体做法可以是y~Uniform[0,1]。

从均匀分布采样u~Uniform[0,1]

如果$u<alpha(x_t,y)=min \left\{ \frac{y ^{\alpha-1}(1-y) ^{\beta-1}}{x_t ^{\alpha-1}(1-x_t) ^{\beta-1}},1 \right\}$,则接受转移$X_t \rightarrow y$,即$X_{t+1}=y$

否则，不接受转移，$X_{t+1}=X_t$

经过循环采样后，就可以得到一组服从beta分布的样本点，我们按照0.4为间隔将这些点放到不同的区间，统计区间数量，就可以得到下图所示的结果。无论是求期望还是曲线拟合，都可以通过采样来做近似。采样点越多，结果正确的可能就越大，但这中采样始终是一种近似，不把所有的点遍历，是无法得到完全正确结果的。

¶Gibbs Sampling

M-H算法已经可以很好的工作了，但是在高维（状态量不是一个而是一组）的情况，由于接受率的存在，状态遍历整个分布的效率还不够高。有没有可能找到一个转移矩阵Q，使接受率始终等于1呢？

右侧相同，所以得到 $$p(x_1,y_1)p(y_2|x_1)=p(x_1,y_2)p(y_1|x_1)$$ 即 $$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)$$,有如下等式

由此，我们可以构造任意两点之间转移矩阵Q

$Q(A \rightarrow B)=p(y_B|x_1) \qquad if \quad x_A=x_B=x_1$
$Q(A \rightarrow C)=p(x_C|y_1) \qquad if \quad y_A=y_C=y_1$
$Q(A \rightarrow D)=0$

由此构造的转移矩阵Q就是满足平稳细致条件的，容易验证

$$p(X)Q(X|Y)=p(Y)p(Y|X)$$

于是这个二维空间上的马氏链将收敛到平稳分布 p(x,y)。而这个算法就称为 Gibbs Sampling 算法。整理出来就是

二维Gibbs Sampling算法

已知分布p(x,y)，要获得满足p(x,y)分布的样本

下面过程中，$X_i,Y_i$表示第i时刻的状态值

1.随机初始化状态量$X_0,Y_0$。

2.对之后的时刻，循环以下过程采样

$y_{t+1}~p(y|x_t)$

$x_{t+1}~p(x|y_{t+1})$

太优雅了，比M-H算法看上去要简洁得多，数学家真TM都是天才。

在上面的过程中，马氏链的转移就是沿着坐标轴并且轮换坐标轴转换的。这里每一个时刻最开始是x，或者y都是可以的，轮换知识一种方便的形式。最终便利整个概率分布空间，不同的区域密度反映了分布的概率密度。

这种只能沿着坐标轴转移的方式就是Gibbs Sampling的特性，很容易把它推广到跟高维的情况。在n为空间中，对于概率分布$p(x_1,x_2,\cdots,x_n)$，可以定义如下矩阵

1.如果当前状态是(x_1,x_2,\cdots,x_n)，马氏链转移的过程中，只能沿着坐标轴做转移。沿着$x_i$这根周转移时，转移概率由条件概率p（x_i|(x_1,x_2,\cdots,x_{i-1},x_{i+1},\cdots,x_n）定义。

2.其他无法沿着坐标轴进行跳转的情况，转移概率都是0。

由此，我们就把Gibbs Sampling算法从二维拓展到高维

Gibbs Sampling算法

已知分布$p(\vec{x})$，要获得满足$p(\vec{x})$分布的样本

其中$\vec{x}={X_i;i=1,2,\cdots,n}$，下面过程中，$X_i^t$表示第i轴第t时刻的状态值

1.随机初始化状态量$\vec{x}={X_i;i=1,2,\cdots,n}$。

2.对之后的时刻，循环以下过程采样

$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})$

[1]Blei D M, Ng A Y, Jordan M I. Latent dirichlet allocation[J]. Journal of machine Learning research, 2003, 3(Jan): 993-1022.]