【机器学习算法推导】隐狄利克雷模型LDA

1.统计学基础

1.1 二项分布

  二项分布是n重伯努利分布，可以看做是将硬币抛n次，出现k次正面向上的概率，每次出现正面向上的概率为p，其概率密度公式为
$$P(X=k)=C_n^k{p}^k(1-p{)}^{n-k}$$

1.2 多项式分布

  多项式分布是二项分布推广到多种结果的情况，在多项式分布下，抛的不是硬币，而是一个骰子，或者是其他多面体。对于一个骰子，每一面出现的结果是$\frac{1}{6}$，每一面出现的概率互斥且和为1，发生其中一个结果X次的概率就是多项式分布，假设投N次骰子，有$x_1,x_2,x_3,x_4,x_5,x_6$六种情况，令$m_1,m_2,m_3,m_4,m_5,m_6$分别表示每种情况出现的次数，且$m_1+m_2+m_3+m_4+m_5+m_6=N$，则其概率为
$$P(X_1=m_1,X_2=m_2,...,X_n=m_n)=\frac{N!}{m_1!m_2!...m_n!}{p}_1^{m_1}{p}_2^{m_2}...p_n^{m_n}$$

1.3 Gamma分布

  $\Gamma$函数的形式如下：
$$\Gamma (x)={\int }_0^{\infty }{t}^{x-1}{e}^{-t}\mathrm{d}t$$
  该函数具有如下的递推性质：
$$\Gamma (x+1)=x\Gamma (x)$$
  容易知道$\Gamma$函数为阶乘在实数集上的延伸，有
$$\Gamma (n)=(n-1)!$$
  对$\Gamma$函数做一下变形，得到：
$${\int }_0^{\infty }\frac{x^{\alpha -1}{e}^{-x}}{\Gamma (\alpha )}\mathrm{d}x=1$$
  取积分中的函数，我们就得到Gamma分布的密度函数：
$$Gamma(x|\alpha )=\frac{x^{\alpha -1}{e}^{-x}}{\Gamma (\alpha )}$$

1.4 Beta分布

  在了解Beta分布前，需要先了解以下概念：

• 先验概率：先验概率就是事情尚未发生前，我们对该事发生概率的估计。
• 后验概率：后验概率是指通过调查或其它方式获取新的附加信息，利用贝叶斯公式对先验概率进行修正，而后得到的概率。
• 共轭分布：在贝叶斯概率理论中，如果后验概率P(θ|x)和先验概率p(θ)满足同样的分布律，那么，先验分布和后验分布被叫做共轭分布，同时，先验分布叫做似然函数的共轭先验分布。
• 极大似然估计(ML，Maximum Likelihood)可以估计模型的参数。其目标是找出一组参数 θ，使得模型产生出观测数据 x 的概率最大，其目标如下：$$\arg \underset{\theta }{\max }P(x|\theta )$$
• 最大后验估计(MAP－max a posterior)可以在要估计的模型的参数存在先验概率的情况下优化后验概率，其目标如下：
  $$\arg \underset{\theta }{\max }P(\theta |x)=\arg \underset{\theta }{\max }\frac{P(x|\theta )P(\theta )}{P(x)}$$
    因为样本x是给定的，所以$P(x)$是定值，可以忽略，则等式变为
  $$\arg \underset{\theta }{\max }P(\theta |x)=\arg \underset{\theta }{\max }P(x|\theta )P(\theta )$$
    即
  $$后验概率\propto 似然函数\times 先验概率$$
    贝塔分布（Beta Distribution) 是一个作为伯努利分布和二项式分布的共轭先验分布的密度函数，在机器学习和数理统计学中有重要应用。其概率密度函数如下：
  $$\begin{array}{rl}f(x;\alpha ,\beta )& =\frac{x^{\alpha -1}(1-x{)}^{\beta -1}}{{\int }_0^1{u}^{\alpha -1}(1-u{)}^{\beta -1}\mathrm{d}u}=\frac{\Gamma (\alpha +\beta )}{\Gamma (\alpha )\Gamma (\beta )}{x}^{\alpha -1}(1-x{)}^{\beta -1}\\ & =\frac{1}{B(\alpha ,\beta )}{x}^{\alpha -1}(1-x{)}^{\beta -1}\end{array}$$
    Beta分布的均值为：
  $$E(p)=\frac{\alpha }{\alpha +\beta }$$
    由于二项分布的似然函数为${\mu }^k(1-\mu {)}^{n-k}$，如果先验概率也是$\mu$和$1-\mu$的次方乘积的关系，那么根据$后验概率\propto 似然函数\times 先验概率$,后验概率分布形式将会和先验概率一样，这样后验概率和先验概率为共轭分布。而Beta分布正好符合这样的要求，因此Beta分布是二项分布的共轭先验。

1.5 狄利克雷分布Dirichlet

  Dirichlet分布则是多项式分布的共轭先验，也可以看做是beta分布推广到多变量的情况。Dirichlet概率密度函数定义如下：
$$Dir(p|\vec{\alpha })=\frac{1}{B(\vec{\alpha })}\prod _{k=1}^K{p}_k^{a_k-1},B(\vec{\alpha })=\frac{\Gamma ({\sum }_{k=1}^K{a}_k)}{{\prod }_{k=1}^K\Gamma (a_k)}$$
  我们要猜测参数$\vec{p}$，其先验分布为$Dir(\vec{p}|\vec{\alpha })$,数据落到不同区间的个数为$\vec{m}$，服从多项分布$Mul(\vec{m}|\vec{p})$，在给定了数据提供的知识$\vec{m}$，p的后验分布变为$Dir(\vec{p}|\vec{\alpha }+\vec{m})$，这个过程层也就是Dirichlet-Multinomial 共轭，用数学形式描述如下：
$$Dir(\vec{p}|\vec{\alpha })+Mul(\vec{m})=Dir(\vec{p}|\vec{\alpha }+\vec{m})$$
  Dirichlet分布的均值为：
$$E(p)=\left(\frac{a_1}{{\sum }_{i=1}^K{a}_i},\frac{a_2}{{\sum }_{i=1}^K{a}_i},...,\frac{a_K}{{\sum }_{i=1}^K{a}_i}\right)$$

2.概率论基础

  统计模拟中有一个重要的问题就是给定一个概率分布p(x)，我们如何在计算机中生成它的样本。常见的概率分布，都可以基于01均匀分布的样本生成，而01均匀分布的样本则可以通过线性同余发生器生成。而如果p(x)是比较复杂的形式，我们就需要一些更加复杂的随机模拟的方法来生成样本。

2.1 马氏链及其收敛性

  在马氏链中，第n+1个状态只依赖于第n个状态。
  在线性代数中，有这样一种矩阵，它的每个元素非负，并且每列之和为1，这样的矩阵被称为马尔科夫矩阵，同时也是一个状态转移矩阵。当我们给定一个初始状态${\pi }_0$时，有${\pi }_1={\pi }_{0}P$,${\pi }_2={\pi }_{1}P={\pi }_0{P}^2$,…,${\pi }_n={\pi }_0{P}^n$,而马尔科夫矩阵有一个性质，就是当任何两个状态联通，并且n足够大时，$P^n$会收敛到某一个概率分布，从而使得我们的分布不再发生变化。
  所以最初的概率分布${\pi }_0$，随着马尔科夫矩阵的转移，最终会收敛到某个分布$\pi (x)$。当$X_n\sim \pi (x)$时，有$X_{n+1}\sim \pi (x),X_{n+2}\sim \pi (x)$，此时$X_n,X_{n+1},X_{n+2},...$属于分布$\pi (x)$的样本。

2.2 MCMC算法

  如果，我们能够构造一个转移矩阵为P的马氏链，使得其平稳分布恰好是$\pi (x)$，那么我们从任意一个初始状态$x_0$出发，如果在n步收敛，那么接下去的样本$x_n,x_{n+1},...$都是$\pi (x)$的样本。
  细致平稳条件：对于达到平稳分布的转移矩阵P，有
$$\pi (i)P_{ij}=\pi (j)P_{ji}$$
  假设我们现在存在一个转移矩阵Q，q(i,j)表示从状态i转移到状态j的概率，通常情况下，细致平稳条件是不成立的。因此我们考虑两边同乘某个值，使得条件成立。令$\alpha (i,j)=p(j)q(j,i),\alpha (j,i)=p(i)q(i,j)$，按照对称性，有
$$p(i)q(i,j)\alpha (i,j)=p(i)q(i,j)\alpha (j,i)$$
  将$q(i,j)\alpha (i,j),q(i,j)\alpha (j,i)$看成是新的转移矩阵$Q^{\prime }(i,j)和Q^{\prime }(j,i)$，而$Q^{\prime }$满足细致平稳条件，于是我们得到了马氏链$Q^{\prime }$的平稳分布就是p(x)。
  而如果站在原来的转移矩阵Q上来考虑，新引进的$\alpha (i,j)$可以看做是接受率，也就是，从状态i以q(i,j)的概率转移到状态j的时候，以$\alpha (i,j)$的概率接受这个转移，从而转移概率为$q(i,j)\alpha (i,j)$。
  总结一下，我们可以得到MCMC算法的流程如下：

1. 初始化马氏链初始状态$X_0=x_0$
2. 进入循环t=0,1,2,…，在时刻t马氏链状态为$X_t=x_t$，采样$y\sim q(x|x_t)$
3. 从01均匀分布中生成随机数u
4. 如果$u<\alpha (x_t,y)=p(y)q(x_t|y)$，则接受转移$x_t\to y$，即X_{t+1}=y,否则$X_{t+1}=x_t$

2.3 Metropolis-Hastings算法

  对于式子$p(i)q(i,j)\times 0.2=p(i)q(i,j)\times 0.1$，虽然细致平稳条件成立，但是接受率太低，收敛得太慢。考虑将左右同乘5，得到
$$p(i)q(i,j)\times 1=p(i)q(i,j)\times 0.5$$
不难发现接受率不仅提高了，而且细致平稳条件依然成立。因此我们稍微改造一下MCMC算法，让
α ( i , j ) = min ⁡ { p ( j ) q ( j , i ) p ( i ) q ( i , j ) , 1 } \alpha(i,j)=\min\{\frac{p(j)q(j,i)}{p(i)q(i,j)},1\} α(i,j)=min{p(i)q(i,j)p(j)q(j,i)​,1}
就得到了效率更高的Metropolis-Hastings算法。

2.4 吉布斯采样(Gibbs Sampling)

  对于高维的情况，由于接受率的存在，Metropolis-Hastings算法的效率不算高，能够找到一个转移矩阵Q使得接受率a=1呢？在二维下，假设有一个概率分布p(x,y)，考虑坐标$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(A)p(y_2|x_1)=p(B)p(y_1|x_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)$$
  因此构造转移矩阵
$$Q=\{\begin{array}{ll}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_C=y_1\\ 0& otherwises\end{array}$$
在上面转移矩阵Q下，二维空间的马氏链将收敛到平稳分布p(x,y)。
  二维Gibbs Sampling的流程总结如下：

1. 随机初始化$X_0=x_0,Y_0=y_0$
2. 对于t=0,1,2,…，循环采样$y_{t+1}\sim p(y|x_t),x_{t+1}\sim p(x|y_{t+1})$

  将二维Gibbs Sampling推广到n维，其算法流程如下：

1. 随机初始化 { x i : i = 1 , 2 , 3 , . . . , n } \{x_i:i=1,2,3,...,n\} {xi​:i=1,2,3,...,n}
2. 对于t=0,1,2,…，循环采样
3. $x_1^{t+1}\sim p(x_1|x_2^t,x_3^t,...,x_n^t)$
4. $x_2^{t+1}\sim p(x_2|x_1^{t+1},x_3^t,...,x_n^t)$
5. $x_3^{t+1}\sim p(x_3|x_1^{t+1},x_2^{t+1},x_4^t,...,x_n^t)$
6. …
7. $x_n^{t+1}\sim p(x_n|x_1^{t+1},x_2^{t+1},x_3^{t+1},...,x_{n-1}^{t+1})$

3. 隐狄利克雷模型LDA

  在pLSA模型中，我们考虑将一篇文档考虑成由不同的主题和主题下的对应的高频词汇组成。但是文档选择某个主题的分布${\vec{\vartheta }}_m$，或者在主题下选择某个词汇的概率分布${\vec{\phi }}_k$，都没有任何的先验分布，而是由EM算法收敛得到其分布。由于${\vec{\vartheta }}_m$和${\vec{\phi }}_k$都对应到多项分布，所以其先验分布的一个好的选择就是Dirichlet分布，于是我们就得到了LDA模型。

  LDA模型分为两个部分：

1. $\vec{\alpha }\to {\vec{\vartheta }}_m\to z_{m,n}$表示，在生成第m篇文档时，先以一定概率抽中某个doc-topic分布${\vec{\vartheta }}_m$，再根据这个分布生成文档中第n个词的主题编号$z_{m,n}$
2. $\vec{\beta }\to {\vec{\phi }}_k\to w_{m,n}|(k=z_{m,n})$表示，在K个topic-word分布${\vec{\phi }}_k$中，挑选编号为k=$z_{m,n}$的分布生成单词$w_{m,n}$
     对第一个过程 $\vec{\alpha }\to {\vec{\vartheta }}_m\to z_{m,n}$来说，${\vec{\vartheta }}_m\to z_{m,n}$对应与多项分布，$\vec{\alpha }\to {\vec{\vartheta }}_m$对应于Dirichlet分布，所以整体是一个Dirichlet-Multinomial共轭结构。

  由于Dirichlet分布则是多项式分布的共轭先验，有
$$Dir({\vec{\theta }}_m|\vec{\alpha })+Mul({\vec{n}}_m)=Dir({\vec{\theta }}_m|\vec{\alpha }+{\vec{n}}_m)$$
  那么在已知先验概率$Dir({\vec{\theta }}_m|\vec{\alpha })$的情况下，我们直接可以得到后验概率
$$p({\vec{\theta }}_m|\vec{z},\vec{\alpha })=Dir({\vec{\theta }}_m|\vec{\alpha }+{\vec{n}}_m)=\frac{1}{\Delta (\vec{\alpha }+{\vec{n}}_m)}\prod _{k=1}^K{\theta }_k^{n_k+a_k-1}d\vec{\theta }$$
  进而我们可得到某个topic的产生概率为
$$\begin{array}{rl}p({\vec{z}}_m|\vec{a})& =\int p({\vec{z}}_m|{\vec{\theta }}_m)p({\vec{\theta }}_m|\vec{a})\mathrm{d}{\vec{\theta }}_m\\ & =\int \prod _{k=1}^K{\theta }_k^{n_k}Dir({\vec{\theta }}_m|\vec{a})\mathrm{d}{\vec{\theta }}_m\\ & =\int \prod _{k=1}^K{\theta }_k^{n_k}\frac{1}{\Delta (\vec{a})}\prod _{k=1}^K{\theta }_k^{a_k-1}\mathrm{d}{\vec{\theta }}_m\\ & =\frac{1}{\Delta (\vec{a})}\int \prod _{k=1}^K{\theta }_k^{n_k+a_k-1}\mathrm{d}{\vec{\theta }}_m\\ & =\frac{\Delta (\vec{a}+{\vec{n}}_m)}{\Delta (\vec{a})}\end{array}$$
  由于M篇文档的topic生成过程相互独立，因此整个topic的生成概率为
$$p(\vec{z}|\vec{\alpha })=\prod _{m=1}^{M}p({\vec{z}}_m|\vec{\alpha })=\prod _{m=1}^M\frac{\Delta ({\vec{n}}_m+\vec{\alpha })}{\Delta (\vec{\alpha })}$$

  令${\vec{w}}^{\prime }=({\vec{w}}_1,...{\vec{w}}_K),{\vec{z}}^{\prime }=({\vec{z}}_1,...{\vec{z}}_K)$，其中${\vec{w}}_i$表示第i类topic的词汇，${\vec{k}}_i$表示第i类topic词汇对应的topic编号。对于LDA模型的第二个部分$\vec{\beta }\to {\vec{\phi }}_k\to w_{m,n}|(k=z_{m,n})$，在$k=z_{m,n}$的限制下，任何两个由主题k生成的词都是可交换的，即使不在同一个文档中。考虑如下过程$\vec{b}eta\to {\vec{\phi }}_k\to {\vec{w}}_k$，其中$\vec{b}eta\to {\vec{\phi }}_k$对应于Dirichlet分布，${\phi }_k\to {\vec{w}}_k$对应于Multinomial分布。于是有
$$p({\vec{w}}_k|\vec{\beta })=\frac{\Delta ({\vec{n}}_k+\vec{\beta })}{\Delta \vec{\beta }}$$
  参数$\vec{\phi }k$的后验分布为
$$Dir({\vec{\phi }}_k|{\vec{n}}_k+\vec{\beta })$$
  整个语料中词生成概率为
$$\begin{array}{rl}p(\vec{w}|\vec{z},\vec{\beta })& =p({\vec{w}}^{\prime }|{\vec{z}}^{\prime },\vec{\beta })\\ & =\prod _{k=1}^{K}p({\vec{w}}_k|{\vec{z}}_k,\vec{\beta })\\ & =\prod _{k=1}^K\frac{\Delta ({\vec{n}}_k+\vec{\beta })}{\Delta (\vec{\beta })}\end{array}$$
c综合式子()和()得到
$$p(\vec{w},\vec{z}|\vec{\alpha },\vec{\beta })=p(\vec{w}|\vec{z},\vec{\beta })p(\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 })}$$
  得到联合分布$p(\vec{w},\vec{z})$后，我们就可以使用Gibbs Sampling对分布进行采样，由于$\vec{w}$已知，所以实际上只需要采样$p(\vec{z}|\vec{w})$ 。将所有的$z_i$看成是坐标轴，令$\neg i$表示去除下标为i的词。根据算法要求，我们需要求得$p(z_i=k|{\vec{z}}_{\neg i},\vec{w})$，假设已经观测到的词$w_i=t$，由贝叶斯法则，容易得到
$$p(z_i=k|{\vec{z}}_{\neg i},\vec{w})\propto p(z_i=k,w_i=t|{\vec{z}}_{\neg i},{\vec{w}}_{\neg i})$$
  去掉某个词并不改变共轭结构，只是对应的计数变少，从而有
$$\begin{gathered} p({\vec{\theta }}_m|{\vec{z}}_{\neg i},{\vec{w}}_{\neg i})=Dir({\vec{\theta }}_m|{\vec{n}}_{m,\neg i}+\vec{\alpha }) \\ p({\vec{\phi }}_k|{\vec{z}}_{\neg i},{\vec{w}}_{\neg i})=Dir({\vec{\phi }}_k|{\vec{n}}_{k,\neg i}+\vec{\beta }) \end{gathered}$$
  整合一下，得到
$$\begin{array}{rl}p(z_i=k|{\vec{z}}_{\neg i},\vec{w})& \propto p(z_i=k,w_i=t|{\vec{z}}_{\neg i},{\vec{w}}_{\neg i})\\ & =\int p(z_i=k,w_i=t,{\vec{\theta }}_m,{\vec{\phi }}_k|{\vec{z}}_{\neg i},{\vec{w}}_{\neg i})d{\vec{\theta }}_{m}d{\vec{\phi }}_k\\ & =\int p(z_i=k,{\vec{\theta }}_m|{\vec{z}}_{\neg i},{\vec{w}}_{\neg i})p(w_i=t,{\vec{\phi }}_k|{\vec{z}}_{\neg i},{\vec{w}}_{\neg i})d{\vec{\theta }}_{m}d{\vec{\phi }}_k\\ & =\int p(z_i=k|{\vec{\theta }}_m)Dir({\vec{\theta }}_m|{\vec{n}}_{m,\neg i}+\vec{\alpha })d{\vec{\theta }}_m\times p(w_i=t|{\vec{\phi }}_k)Dir({\vec{\phi }}_k|{\vec{n}}_{k,\neg i}+\vec{\beta })d{\vec{\phi }}_k\\ & =\int {\theta }_{mk}Dir({\vec{\theta }}_m|{\vec{n}}_{m,\neg i}+\vec{\alpha })d{\vec{\theta }}_m\times \int {\phi }_{kt}Dir({\vec{\phi }}_k|{\vec{n}}_{k,\neg i}+\vec{\beta })d{\vec{\phi }}_k\\ & =E({\theta }_{mk})E({\phi }_{kt})\\ & ={\hat{\theta }}_{mk}{\hat{\phi }}_{kt}\\ & =\frac{n_{m.\neg i}^{(k)}+{\alpha }_k}{{\sum }_{x=1}^K{n}_{m.\neg i}^{(x)}+{\alpha }_x}\frac{n_{k.\neg i}^{(t)}+{\beta }_t}{{\sum }_{x=1}^V{n}_{k.\neg i}^{(x)}+{\beta }_x}\end{array}$$
  上式便是LDA模型的Gibbs Samppling 公式，基于该公式我们可以训练LDA模型并应用到新的文档进行topic语义分析。训练的过程如下：

1. 对文档中每个词w，随机赋予一个topic编号z
2. 扫描整个语料库，对每个词，按照Gibbs Sampling公式重新采样它的topic，在语料中进行更新。重复直至Gibbs Sampling收敛。
3. 统计topic-word共现频率矩阵，该矩阵就是LDA的模型。

  利用LDA模型，我们可以对新来的文档进行主题的分类，分类的过程如下：

1. 对新来的文档每个词w随机赋予一个topic编号z
2. 重新扫描当前文档，按照Gibbs Sampling公式，对每个词w，重新采样它的topic。重复直至Gibbs Sampling收敛。
3. 统计文档中的topic分布得到${\vec{\theta }}_{new}$