Latent Dirichlet Allocation（LDA）学习笔记

1，Gamma函数

Gamma函数

[Gamma (x) = int_0^infty  {{e^{ - t}}{t^{x - 1}}dt} ]

是阶乘的从整数域到实数域的扩展

[Gamma (n) = (n - 1)!,n in { 0,1,2,3...} ]

函数递推推导如下，根据分布积分公式

[uv = int {(uv} )'dt = int {uv'dt + int {u'vdt} } ]

令

[u = frac{{{t^x}}}{x},v = {e^{ - t}}]

得

[Gamma (x) = int_0^infty  {{e^{ - t}}{t^{x - 1}}dt}  = int_0^infty  {u'vdt}  = uv|_0^infty  - int_0^infty  {uv'dt}  = frac{{{t^x}}}{{x{e^t}}}|_0^infty  - int_0^infty  { - {e^{ - t}}frac{{{t^x}}}{x}dt}  = frac{1}{x}int_0^infty  {{e^{ - t}}{t^x}dt}  = frac{1}{x}Gamma (x + 1)]

Gamma函数以及变形得到的Gamma分布在数学上应用很广，其发现过程可以参考《LDA数据八卦》

2，Bernoulli（伯努利）分布

伯努利分布又称为两点分布或者0-1分布，是一种离散型概率分布，指事件可能性存在两种，1或0，发生概率分别为p与1-p，即

[P(x) = left{ {egin{array}{*{20}{c}}
{p,x = 1}\
{1 - p,x = 0}
end{array}} ight.]

可知

[egin{array}{l}
E[x] = p*1 + (1 - p)*0 = p\
D[x] = E[{x^2}] - {E^2}[x] = p(1 - p)
end{array}]

3，二项分布

二项分布指重复n次的概率为p的伯努利分布的试验，记为X~B(n,p)，也是一种离散型概率分布，此时事件a发生k次（即事件b发生n-k次）的概率为

[P(k) = frac{{n!}}{{k!(n - k)!}}{p^k}{(1 - p)^{n - k}}]

可知

[E[x] = sumlimits_0^n {kfrac{{n!}}{{k!(n - k)!}}{p^k}} {(1 - p)^{n - k}} = sumlimits_1^n {frac{{n*(n - 1)!}}{{(k - 1)!(n - k)!}}p*{p^{k - 1}}} {(1 - p)^{n - k}} = npsumlimits_1^n {frac{{(n - 1)!}}{{(k - 1)!(n - k)!}}{p^{k - 1}}} {(1 - p)^{n - k}} ]

可知，式子后半部分为X~B(n-1,p)的全部可能取值之和

[E[x] = npsumlimits_0^t {frac{{t!}}{{s!(t - s)!}}{p^s}} {(1 - p)^{t - s}}{|_{s = k - 1,t = n - 1}} = np]

另，由

[E[{x^2}] = sumlimits_0^n {{k^2}frac{{n!}}{{k!(n - k)!}}{p^k}} {(1 - p)^{n - k}} = sumlimits_0^n {kfrac{{n!}}{{(k - 1)!(n - k)!}}{p^k}} {(1 - p)^{n - k}}]

想办法拆项消去k

[egin{array}{l}
E[{x^2}] = sumlimits_0^n {((k - 1) + 1)frac{{n!}}{{(k - 1)!(n - k)!}}{p^k}} {(1 - p)^{n - k}}\
= kfrac{{n!}}{{(k - 1)!(n - k)!}}{p^k}{(1 - p)^{n - k}}{|_{k = 1}} + sumlimits_2^n {((k - 1) + 1)frac{{n!}}{{(k - 1)!(n - k)!}}{p^k}} {(1 - p)^{n - k}}\
= np{(1 - p)^{n - 1}} + sumlimits_2^n {frac{{n!}}{{(k - 2)!(n - k)!}}{p^k}} {(1 - p)^{n - k}} + sumlimits_1^n {frac{{n!}}{{(k - 1)!(n - k)!}}{p^k}} {(1 - p)^{n - k}} - frac{{n!}}{{(k - 1)!(n - k)!}}{p^k}{(1 - p)^{n - k}}{|_{k = 1}}\
= np{(1 - p)^{n - 1}} + n(n - 1){p^2}sumlimits_2^n {frac{{(n - 2)!}}{{(k - 2)!(n - k)!}}{p^{k - 2}}} {(1 - p)^{n - k}} + npsumlimits_1^n {frac{{(n - 1)!}}{{(k - 1)!(n - k)!}}{p^{k - 1}}} {(1 - p)^{n - k}} - np{(1 - p)^{n - 1}}\
= {n^2}{p^2} - {n^2}{p} + np
end{array}]

可得

[D[x] = E[{x^2}] - {E^2}[x] = {n^2}{p^2} - {n^2}{p} + np - {(np)^2} = np(1 - p)]

4，Beta函数

Beta函数跟Gamma函数关系十分密切

[B(alpha ,eta ) = int_0^1 {{x^{alpha  - 1}}{{(1 - x)}^{eta  - 1}}dx = frac{{Gamma (alpha )Gamma (eta )}}{{Gamma (alpha  + eta )}}} ]

这个等式的证明过于复杂，详见参考文献

5，Beta分布

将Beta函数作为归一化分母，将Beta函数积分内部的式子作为分子，可以得到一个自变量介于0-1之间的、积分结果为1的函数

[f(x;alpha ,eta ) = frac{{{x^{alpha  - 1}}{{(1 - x)}^{eta  - 1}}}}{{int_0^1 {{x^{alpha  - 1}}{{(1 - x)}^{eta  - 1}}dx} }} = frac{{Gamma (alpha  + eta )}}{{Gamma (alpha )Gamma (eta )}}{x^{alpha  - 1}}{(1 - x)^{eta  - 1}}]

将这个函数定义为Beta分布的概率密度函数，由于自变量介于0-1之间，可以认为这是一个概率的概率分布，记为X~Beta(α,β)

[E[x] = int_0^1 {xfrac{{Gamma (alpha  + eta )}}{{Gamma (alpha )Gamma (eta )}}{x^{alpha  - 1}}{{(1 - x)}^{eta  - 1}}dx}  = int_0^1 {frac{{Gamma (alpha  + eta )}}{{Gamma (alpha )Gamma (eta )}}{x^alpha }{{(1 - x)}^{eta  - 1}}dx}  = int_0^1 {frac{alpha }{{(alpha  + eta )}}frac{{Gamma (alpha  + eta  + 1)}}{{Gamma (alpha  + 1)Gamma (eta )}}{x^alpha }{{(1 - x)}^{eta  - 1}}dx}  = frac{alpha }{{(alpha  + eta )}}]

可以看到Beta分布的概率密度函数与二项分布的概率密度函数形式非常接近，其实两者的确是同根同源的，关于两者的关系，可以参考《LDA数学八卦》中关于撒旦的游戏的描述

6，多项分布

将二项分布的事件可能结果扩展为k个，则变形为多项分布，概率分布函数为

[P({x_1},{x_2},...,{x_k}) = frac{{n!}}{{sumlimits_1^k {({x_i}!)} }}prodlimits_1^k {{p_i}^k} ,{x_i} ge 0,sumlimits_1^k {{x_i} = 1} ]

7，Dirichlet（狄利克雷）分布

将Beta分布的概率分布扩展到多维，就可以得到狄利克雷分布，记为X~Dir(α)

[f({x_i};{alpha _i}) = frac{{prodlimits_1^k {{x_i}^{{alpha _i} - 1}} }}{{int_0^1 {prodlimits_1^k {{x_i}^{{alpha _i} - 1}} dvec x} }} = frac{{prodlimits_1^k {{x_i}^{{alpha _i} - 1}} }}{{Delta (vec alpha )}} = frac{{Gamma (sumlimits_1^k {{alpha _i}} )}}{{prodlimits_1^k {Gamma ({alpha _i})} }}prodlimits_1^k {{x_i}^{{alpha _i} - 1}} ]

[Delta (vec alpha ) = int_0^1 {prodlimits_1^k {{x_i}^{{alpha _i} - 1}} dx}  = frac{{prodlimits_1^k {Gamma ({alpha _i})} }}{{Gamma (sumlimits_1^k {{alpha _i}} )}}]

这个式子的证明更加复杂

同Beta分布，Dirichlet分布的期望为

[E[{x_i}] = frac{{{alpha _i}}}{{sumlimits_1^k {{alpha _i}} }}]

8，先验分布、后验分布、似然函数与共轭先验分布

根据贝叶斯公式，先验概率、似然函数与后验概率的关系如下

[P( heta |x) = frac{{P(x| heta )P( heta )}}{{P(x)}}]

如果P(θ)不是一个确定的概率，而是一个概率分布，则可得先验分布、似然函数与后验分布的关系

[f(x;alpha ,k) = frac{{P(k|x)f(x;alpha )}}{{int {P(k|x)f(x;a)dx} }}]

假设先验分布满足Beta分布，X~Beta(α,β)，似然函数满足二项分布，则后验分布为

[f(x;alpha ,eta ,k) = frac{{{x^{{k_1} - 1}}{{(1 - x)}^{{k_2} - 1}}frac{{Gamma (alpha  + eta )}}{{Gamma (alpha )Gamma (eta )}}{x^{alpha  - 1}}{{(1 - x)}^{eta  - 1}}}}{{int {{x^{{k_1} - 1}}{{(1 - x)}^{{k_2} - 1}}frac{{Gamma (alpha  + eta )}}{{Gamma (alpha )Gamma (eta )}}{x^{alpha  - 1}}{{(1 - x)}^{eta  - 1}}dx} }} =  frac{{frac{{Gamma (alpha  + eta  + {k_1} + {k_2})}}{{Gamma (alpha  + {k_1})Gamma (eta  + {k_2})}}{x^{alpha  + {k_1} - 1}}{{(1 - x)}^{eta  + {k_2} - 1}}}}{{int {frac{{Gamma (alpha  + eta  + {k_1} + {k_2})}}{{Gamma (alpha  + {k_1})Gamma (eta  + {k_2})}}{x^{alpha  + {k_1} - 1}}{{(1 - x)}^{eta  + {k_2} - 1}}dx} }}]

可见后验分布还是Beta分布，X~Beta(α+k₁,β+k₂)

因此，称Beta分布为二项分布的共轭先验分布

扩展到多维，同样可证，Dirichlet分布为多项分布的共轭先验分布

8，随机模拟（Monte Carlo Method）

给定一个概率分布P(x)，利用计算机随机生成样本，称为随机模拟

均匀分布是最容易生成的样本的分布，一些常见分布如Gamma分布、正态分布，都可以通过均匀分布变换得到

但有些复杂概率分布的样本不能通过均匀分布变换得到

9，马尔科夫链（Markov Chain）与平稳分布

马尔科夫链的定义为

[P({s_n}|{s_{n - 1}},{s_{n - 2}},{s_{n - 3}},...) = P({s_n}|{s_{n - 1}})]

根据定义，马尔科夫链中各个状态相互转移只取决于当前状态，将状态相互转移的概率P_i,j构成的矩阵称为转移矩阵

根据马尔科夫链收敛定理，如果一个马尔科夫链的转移概率矩阵存在，并且任意两个状态都是联通的（即有限步内可达），那么这个马尔科夫链存在平稳分布，即经过若干次转移后，处于每个状态的概率稳定，这个稳定的概率分布即平稳分布

马尔科夫链处于平稳状态下的充分必要条件是满足细致平稳条件

[{s_i}{P_{i,j}} = {s_j}{P_{j,i}}]

10，MCMC（Monte Carlo Markov Chain）方法

如果能构造一个马尔科夫链，此链的平稳分布为待采样的概率分布p(x)，则可以通过在已经收敛到平稳状态的马尔科夫链上采样，作为p(x)的采样样本

构造一个转移矩阵为P的马尔科夫链，令其平稳分布为s(x)，则一般情况下

[{s_i}{P_{i,j}} e {s_j}{P_{j,i}}]

此时引入

[{alpha _{i,j}} = {s_j}{P_{j,i}},{alpha _{j,i}} = {s_i}{P_{i,j}}]

则可得

[{s_i}{P_{i,j}}{alpha _{i,j}} = {s_j}{P_{j,i}}{alpha _{j,i}}]

称α为接收矩阵，意义为在原转移概率基础上，有一定几率接收转移，有一定几率拒绝转移，得到平稳分布

11，M-H（Metropolis-Hastings）采样方法

Metropolis-Hastings方法是改造的MCMC方法，由于在一般MCMC方法中，转移概率还会有一定接收率，有可能很小，则采样时会大概率拒绝转移，为了提高转移概率，Metropolis-Hastings方法将接收概率在满足细致平稳条件的前提下尽量增大

[{alpha _{i,j}} = min { frac{{{s_j}{P_{j,i}}}}{{{s_i}{P_{i,j}}}},1} ]

12，Gibbs采样方法

Gibbs采样方法是一种特殊的Metropolis-Hastings方法，将接受率增大为1,

考虑二维的情况，在概率分布上取两点A(x₁,y₁)，B(x₁,y₂)，并将概率分布p(y₁|x₁)，p(y₂|x₁)作为转移概率，可得

[egin{array}{l}
p({x_1},{y_1})p({y_2}|{x_1}) = p({y_1}|{x_1})p({x_1})p({y_2}|{x_1})\
p({x_1},{y_2})p({y_1}|{x_1}) = p({y_2}|{x_1})p({x_1})p({y_1}|{x_1})
end{array}]

可见，则平稳分布条件是成立的

[p({x_1},{y_2})p({y_1}|{x_1}) = p({x_1},{y_2})p({y_1}|{x_1})]

Gibbs采样即基于此等式，为此马尔科夫链构造转移矩阵P

[{P_{AB}} = left{ {egin{array}{*{20}{c}}
{p({y_B}|{x_A}),{x_A} = {x_B}}\
{p({x_B}|{y_A}),{y_A} = {y_B}}\
{0,{x_A} e {x_B},{y_A} e {y_B}}
end{array}} ight.]

在采样时，选取一个维度，利用现有所有样本点估计这个维度的概率分布，并随机选取采样值得到新的样本点位置，选取到足够多的样本点即认为可以代表真实分布

13，文本建模

有若干文档，假设只考虑文档中各个单词的个数，不考虑单词出现次序，即Bag-Of-Words模型

文本建模问题即将文档中单词的生成看做上帝抛掷骰子的结果，模型预测的参数有两个：1）有哪些样式的骰子，2）怎样抛掷这些骰子

14，简单Unigram模型

简单Unigram Model认为：

1. 只有一枚骰子，每面代表一个单词，各个单词概率不尽相同
2. 不断抛掷这枚骰子，生成所有文档的所有单词

可见所有文档都是统一的，即所有单词都在一个袋子里

令V为单词集合（字典），M为单词总数（字典大小），w_m为所有文档中单词v_m的出现次数

令p_m为模型中单词v_m的出现概率，文档生成的过程满足多项分布Multi(w|p)

文档生成的概率为

[P(W) = prodlimits_m^M {p_m^{{w_m}}} ]

根据最大似然估计

[{{hat p}_m} = frac{{{w_m}}}{{sumlimits_m^M {{w_m}} }}]

15，贝叶斯Unigram模型

贝叶斯认为骰子也是有概率分布的，所以等同于：

1. 有若干枚骰子，每个骰子各面概率不尽相同，每面代表一个单词，各个单词概率不尽相同
2. 随机选出一枚骰子，不断抛掷这枚骰子，生成所有文档的所有单词

所以贝叶斯Unigram模型也是将所有单词都放在一个袋子里

令V为单词集合（字典），M为单词总数（字典大小），w_m为所有文档中单词v_m的出现次数

令p_m为模型中单词v_m的出现概率，由于文档生成的过程满足多项分布Multi(w|p)

[P(vec w|vec p) = prodlimits_m^M {p_m^{{w_m}}} ]

所以采用多项分布的共轭先验分布Dirichlet分布作为先验分布Dir(p|α)，此处α为先验分布的超参数

[f(vec p;vec alpha ) = frac{{prodlimits_m^M {p_m^{{alpha _m} - 1}} }}{{int {prodlimits_m^M {p_m^{{alpha _m} - 1}} dvec p} }}]

根据多项分布于Dirichlet分布的关系，后验分布也应是Dirichlet分布，后验分布满足

[f(vec p;vec alpha ,vec w) = frac{{prodlimits_m^M {p_m^{{alpha _m} + {w_m} - 1}} }}{{int {prodlimits_m^M {p_m^{{alpha _m} + {w_m} - 1}} dvec p} }}]

此时文档生成概率为

[P(W) = int {P(vec w|vec p)f(vec p;vec alpha )dvec p}  = int {f(vec p;vec alpha ,vec w)dvec p}  = frac{{int {prodlimits_m^M {p_m^{{alpha _m} + {w_m} - 1}} dvec p} }}{{int {prodlimits_m^M {p_m^{{alpha _m} - 1}} dvec p} }}]

将后验分布的平均值作为参数估计值

[{{hat p}_m} = frac{{{alpha _m} + {w_m}}}{{sumlimits_m^M {({alpha _m} + {w_m})} }}]

16，PLSA（Probabilistic Latent Semantic Analyse）

Unigram模型认为所有文档都是独立同分布的，PLSA在Unigram模型基础上，认为每个文档都由若干个不尽相同的主题构成，每个单词都由主题以一定概率生成，即：

1. 每篇文档有一枚特定的doc-topic骰子
2. 有若干枚topic-word骰子
3. 生成每个词前，先抛掷doc-topic骰子，再根据doc-topic骰子结果选择topic-word骰子，抛掷topic-word骰子，生成单词

令V为单词集合（字典），M为单词总数（字典大小），N为文档个数，K为主题个数

令d_n代表文档n，z_k代表主题k，v_m代表单词m，则p(z_k|d_n)为文档d_n生成主题t_k的概率，p(v_m|z_k)为主题t_k生成单词v_m的概率，并且模型假设p(z_k|d_n)、p(v_m|z_k)相互独立

令w_nm为文档d_n中单词v_m的个数，则文档生成概率为

[egin{array}{*{20}{l}}
{P(W) = prodlimits_n^N {prodlimits_m^M {p{{({v_m},{d_n})}^{{w_{nm}}}}} } }\
{p({v_m},{d_n}) = p({d_n})sumlimits_k^K {p({v_m}|{z_k})p({z_k}|{d_n})} }
end{array}]

可得log似然函数如下

[L(W) = sumlimits_n^N {sumlimits_m^M {{w_{nm}}log (p({v_m},{d_n}))} }  = sumlimits_n^N {sumlimits_m^M {{w_{nm}}(log (p({d_n})) + log (sumlimits_k^K {p({v_m}|{z_k})p({z_k}|{d_n})))} } } ]

去除常数项

[L'(W) = sumlimits_n^N {sumlimits_m^M {{w_{nm}}log (sumlimits_k^K {p({v_m}|{z_k})p({z_k}|{d_n})))} } } ]

使用EM算法估计似然函数最大值

首先假设p(z_k|d_n)与p(v_m|z_k)都已知，则可求后验概率

[p({z_k}|{v_m},{d_n}) = frac{{p({z_k},{v_m},{d_n})}}{{p({v_m},{d_n})}} = frac{{p({v_m}|{z_k})p({z_k}|{d_n})p({d_n})}}{{p({v_m}|{d_n})p({d_n})}}]

可知在后验概率中m、n的相关项都是一定的，则后验概率可化为下式求解

[p({z_k}|{v_m},{d_n}) = frac{{p({v_m}|{z_k})p({z_k}|{d_n})}}{{sumlimits_k^K {p({v_m}|{z_k})p({z_k}|{d_n})} }}]

然后根据当该后验概率求令似然函数期望最大化的p(z_k|d_n)与p(v_m|z_k)，注意此处log函数中的求和消失了，改到了函数外面，因为对于每一个特定的主题k，p(z_k|d_n)与p(v_m|z_k)也是特定的，不用求和

[E[L'(W)] = sumlimits_n^N {sumlimits_m^M {{w_{nm}}sumlimits_k^K {p({z_k}|{v_m},{d_n})log (p({v_m}|{z_k})p({z_k}|{d_n})))} } } ,sumlimits_m^M {p({v_m}|{z_k}) = 1,sumlimits_k^K {p({z_k}|{d_n}) = 1} } ]

期望函数有约束，则利用Lagrange Multiplier拉格朗日乘子法构造Lagrange函数

[f(x) = E[L'(W)] + alpha (sumlimits_m^M {p({v_m}|{z_k})}  - 1) + eta (sumlimits_k^K {p({z_k}|{d_n})}  - 1)]

令函数偏导数为0

[egin{array}{l}
frac{{partial f(x)}}{{partial p({v_m}|{z_k})}} = sumlimits_n^N {frac{{{w_{nm}}p({z_k}|{v_m},{d_n})}}{{p({v_m}|{z_k})}}} + alpha = 0\
frac{{partial f(x)}}{{partial p({z_k}|{d_n})}} = sumlimits_m^M {frac{{{w_{nm}}p({z_k}|{v_m},{d_n})}}{{p({z_k}|{d_n})}}} + eta = 0
end{array}]

根据p(v_m|z_k)和为1，代入上式消去α

[egin{array}{l}
sumlimits_m^M {p({v_m}|{z_k})} = sumlimits_m^M {frac{{sumlimits_n^N {{w_{nm}}p({z_k}|{v_m},{d_n})} }}{{ - alpha }}} = 1\
alpha = - sumlimits_n^N {sumlimits_m^M {{w_{nm}}p({z_k}|{v_m},{d_n})} }
end{array}]

可得p(v_m|z_k)的迭代公式

[p({v_m}|{z_k}) = frac{{sumlimits_n^N {{w_{nm}}p({z_k}|{v_m},{d_n})} }}{{sumlimits_n^N {sumlimits_m^M {{w_{nm}}p({z_k}|{v_m},{d_n})} } }}]

同样可得p(z_k|d_n)的迭代公式

[egin{array}{l}
sumlimits_k^K {p({z_k}|{d_n})} = sumlimits_k^K {frac{{sumlimits_m^M {{w_{nm}}p({z_k}|{v_m},{d_n})} }}{{ - eta }}} = 1\
eta = - sumlimits_m^M {sumlimits_k^K {{w_{nm}}p({z_k}|{v_m},{d_n})} } \
p({z_k}|{d_n}) = frac{{sumlimits_m^M {{w_{nm}}p({z_k}|{v_m},{d_n})} }}{{sumlimits_m^M {sumlimits_k^K {{w_{nm}}p({z_k}|{v_m},{d_n})} } }}
end{array}]

依样迭代，直至收敛

17，LDA（Latent Dirichlet Allocation）

LDA将PLSA做了贝叶斯改造，即认为doc-topic的概率与topic-word的概率都是有分布的，对应于骰子模型即：

1. 先从topic-word分布中选择K个topic-word骰子
2. 对每一篇文档，从doc-topic分布中选择1个doc-topic骰子
3. 对每一个单词，先抛掷doc-topic骰子得到topic编号k，再抛掷第k个topic-word骰子得到word

假设doc-topic先验分布为Dirichlet分布Dir(θ|α)，抛掷doc-topic骰子符合多项分布Multi(z|θ)，则可知doc-topic后验分布也是满足Dirichlet分布，设z_(d)为文档d的doc-topic分布，θ_(d),i为文档d中主题i的概率，n_(d),i为文档d中主题i出现的次数，K为主题总数

[P({{vec z}_{(d)}}|vec alpha ) = int {p({{vec z}_{(d)}}|{{vec heta }_{(d)}})f({{vec heta }_{(d)}};vec alpha )dvec heta }  = int {prodlimits_i^K {{{({ heta _{(d),i}})}^{{n_{(d),i}}}}} frac{{prodlimits_i^K {{{({ heta _{(d),i}})}^{{alpha _i} - 1}}} }}{{int {prodlimits_i^K {{{({ heta _{(d),i}})}^{{alpha _i} - 1}}dvec heta } } }}dvec heta }  = frac{{Delta ({{vec n}_{(d)}} + vec alpha )}}{{Delta (vec alpha )}}]

再假设topic-word先验分布为Dirichlet分布Dir(φ|β)，抛掷topic-word骰子符合多项分别Multi(v|φ)，则可知topic-word后验分布也是满足Dirichlet分布，设v_(k)为主题k的topic-word分布，φ_(k),i为主题k生成词语i的概率，m_(k),i为主题k生成词语i的次数，M为词库总数

[P({{vec v}_{(k)}}|vec eta ) = int {p({{vec v}_{(k)}}|{{vec varphi }_{(k)}})f({{vec varphi }_{(k)}};vec eta )dvec varphi }  = int {prodlimits_i^M {{{({varphi _{(k),i}})}^{{m_{(k),i}}}}} frac{{prodlimits_i^M {{{({varphi _{(k),i}})}^{{eta _i} - 1}}} }}{{int {prodlimits_i^M {{{({varphi _{(k),i}})}^{{eta _i} - 1}}dvec varphi } } }}dvec varphi }  = frac{{Delta ({{vec m}_{(k)}} + vec eta )}}{{Delta (vec eta )}}]

令D为文档总数，综合上两式可得，所有文档生成的后验概率为

[egin{array}{*{20}{l}}
{P(vec v,vec z|vec alpha ,vec eta )}\
{ = P(vec v|vec z,vec eta )P(vec z|vec alpha )}\
{ = prodlimits_i^K {P({{vec v}_{(i)}}|{{vec z}_{(i)}},vec eta )} prodlimits_i^D {P({{vec z}_{(i)}}|vec alpha )} }\
{ = prodlimits_i^K {frac{{Delta ({{vec m}_{(i)}} + vec eta )}}{{Delta (vec eta )}}} prodlimits_i^D {frac{{Delta ({{vec n}_{(i)}} + vec alpha )}}{{Delta (vec alpha )}}} }
end{array}]

令v_î为当前选取的样本的词语，z_î为当前选取的样本的主题，v_¬î为除了当前样本其他所有的样本的词语，z_¬î为除了当前样本其他所有的样本的主题

i为当前样本词语，k为所求主题，d为当前样本所在文档。注意“î”代表当前样本，只代表一个样本，“i”则代表当前样本的词语，代表该词语的若干个样本

则n_(d),k为文档d中主题k出现的次数，n_(d),k,¬î为文档d中除样本î外主题k出现的次数，m_(k),i为主题k生成词语i的次数，m_(k),i,¬î为主题k中除样本î外生成词语i的次数

Gibbs采样时使用的概率分布如下

[egin{array}{l}
P({z_{hat i}} = k)\
= P({v_{hat i}},{z_{hat i}}|{{vec v}_{ eg hat i}},{{vec z}_{ eg hat i}},vec alpha ,vec eta )\
= frac{{P({v_{hat i}},{z_{hat i}},{{vec v}_{ eg hat i}},{{vec z}_{ eg hat i}}|vec alpha ,vec eta )}}{{P({{vec v}_{ eg hat i}},{{vec z}_{ eg hat i}}|vec alpha ,vec eta )}}\
= frac{{P(vec v,vec z|vec alpha ,vec eta )}}{{P({{vec v}_{ eg hat i}},{{vec z}_{ eg hat i}}|vec alpha ,vec eta )}}\
= frac{{prodlimits_j^K {frac{{Delta ({{vec m}_{(j)}} + vec eta )}}{{Delta (vec eta )}}} cdot prodlimits_j^D {frac{{Delta ({{vec n}_{(j)}} + vec alpha )}}{{Delta (vec alpha )}}} }}{{prodlimits_j^{K, eg k} {frac{{Delta ({{vec m}_{(j)}} + vec eta )}}{{Delta (vec eta )}}} cdot frac{{Delta ({{vec m}_{(k), eg hat i}} + vec eta )}}{{Delta (vec eta )}} cdot prodlimits_j^{D, eg d} {frac{{Delta ({{vec n}_{(j)}} + vec alpha )}}{{Delta (vec alpha )}}} cdot frac{{Delta ({{vec n}_{(d), eg hat i}} + vec alpha )}}{{Delta (vec alpha )}}}}\
= frac{{Delta ({{vec m}_{(k)}} + vec eta )Delta ({{vec n}_{(d)}} + vec alpha )}}{{Delta ({{vec m}_{(k), eg hat i}} + vec eta )Delta ({{vec n}_{(d), eg hat i}} + vec alpha )}}\
= frac{{frac{{prodlimits_j^M {Gamma ({m_{(k),j}} + {eta _j})} }}{{Gamma (sumlimits_j^M {({m_{(k),j}} + {eta _j})} )}} cdot frac{{prodlimits_j^K {Gamma ({n_{(d),j}} + {alpha _j})} }}{{Gamma (sumlimits_j^K {({n_{(d),j}} + {alpha _j}} ))}}}}{{frac{{prodlimits_j^M {Gamma ({m_{(k),j, eg hat i}} + {eta _j})} }}{{Gamma (sumlimits_j^M {({m_{(k),j, eg hat i}} + {eta _j})} )}} cdot frac{{prodlimits_j^K {Gamma ({n_{(d),j, eg hat i}} + {alpha _j})} }}{{Gamma (sumlimits_j^K {({n_{(d),j, eg hat i}} + {alpha _j})} )}}}}\
= frac{{prodlimits_j^M {Gamma ({m_{(k),j}} + {eta _j})} }}{{prodlimits_j^M {Gamma ({m_{(k),j, eg hat i}} + {eta _j})} }} cdot frac{{prodlimits_j^K {Gamma ({n_{(d),j}} + {alpha _j})} }}{{prodlimits_j^K {Gamma ({n_{(d),j, eg hat i}} + {alpha _j})} }} cdot frac{{Gamma (sumlimits_j^M {({m_{(k),j, eg hat i}} + {eta _j})} )}}{{Gamma (sumlimits_j^M {({m_{(k),j}} + {eta _j})} )}} cdot frac{{Gamma (sumlimits_j^K {({n_{(d),j, eg hat i}} + {alpha _j})} )}}{{Gamma (sumlimits_j^K {({n_{(d),j}} + {alpha _j}} ))}}\
= frac{{Gamma ({m_{(k),i}} + {eta _i})}}{{Gamma ({m_{(k),i, eg hat i}} + {eta _i})}} cdot frac{{Gamma ({n_{(d),k}} + {alpha _k})}}{{Gamma ({n_{(d),k, eg hat i}} + {alpha _k})}} cdot frac{{Gamma (sumlimits_j^M {({m_{(k),j, eg hat i}} + {eta _j})} )}}{{Gamma (sumlimits_j^M {({m_{(k),j, eg hat i}} + {eta _j})} + 1)}} cdot frac{{Gamma (sumlimits_j^K {({n_{(d),j, eg hat i}} + {alpha _j})} )}}{{Gamma (sumlimits_j^K {({n_{(d),j, eg hat i}} + {alpha _j}} ) + 1)}}\
= ({m_{(k),i, eg hat i}} + {eta _i}) cdot ({n_{(d),k, eg hat i}} + {alpha _k}) cdot frac{1}{{sumlimits_j^M {({m_{(k),j, eg hat i}} + {eta _j})} }} cdot frac{1}{{sumlimits_j^K {({n_{(d),j, eg hat i}} + {alpha _j})} }}\
propto ({n_{(d),k, eg hat i}} + {alpha _k}) cdot frac{{{m_{(k),i, eg hat i}} + {eta _i}}}{{sumlimits_j^M {({m_{(k),j, eg hat i}} + {eta _j})} }}
end{array}]

采样时首先给每个词语随机分配一个主题，作为初始采样点，然后迭代每篇文档的每个词语，利用上述概率分布为这个词语随机选取一个新主题作为新采样点，不断迭代直至采样点达到一定数量，令现有采样点分布作为真实分布

由于doc-topic和topic-word后验概率都满足Dirichelet分布，根据Dirichlet分布的期望，参数估计方式如下，n_(d),i为文档d中主题i出现的次数，m_(k),i为主题k生成词语i的次数

[egin{array}{l}
P({{hat heta }_{(d)}} = k|alpha ) = frac{{{n_{(d),k}} + {alpha _k}}}{{sumlimits_j^K {({n_{(d),j}} + {alpha _j})} }}\
P({{hat varphi }_{(k)}} = i|eta ) = frac{{{m_{(k),i}} + {alpha _i}}}{{sumlimits_j^K {({m_{(k),j}} + {alpha _j})} }}
end{array}]

参考文献：

数学基础：

    《LDA数学八卦》

LDA漫游系列：

    https://www.jianshu.com/p/d8485c623669

    https://www.jianshu.com/p/8fb2fcb52a3a

    https://www.jianshu.com/p/e7fbd3a2b786

    https://www.jianshu.com/p/1511c94b2ac3

文本主题模型之LDA系列：

    http://www.cnblogs.com/pinard/p/6831308.html

    http://www.cnblogs.com/pinard/p/6867828.html

    http://www.cnblogs.com/pinard/p/6873703.html

LSA，pLSA原理及其代码实现：

    http://www.cnblogs.com/bentuwuying/p/6219970.html

Beta函数与Gamma 函数关系推导：

    https://blog.csdn.net/xhf0374/article/details/53946146

LDA Gibbs采样公式推导：

    https://blog.csdn.net/z_q3322183/article/details/80969570

扩展模型ATM：

    《The Author-Topic Model for Authors and Documents》

集群并行LDA：

    《PLDA: parallel latent Dirichlet allocation for large-scale applications》

LDA模型变种：

    https://blog.csdn.net/pirage/article/details/9611885