LDA

Latent Dirichlet

Dirichlet Distribution & Dirichlet Process:狄利克雷分布于狄利克雷过程

作者:肉很多链接:https://www.zhihu.com/question/26751755/answer/80931791来源:知乎著作权归作者所有,转载请联系作者获得授权。要想易懂地理解dirichlet distribution,首先先得知道它的特殊版本 beta distribution 干了什么。而要理解 beta distribution 有什么用,还得了解 Bernoulli process。

首先先看Bernoulli process。要理解什么是 Bernoulli process,首先先看什么 Bernoulli trial。Bernoulli trial 简单地说就是一个只有两个结果的简单 trial,比如*抛硬币*。 那我们就用抛一个(不均匀)硬币来说好了,X = 1 就是头,X = 0 就是字,我们设定 q 是抛出字的概率。 那什么是 bernoulli process?就是从 Bernoulli population 里随机抽样,或者说就是重复的独立 Bernoulli trials,再或者说就是狂抛这枚硬币 n 次记结果吧(汗=_=)。好吧,我们就一直抛吧,我们记下 X=0 的次数 k.

现在问题来了。 Q:我们如何知道这枚硬币抛出字的概率?我们知道,如果可以一直抛下去,最后 k/n 一定会趋近于 q;可是现实中有很多场合不允许我们总抛硬币,比如我只允许你抛 4 次。你该怎么回答这个问题?显然你在只抛 4 次的情况下,k/n 基本不靠谱;那你只能"猜一下 q 大致分布在[0,1]中间的哪些值里会比较合理",但绝不可能得到一个准确的结果比如 q 就是等于 k/n。

举个例子,比如:4 次抛掷出现“头头字字”,你肯定觉得 q 在 0.5 附近比较合理,q 在 0.2 和 0.8 附近的硬币抛出这个结果应该有点不太可能,q = 0.05 和 0.95 那是有点扯淡了。 你如果把这些值画出来,你会发现 q 在[0,1]区间内呈现的就是一个中间最高,两边低的情况。从感性上说,这样应当是比较符合常理的。

那我们如果有个什么工具能描述一下这个 q 可能的分布就好了,比如用一个概率密度函数来描述一下? 这当然可以,可是我们还需要注意另一个问题,那就是随着 n 增长观测变多,你每次的概率密度函数该怎么计算?该怎么利用以前的结果更新(这个在形式上和计算上都很重要)?

到这里,其实很自然地会想到把 bayes theorem 引进来,因为 Bayes 能随着不断的观测而更新概率;而且每次只需要前一次的 prior 等等…在这先不多说 bayes 有什么好,接下来用更形式化语言来讲其实说得更清楚。

我们现在用更正规的语言重新整理一下思路。现在有个硬币得到 random sample X = (x1,x2,...xn),我们需要基于这 n 次观察的结果来估算一下q 在[0,1]中取哪个值比较靠谱,由于我们不能再用单一一个确定的值描述 q,所以我们用一个分布函数来描述:有关 q 的概率密度函数(说得再简单点,即是 q 在[0,1]“分布律”)。当然,这应当写成一个条件密度:f(q|X),因为我们总是观测到 X 的情况下,来猜的 q。

现在我们来看看 Bayes theorem,看看它能带来什么不同: P(q|x) P(x) = P(X=x|q)P(q)

在这里 P(q)就是关于 q 的先验概率(所谓先验,就是在得到观察 X 之前,我们设定的关于 q 的概率密度函数)。P(q|x)是观测到 x 之后得到的关于 q 的后验概率。注意,到这里公式里出现的都是"概率",并没有在[0,1]上的概率密度函数出现。为了让贝叶斯定理和密度函数结合到一块。我们可以从方程两边由 P(q)得到 f(q),而由 P(q|x)得到 f(q|x)。 又注意到 P(x)可以认定为是个常量(Q:why?),可以在分析这类问题时不用管。那么,这里就有个简单的结论——**关于 q 的后验概率密度 f(q|x)就和“关于 q 的**先验概率密度乘以一个条件概率"成比例,即: f(q|x)\sim P(X=x|q)f(q)

带着以上这个结论,我们再来看这个抛硬币问题: 连续抛 n 次,即为一个 bernoulli process,则在 q 确定时,n 次抛掷结果确定时,又观察得到 k 次字的概率可以描述为:P(X=x|p) = q^{k}(1-q)^{n-k} 那么 f(q|x)就和先验概率密度乘以以上的条件概率是成比例的: f(q|x) \sim q^{k}(1-q)^{n-k}f(q) 虽然我们不知道,也求不出那个 P(x),但我们知道它是固定的,我们这时其实已经得到了一个求 f(q|x)的公式(只要在 n 次观测下确定了,f(q)确定了,那么 f(q|x)也确定了)。

现在在来看 f(q)。显然,在我们对硬币一无所知的时候,我们应当认为硬币抛出字的概率 q 有可能在[0,1]上任意处取值。f(q)在这里取个均匀分布的密度函数是比较合适的,即 f(q) = 1 (for q in [0,1]) 。 有些同学可能发现了,这里面f(q|x) \sim q^{k}(1-q)^{n-k}那个q^{k}(1-q)^{n-k}乘上[0,1]的均匀分布不就是一个 Beta distribution 么? 对,它就是一个 Beta distribution。Beta distribution 由两个参数 alpha、beta 确定;在这里对应的 alpha 等于 k+1,beta 等于 n+1-k。而均匀分布的先验密度函数,就是那个 f(q)也可以被 beta distribution 描述,这时 alpha 等于 1,beta 也等于 1。

更有意思的是,当我们每多抛一次硬币,出现字时,我们只需要 alpha = alpha + 1;出现头只需要 beta = beta + 1。这样就能得到需要估计的概率密度 f(q|x)…

其实之所以计算会变得这么简单,是因为被 beta distribution 描述的 prior 经过 bayes formula 前后还是一个 beta distribution;这种不改变函数本身所属 family 的特性,叫共轭(conjugate)

ok。讲到这你应该明白,对于有两个结果的重复 Bernoulli trial,我们用 beta prior/distribution 就能解决。那么加入我们有 n 个结果呢?比如抛的是骰子? 这时候上面的 Bernoulli trial 就要变成有一次 trial 有 k 个可能的结果; Bernoulli distribution 就变成 multinomial distribution。而 beta distribution 所表述的先验分布,也要改写成一个多结果版本的先验分布。那就是 dirichlet distribution。 均匀的先验分布 Beta(1,1)也要变成 k 个结果的 Dir(alpha/K)。dirichlet prior 也有共轭的性质,所以也是非常好计算的。 简而言之,就是由 2 种外推到 k 种,而看待它们的视角并没有什么不同。 他们有着非常非常非常相似的形式。

结论 1:dirichlet distribution 就是由 2 种结果 bernoulli trial 导出的 beta distribution 外推到 k 种的 generalization

from scipy.stats import dirichlet, poisson
from numpy.random import choice
from collections import defaultdict
num_documents = 5
num_topics = 2
topic_dirichlet_parameter = 1 # beta
term_dirichlet_parameter = 1 # alpha
vocabulary = ["see", "spot", "run"]
num_terms = len(vocabulary)
length_param = 10 # xi
term_distribution_by_topic = {} # Phi
topic_distribution_by_document = {} # Theta
document_length = {}
topic_index = defaultdict(list)
word_index = defaultdict(list)
term_distribution = dirichlet(num_terms * [term_dirichlet_parameter])
topic_distribution = dirichlet(num_topics * [topic_dirichlet_parameter])
# 遍历每个主题
for topic in range(num_topics):
# 采样得出每个主题对应的词分布
term_distribution_by_topic[topic] = term_distribution.rvs()[0]
# 遍历所有的文档
for document in range(num_documents):
# 采样出该文档对应的主题分布
topic_distribution_by_document[document] = topic_distribution.rvs()[0]
topic_distribution_param = topic_distribution_by_document[document]
# 从泊松分布中采样出文档长度
document_length[document] = poisson(length_param).rvs()
# 遍历整个文档中的所有词
for word in range(document_length[document]):
topics = range(num_topics)
# 采样出某个生成主题
topic = choice(topics, p=topic_distribution_param)
topic_index[document].append(topic)
# 采样出某个生成词
term_distribution_param = term_distribution_by_topic[topic]
word_index[document].append(choice(vocabulary, p=term_distribution_param))

如果还有困惑的同学可以参考如下 Python 代码:

def perplexity(self, docs=None):
if docs == None: docs = self.docs
# 单词在主题上的分布矩阵
phi = self.worddist()
log_per = 0
N = 0
Kalpha = self.K * self.alpha
//遍历语料集中的所有文档
for m, doc in enumerate(docs):
// n_m_z 为每个文档中每个主题的单词数,theta 即是每个单词出现的频次占比
theta = self.n_m_z[m] / (len(self.docs[m]) + Kalpha)
for w in doc:
// numpy.inner(phi[:,w], theta) 即是某个出现的概率统计值
log_per -= numpy.log(numpy.inner(phi[:,w], theta))
N += len(doc)
return numpy.exp(log_per / N)

Introduction

LDA has been widely used in textual analysis,

LDA 是标准的词袋模型。

通俗理解 LDA 主题模型

LDA 主要涉及的问题包括共轭先验分布、Dirichlet 分布以及 Gibbs 采样算法学习参数。LDA 的输入为文档数目$M$,词数目$V$(非重复的 term),主题数目$K$。

Mathematics

Beta 分布:Dirichlet 分布的基础

Beta 分布的概率密度为:

f(x)={1B(α,β)xα1(1x)β1,x(0,1)0,f(x) = \left \{ \begin{aligned} \frac{1}{B(\alpha,\beta)}x^{\alpha-1}(1-x)^{\beta-1}, x \in (0,1) \\ 0,其他 \end{aligned} \right.

其中B(α,β)=01xα1(1x)β1dx=Γ(α+β)Γ(α)Γ(β)B(\alpha,\beta) = \int_0^1 x^{\alpha - 1}(1-x)^{\beta-1}dx=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} 其中 Gamma 函数可以看做阶乘的实数域的推广:

Γ(x)=0tx1etdtΓ(n)=(n1)!B(α,β)=Γ(α+β)Γ(α)Γ(β)\Gamma(x) = \int_0^{\infty}t^{x-1}e^{-t}dt \Rightarrow \Gamma(n) = (n-1)! \Rightarrow B(\alpha,\beta)=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}

Beta 分布的期望为: E(X)=α+βαE(X) = \frac{\alpha + \beta}{\alpha}

Dirichlet 分布:多项分布的共轭分布

Dirichlet 分布实际上就是把:

α=α1,β=α2,x=x1,x1=x2\alpha = \alpha_1 , \beta = \alpha_2 , x = x_1 , x - 1 = x_2
f(pα)={1Δ(α)k=1Kpkαk1,pk(0,1)0,f(\vec{p} | \vec{\alpha}) = \left \{ \begin{aligned} \frac{1}{\Delta(\vec{\alpha})} \prod_{k=1}^{K} p_k^{\alpha_k - 1} ,p_k \in (0,1) \\ 0,其他 \end{aligned} \right.

可以简记为:

Dir(pα)=1Δ(α)k=1Kpkαk1Dir(\vec{p} | \vec{\alpha}) = \frac{1}{\Delta(\vec{\alpha})} \prod_{k=1}^{K} p_k^{\alpha_k - 1}

其中

Δ(α)=k=1KΓ(αk)Γ(k=1Kαk)\Delta(\vec{\alpha}) = \frac{ \prod_{k=1}^K \Gamma(\alpha_k)}{ \Gamma(\sum_{k=1}^{K}\alpha_k)}

该部分在给定的$\vec{\alpha}$情况下是可以计算出来值的。假设给定的一篇文档有 50 个主题,那么$\vec{\alpha}$就是维度为 50 的向量。在没有任何先验知识的情况下,最方便的也是最稳妥的初始化就是将这个 50 个值设置为同一个值。

Symmetric Dirichlet Distribution(对称 Dirichlet 分布)

一旦采取了对称的 Dirichlet 分布,因为参数向量中的所有值都一样,公式可以改变为:

Dir(pα,K)=1ΔK(α)k=1Kpkα1ΔK(α)=ΓK(α)Γ(Kα)Dir(\vec{p} | \alpha,K) = \frac{1}{\Delta_K(\alpha)} \prod_{k=1}^{K} p_k^{\alpha - 1} \\ \Delta_K(\vec{\alpha}) = \Gamma^K(\alpha){ \Gamma(K * \alpha)}

而不同的$\alpha$取值,当$\alpha=1$时候,退化为均匀分布。当$\alpha>1$时候,$p1 = p_2 = \dots = p_k$的概率增大。当$\alpha<1$时候,$p_1 = 1 , p{非i} = 0$的概率增大。映射到具体的文档分类中,$\alpha$取值越小,说明各个主题之间的离差越大。而$\alpha$值越大,说明该文档中各个主题出现的概率约接近。

在实际的应用中,一般会选用$1/K$作为$\alpha$的初始值。

模型解释

上图的箭头指向即是条件依赖。

Terminology

  • 字典中共有$V$个不可重复的 term,如果这些 term 出现在了具体的文章中,就是 word。在具体的某文章中的 word 当然是可能重复的。

  • 语料库(Corpus)中共有$m$篇文档,分别是$d_1,d_2,\dots,d_m$,每篇文章长度为$N_m$,即由$N_i$个 word 组成。每篇文章都有各自的主题分布,主题分布服从多项式分布,该多项式分布的参数服从 Dirichlet 分布,该 Dirichlet 分布的参数为$\vec{ \alpha }$。注意,多项分布的共轭先验分布为 Dirichlet 分布。

    怎么来看待所谓的文章主题服从多项分布呢。你每一个文章等于多一次实验,$m$篇文档就等于做了$m$次实验。而每次实验中有$K$个结果,每个结果以一定概率出现。

  • 一共涉及到$K$(值给定)个主题,$T_1,T_2,\dots,T_k$。每个主题都有各自的词分布,词分布为多项式分布,该多项式分布的参数服从 Dirichlet 分布,该 Diriclet 分布的参数为$\vec{\beta}$。注意,一个词可能从属于多个主题。

模型过程

$\vec{\alpha}$与$\vec{\beta}$为先验分布的参数,一般会实现给定。如取 0.1 的对称 Dirichlet 分布,表示在参数学习结束后,期望每个文档的主题不会十分集中。

(1)选定文档主题

(2)根据主题选定词

参数学习

给定一个文档集合,$w{mn}$是可以观察到的已知变量,$\vec{\alpha}$与$\vec{\beta}$是根据经验给定的先验参数,其他的变量$z{mn}$、$\vec{\theta}$、$\vec{\varphi}$都是未知的隐含变量,需要根据观察到的变量来学习估计的。根据上图,可以写出所有变量的联合分布:

似然概率

一个词$w_{mn}$(即 word,可重复的词)初始化为一个词$t$(term/token,不重复的词汇)的概率是:

p(wm,n=tθm,Φ)=k=1Kp(wm,n=tϕk)p(zm,n=kθm)p(w_{m,n}=t | \vec{\theta_m},\Phi) = \sum_{k=1}^K p(w_{m,n}=t | \vec{\phi_k})p(z_{m,n}=k|\vec{\theta}_m)

上式即给定某个主题的情况下能够看到某个词的概率的总和。每个文档中出现主题$k$的概率乘以主题$k$下出现词$t$的概率,然后枚举所有主题求和得到。整个文档集合的似然函数为:

p(WΘ,Φ)=m=1Mp(wmθm,Φ)=m=1Mn=1Nmp(wm,nθm,Φ)p(W | \Theta,\Phi) = \prod_{m=1}^{M}p(\vec{w_m} | \vec{\theta_m},\Phi) = \prod_{m=1}^M \prod_{n=1}^{N_m}p(w_{m,n}|\vec{\theta_m},\Phi)

Gibbs Sampling

首先通俗理解一下,在某篇文档中存在着$N_m$个词,依次根据其他的词推算某个词来自于某个主题的概率,从而达到收敛。最开始的时候,某个词属于某个主题是随机分配的。Gibbs Sampling 的核心在于找出某个词到底属于哪个主题。

Gibbs Sampling 算法的运行方式是每次选取概率向量的一个维度,给定其他维度的变量值采样当前度的值,不断迭代直到收敛输出待估计的参数。初始时随机给文本中的每个词分配主题$z^{(0)}$,然后统计每个主题$z$下出现词$t$的数量以及每个文档$m$下出现主题$z$的数量,每一轮计算$p(zi|z{\neq i},d,w)$,即排除当前词的主题分布。 这里的联合分布:

p(w,zα,β)=p(wz,β)p(zα)p(\vec{w},\vec{z} | \vec{\alpha},\vec{\beta}) = p(\vec{w} | \vec{z},\vec{\beta})p(\vec{z} | \vec{\alpha})

第一项因子是给定主题采样词的过程。后面的因此计算,$n_z^{(t)}$表示词$t$被观察到分配给主题$z$的次数,$n_m^{(k)}$表示主题$k$分配给文档$m$的次数。

p(w,z,β)=p(wz,Φ)p(Φβ)dΦ=z=1K1Δ(β)t=1Vϕz,tnz(t)+βt1dϕz=z=1KΔ(nz+β)Δ(β),nz={nz(t)}t=1Vp(\vec{w} | ,\vec{z},\vec{\beta}) = \int p(\vec{w} | \vec{z},\vec{\Phi})p(\Phi | \vec{\beta})d \Phi \\ = \int \prod_{z=1}^{K} \frac{1}{\Delta(\vec{\beta})}\prod_{t=1}^V \phi_{z,t}^{n_z^{(t)} + \beta_t - 1}d\vec{\phi_z} \\ = \prod_{z=1}^{K}\frac{\Delta(\vec{n_z} + \vec{\beta})}{\Delta(\vec{ \beta })} , \vec{n_z} = \{ n_z^{(t)} \}_{t=1}^V
p(zα)=p(zΘ)p(Θα)dΘ=m=1M1Δ(α)k=1Kθm,knm(k)+αk1dθm=m=1MΔ(nm+α)Δ(α),nm={nm(k)}k=1Kp(\vec{z} | \vec{\alpha}) \\ = \int p(\vec{z} | \Theta) p(\Theta|\vec{\alpha}) d\Theta \\ = \int \prod_{m=1}^{M} \frac{1}{\Delta(\vec\alpha)} \prod_{k=1}^K\theta_{m,k}^{ n_m^{(k)} + \alpha_k - 1 }d\vec{\theta_m} \\ = \prod_{m=1}^M \frac{ \Delta(\vec{n_m} + \vec\alpha) }{ \Delta(\vec\alpha) }, \vec{n_m}=\{ n_m^{(k)} \}_{k=1}^K

Gibbs Updating Rule

词分布和主题分布总结

经过上面的 Gibbs 采样,各个词所被分配到的主题已经完成了收敛,在这里就可以计算出文档属于主题的概率以及词属于文档的概率了。

ϕk,t=nk(t)+βtt=1Vnk(t)+βtθm,k=nm(k)+αkk=1Knm(k)+αk\phi_{k,t} = \frac{ n_k^{(t)} + \beta_t }{ \sum^V_{t=1}n_k^{(t)} + \beta_t } \\ \theta_{m,k} = \frac{ n_m^{(k)} + \alpha_k }{ \sum^K_{k=1}n_m^{(k)} + \alpha_k } \\
p(θmzm,α)=1Zθmn=1Nmp(zm,nθmp(θmalpha))=Dir(θmnm+α)p(ϕkz,w,β)=1Zϕki:zi=kp(wiϕk)p(ϕkβ)=Dir(ϕknk+β)p(\vec{\theta_m} | \vec{z_m}, \vec{\alpha} ) = \frac{1}{Z_{\theta_m}} \prod_{n=1}^{N_m} p(z_{m,n} | \vec{\theta_m} * p(\vec{\theta_m} | \vec{alpha} )) = Dir(\vec{\theta_m} | \vec{n_m} + \vec{\alpha}) \\ p(\vec{\phi_k} | \vec{z}, \vec{w}, \vec{\beta} ) = \frac{1}{Z_{\phi_k}} \prod_{i:z_i=k} p(w_i | \vec{\phi_k}) * p(\vec{\phi_k} | \vec{\beta}) = Dir(\vec{\phi_k} | \vec{n_k} + \vec{\beta})

代码实现

代码的输入有文档数目$M$、词的数目$V$(非重复的 term)、主题数目$K$,以及用$d$表示第几个文档,$k$表示主题,$w$表示词汇(term),$n$表示词(word)。 $z[d][w]$:第$d$篇文档的第$w$个词来自哪个主题。$M$行,$X$列,$X$为对应的文档长度:即词(可重复)的数目。 $nw[w][t]$:第 w 个词是第 t 个主题的次数。word-topic 矩阵,列向量$nw[][t]$表示主题 t 的词频数分布;V 行 K 列。 $nd[d][t]$:第 d 篇文档中第 t 个主题出现的次数,doc-topic 矩阵,行向量$nd[d]$表示文档$d$的主题频数分布。M 行,K 列。 辅助向量: $ntSum[t]$:第 t 个主题在所有语料出现的次数,K 维 $ndSum[d]$:第 d 篇文档中词的数目(可重复),M 维 $P[t]$:对于当前计算的某词属于主题 t 的概率,K 维

超参数的确定

  • 交叉验证

  • $\alpha$表达了不同文档间主题是否鲜明,$\beta$度量了有多少近义词能够属于同一个类别。

  • 给定主题数目$K$,可以使用:

α=50/Kβ=0.01\alpha = 50 / K \\ \beta = 0.01