主题模型（4）——LDA模型及其Gibbs Sample求解_k 代表主题数量,m 代表语料中文档的数量, n 代表文档的单词数,nm表示第 m 篇文档

之前关于主题模型整理了《文本建模之Unigram Model，PLSA与LDA》与《再看LDA主题模型》两篇博客，以及针对PLSA的求解整理了博客《主题模型（3）——PLSA模型及其EM算法求解》，这一篇博客将继续整理LDA（Latent Dirichlet Allocation）模型的Gibbs Sample求解方法。

LDA回顾

同样，首先回归下LDA模型的文档生成过程。我们知道，LDA在PLSA的基础上引入了贝叶斯框架，即参数不仅未知，而且取值不固定。 按照贝叶斯学派的思想，我们先给参数指定一个先验分布（prior），然后我们观察到一些样本，进而对参数的分布进行调整，得到参数的后验分布（posterior），参数先验分布与后验分布通过贝叶斯定理连接：
$$p(x|y)=\frac{p(x,y)}{p(y)}=\frac{p(y|x)p(x)}{p(y)}$$

LDA模型将Dirichlet分布作为参数的先验分布，至于原因，是因为Dirichlet分布于多项式分布是一对共轭分布，保证了计算后验概率的便利。在给出LDA模型文档生成概率之前，我们先规定一些记号：

• $V$：语料中词汇构成的字典大小
• $K$：人工定义的主题个数
• $M$：语料中文档数目
• $N_m$：语料中第$m$篇文档的单词数目
• $w_{m,n}$：语料中第$m$篇文档第$n$个单词
• $z_{m,n}$：语料中第$m$篇文档第$n$个单词$w_{m,n}$对应的主题
• ${\vec{\theta }}_m$：第$m$篇文档的主题分布参数，长度为$K$，因此$M$篇文档的主题分布构成了一个$M\ast K$的矩阵，记为$\Theta$，每一行代表一篇文档的主题分布。
• ${\vec{\varphi }}_k$：第$k$个主题的词分布参数，长度为$V$，因此$K$个主题的词分布构成了一个$K\ast V$的矩阵，记为$\Phi$，每一行代表一个主题的词分布。

LDA通过模拟文档生成过程，找出文档最有可能的主题参数。 在生成文档时，首先根据主题参数的先验分布随机选取一个文档主题分布${\vec{\theta }}_m$与主题词分布$\Phi$，然后从主题分布下生成主题，从主题对应的词分布下生成单词，因此第$m$篇文档第$n$个单词生成概率是
$$p(w_{m,n},z_{m,n},{\vec{\theta }}_m,{\vec{\varphi }}_{z_{m,n}}|\vec{\alpha },\vec{\beta })=p(w_{m,n}|{\vec{\varphi }}_{z_{m,n}})p(z_{m,n}|{\vec{\theta }}_m)p({\vec{\theta }}_m|\vec{\alpha })p({\vec{\varphi }}_{z_{m,n}}|\vec{\beta })\text{\hspace{1em}(1)}$$

上式中，${\vec{\theta }}_m$、$z_{m,n}$和${\vec{\varphi }}_{z_{m,n}}$都是隐变量，将${\vec{\theta }}_m$、$z_{m,n}$和${\vec{\varphi }}_{z_{m,n}}$积分掉，将得到全概率$p(w_{m,n}|\vec{\alpha },\vec{\beta })$
$$p(w_{m,n}|\vec{\alpha },\vec{\beta })=\sum _{z_{m,n}}\int \int p(w_{m,n}|{\vec{\varphi }}_{z_{m,n}})p(z_{m,n}|{\vec{\theta }}_m)p({\vec{\theta }}_m|\vec{\alpha })p({\vec{\varphi }}_{z_{m,n}}|\vec{\beta })d{\vec{\varphi }}_{z_{m,n}}d{\vec{\theta }}_m\text{\hspace{1em}(2)}$$

重复上述单词生成过程，生成一篇文档；重复文档生成过程，产生整个文档集，文档集的生成概率为：
(3) p ( D ) = ∏ m = 1 M ∏ n = 1 N m p ( w m , n ∣ α ⃗ , β ⃗ ) = ∏ m = 1 M ∏ n = 1 N m ∑ z m , n ∫ ∫ p ( w m , n ∣ ϕ ⃗ z m , n ) p ( z m , n ∣ θ ⃗ m ) p ( θ ⃗ m ∣ α ⃗ ) p ( ϕ ⃗ z m , n ∣ β ⃗ ) d ϕ ⃗ z m , n d θ ⃗ m

$$\begin{aligned} p(D)&=\prod_{m=1}^{M}\prod_{n=1}^{N_m}p(w_{m,n}|\vec{\alpha},\vec{\beta})\\ &= \prod_{m=1}^{M}\prod_{n=1}^{N_m}\sum_{z_{m,n}}\int\int p(w_{m,n}|\vec{\phi}_{z_{m,n}})p(z_{m,n}|\vec\theta_m)p(\vec{\theta}_m|\vec{\alpha})p(\vec{\phi}_{z_{m,n}}|\vec{\beta})d\vec{\phi}_{z_{m,n}}d\vec{\theta}_m \tag 3 \end{aligned}$$ p(D)​=m=1∏M​n=1∏Nm​​p(wm,n​∣α ,β ​)=m=1∏M​n=1∏Nm​​zm,n​∑​∫∫p(wm,n​∣ϕ ​zm,n​​)p(zm,n​∣θ m​)p(θ m​∣α )p(ϕ ​zm,n​​∣β ​)dϕ ​zm,n​​dθ m​​(3)

按照常理，最大化$p(D)$，我们就能求出主题参数，但是，上式包含无数隐变量的累加，即使给定一组主题参数，直接求出$p(D)$都几乎不可能。好在计算机科学家们想到了通过随机模拟的方法（又称蒙特卡洛方法，Monte Carlo Simulation），生成分布的无数样本，用频率近似代替概率的方式解决概率求解的问题，Gibbs Sample便是其中的代表方法。

Gibbs Sample算法

采样是从特定概率分布中抽取样本的过程，采样方法用于获得服从指定概率分布的样本。 因此要采样，我们就必须知道概率分布，但是，对于一个概率分布$p(\vec{x})$，其分布可能非常复杂，无法直接求出其概率值，将导致采样遇到困难。此时就需要Gibbs Sample等更加复杂的能够work的采样方法。那Gibbs Sample是怎么操作的呢？

Gibbs Sample是马尔可夫链蒙特卡罗理论中用来近似求解多维概率分布（2个或多个随机变量的联合分布）的算法。Gibbs Sample形式化描述如下：对于一个复杂的概率分布$p(\vec{x})$，我们没法求出其对应的概率分布，但如果我们能够求出在给定其他分量$x_{\neg i}$情况下$x_i$的概率$p(x_i|x_{\neg i})$，那么我们可以根据$p(x_i|x_{\neg i})$采样每一个分量$x_i$组成一个样本$\vec{x}$，重复上述采样过程$T$次，最后统计频率作为概率。 执行过程如下：

• 随机初始化 x ⃗ 0 = { x 1 0 , x 2 0 , ⋯   , x N 0 } \vec{x}^0=\{x_1^0,x_2^0,\cdots,x_N^0\} x 0={x10​,x20​,⋯,xN0​}
• for $t=1,2,3,\cdots ,T:$
  • $x_1^t\sim p(x_1|x_2^{t-1},x_3^{t-1},\cdots ,x_N^{t-1})$
  • $x_2^t\sim p(x_2|x_1^t,x_3^{t-1},\cdots ,x_N^{t-1})$
  • $x_3^t\sim p(x_2|x_1^t,x_2^t,\cdots ,x_N^{t-1})$
  • $\cdots$
  • $x_N^t\sim p(x_2|x_1^t,x_2^t,\cdots ,x_{N-1}^t)$
  • 得到一个样本 x ⃗ t = { x 1 t , x 2 t , ⋯   , x N t } \vec{x}^t=\{x_1^t,x_2^t,\cdots,x_N^t\} x t={x1t​,x2t​,⋯,xNt​}
• 直至采样过程收敛，得到的样本$[{\vec{x}}^{n+1},{\vec{x}}^{n+2},\cdots ,{\vec{x}}^T]$就是真实的$p(\vec{x})$样本。

马尔可夫链及其性质

在引出Gibbs Sample前，需要先简要介绍下马尔可夫链及其平稳分布，马尔可夫链数学定义如下：
$$p(X_{t+1}|X_t,X_{t-1},\cdots )=p(X_{t+1}|X_t)$$

即当前状态只受到前一个状态的影响，而与其他状态无关。通常状态之间的转移概率用转移矩阵描述，如下图，展示了状态1，2，3之间的概率转移

其对应的状态转移矩阵如下：
P = [ 0.65 0.28 0.07 0.15 0.67 0.18 0.12 0.36 0.52 ] P=

$$\begin{bmatrix} 0.65 & 0.28 & 0.07 \\ 0.15 & 0.67 & 0.18 \\ 0.12 & 0.36 & 0.52 \\ \end{bmatrix}$$ P=⎣⎡​0.650.150.12​0.280.670.36​0.070.180.52​⎦⎤​

对于马尔可夫链状态转移矩阵，存在马氏链定理： 如果一个非周期马氏链具有状态转移矩阵$P$，且它的任意两个状态是联通的（从任意一个状态能到其他任意状态），那么$\underset{n\to \infty }{\lim }{P}_{ij}^n$存在且与$i$无关，记$\underset{n\to \infty }{\lim }{P}_{ij}^n=\pi (j)$，我们有

1. $\underset{n\to \infty }{\lim }{P}^n=\left[\begin{array}{ccccc}\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 \end{array}\right]$
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\pi (i)=1$$

称为马氏链的平稳分布。

上面1，2，3分别说明

1. 状态转移矩阵经过多次转移后，会达到稳定。
2. 马氏链稳定后，转移到任意状态的概率是稳定的。
3. 一个状态转移矩阵只对应唯一一个平稳分布。

如果我们从一个具体地初始状态$x_0$出发，时刻1所处状态$x_1$满足概率分布$P_i$（$P_i$是状态转移矩阵中状态$x_0$对应的转移概率），将其记为${\pi }_1(x)$。时刻2所处状态$x_2$满足概率分布${\pi }_1(x)P$，将其记为${\pi }_2(x)$。以此类推，时刻n所处状态$x_n$满足概率分布${\pi }_{n-1}(x)P={\pi }_1(x)P^{n-1}$，即
x 0 x 1 ∼ π 1 ( x ) = P i ⋯ x n ∼ π n ( x ) = π n − 1 ( x ) P = P i P n − 1 x n + 1 ∼ π n + 1 ( x ) = π n ( x ) P = P i P n x n + 2 ∼ π n + 2 ( x ) = π n + 1 ( x ) P = P i P n + 1 ⋯

$$\begin{aligned} x_0&\\ x_1&\thicksim \pi_1(x)=P_i\\ \cdots \\ x_n&\thicksim \pi_n(x)=\pi_{n-1}(x)P=P_iP^{n-1} \\ x_{n+1}&\thicksim \pi_{n+1}(x)=\pi_{n}(x)P=P_iP^{n} \\ x_{n+2}&\thicksim \pi_{n+2}(x)=\pi_{n+1}(x)P=P_iP^{n+1} \\ \cdots \end{aligned}$$ x0​x1​⋯xn​xn+1​xn+2​⋯​∼π1​(x)=Pi​∼πn​(x)=πn−1​(x)P=Pi​Pn−1∼πn+1​(x)=πn​(x)P=Pi​Pn∼πn+2​(x)=πn+1​(x)P=Pi​Pn+1​

根据马氏链的收敛定理，假设$n$时刻马氏链已经收敛到平稳分布$\pi (x)$，那么$x_{n+1},x_{n+2},\cdots$都将是平稳分布$\pi (x)$的样本。这不正是我们想要实现的从某个概率分布中采样的过程吗？

Markov Chain Monte Carlo

根据上面的推导，对于一个概率分布$p(x)$，如果我们能够构造一个转移矩阵为$P$的马氏链，并且该马氏链的平稳分布恰好是$p(x)$，那么当马氏链收敛后，我们就得到了概率分布$p(x)$的样本。 因为上述方法属于随机模拟方法的扩展，并且用到了马尔可夫链，因此这类方法被称为马尔可夫链蒙特卡罗方法（Markov Chain Monte Carlo，MCMC）。

基于马氏链做采样的关键是如何构造概率分布$p(x)$对应的状态转移矩阵$P$，直接找到这个矩阵$P$通常很难做到，需要借助细致平稳条件（detailed balance condition）实现。

细致平稳条件的定义如下：如果非周期马氏链的转移矩阵$P$和概率分布$p(x)$满足
$$p(i)P(i,j)=p(j)P(j,i)\text{for all}i,j$$

则概率分布$p(x)$是马氏链转移矩阵$P$的平稳分布。

也就是说我们只需要找到使概率分布满足细致平稳条件的矩阵$P$。通常情况下，对于一个状态转移矩阵$Q$，$p(i)q(i,j)/=p(j)q(j,i)$，因此我们需要对转移矩阵$Q$进行改造，使得$p(i)q^{\prime }(i,j)=p(j)q^{\prime }(j,i)$。最直观的方法是引入一个$\alpha (i,j)$，使得
$$p(i)q(i,j)\alpha (i,j)=p(j)q(j,i)\alpha (j,i)$$

按照对称性，我们可以取
$$\begin{gathered} \alpha (i,j)=p(j)q(j,i) \\ \alpha (j,i)=p(i)q(i,j) \end{gathered}$$

此时，新的转移矩阵$q^{\prime }(i,j)=q(i,j)\alpha (i,j)$，即$Q^{\prime }=Q\odot \alpha$，此处是矩阵的Hadamard乘积。直接从$Q^{\prime }$中采样还是难以实现，可以将$\alpha$看作接受率，应用接受-拒绝采样得到分布的一系列样本。

Gibbs Sample原理

分析上面的MCMC采样，由于接受率$\alpha$的存在，导致算法效率不高；另一方面，更为致命的是，很多时候我们很难求出多维联合概率分布$p(\vec{x})$，导致上面的算法无法work，这时就需要使用Gibbs Sample。

考虑由二维特征$(x,y)$描述的状态，假设其概率分布为$p(x,y)$，现在固定其中的$x$，对于两个状态$A(x_1,y_1)$和$B(x_1,y_2)$，根据概率公式：
$$\begin{gathered} 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{gathered}$$

可得
$$p(x_1,y_1)p(y_2|x_1)=p(x_1,y_2)p(y_1|x_1)\text{\hspace{1em}(4)}$$

其中，因为$A,B$两状态的$x$特征相同，因此从$A$状态转移到$B$状态等价于在给定$x_1$时$y_2$的取值概率，即$p(y_2|x_1)$，所以$p(y_2|x_1)$表示了转移概率$p(A\to B)$，$p(y_1|x_1)$表示了转移概率$p(B\to A)$，所以上面等式（4）表达了
$$p(A)p(A\to B)=p(B)p(B\to A)\text{\hspace{1em}(5)}$$

同理，固定特征$y$，我们也会得到同样的等式，这表明在二维概率分布中，固定其中一维，使用条件分布$p(y|x)$作为状态之间的转移概率，那么任意两个状态之间的转移满足细致平稳条件。即构造如下状态转移矩阵
Q ( A → B ) = p ( y B ∣ x 1 ) if x A = x B = x 1 Q ( A → C ) = p ( x C ∣ y 1 ) if y A = y B = y 1 Q ( A → D ) = 0 else

$$\begin{aligned} Q(A\to B)=&p(y_B|x_1) \qquad \text{if}\enspace x_A=x_B=x_1\\ Q(A\to C)=&p(x_C|y_1) \qquad \text{if}\enspace y_A=y_B=y_1\\ Q(A\to D)=&0\qquad\qquad\qquad\qquad\qquad\enspace \text{else} \end{aligned}$$ Q(A→B)=Q(A→C)=Q(A→D)=​p(yB​∣x1​)ifxA​=xB​=x1​p(xC​∣y1​)ifyA​=yB​=y1​0else​

如下图所示：

根据上面的状态转移矩阵，我们可以很容易的验证任意两个状态$X,Y$，满足细致平稳条件：
$$p(X)Q(X\to Y)=p(Y)Q(y\to X)$$

因此，在二维平面中，马氏链轮换的沿着$x$轴和$y$轴转移，得到样本$(x_0,y_0),(x_1,y_1),\cdots$，等到马氏链收敛后，最终得到的样本就是概率分布$p(x,y)$的样本，上述过程就是二维Gibbs Sample算法。推广到多维分布，我们就得到了Gibbs Sample的一般形式。

LDA模型的Gibbs Sample求解步骤

回到我们的LDA模型，我们已经写出文档集的生成概率为：
p ( D ) = ∏ m = 1 M ∏ n = 1 N m p ( w m , n ∣ α ⃗ , β ⃗ ) = ∏ m = 1 M ∏ n = 1 N m ∑ z m , n ∫ ∫ p ( w m , n ∣ ϕ ⃗ z m , n ) p ( z m , n ∣ θ ⃗ m ) p ( θ ⃗ m ∣ α ⃗ ) p ( ϕ ⃗ z m , n ∣ β ⃗ ) d ϕ ⃗ z m , n d θ ⃗ m

$$\begin{aligned} p(D)&=\prod_{m=1}^{M}\prod_{n=1}^{N_m}p(w_{m,n}|\vec{\alpha},\vec{\beta})\\ &= \prod_{m=1}^{M}\prod_{n=1}^{N_m}\sum_{z_{m,n}}\int\int p(w_{m,n}|\vec{\phi}_{z_{m,n}})p(z_{m,n}|\vec\theta_m)p(\vec{\theta}_m|\vec{\alpha})p(\vec{\phi}_{z_{m,n}}|\vec{\beta})d\vec{\phi}_{z_{m,n}}d\vec{\theta}_m \end{aligned}$$ p(D)​=m=1∏M​n=1∏Nm​​p(wm,n​∣α ,β ​)=m=1∏M​n=1∏Nm​​zm,n​∑​∫∫p(wm,n​∣ϕ ​zm,n​​)p(zm,n​∣θ m​)p(θ m​∣α )p(ϕ ​zm,n​​∣β ​)dϕ ​zm,n​​dθ m​​

这个概率直接求解几乎不可能。遇到困难，再看一下目标，能保证我们不迷失方向。

单词主题联合概率

LDA模型的目标是：找出文档的主题参数，也就是文档中每个单词对应的主题，那我们是不是可以对上面的$p(D)$放松一下，求出$p(\vec{w},\vec{z})$的概率呢？
p ( w ⃗ , z ⃗ ) = p ( w ⃗ ∣ z ⃗ , β ⃗ ) p ( z ⃗ ∣ α ⃗ )

$$\begin{aligned} p(\vec{w},\vec{z})&=p(\vec{w}|\vec{z},\vec{\beta})p(\vec{z}|\vec{\alpha})\\ \end{aligned}$$ p(w ,z )​=p(w ∣z ,β ​)p(z ∣α )​

上式在博客《再看LDA主题模型》中已经分析过，其包含两个独立的过程：先生成文档每个位置对应的主题，然后根据主题生成每个位置的单词。

文档$m$的主题生成概率$p(\overrightarrow{z_m}|\vec{\alpha })$计算如下：
p ( z m ⃗ ∣ α ⃗ ) = ∫ p ( z m ⃗ ∣ θ m ⃗ ) p ( θ m ⃗ ∣ α ⃗ ) d θ m ⃗ = ∫ p ( z m ⃗ ∣ θ m ⃗ ) D i r ( θ m ⃗ ∣ α ⃗ ) d θ m ⃗ = ∫ ∏ k = 1 K θ k , m n k , m 1 Δ ( α ⃗ ) ∏ k = 1 K θ k , m α k , m − 1 d θ m ⃗ = 1 Δ ( α ⃗ ) ∫ ∏ k = 1 K θ k , m n k , m + α k − 1 d θ m ⃗ = Δ ( n m ⃗ + α ⃗ ) Δ ( α ⃗ )

$$\begin{aligned} p(\vec{z_m}|\vec{\alpha}) =& \int p(\vec{z_m}|\vec{\theta_m})p(\vec{\theta_m}|\vec{\alpha})d\vec{\theta_m} \\ =& \int p(\vec{z_m}|\vec{\theta_m})Dir(\vec{\theta_m}|\vec{\alpha})d\vec{\theta_m} \\ =& \int \prod_{k=1}^K\theta_{k,m}^{n_{k,m}}\frac{1}{\Delta(\vec{\alpha})}\prod_{k=1}^K\theta_{k,m}^{\alpha_{k,m}-1}d\vec{\theta_m} \\ =& \frac{1}{\Delta(\vec\alpha)}\int\prod_{k=1}^K\theta_{k,m}^{n_{k,m}+\alpha_{k}-1}d\vec{\theta_m} \\ =& \frac{\Delta(\vec{n_m}+\vec{\alpha})}{\Delta(\vec{\alpha})} \end{aligned}$$ p(zm​ ​∣α )=====​∫p(zm​ ​∣θm​ ​)p(θm​ ​∣α )dθm​ ​∫p(zm​ ​∣θm​ ​)Dir(θm​ ​∣α )dθm​ ​∫k=1∏K​θk,mnk,m​​Δ(α )1​k=1∏K​θk,mαk,m​−1​dθm​ ​Δ(α )1​∫k=1∏K​θk,mnk,m​+αk​−1​dθm​ ​Δ(α )Δ(nm​ ​+α )​​

主题$k$下单词生成概率$p(\overrightarrow{w_k}|z_k,\vec{\beta })$计算如下：
p ( w k ⃗ ∣ z k , β ⃗ ) = ∫ p ( w k ⃗ ∣ φ k ⃗ ) p ( φ k ⃗ ∣ β ⃗ ) d φ k ⃗ = ∫ p ( w k ⃗ ∣ φ k ⃗ ) D i r ( φ k ⃗ ∣ z k , β ⃗ ) d φ k ⃗ = ∫ ∏ v = 1 V φ v , k n v , k 1 Δ ( β ⃗ ) ∏ v = 1 V φ v , k n v , k d φ k ⃗ = 1 Δ ( β ⃗ ) ∫ ∏ v = 1 V φ v , k n v , k + β v − 1 d φ k ⃗ = Δ ( n k ⃗ + β ⃗ ) Δ ( β ⃗ )

$$\begin{aligned} p(\vec{w_k}|z_k,\vec{\beta}) =& \int p(\vec{w_k}|\vec{\varphi_k})p(\vec{\varphi_k}|\vec{\beta})d\vec{\varphi_k} \\ =& \int p(\vec{w_k}|\vec{\varphi_k})Dir(\vec{\varphi_k}|z_k,\vec{\beta})d\vec{\varphi_k} \\ =& \int \prod_{v=1}^V\varphi_{v,k}^{n_{v,k}}\frac{1}{\Delta(\vec{\beta})}\prod_{v=1}^V\varphi_{v,k}^{n_{v,k}}d\vec{\varphi_k} \\ =& \frac{1}{\Delta(\vec{\beta})}\int\prod_{v=1}^V\varphi_{v,k}^{n_{v,k}+\beta_{v}-1}d\vec{\varphi_k}\\ =& \frac{\Delta(\vec{n_k}+\vec{\beta})}{\Delta(\vec{\beta})} \end{aligned}$$ p(wk​ ​∣zk​,β ​)=====​∫p(wk​ ​∣φk​ ​)p(φk​ ​∣β ​)dφk​ ​∫p(wk​ ​∣φk​ ​)Dir(φk​ ​∣zk​,β ​)dφk​ ​∫v=1∏V​φv,knv,k​​Δ(β ​)1​v=1∏V​φv,knv,k​​dφk​ ​Δ(β ​)1​∫v=1∏V​φv,knv,k​+βv​−1​dφk​ ​Δ(β ​)Δ(nk​ ​+β ​)​​

因此
p ( w ⃗ , z ⃗ ) = p ( w ⃗ ∣ z ⃗ , β ⃗ ) p ( z ⃗ ∣ α ⃗ ) = ∏ k = 1 K p ( w ⃗ ∣ z ⃗ , β ⃗ ) ∏ m = 1 M p ( z ⃗ ∣ α ⃗ ) = ∏ k = 1 K Δ ( n k ⃗ + β ⃗ ) Δ ( β ⃗ ) ∏ m = 1 M Δ ( n m ⃗ + α ⃗ ) Δ ( α ⃗ )

$$\begin{aligned} p(\vec{w},\vec{z}) =& p(\vec{w}|\vec{z},\vec{\beta})p(\vec{z}|\vec{\alpha}) \\ =& \prod_{k=1}^K p(\vec{w}|\vec{z},\vec{\beta})\prod_{m=1}^Mp(\vec{z}|\vec{\alpha})\\ =& \prod_{k=1}^K\frac{\Delta(\vec{n_k}+\vec{\beta})}{\Delta(\vec{\beta})}\prod_{m=1}^M\frac{\Delta(\vec{n_m}+\vec{\alpha})}{\Delta(\vec{\alpha})} \end{aligned}$$ p(w ,z )===​p(w ∣z ,β ​)p(z ∣α )k=1∏K​p(w ∣z ,β ​)m=1∏M​p(z ∣α )k=1∏K​Δ(β ​)Δ(nk​ ​+β ​)​m=1∏M​Δ(α )Δ(nm​ ​+α )​​

主题条件概率

有了联合概率，下面就需要求解条件概率$p(z_i=k|\vec{w},{\vec{z}}_{\neg i})$，这里$i$是指第$m$篇文档的第$n$个词。因为词$w_i$是可以观察到的，因此我们有：
p ( z i = k ∣ w ⃗ , z ⃗ ¬ i ) ∝ p ( z i = k , w i = t ∣ w ⃗ ¬ i , z ⃗ ¬ i ) = ∫ ∫ p ( z i = k , w i = t , θ ⃗ m , φ ⃗ k ∣ w ⃗ ¬ i , z ⃗ ¬ i ) d θ ⃗ m d φ ⃗ k = ∫ ∫ p ( z i = k , θ ⃗ m ∣ w ⃗ ¬ i , z ⃗ ¬ i ) p ( w i = t , φ ⃗ k ∣ w ⃗ ¬ i , z ⃗ ¬ i ) d θ ⃗ m d φ ⃗ k = ∫ p ( z i = k , θ ⃗ m ∣ w ⃗ ¬ i , z ⃗ ¬ i ) d θ ⃗ m ∗ ∫ p ( w i = t , φ ⃗ k ∣ w ⃗ ¬ i , z ⃗ ¬ i ) d φ ⃗ k = ∫ p ( z i = k ∣ θ ⃗ m ) p ( θ ⃗ m ∣ w ⃗ ¬ i , z ⃗ ¬ i ) d θ ⃗ m ∗ ∫ p ( w i = t ∣ φ ⃗ k ) p ( φ ⃗ k ∣ w ⃗ ¬ i , z ⃗ ¬ i ) d φ ⃗ k = ∫ p ( z i = k ∣ θ ⃗ m ) D i r ( θ ⃗ m ∣ n ⃗ m , ¬ i + α ⃗ ) d θ ⃗ m ∗ ∫ p ( w i = t ∣ φ ⃗ k ) D i r ( φ ⃗ k ∣ n ⃗ k , ¬ i + β ⃗ ) d φ ⃗ k = ∫ θ m , k D i r ( θ ⃗ m ∣ n ⃗ m , ¬ i + α ⃗ ) d θ ⃗ m ∗ ∫ φ k , t D i r ( φ ⃗ k ∣ n ⃗ k , ¬ i + β ⃗ ) d φ ⃗ k = E ( θ m , k ) ∗ E ( φ k , t )

$$\begin{aligned} p(z_i=k|\vec{w},\vec{z}_{¬i})&\propto p(z_i=k,w_i=t|\vec{w}_{¬i},\vec{z}_{¬i})\\ &=\int\int p(z_i=k,w_i=t,\vec{\theta}_m,\vec{\varphi}_k|\vec{w}_{¬i},\vec{z}_{¬i})d\vec{\theta}_md\vec{\varphi}_k\\ &=\int\int p(z_i=k,\vec{\theta}_m|\vec{w}_{¬i},\vec{z}_{¬i})p(w_i=t,\vec{\varphi}_k|\vec{w}_{¬i},\vec{z}_{¬i})d\vec{\theta}_md\vec{\varphi}_k\\ &=\int p(z_i=k,\vec{\theta}_m|\vec{w}_{¬i},\vec{z}_{¬i})d\vec{\theta}_m*\int p(w_i=t,\vec{\varphi}_k|\vec{w}_{¬i},\vec{z}_{¬i})d\vec{\varphi}_k\\ &=\int p(z_i=k|\vec{\theta}_m)p(\vec{\theta}_m|\vec{w}_{¬i},\vec{z}_{¬i})d\vec{\theta}_m*\int p(w_i=t|\vec{\varphi}_k)p(\vec{\varphi}_k|\vec{w}_{¬i},\vec{z}_{¬i})d\vec{\varphi}_k\\ &=\int p(z_i=k|\vec{\theta}_m)Dir(\vec{\theta}_m|\vec{n}_{m,¬i}+\vec{\alpha})d\vec{\theta}_m*\int p(w_i=t|\vec{\varphi}_k)Dir(\vec{\varphi}_k|\vec{n}_{k,¬i}+\vec{\beta})d\vec{\varphi}_k\\ &=\int \theta_{m,k}Dir(\vec{\theta}_m|\vec{n}_{m,¬i}+\vec{\alpha})d\vec{\theta}_m*\int \varphi_{k,t}Dir(\vec{\varphi}_k|\vec{n}_{k,¬i}+\vec{\beta})d\vec{\varphi}_k\\ &=E(\theta_{m,k})*E(\varphi_{k,t}) \end{aligned}$$ p(zi​=k∣w ,z ¬i​)​∝p(zi​=k,wi​=t∣w ¬i​,z ¬i​)=∫∫p(zi​=k,wi​=t,θ m​,φ ​k​∣w ¬i​,z ¬i​)dθ m​dφ ​k​=∫∫p(zi​=k,θ m​∣w ¬i​,z ¬i​)p(wi​=t,φ ​k​∣w ¬i​,z ¬i​)dθ m​dφ ​k​=∫p(zi​=k,θ m​∣w ¬i​,z ¬i​)dθ m​∗∫p(wi​=t,φ ​k​∣w ¬i​,z ¬i​)dφ ​k​=∫p(zi​=k∣θ m​)p(θ m​∣w ¬i​,z ¬i​)dθ m​∗∫p(wi​=t∣φ ​k​)p(φ ​k​∣w ¬i​,z ¬i​)dφ ​k​=∫p(zi​=k∣θ m​)Dir(θ m​∣n m,¬i​+α )dθ m​∗∫p(wi​=t∣φ ​k​)Dir(φ ​k​∣n k,¬i​+β ​)dφ ​k​=∫θm,k​Dir(θ m​∣n m,¬i​+α )dθ m​∗∫φk,t​Dir(φ ​k​∣n k,¬i​+β ​)dφ ​k​=E(θm,k​)∗E(φk,t​)​

上式正是Dirichlet分布的期望，对于其中任意一项$i$，我们有：
E ( p k ) = ∫ 0 1 p k ∗ Γ ( ∑ i = 1 K α i ) ∏ i = 1 K Γ ( α i ) ∏ i = 1 K p i α i − 1 d p ⃗ = Γ ( ∑ i = 1 K α i ) ∏ i = 1 K Γ ( α i ) ∫ 0 1 ∏ i = 1 k − 1 p i α i − 1 ∗ p k α k ∗ ∏ i = k + 1 K p i α i − 1 d p ⃗ = Γ ( ∑ i = 1 K α i ) ∏ i = 1 K Γ ( α i ) ∗ ∏ i = 1 k − 1 Γ ( α i ) ∗ Γ ( α k + 1 ) ∗ ∏ i = k + 1 K Γ ( α i ) Γ ( ∑ i = 1 k − 1 α i + ( α k + 1 ) + ∑ i = k + 1 K α i ) = α k ∑ k = 1 K α k

$$\begin{aligned} E(p_k)&=\int_0^1p_k*\frac{\varGamma(\sum_{i=1}^K\alpha_i)}{\prod_{i=1}^K\varGamma(\alpha_i)}\prod_{i=1}^Kp_i^{\alpha_i-1}d\vec{p}\\ &=\frac{\varGamma(\sum_{i=1}^K\alpha_i)}{\prod_{i=1}^K\varGamma(\alpha_i)}\int_0^1\prod_{i=1}^{k-1}p_i^{\alpha_i-1}*p_k^{\alpha_k}*\prod_{i=k+1}^{K}p_i^{\alpha_i-1}d\vec{p}\\ &=\frac{\varGamma(\sum_{i=1}^K\alpha_i)}{\prod_{i=1}^K\varGamma(\alpha_i)}*\frac{\prod_{i=1}^{k-1}\varGamma(\alpha_i)*\varGamma(\alpha_k+1)*\prod_{i=k+1}^{K}\varGamma(\alpha_i)}{\varGamma(\sum_{i=1}^{k-1}\alpha_i+(\alpha_k+1)+\sum_{i=k+1}^K\alpha_i)}\\ &=\frac{\alpha_k}{\sum_{k=1}^K\alpha_k} \end{aligned}$$ E(pk​)​=∫01​pk​∗∏i=1K​Γ(αi​)Γ(∑i=1K​αi​)​i=1∏K​piαi​−1​dp ​=∏i=1K​Γ(αi​)Γ(∑i=1K​αi​)​∫01​i=1∏k−1​piαi​−1​∗pkαk​​∗i=k+1∏K​piαi​−1​dp ​=∏i=1K​Γ(αi​)Γ(∑i=1K​αi​)​∗Γ(∑i=1k−1​αi​+(αk​+1)+∑i=k+1K​αi​)∏i=1k−1​Γ(αi​)∗Γ(αk​+1)∗∏i=k+1K​Γ(αi​)​=∑k=1K​αk​αk​​​

因此
$$p(z_i=k|\vec{w},{\vec{z}}_{\neg i})\propto E({\theta }_{m,k})\ast E({\phi }_{k,t})=\frac{n_{m,\neg i}^{(k)}+{\alpha }_k}{{\sum }_{k=1}^K(n_{m,\neg i}^{(k)}+{\alpha }_k)}\ast \frac{n_{k,\neg i}^{(t)}+{\beta }_t}{{\sum }_{t=1}^V(n_{k,\neg i}^{(t)}+{\beta }_t)}$$

有了上式，我们就可以借助Gibbs Sample算法对概率分布进行采样，待算法收敛后，采样出的样本就是概率分布$p(\vec{w},\vec{z})$的样本。通过采样我们得到了每个词的主题，统计各个主题下单词出现频率，就可以得到各个主题的词分布。同样的，统计每篇文档的主题频率，我们就得到了各个文档的主题分布。通常，我们取Gibbs Sample收敛后的$n$个样本进行平均做参数估计。

总结

LDA参数估计是Gibbs Sample的一个非常好的应用实例，本文主要总结了Gibbs Sample的工作原理以及LDA模型的Gibbs Sample求解步骤。
关于采样：

1. 采样可以用来计算一些复杂的代数式子的值，这些式子可能在机器学习的优化过程中遇到，比如一些很复杂的积分，级数，没有办法去直接求解，那么Gibbs采样后可以用采样样本去近似求值。
2. 对于简单的概率分布，如果可以计算机程序实现，就直接采样。
3. 如果概率分布的形式比较复杂或者概率分布根本不知道，我们就需要借助接受-拒绝采样，MCMC等算法，避开直接对原始分布进行采样。
4. 在使用MCMC方法时，如果概率分布已知，一般可以考虑M-H采样（M-H采样必须知道要采样的概率分布）。但是如果概率分布未知，或者特征维度特别大时，M-H采样就不那么work了。
5. 在M-H采样中，即使我们构造出使得概率分布满足平稳分布的转移矩阵$P$，$p(i,j)=q(i,j)\pi (j)q(j,i)$，因为$\pi (j)$的存在导致我们直接从P中采样无法程序实现，所以我们需要用接受-拒绝采样。
6. 即使不知道联合概率，只知道特征之间转换的条件概率，我们也可以采用Gibbs Sample获得分布的样本。

个人也是最近才开始深入了采样方法，欢迎指正，多多交流。

参考文献

MCMC(一)蒙特卡罗方法
一文了解采样方法
深度学习中的采样以及采样算法
各领域中采样方式研究 (持续更新)
从随机过程到马尔科夫链蒙特卡洛方法
MCMC(三)MCMC采样和M-H采样
MCMC(四)Gibbs采样
文本主题模型之LDA(二) LDA求解之Gibbs采样算法
LDA数学八卦