アカリの部屋

主题模型(LDA算法)

主题模型概念

直观的说,一个企业要招聘工程师,想通过看一大把简历来分出小白和大牛,那么文本里就能找到很多特征,比如“百度实习”、“XX大学第一名”、“有GitHub账号”等,那么HR做判断的依据就是,如果很多好的码农都有某一些特征,而你的简历里正好也有,那你是个好码农的概率就高。

即:特征组成简历,简历可以通过HR通过经验学习分出好坏。

一个HR如果没有先验知识、不知道这些特征,那么他只能遍历所有的简历、遍历每份简历里所有的特征,去猜每个特征是不是一个好工程师的特征,对每个特征他都会有一个自信值来判断这个人怎么样,随着时间增加,他会通过经验统计来修正他的猜测,最终他会学到一个好的判断模型。

为什么他这样做就是对的呢?因为:

1
2
一份简历中 每个特征的存在 都是因为 这个人本身 有一定概率 落在好/坏分类当中
并且 这个人写简历的过程 是从好/坏这个分类 各自对应的特征中 以一定的概率 选择了某些特征

这里我们的“特征”组成了“简历”,而通过“简历”可以分辨出“好坏”,这就是一种主题模型(LDA)
主题模型是指通过“单词”可以组成“文档”,通过“文档”可以分辨出“主题”。
这是一种无监督的贝叶斯模型,它的特点是,对于每一个“主题”都能找到一个“单词”来描述它。
LDA是一种词袋模型,即不考虑词的前后关系,只考虑出现频次。

LDA概念

我们之所以可以给一篇文章分主题是因为:

1
2
一篇文章的每个词 都是以一定概率 对应到了某个主题
并且从这个主题对应的词的范围中 以一定的概率 选择了这些词

LDA模型要做的事,就是要得到这样的分布,它要找的是:
1.当文章的主题一定时,能找到词的分布,即哪些词会冒出来,能以高概率判定文章属于这个主题
2.当这篇文章一定时,能找到主题的分布,即各个单词,对应到哪些主题的概率高

使用公式来表示:

$ P(单词|文档) = P(单词|主题) * P(主题|文档) $

LDA学习

ps.来自:http://blog.sina.com.cn/s/blog_8eee7fb60101czhx.html

在算法开始时,对所有 $ d $ 和 $ t $ 随机给定 $ θd $ 和 $ φt $,然后对一个特定文档 $ ds $ 中的某个单词 $ wi $ 假定为某一个主题 $ tj $,有 $ pj(wi|ds)=p(wi|tj)*p(tj|ds) $,枚举所有主题,可以得到所有 $ pj(wi|ds) $

现在我们可以枚举T中的t,得到所有的 $ pj(wi|ds) $。然后可以根据这些概率值结果为 $ ds $ 中的第 $ i $ 个单词 $ wi $ 选择一个 $ t $。最简单的想法是取令 $ pj(wi|ds) $ 最大的 $ tj $。

然后,如果 $ ds $ 中的第 $ i $ 个单词 $ wi $ 在这里选择了一个与原先不同的 $ t $,就会对$ θd $ 和 $ φt $有影响了。它们的影响又会反过来影响对上面提到的 $ p(w|d) $ 的计算。对 $ D $ 中所有的 $ d $ 中的所有 $ w $ 进行一次 $ p(w|d) $ 的计算并重新选择 $ t $ 看作一次迭代。这样进行n次循环迭代之后,就会收敛到LDA所需要的结果了。

马尔可夫链

在之前的概率图模型笔记里有对马尔可夫模型的描述,其中有一个假设说状态转移矩阵是不会变化的,举个例子,对于天气在晴和阴两种情况中的转化,假设晴->晴是0.9,晴->阴是0.1,晴->晴和阴->阴都是0.5,那么就能得到一个矩阵

$ P = ( \begin{matrix} 0.9 & 0.1 \\ 0.5 & 0.5 \\ \end{matrix} ) $

假设今天是晴,那么明天就是

$ ( \begin{matrix} 1 & 0 \end{matrix} ) * P $

那么,由于 $ P $ 是稳定不变的,n天后的概率就会收敛到

$ q = lim (\begin{matrix} 1 & 0 \end{matrix}) P^n = (\begin{matrix} 0.833 & 0.167 \end{matrix})$

这样一来,这个状态转移矩阵就成为了一个稳定的向量,这个固定分布叫做稳态分布

不论初始状态如何,维度多大,它的每一行相同,即

$ q = p^∞ = ( \begin{matrix} 0.833 & 0.167 \\ 0.833 & 0.167 \\ \end{matrix} ) $

且有

$ ( \begin{matrix} 1 & 0 \end{matrix} ) P^∞ = ( \begin{matrix} 0 & 1 \end{matrix} ) P^∞ $

这样一来,每个状态就都是 $ q $ 的一个样本了。

ps.统计知识

贝叶斯模型就是概率的交换律:$ P(θ|y)P(y) = P(y|θ)P(θ) $

在简历的例子里,可以说,$ P(好工程师|简历) ∝ P(好工程师)P(简历|好工程师) $,其中 $ P(好工程师) $ 叫先验, $ P(好工程师|简历) $ 叫后验, $ P(简历|好工程师) $ 叫似然。

世界上搞概率的有频率学派和贝叶斯学派两种,频率学派就单纯的认为P=A/B,但是贝叶斯学派认为概率是有先验和后验区别的,先验就是事情没发生的时候我猜它是什么样,后验就是事情发生后当我观测到事情什么样之后我再说是什么样。

比如,频率学派认为投硬币正面的概率是1/2,而贝叶斯学派认为猜出来的先验是1/2,而后验则是做实验抛了100观察到100次正面后,把后验概率修正得更大。

假设D是一组观测实验,x是个随机变量,$ P(x) $ 代表先验概率, $ μ $ 表示实际实验当中的结果到底如何,那么 $ P(x|D) ∝ P(μ|D) P(x) $, 为了让正比符号变等号,通过贝叶斯公式归一化,有 $ P(x|D) = P(μ|D) P(x)/P(D) $, 其中 $ P(μ|D) $ 叫似然函数,贝叶斯学派认为因为环境影响,抛硬币的结果可能不是1/2,可以通过似然函数调整先验概率。

MCMC算法(蒙特卡洛马尔可夫链)和Gibbs采样

在实际计算当中,先验好猜,似然可以观察,但是后验难算,为了解决这个问题可以构造一个无限状态马尔可夫链,让它的稳态为 $ q(X{_1},…,X{_m}) $ ,这样在足够多步数之后,就可以采几步估出分布$ q $的样子。

这个算法的一个直观的解释就是,拿中国地图来说,因为东西海拔不一样,那么在高的地方扔一个球,这个球的下落自然会被牵引(收敛)到低的地方,经过无数多次之后,就能拿到海拔的分布了。

MCMC的含义就是,第二个MC表示马尔科可链,暗示含有稳态分布,第一个MC是蒙特卡洛,即随机过程,可以通过随机过程搞到稳态分布。

假设我们知道多个条件分布,但是不知道联合分布,就可以通过Gibbs采样的方式逼近真实的联合分布。做法是,假定一个初始组合,然后根据条件概率改变其中的变量,形成一个比较长的马尔可夫链,之后,跳过一定条数(比如前100条),然后每间隔几个点(比如20个点)采样一次,这样得到的分布是逼近联合分布的。

为什么LDA算法用于主题分析可行呢?因为我们认定了文档中的单词和主题的关系是有客观的分布规律的,也就是一个马尔可夫链的隐藏关系,即存在一个能够收敛到稳定分布的矩阵。所以才能轮流去找两个分布,不断的去迭代来估计这个联合分布矩阵。

所以,可以通过抽样的方式计算两个条件概率,之后更新分布 $ q $,拿这个分布再采样,再计算,反复迭代,最终这个分布就会收敛逼近实际的文档到话题的分布 $ p $。

变分(KL-divergence)

这个过程相当于把自己的 $ Q $ 拟合到 $ P $,那么就需要一个所谓损失函数定义二者相差多少,这个度量衡就叫变分,也叫相对熵,它是用来衡量两个分布的相似度的度量。

对于这两个分布,变分的表述形式为:

离散型:$ KL(p{_A}||p{_B}) = \sum{_i}p{_A}^{(i)} log \frac{p{_A}^{(i)}}{p{_B}^{(i)}} $
连续型:$ KL(P||Q) = \int p(x)log \frac{p(x)}{q(x)}dx $

因为log前面乘上了用于去拟合到某个标准的那个分布,所以这个公式不太care这个分布中密度比较低的地方,也就是说,比如是两个正态分布,两边的部分权重就低。

另外这个函数是非对称的,$ KL(P||Q) $ 表示用P逼近Q,而 $ KL(Q||P) $ 表示用Q逼近P,二者不等。如果 $ P $ 是个双峰,那么前者会去把 $ Q $ 拉成一个比较宽的单峰,而后者只会去很好的拟合其中一个峰而忽略另一个峰。在机器学习中通常采用后者 $ KL(P||Q) $,这是ML的一个局限,因为更多的时候想要的就是局部最优解。

LDA的变分推断就是要用变分分布 $ p $ 去逼近 $ q $,这就是Gensim里的实现方式。

手写GibbsLDA算法

设定真理

首先我们要作为上帝,造一些东西出来:

假设词库里只有5个单词:monkey,loan,bank,river,stream,其中bank又表示银行又表示岸滩。假设只有两个话题T::T1,T2 一个谈钱一个谈水。设置一个真理,前三个单词在话题1的可能性各是1/3,后三个单词在话题2的可能性也各是1/3。最终我们要拟合出来的单词到话题的分布一定是逼近上面两个分布的,需要用一个phi矩阵保存这个结果。然后要生成语料,假设要做16个文档,每个文档平均单词个数是10,用泊松分布给每个不同的文档设置句子有多长,然后从概率分布里抽单词出来组成句子。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#coding=utf-8
import numpy as np
np.random.seed(0)
# 假设词库里只有5个单词 money,loan,bank,river,stream 其中bank又表示银行又表示岸滩
vocab = ["monkey", "loan", "bank", "river", "stream"]
# 假设只有两个话题T: T1,T2 一个谈钱一个谈水
T = 2
# 设置一个真理,前三个单词在话题1的可能性各是1/3,后三个单词在话题2的可能性也各是1/3
z_1 = np.array([1.0/3, 1.0/3, 1.0/3, 0.0, 0.0])
z_2 = np.array([0.0, 0.0, 1.0/3, 1.0/3, 1.0/3])
# 最终我们要拟合出来的单词到话题的分布,一定是逼近上面两个分布的
# 把两个话题分布转化成phi矩阵
phi_actual = np.array([z_1, z_2]).T.reshape(len(z_2), 2)
# 生成语料
# 假设要做16个文档,每个文档平均单词个数是10
# 用泊松分布给每个不同的文档设置句子有多长
# 然后从概率分布里抽单词出来组成句子
D = 16
mean_length = 10
len_doc = np.random.poisson(mean_length, size=D)
docs = []
orig_topics = []
for i in range(D):
z = np.random.randint(0, 2)
if z == 0:
words = np.random.choice(vocab, size=(len_doc[i]), p=z_1).tolist()
else:
words = np.random.choice(vocab, size=(len_doc[i]), p=z_2).tolist()
orig_topics.append(z)
docs.append(words)
# 如果打印出这些句子,我们就能看出每个句子长什么样,和它应该属于哪个话题
print docs[3]
# out: ['bank', 'bank', 'bank', 'river', 'bank', 'stream', 'river', 'bank', 'river']
print orig_topics[3]
# out: 1

统计表现

做完了上帝的事情就要做GibbsLDA了,但是在这之前先要做一些初始化:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# 初始化变量
w_i = [] # 单词在V中的位置
d_i = [] # 文档位置
z_i = [] # topic分布
i = []
counter = 0 # 循环计数
# 走一遍每个文档
for doc_idx, doc in enumerate(docs):
# 当中的每个单词
for word_idx, word in enumerate(doc):
# 找到目前这个单词在V中的位置
w_i.append(np.where(np.array(vocab) == word)[0][0])
# 记录i和文档位置
i.append(counter)
d_i.append(doc_idx)
# 最后初始化随机的它属于哪个topic的分布
z_i.append(np.random.randint(0, T))
counter += 1
# 为了统一,把数组转成numpy数组
w_i = np.array(w_i)
d_i = np.array(d_i)
z_i = np.array(z_i)

之后就要初始化WT和DT两个分布了,这两个分布可以通过计数拿到:

1
2
3
4
5
6
7
8
9
10
11
12
13
# WT记录每个单词出现在某个topic里多少次
WT = np.zeros((len(vocab), T))
for idx, word_ in enumerate(vocab):
topics = z_i[np.where(w_i == idx)]
for t in range(T):
WT[idx, t] = sum(topics == t)
# DT记录每个文档属于某个topic多少次
DT = np.zeros((D, T))
for idx, doc_ in enumerate(range(D)):
topics = z_i[np.where(d_i == idx)]
for t in range(T):
DT[idx, t] = sum(topics == t)

以DT来说,它的结果是:

1
2
3
4
5
6
[[ 8. 2.]
[ 6. 5.]
[ 3. 6.]
[ 5. 4.]
[ 8. 10.]
...

也就是说,对第一篇文档,其中有8个单词属于话题1,2个单词属于话题2。

打印出WT可以看到

1
2
3
4
5
[[ 11. 14.]
[ 20. 7.]
[ 29. 28.]
[ 11. 17.]
[ 18. 17.]]

也就是说第一个单词money出现在第一个话题中11次,第二个话题中14次。

Gibbs采样

假设 $ w $ 代表单词,$ d $ 代表文档,$ t $ 代表主题,则 $ D $ 中每一个文档 $ d $ 都能看做一个 $ $ 的序列,$ D $ 中的所有不同单词会组成一个词汇集合 $ V $。LDA以 $ D $ 作为输入,希望训练出两个向量:

1.对于 $ D $ 中的每个 $ d $,对应到 $ T $ 中不同 $ t $ 的概率 $ θd $,其中 $ pti $ 表示 $ d $ 对应 $ T $ 中第 $ i $ 个 $ t $ 的概率。这里 $ pti = nti/n $,$ nti $ 表示 $ d $ 中对应第 $ i $ 个 $ t $ 的词的数目, $ n $ 是 $ d $ 中所有词的总数。

2.对于 $ T $ 中的每个 $ t $,生成不同单词的概率 $ φt $,其中 $ pwi $ 表示 $ t $ 生成 $ V $ 中第 $ i $ 个单词的概率。$ pwi = Nwi/N $,其中 $ Nwi $ 表示对应到 $ t $ 的 $ v $ 中第 $ i $ 个单词的数目,$ N $ 表示所有对应高 $ t $ 的单词总数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
# 记录每个phi的变化结果
phi_1 = np.zeros((len(vocab), 100))
phi_2 = np.zeros((len(vocab), 100))
# 跑100次
iters = 100
# 定义先验分布参数
alpha = 1.0
beta = 1.0
for step in range(iters):
for current in i:
# 把D和W分别拿出来
doc_idx = d_i[current]
w_idx = w_i[current]
# 并把这两个从总体集合中减去
DT[doc_idx, z_i[current]] -= 1
WT[w_idx, z_i[current]] -= 1
# 计算新的W和D的分布
prob_word = (WT[w_idx, :] + beta) / (WT[:, :].sum(axis=0) + len(vocab) * beta)
prob_document = (DT[doc_idx, :] + alpha) / (DT.sum(axis=0) + D * alpha)
# 这其实就是对于每个topic的概率
prob = prob_word * prob_document
# 把Z更新(根据刚刚求得的prob)
z_i[current] = np.random.choice([0, 1], 1, p=prob / prob.sum())[0]
# 更新计数器
DT[doc_idx, z_i[current]] += 1
WT[w_idx, z_i[current]] += 1
# 记录下Phi的变化
phi = WT / (WT.sum(axis=0))
phi_1[:, step] = phi[:, 0]
phi_2[:, step] = phi[:, 1]
# 得到后验值
phi = WT/(WT.sum(axis=0))
print phi

可以得到

1
2
3
4
5
[[ 0.3164557 0. ]
[ 0.34177215 0. ]
[ 0.34177215 0.32258065]
[ 0. 0.30107527]
[ 0. 0.37634409]]

会发现这个学习结果和预先设定的值几乎是一样的。

归一化后输出预测结果:

1
2
3
theta = DT/DT.sum(axis=0)
theta = theta/np.sum(theta, axis=1).reshape(16,1)
print np.argmax(theta, axis=1) == orig_topics

得到

1
2
[ True True True True True True True True True True True True
True True True True]