0. latent dirichlet allocation 前言
最近公司分享了这个topic， 我自己钻研了一下写一下自己对这个模型的理解。
因为本人在阅读网上各种资料时，发现写的比较不适合新手。所以写下这篇博文。
一是，对自己的学习收获做一个总结。另外，是希望可以帮助新手在学习这个模型的过程中能够获得比较好的理解。
我觉得理解一个事物的最好方式是，是从最高层级开始，知道它的全貌，忽略它的细枝末节。随后看到一个高粒度抽象的架构，去理解它的工作模式。最后才是深入细节底层。去逐个击破各个抽象的背后实现原理。
那么我就用这种方式来写这篇文章。
所以这篇文章的写法会和网上大多数的不同。他们都是从最底层需要的数学知识讲起，做了大量公式推导。开始看的云里雾里，之后串的时候就会看了后面忘了前面，思维跳跃跨度很大。下面是这篇文章是怎么展开的，好让读者心中有数。
- 把这个模型看成黑盒，它需要什么，它能给你什么（最顶层）
- 理解这个模型有哪些组件构成，这些组件怎么合作的 （看到高粒度的组件）
- 这些组件背后的数学和算法是什么，这些原理间的联系是什么 （高粒度组件下用了哪些数学算法工具）
- 我们如何基于这些原理去实现自己的LDA （不用知道推导，已经可以实现这个模型了）
- 背后的数学详细推导 （每个数学结论是怎么得来的）
1. 这个模型可以解决什么问题(顶层）
当我们有一组文章（如文章A， 文章B， 文章C）； 我们希望知道每个文章有哪些主题。（如主题A， 主题B， 主题C）我们就可以使用这个模型。与此同时，又来了一篇新的文章，我们也可以丢进这个模型里，它会告诉我们新文章里有哪些主题。
解决的问题分为2个阶段
- 训练阶段：
输入是一组文章（article = list <word>, 所以输入是list<article>）和 一个数字K（表示这组文章里你觉得大概有多少个TOPIC）。 这个K是个需要人们给的参数，会决定这个模型的表现，所以需要调参侠，去对K做调整，得到最好的模型效果。
输出 输入的文章，每一篇都会给一个文章到TOPIC的概率分布。 每一个TOPIC，都会告诉你这个TOPIC下每一个单词的概率分布。
具体例子：
假设我们有4篇文章, 并且指定K为3（我认为这组文章应该用3个TOPIC就可以囊括了）
文章a: ball ball ball planet galaxy
文章b: referendum planet planet referendum referendum
文章c: planet planet galaxy planet ball
文章d: planet galaxy referendum planet ball
可能的输出如下：
TOPIC下每一个单词的概率分布
topic 1: {ball:70%, galaxy:25%, referendum :1%, planet : 4%}
topic 2: {ball:9%, galaxy:35%, referendum :1%, planet : 55%}
topic 3: {ball:2%, galaxy:3%, referendum :80%, planet : 15%}
根据观察我们可以根据人类的知识，总结出topic1可能是体育。因为有球，而GALAXY可能是一只球队名。 topic2可能是科学，因为讲行星和银河系。 topic3可能是讲社会科学
每一个文章到TOPIC的概率分布
文章a: {topic1: 80%, topic2: 10%, topic3: 10%} 文章极大概率是体育类
文章b: {topic1: 5%, topic2: 85%, topic3: 10%} 文章极大概率是社科类
文章c: {topic1: 10%, topic2: 2%, topic3: 88%} 文章极大概率是科学类
文章d: {topic1: 25%, topic2: 15%, topic3: 60%} 文章大概率是科学类
- 预测阶段：
输入一篇新文章
galaxy ball ball ball ball
输出这篇文章的概率分布
{topic1: 90%, topic2: 1%, topic3: 9%}
2. 理解这个模型有哪些组件构成，这些组件怎么合作的 （高粒度的组件）
这个模型有4个组件，2个旋钮，2组齿轮。
当我们决定了K是多少之后。整个LDA模型就可以想象成一个机器如下图。它有2个可以调控的旋钮，旋钮调整的好坏决定了模型的好坏。调整这2个旋钮的作用就是决定内部齿轮的工作方式。
我们按下红色按钮，这个机器就会输出一篇文章。这个文章是一个装着很多单词的袋子。（就是说单词之间的顺序不重要）
那么模型输出的文章 越 接近 我们预先给它的文集，我们就认为这个模型质量越好。
那么如何定量的计算这个模型的好坏呢？
我们假设左边的模型生成的那篇REAL DOCUMENT的概率为P(left), 右边的为P（right） 。那么我们就可以比较谁概率大就知道哪个模型好了。
下面根据LDA的设计思想，我们就来详细说下这4个组件是怎么搭配在一起工作，产生出一篇文章，以及它的概率是什么？
首先第一个旋钮，会按照概率分布D1，随机出一个多项式分布值M1（eg. topic1:50%，topic2：30%, topic3:20%）我们告诉第一个齿轮。然后第二个旋钮按照概率分布D2，随机出K(例子中是3)个多项式分布值
如：
topic 1: {ball:70%, galaxy:25%, referendum :1%, planet : 4%}
topic 2: {ball:9%, galaxy:35%, referendum :1%, planet : 55%}
topic 3: {ball:2%, galaxy:3%, referendum :80%, planet : 15%}
随后第一个齿轮拿到M1（eg. topic1:50%，topic2：30%, topic3:20%）， 根据这个概率随机出一个topic 告诉第二个齿轮（eg. topic1）。 第二个齿轮根据这个TOPIC 1和第二个旋钮给他的这组分布，用topic 1: {ball:70%, galaxy:25%, referendum :1%, planet : 4%} 随机出一个词=》比如ball
循环上述这个过程，又可以得到一个词。这样一直输出到文章要求的词数N，那么就产生了一篇文章。
这2个旋钮控制的是一个叫狄里克雷分布的函数。而2个齿轮则是多项式分布。这里你只要知道这2个名词。也不妨碍理解第二层的结构。我会在第3层讲这些组件背后的数学和算法是什么的时候，展开介绍你需要理解第3层要知道的这2个分布的意义。那么我们就有了下图。
一个参数为alpha的狄里克雷分布 随机出 一个多项式分布叫theta，theta会随机出一个数字z. 一个参数为beta的狄里克雷分布 随机k次 拿到k个多项式分布叫phi, 单位w 会根据 第z个phi 随机得到。w的颜色是灰色，表示这个单词，我们最后是可以在文章里看到的。而别的字母的值都是在这个文章中看不到的。所以我们称它为隐式狄里克雷分配模型。我们可以看到决定这个单词概率的源头是alpha和beta。所以我们可以做的是就是转动这2个旋钮，让我们得到的文章里的词和真实的文集越接近越好。
注意之前是文章，现在是文集了。因为我们输入的是一组文章，所以我们希望也输出一组文章去比较和原来的那组文章比较。（文章间顺序无关，每篇文章内的词语顺序无关）。也就是说我们的输入是（（a,c）,(a,b)）. 模型输出((b,a),(c,a)). 那么就是完全相等的结果。
那么不考虑的顺序的概率就可以写成如下形式。
想生成和原来的文集一样的概率是。
在知道alpha下取到theta的概率 乘以 在知道theta下取到z的概率。随后是在知道beta下取到phi的概率，随后是在phi (z)下得到这个w1的概率。一篇文章有N个单词，所以每个概率连乘起来。因为phi有K个 所以K个概率练成起来。最后有M篇文章，所以这些概率也乘起来。（乘法表示彼此独立，即产生下一个词语和上一个词语是什么无关，topic和文章同理）
下面我们直观的来感受下这个公式，用可视化的方式知道这个模型的运作机制。
在这里我们把一个2维的狄里克雷分布（又叫beta分布）表示为一条线。一个样本点可以落在这条线的任何一个地方。这条线的范围区间是【0,1】。
小伙伴这边可以打开这个网站和我一起动手玩一些，一旦你理解了beta的分布和二项分布的联系。那么狄里克雷只不过是把beta分布的2维变到K维。二项分布的2项变成多项分布的K项。
2.1 理解beta分布和二项分布
我们从上面的模型介绍可以知道，狄里克雷分布可以按一定的概率分布来随机出一个多项分布。那么对应到2维的情况。等价于一个beta分布 按一定的概率分布可以随机出一个二项分布。
我们控制beta分布的旋钮有2个变量（a, b）代表2个维度。比如丢硬币的话，a代表正面向上，b代表背面向上。
beta分布，试玩https://keisan.casio.com/exec/system/1180573226； 我们先把A=1, B=1; 观察一下beta分布会长成什么样？(下面蓝色的区域是对X坐标轴缩放用的，我们暂时用不到。我们只要关于调整a, b就可以了。)
我们发现图像变成了一个直线，并且围成的面积是1. 上面图中的面积是1，所以要求一个二项分布产生的概率，我们就要对一个X的区间做积分。比如a=1,b=1的beta分布产生 二项分布在0.5~0.6的概率 为下图这个面积的大小。
如果我们想更精确到0.5的概率。我们就让delta X比较小即可。比如0.5~0.500001
二项分布在0.5~0.6的概率 的意思是 我们从这个beta分布随机出一个2项分布。
这个二项分布，（正面向上的概率为(0.5~0.6)中的一个数。背面向上的概率为1-正面向上的概率）的概率
比如我们这次用这个beta分布随机出0.55这个值。 就等价我们随机出2个topic下，每个topic的概率是多少的多项分布。为（topic1:55%, topic2:45%）
接下来会引入，共轭的概念，beta分布是二项分布的共轭分布。
大白话的意思是，我们有一个beta分布的图像。然后我们观察到一批二项分布的样本。我们把样本的数据丢进这个 beta分布，他会根据这组数据重新生成一个同时贴近数据和原来beta分布的新beta分布。
比如这个硬币，我们丢了50次（在二项分布上采样了50个样本）；那么我们就得到50个数据。其中40次正面向上，10次反面向上。
按照beta分布和二项分布的公式（式2， 数学推导在最后，现在只要记住结论理解过程即可）
就是把M1=40， M2=10， 直接作用在beta分布的2个参数上。我们得到新的图像为：
我们发现这个beta分布，就向0.8这个概率聚拢。所以他下次丢出的二项分布就的概率值更加接近（80%，20%）
随着样本增多，这个beta分布能学习的更加精确。比如我们丢了10W次，有8W次向上。beta图像会进一步收敛
值得注意的是，因为围起来的总面积始终要等于1（因为概率和为1）.因为X的可能范围不断缩小。所以Y会增大。所以我们可以看纵坐标轴的大小是在变大的。
那么总结一下：
1. 把一个二项分布产生出的数据结果 丢给一个beta分布，我们会得到一个新的beta分布。在专业术语里，就是把一个似然函数 丢给一个 先验分布，会得到一个后验分布，并且后验分布和先验分布 是 同一种分布。
2. 数据量越大，模型就越可以聚拢。当数量趋近于无穷。beta分布退化至二项分布。比如上图，就会成为在0.8上的一条线，Y在【0.8,0.8+无穷小】为无穷大。
我们理解了2维的情况。我们可以以此类推到3维。就引出了狄里克雷分布和多项式分布。其实本质是一样的。为了更好的用图像表示，我们把一个Y值交高的区域，用更深的颜色去表示。就是概率密度高的地方颜色深，然后低的地方浅。
这样我们就可以用画出一个3维的狄里克雷分布的分布情况。不过要强调的是，如果颜色的深度用一个值可以表示的话，那么区域里颜色的总和是个恒定值。就和上面的1一样。也就是颜色总量不变。我们可以让他均匀分配。也可以有些地方深（出现的概率高），有些地方浅（出现的概率低）
因为是3维了，所以我们由原来的beta分布里的a, b 2个参数。变为a1, a2, a3 3个参数。 那么我们就可以从下面的等边三角形里，随机出一个点。这个点到3个顶点的距离的远近分别代表对应3种可能的概率。越近代表概率越高。
如下图，
可以看到有个点在体育和科学的距离这里5,5开了。
当然我们的分布也可以是，
你们就会发现，3个角的颜色很深，就代表点出现在角落的概率就很高。这个就好比一个派对。每个角落都有一个正向的刺激，那么人们就会倾向于聚集到各个角落。
如果是危险，人们就倾向与原理角落。
当然上述2个例子，都是中心对称的。也有可能有一群热爱火的民族。会让这个黑色区域偏向某一个角落。
我们可以看到三角形中黄色的点，就是随机出来的一个3项分布。这里比如某个点可能是（{topic1: 90%, topic2: 1%, topic3: 9%} ）， 在图中具体的位置，就是离角落1很近。离角落2和角落3都很远。
如果我们有了大量样本数据：如这组文集里。有1000个TOPIC1， 10个TOPIC2， 10个TOPIC3. 我们可以让原本颜色均匀的三角形，会根据新到来的数据调整自己的颜色分布。
这个调整的目的就是为了这个三角形下次随机出来的结果，能够更加接近于样本集。
这里就和上面串起来了，训练模型，其实就是把一组文章当成输入，而获得一个更接近样本的三角形（狄里克雷分布）
根据我们之前的例子，我们设定了3个主题（体育，社科，科学）， 我们有4个词（planet, galaxy, ball, referendum）那么从主题到词的那个狄里克雷分布，我们就可以用下图来表示。这张图恒定只有3个点，因为我们有3个主题。
最后我们把那个公式可视化就变成了下图。后面2个方盒就是多项分布。3个颜色的球各1个再箱子里，代表{1/3， 1/3, 1/3}的多项分布
所以整个模型的运作，就是现在第一个三角形李根据颜色分布的概率随机出一个球的位置，得到topic一个分布数值
第二个正四面体，根据颜色密度，随机出3个球，得到TOPIC到词的概率。
上面这2个概率值，就能构建出对应的BOX。然后针对每次要随机的词，我们先选TOPIC，然后在到对应的TOPIC->词的BOX里选词。
最后构建出了整个文章。
上面就是LDA模型的全过程。
下面我们就是要了解我们如何训练模型，找到那个best setting？
在深入介绍训练之前我们先
2.2 暂停回顾一下
到目前为止，我们介绍了LDA的模型。
希望你可以
- 自己画出我们上面贴过的那张概率图模型的样子。也就是从alpha, beta为源头 最后生成到w的过程。并且可以自己写一下每个词的概率的计算公式。
- 在思考一下， beta分布和二项分布的关系。beta分布如何随机出一个二项分布。二项分布的数据是怎么反向更新beta分布的。那么扩展到K项，本质是一样的。
-
结合我们的第一节的顶层设计，再理解一下，我们的输入其实扮演的是多项分布的数据作用。但是他不是直接能反向更新2个狄里克雷分布的。因为他的比例，反应的是文章到词的多项分布。我们需要的是文章到主题的多项分布，和 主题到词的多项分布。这中间有一个主题在这里。那么我们应该怎么根据我们先有的（文章到词的样本） 转换 为我们需要的 （文章到主题的样本，主题到词的样本）才是我们真正的目标。
而根据狄里克雷分布和多项分布的公式。有了样本，我们只要加回去，就能得到新的参数了。这不是难点。
image.png
下一章我就会介绍LDA训练里最核心的吉布斯采样，他是怎么根据
（文章到词的样本），来达成对（文章到主题的样本，主题到词的样本）。 而有了这些样本，我们又是如何得到我们顶层设计里提到的输出。
输出：（每一篇都会给一个文章到TOPIC的概率分布。 每一个TOPIC，都会告诉你这个TOPIC下每一个单词的概率分布。）
3. 这些组件背后的数学和算法是什么，这些原理间的联系是什么 （高粒度组件下用了哪些数学算法工具）
我们先来看一下我们手上有哪些组件了。
- 输入（一组文章，一个K）
- 输出 (每篇文章到主题的分布， 每个主题到单词的分布）
- 文章到主题的多项分布生成器-狄里克雷分布alpha（旋钮1）
- 主题到词汇的多项分布生成器-狄里克雷分布beta (旋钮2）
- m个文章到主题的多项式分布（齿轮1）
- k个主题到词的多项式分布（齿轮2）
- 吉布斯采样算法（内置算法，用于获得文章到主题的样本，主题到词的样本）
3，4号组件背后的数学： 狄里克雷分布 可以 根据概率密度 生成 多项式分布
5,6 号组件背后的数学： 如果我们有文章下主题的样本数据，我们可以优化3号组件。如果我们有主题下词的样本数据，我们可以优化4号组件。
如果3，4号组件优化的越好，那么根据概率公式，我们产生的文档所构成的词袋（顺序无关的词的集合）和原始输入本章集合的词袋越接近。
因为topic的定义是由原始输入集定义的。所以来了一篇新的文章，我们的模型参数更加符合原始的词袋，也就意味着预测它基于文章里的词，属于的topic有哪些就更准确。
3.1下面就是如何训练的问题
如果没有吉布斯采样，我们会怎么解决这个问题。
我们现在手里有的是，文章-》单词的样本数据。我们希望有的是文章-》主题，和 主题-》单词的样本数据。因为我们希望反向更新的旋钮的定义是文章-》主题，主题-》单词。
那么我们的目标就是可以估算出
p(z | w) 的概率。 这个式子的定义是，已知词是什么的情况下，它是z topic 的条件概率是多少。
比如我们有一个文章 （ball, galaxy）。 我们希望知道 第一个ball 它是（topic1: x%, topic2: y%, topic3: d%）
p(z | w) 如果word 是这篇文章的第一个单词,z=2，他的值就是y%
随后我们在W和Z上加一个剪头，表示W是向量，那么Z也是一个向量了。也就是已知多个单词，那么每个单词分别是Z向量里对应的z的概率是多少。
我们希望得到这个是因为 有了这个概率分布，我们就可以用这个分布 对每个单词去采样。那么每篇文章的每个单词 就会得到 一个他是属于哪个TOPIC的样本数据。比如上面的例子文章 （ball, galaxy）， 我们可以运用p(z|w) 采样得到 (ball:topic1, galaxy: topic 3)
那么我们有M篇文章，假设每篇文章N个词。我们就可以统计出M*N个词，哪些属于TOPIC1， 哪些属于TOPIC2， 哪些属于TOPIC3。然后看属于TOPIC1的词，分别BALL占几个，GALAXY占几个。 这样我们就有了TOPIC 到 词的多项分布采样数据了。
同时我们可以分别统计每篇文章下，这个词是TOPIC1，那么就TOPIC1 的计数+1. 下一个词是TOPIC2，那么TOPIC2的计数+1. 最后我们可以得到这篇文章如果100个单词，那么会统计得到（topic1: 90, topic2: 5, topic3: 5）的文章到主题的单词分布。
而这2个采样数据正是我们需要的。
所以核心是要知道p(z | w) 。
如果我们直接去求 我们会怎么做？
首先我们把前面的公式
用向量来表示得到如下公式
p(z) 表示了第一个旋钮和第一个齿轮做的事情，也就是根据文章的aplha 狄里克雷分布，得到theta的多项分布，在根据这个多项分布对每一个词选一个主题Z的过程。
p(w|z) 表示了我们前一个公式的第二个旋钮和第二个齿轮做的事情。就是根据beta的狄里克雷分布，随机出phi的多项分布，并且在知道Z是哪个的条件下，随机出word是哪个的过程。
那么现在我们的式子里有了p(w|z) 我们想要的p(z | w). 那么我就想到了贝叶斯公式，可以把已知的和要求的做一个交换。
套公式得到：
p(z|w) = p(w|z) * p(z) / p(w)
根据上上图的公式化简，也就是
p(z|w) = p(w,z) / p(w)
p(w,z) 我们可以根据定义把它给求出来。（式1， 数学推导在最后，现在只要记住结论理解过程即可）
但是p(w) 非常不好求。原因是
上面应该还是比较好理解的，文章中某一个单词的概率，因为K个TOPIC 都有一定概率会生成这个单词。所以需要分别求每个TOPIC下生成这个的概率，然后求和。
我们假设第一个单词的概率（p1k1 + p1k2 + p1k3）。 那么如果一共这篇文章有2个单词，就会有
（p1k1 + p1k2 + p1k3）* (p2k1 + p2k2 + p2k3)
我们展开就会有9个项。3^2
如果有N个单词的话，K个TOPIC的话。就可以有K^N 的项。这个离散状态空间太大，无法列举。
我们试了正统方法去算p（z|w）遇到困难。
所以我们需要用一个别的方法来采样到p(z|w)
下面就是MCMC方法的介绍，以及MCMC算法家族中吉布斯采样具体解决这个问题的思路。
3.2 MCMC 的工作原理
这块部分其他介绍LDA的文章，都有比较好的说明。我就简要讲解最核心的部分，不详细展开。如果需要深入掌握这块的，可以参考别的文章。
整个逻辑思路是这样的。首先在计算机中我们用线性同余法可以很容易生成均匀分布的随机数。可是有一些分布，我们没有方法可以很好的采样（从均匀分布变换过去比较困难）。这时就需要用复杂的随机模拟的方法来生成样本。
其核心1就是马尔科夫链，这种链的特性就是基于当前状态，就可以推出下一状态。并且大多数马氏链在满足平稳条件下，最终是会收敛的。
核心2就是蒙特卡罗方法，我们通过自己给出一个建议概率转移矩阵，使得这个建议概率产生的新状态在根据接受率有些采纳，有些放弃。来获得符合原来分布的样本集。
这就是MCMC的核心2个概念。
3.2.1 理解马氏链会收敛的性质（理解可跳过）
上面的矩阵就是概率转移矩阵。这个矩阵的作用就是我们当前在某一个状态。比如我现在在I状态，我们可以找到矩阵中对应的那一行，就能知道我下一个状态是I,J,K,L的概率分别是多少。根据这个概率，我就可以去到下一个状态。
有这么一种情况，我的位置向量在若干次转移后会达到一个稳定状态。
举个例子，我最开始在一个（25%是i, 25% 是j, 25%是k, 25%是l）经过前面的矩阵转换。我下一个状态可能是（30%， 15%，25%， 30%）。在经过前面的矩阵，会继续变化自己的状态。可是经过了足够多次之后，我的状态会收敛到一个固定的值上。
下面给一段python代码，小伙伴可以自己做下实验，如果转移矩阵不变。你们随便初始化初始状态，最后的结果都会是个定值。
import numpy as np
input = np.array([[0.05, 0.70, 0.25]])
matrix = np.matrix([[0.9, 0.075, 0.025],
[0.15, 0.8, 0.05],
[0.25, 0.25, 0.5]])
for i in range(100):
input = np.dot(input, matrix)
print(input)
....
[[ 0.625 0.3125 0.0625]]
[[ 0.625 0.3125 0.0625]]
[[ 0.625 0.3125 0.0625]]
[[ 0.625 0.3125 0.0625]]
[[ 0.625 0.3125 0.0625]]
[[ 0.625 0.3125 0.0625]]
[[ 0.625 0.3125 0.0625]]
那么insight就是，转移矩阵决定最后的概率分布和初始值无关。那么如果一个概率分布不好求，我们可以通过构造一个马尔科夫链的转移矩阵，使得最后收敛的这个概率分布正好等于我们想求的这个概率分布，即可。
重点，整体概率分布不好求，但是具体给定某个状态，我们还是要能够求这个状态的产生的概率
其实到这里，我们依然会问，原始概率分布不好求，那收敛到这个概率分布的转移矩阵就好求了吗？
3.2.2 理解蒙特卡罗思想（理解可跳过）
这里需要用到1个技巧，这就是第二个MC，也就是蒙特卡罗算法。我们通过自己给出一个建议转移矩阵。再根据原来的当前状态的概率值和自己的建议的转移概率值，可以算出一个接受率。就可以知道我们随机出来的新状态，是保留还是放弃。
而蒙特卡罗的思想我们可以举个例子，我们想求圆周率pi的值是多少。直接求困难的话，我们可以用均匀分布随机出很多在范围【0,1】【0,1】的（x,y）。只对满足（x^2 + y^2 <= 1）的点留下来。最后比如我们丢了1W个点。计数的只有X个
那么4 * X / 10000就是PI的值。因为计数的是半径为1的1/4圆的面积 ，总量是面积变长为1的正方形面积
python代码如下：
import random
cnt = 0
samples = 1000000
for i in range(samples):
x = random.random()
y = random.random()
if x*x + y*y <= 1:
cnt += 1
print(4 * cnt / samples)
3.2.3 理解MH算法（理解可跳过）
那么我们要构造的就是一个采样算法。那么构造这样的算法的洞察在于只有满足一个叫detailed balance 的马尔科夫链才会收敛到一个平稳分布。而我们需要这个平稳分布。所以我们需要去构造。MH是构造平稳分布的其中一种算法
这里的pi 就是最后已经收敛的分布。在我们上文给出的程序代码里是
[[ 0.625 0.3125 0.0625]]
所以pi1 = 0.625, pi2 = 0.3125, pi3 = 0.0625
转移矩阵是
[[0.9, 0.075, 0.025],
[0.15, 0.8, 0.05],
[0.25, 0.25, 0.5]]
那么假设i = 1, j = 2（1开始计数）
对于任意的I,J。如果detailed balance 都满足。我们就可以得到一个式子
我觉得上图把I 理解为前一个状态，J理解为后一个状态，还是比较好理解的。
pi * P = pi ， 也向我们说明了，转移矩阵和最后我们要逼近的平稳分布间的关系是什么
那么我们先从一个一般的算法讲起叫MH算法（吉布斯采样是它的特例）
MH 算法的目的：是根据一个需求的(desired distribution)概率分布 P(x)生成一系列样本状态点（因此，这个算法可以生成任意的概率分布）。为了达到这个目的，该算法使用马尔可夫过程去到达一个平稳分布pi(x), 以使pi(x)=P(x)。
根据下面的做法，可以满足平稳分布（式3， 数学推导在最后，现在只要记住结论理解过程即可）
第一个红圈我们需要根据我们提供的建议概率g， 随机出一个新状态。
后面那个红圈，我们是要根据建议概率g，和原来的我们要逼近的概率分布p，去算出一个概率值。这也是上文提到的重点。
我们没法在该分布下采样。但是我们可以根据某个样本算出它的概率。
这就是MH算法。
还有一个容易混淆的点，上文提及的丢弃。不是说完全放弃这次迭代。而是如果没有在接受率范围里。就使用旧的样本作为当前这次采样的样本。所以迭代1000次，接受了600次，还是会有1000个样本
3.2.4 理解吉布斯采样
而吉布斯算法，就是这个接受率变为了永远的1. 也就是不再会放弃一些新样本。当然这就要求这个算法的建议概率为原概率分布的完全条件概率。
其实用大白话说很容易，如果一个状态是N维的，这个完全条件概率，就是只看其中1维，在别的维度都已知的条件下，这个我看中的这维的概率分布会是怎么样。
然后再这个基础上对这维的数据进行采样，并且更新。
把所有维度都求一遍，才得到一个新的样本
FOR里面所有行即采样一个样本要做的事情。N为采样N个样本
算法如下
因为接受概率在这个情况下为1.并且上述做法依然可以保证detailed condition （式4， 数学推导在最后，现在只要记住结论理解过程即可）
对我们的影响就是：
- 不再需要自己给一个建议概率g
- 不再需要计算原有的概率分布下的概率值P（也是平稳分布pi）
- 需要我们可以计算完全条件概率分布值，并采样。
3.2.5 吉布斯采样完全条件概率如何计算
你理解完这段，你就把前面所有知识点串起来了 这里面涉及很多数学推导，我都只给结论。我们先抛开细节，关于结论和结论之间的逻辑关系。最后在回到每个结论的推导证明上。
步骤1. 找到LDA下的要求的完全条件概率。
我们的初心是得到 p(z | w) 也就是在已知单词的情况下知道每个单词的主题分布是什么。
那么我们继续思考，LDA下的维度是什么。因为完全条件概率，就是在维度上推演的。
我们可以发现M篇文章，每篇文章的N个单词。都是可以变换TOPIC的。所以这里的维度是M * N维。
那么基于上述2点， 我们
我们知道p(w, z) 我们根据LDA模型的定义是可以算出来的。那么这个概率应该怎么算呢？
步骤二， 想办法计算这个概率
（式5， 数学推导在最后， 只给结论之间的关系，具体证明在最后）
给出狄里克雷的概率密度函数
基于上述概率分布。我们可以得到p(w,z)形式如下：（式6， 数学推导在最后）
有了p(w, z)的公式，我们就可以 扩展到 p(w, z（除去某一维度））
维度是某一篇文章下的某一个单词。
最后可以根据p(w,z)推得
那么这个概率已经变成我们可以计算的形式了。
因为先验狄里克雷分布的alpha和beta不影响最后模型的结果，所以我们可以把它们全部取为同一个值，这样我们就可以把上面的式子转换为下面的式子
拿beta分布来举例，其实就是alpha 和 beta取等值。在k维狄里克雷分布下，就是alpha 1 ~k 都是一个值。这种所有超参数都为同一个值称为对称超参数
步骤三， 有了这个概率值，我们就可以采样了。
因为马尔科夫链和初始状态是什么无关，所以我们可以用均匀分布随机给所有维度一个初始的状态（topic 号）
那么现在采样的步骤就是，比如我们在处理文章1的第一个单词，假设是ball。
我们先要计算他在topic1下的概率值是多少。
那就是看除了它之外的其他所有的单词的4个统计量。
分别是乘号左边的2个。
就是在全文本集中，topic1有多少个ball词是分子。 和 topic1下一共有多少个词是分母。
乘号右边的2个。
在第一篇文章中，topic1的词有多少个。 第一篇文章一共有多少个词。
那么整个采样过程如下(java代码）：
private int multiSamplingAcc(int m, int n) {
// do multinomial sampling via cumulative method:
double[] p = new double[K];
for (int k = 0; k < K; k++) {
// nw[w][k] 第k个topic下有多少个w号单词， nwsum[k] 第k个topic下一共有多少个单词
// nd[m][k] 第m篇文章里有多少个属于topic k的单词， ndsum[m] 第m篇文章里有多少个单词
p[k] = (nw[docs[m][n]][k] + beta) / (nwsum[k] + V * beta)
* (nd[m][k] + alpha) / (ndsum[m] + K * alpha);
}
for (int k = 1; k < p.length; k++) {
p[k] += p[k - 1];
}
double u = Math.random() * p[K - 1];
int k_new;
for (k_new = 0; k_new < p.length; k_new++) {
if (u <= p[k_new])
break;
}
z[m][n] = k_new;
return k_new;
}
步骤四，那么我们就可以一直采样，直到收敛。得到我们想要的输出
还记得顶层我们程序需要的输出吗？
每一篇都会给一个文章到TOPIC的概率分布。 每一个TOPIC，都会告诉你这个TOPIC下每一个单词的概率分布。
那么吉布斯采样在收敛完以后，我们手里有的是。一套数据样本。这个样本是每个词分别属于哪个topic。
我们在来看一下，一开始我们有的是文章-》单词的样本数据。
我们希望有的是文章-》主题，和 主题-》单词的样本数据。因为我们希望反向更新的旋钮的定义是文章-》主题，主题-》单词。
那么现在每个词都有一个主题。所以我们可以统计这个文章下，每个主题所占的比例（根据文章内的词所占的比例）； 已经在所有文集里，这个主题下的所有单词（每个单词所占的比例）那么我们就得到了主题到单词的多项分布。
这样我们就得到了我们的输出。没错此时我们已经通过这2个输出，作为似然函数的数据值，把我们最开始的LDA模型的2个狄里克雷分布的旋钮给训练好了。我们不需要更新alpha和beta的数值，因为有先验和数据也是等价的，并且我们不需要模型真的去随机产生文章。而是要对新的文章去做预测。
上述是一个直观的理解。其实这背后也有数学上的推导和证明。我们把这个称为（式7，在最后统一处理）
那么我们来看一下这块对应的代码
// nw[w][k] 第k个topic下有多少个w号单词， nwsum[k] 第k个topic下一共有多少个单词
// nd[m][k] 第m篇文章里有多少个属于topic k的单词， ndsum[m] 第m篇文章里有多少个单词
// V为distinct的文集单词集合数量
// K为指定的主题数量
// M为文章数量
// 得到每个主题下，每个词的概率
public double[][] getPhi() {
double[][] phi = new double[K][V];
for (int k = 0; k < K; k++) {
for (int t = 0; t < V; t++) {
phi[k][t] = (nw[t][k] + beta) / (nwsum[k] + V * beta);
}
}
return phi;
}
// 得到每个文章下，每个主题的概率分布
public double[][] getTheta() {
double[][] theta = new double[M][K];
for (int m = 0; m < M; m++) {
for (int k = 0; k < K; k++) {
theta[m][k] = (nd[m][k] + alpha) / (ndsum[m] + K * alpha);
}
}
return theta;
}
3.2.6 回顾一下全貌
我们现在已经掌握可以编写LDA代码的全部逻辑链条和组件。我们来回顾一下。
- 最顶层上，我们明确了输入是list<list<words>>, K; 输出是double[K][V] phi (代表主题到词的概率分布）
double[M][K] theta(代表文章到主题的分布 ） - 模型的好坏 由2个旋钮决定，分别是文章到主题的alpha的狄里克雷分布，和主题到词的beta狄里克雷分布；随后会由2组齿轮去产生单词，分别是文章到主题的多项分布，主题到词的多项分布）
- 根据上述4个组件，我们可以算出在已知p(w, z|a,b) 的概率；我们需要的是p(z |w,a,b) 的概率， 因为我们要知道每个词具体的主题才可以采样得到概率期望值作为顶层模型的输出
- 在直接计算的时候，遇到了p(w) 项数指数爆炸不好计算。所以我们曲线救国，引入了MCMC的gibbs 采样算法
- gibbs采样算法，需要用到完全条件概率。这个可以由p(w,z) 去推导出来。
- 有了完全条件概率的值的计算方法，我们就可以知道当前每一个维度，转换到下一个状态，这个维度的新的值是基于什么概率转移方程去做的。那么我们就能采样到一个新的值对这个维度来说。
- 这样依次对每一个维度，基于其他维度的信息去采样。做完之后，我们就得到一个全新的采样样本。这样做直到gibbs收敛。
- gibbs收敛后，那么现在的每个词上的TOPIC分布，就满足了和文章输入最密切相关的联合概率分布 p(w, z). 那么我们可以根据狄里克雷概率的期望公式拿到顶层设计的2个输出。
4. 我们如何基于这些原理去实现自己的LDA
4.1 根据上面的算法， 定义出我们需要的数据结构。
public class GibbsSampler {
int K, V, M;
double beta = 1, alpha = 1;
// nw[:][k] that of term observation counts for topic k
// topic–term count
int[][] nw;
// topic–term sum
int[] nwsum;
// document–topic count
int[][] nd;
// document–topic sum
int[] ndsum;
// Z(m,n) => m doc, n word 's topic idx
int[][] z;
// docs[m][n]
int[][] docs;
Random random = new Random();
4.2 写构造函数处理输入
这个时候我们要传进来输入了。这里需要对每个单词做一个唯一编号。也就是给每个word一个递增id，来加速运算。（否则要用map<string>）， 在外层，我们需要得到这个编号ID的最大值+1之后，就是distinct 单词集合的数量，作为V传进来
public GibbsSampler(int[][] docs, int k, int v) {
this.docs = docs;
K = k;
M = docs.length;
V = v;
}
4.3 gibbs 采样框架书写
init 是随机初始化，因为我们的最后的结果和初始的概率分布无关。所以随机初始化，把几个值统计一下就好了。
在gibbs采样的每一个维度下，都要先排除这个维度的数据统计值，这个时候需要调用
decrementCntAndSum
multiSamplingAcc 就是上文提到了完全条件概率具体值的计算和运用。返回出一个这个维度下的新topic
在得到一个这个维度下的新的topic的时候，要把统计值更新回来。
伪代码为
init()
for 一个足够大使得gibbs可以收敛的迭代次数
for 每一个维度
排除当前维度的统计信息
用完全条件概率公式的计算方式去采样得到一个新的topic
更新当前维度的topic和统计信息
把模型写到磁盘里
public void gibbs() {
init();
for (int i = 0; i < ITERATIONS; i++) {
for (int m = 0; m < M; m++) {
int N = docs[m].length;
for (int n = 0; n < N; n++) {
int k = z[m][n];
decrementCntAndSum(m, k, n);
int k_new = multiSamplingAcc(m, n);
incrementCntAndSum(m, k_new, n);
}
}
}
}
4.4 实现辅助函数
private void init() {
nw = new int[V][K];
nwsum = new int[K];
nd = new int[M][K];
ndsum = new int[M];
z = new int[M][];
for (int m = 0; m < M; m++) {
int N = docs[m].length;
z[m] = new int[N];
for (int n = 0; n < N; n++) {
int k = random.nextInt(K);
z[m][n] = k;
incrementCntAndSum(m, k, n);
}
}
}
private void incrementCntAndSum(int m, int k, int n) {
nd[m][k]++;
ndsum[m]++;
int t = docs[m][n];
nw[t][k]++;
nwsum[k]++;
}
private void decrementCntAndSum(int m, int k, int n) {
nd[m][k]--;
ndsum[m]--;
int t = docs[m][n];
nw[t][k]--;
nwsum[k]--;
}
4.5 实现核心完全条件概率函数采样算法
private int multiSamplingAcc(int m, int n) {
// do multinomial sampling via cumulative method:
double[] p = new double[K];
for (int k = 0; k < K; k++) {
// nw[w][k] 第k个topic下有多少个w号单词， nwsum[k] 第k个topic下一共有多少个单词
// nd[m][k] 第m篇文章里有多少个属于topic k的单词， ndsum[m] 第m篇文章里有多少个单词
p[k] = (nw[docs[m][n]][k] + beta) / (nwsum[k] + V * beta)
* (nd[m][k] + alpha) / (ndsum[m] + K * alpha);
}
for (int k = 1; k < p.length; k++) {
p[k] += p[k - 1];
}
double u = Math.random() * p[K - 1];
int k_new;
for (k_new = 0; k_new < p.length; k_new++) {
if (u <= p[k_new])
break;
}
z[m][n] = k_new;
return k_new;
}
4.6 计算输出的2个数组
// 得到每个主题下，每个词的概率
public double[][] getPhi() {
double[][] phi = new double[K][V];
for (int k = 0; k < K; k++) {
for (int t = 0; t < V; t++) {
phi[k][t] = (nw[t][k] + beta) / (nwsum[k] + V * beta);
}
}
return phi;
}
// 得到每个文章下，每个主题的概率分布
public double[][] getTheta() {
double[][] theta = new double[M][K];
for (int m = 0; m < M; m++) {
for (int k = 0; k < K; k++) {
theta[m][k] = (nd[m][k] + alpha) / (ndsum[m] + K * alpha);
}
}
return theta;
}
所以整个训练的代码 也就100多行，代码量其实不大。
4.7 写预测文章的主题分布代码
预测其实就是还是在同样的LDA模型下做预测。那么我们就可以假定原来训练出来的参数是最佳的。也就是主题到词的概率。所以我们需要的输出是针对这篇文章的文章到主题的概率分布。
而我们有的输入，是这篇文章的文章到词。
那么我们就单独针对这篇文章用gibbs采样，来获得这篇文章的p(z|w). 其中算法的一个乘法因子（文章到主题的）是根据文章自己的词语分布特性重新计算的。另一个因子（主题到词的），我们就用模型算好的phi数组即可。
那么模型收敛以后，我们就可以用样本数据求期望获得对文章到主题的分布了。
也就是在最上面的我们画出来的狄里克雷三角形中，再找到一个最合适的黄色小球在三角形中的位置。
public class LdaPredicator {
public static final int ITERATIONS = 1000;
public static double[] predict(double[][] phi, double alpha, double beta, int[] doc) {
int K = phi.length;
int V = phi[0].length;
int[][] nw = new int[V][K];
int[] nd = new int[K];
int[] nwsum = new int[K];
int ndsum = 0;
// init
int N = doc.length;
int[] z = new int[N]; // z_i := 1到K之间的值，表示马氏链的初始状态
for (int n = 0; n < N; n++) {
int k = (int) (Math.random() * K);
z[n] = k;
nw[doc[n]][k]++;
nd[k]++;
nwsum[k]++;
}
ndsum = N;
for (int i = 0; i < ITERATIONS; i++) {
for (int n = 0; n < N; n++) {
int topic = z[n];
nw[doc[n]][topic]--;
nd[topic]--;
nwsum[topic]--;
ndsum--;
// do multinomial sampling via cumulative method: 通过多项式方法采样多项式分布
double[] p = new double[K];
for (int k = 0; k < K; k++) {
p[k] = phi[k][doc[n]]
* (nd[k] + alpha) / (ndsum + K * alpha);
}
// cumulate multinomial parameters 累加多项式分布的参数
for (int k = 1; k < p.length; k++) {
p[k] += p[k - 1];
}
// scaled sample because of unnormalised p[] 正则化
double u = Math.random() * p[K - 1];
for (topic = 0; topic < p.length; topic++) {
if (u < p[topic])
break;
}
// add newly estimated z_i to count variables 将重新估计的该词语加入计数器
nw[doc[n]][topic]++;
nd[topic]++;
nwsum[topic]++;
ndsum++;
z[n] = topic;
}
}
double[] theta = new double[K];
for (int k = 0; k < K; k++) {
theta[k] = (nd[k] + alpha) / (ndsum + K * alpha);
}
return theta;
}
}
最后写一个main函数，大家自己来玩一下吧
import java.util.Arrays;
import java.util.HashMap;
import java.util.Map;
public class LdaMain {
private int[][] docs;
private Map<String, Integer> w2id = new HashMap<>();
private Map<Integer,String> id2w = new HashMap<>();
private GibbsSampler gibbsSampler;
public LdaMain(String[] oriDocs, int k) {
int m = oriDocs.length;
docs = new int[m][];
for (int i = 0; i < m; i++) {
String[] ws = oriDocs[i].split(" ");
docs[i] = new int[ws.length];
int j = 0;
for (String w : ws) {
w2id.putIfAbsent(w, w2id.size());
id2w.putIfAbsent(w2id.get(w), w);
docs[i][j++] = w2id.get(w);
}
}
gibbsSampler = new GibbsSampler(docs, k, w2id.size());
}
public static void main(String[] args) {
String[] docs = {
"ball ball ball planet galaxy ball ball ball ball ball ball ball ball ball ball",
"referendum planet planet referendum referendum referendum referendum referendum referendum referendum referendum referendum referendum referendum referendum",
"planet planet galaxy planet ball planet galaxy planet galaxy planet galaxy planet galaxy planet planet planet planet",
"planet galaxy referendum planet ball planet galaxy planet galaxy planet galaxy planet galaxy planet galaxy planet planet planet planet"};
LdaMain main = new LdaMain(docs, 3);
main.train();
main.predict("planet galaxy planet galaxy ball planet galaxy planet galaxy planet galaxy planet galaxy planet planet planet planet planet planet");
}
private void predict(String doc) {
double[][] phi = gibbsSampler.getPhi();
String[] ws = doc.split(" ");
int[] ddoc = new int[ws.length];
for (int i = 0; i < ws.length; i++) {
ddoc[i] = w2id.get(ws[i]);
}
printTopic2Word(phi);
double[] theta = LdaPredicator.predict(phi, gibbsSampler.alpha, gibbsSampler.beta, ddoc);
System.out.println(Arrays.toString(theta));
}
private void printTopic2Word(double[][] phi) {
for (int i = 0; i < phi.length; i++) {
System.out.println("topic " + i + ": ");
for (int j = 0; j < w2id.size(); j++) {
System.out.print(id2w.get(j) + ":" + phi[i][j] + " ");
}
System.out.println();
}
}
private void train() {
gibbsSampler.gibbs();
}
}
5. 背后的数学详细推导 （每个数学结论是怎么得来的）
这篇文章的重点我已经讲完了。我们需要数学推导的地方，我都留了（式X）的字样， 应该是有7个，需要数学推导一下。
我们把这块，当做读后作业，交给读者自己去寻找。所有证明都可以在LDA漫游指南里找到。
当然有余力的小伙伴，还可以看看LDA数学八卦。
看懂数学推导需要一定的微积分基础和概率统计基础。小伙伴们加油！
感谢
看到最后的小伙伴，非常不容易。这个模型的原理还是有不小的门槛。我也不确定哪里有笔误，或者理解不准确的地方。
不过把整个逻辑链条理清楚了还是很畅快的。
你们投入了时间希望你们也可以收获到了这个LDA模型的魅力。如果有任何疑问，可以在文章下留言。thank you for your time and suggestion~
