这一节我们的目的是利用Gibbs采样去估计LDA的参数$\theta$和$\varphi$,如果直接对联合概率分布$p(W,Z,\theta,\varphi\mid\alpha,\beta)$进行采样估计,算法复杂且没必要,所以LDA通常采用一种收缩的Gibbs采样方法,基本想法是,通过对参数$\theta$和$\varphi$积分,得到边缘概率分布$p(W,Z\mid\alpha,\beta)$,其中$W$是我们所有文本的观测序列,$Z$是其对应的隐变量;然后我们对后验概率分布$p(Z\mid W,\alpha,\beta)$进行Gibbs采样,得到分布$p(Z\mid W,\alpha,\beta)$的样本集合,最后再利用这个样本集合对参数$\theta$和$\varphi$进行估计
首先:
$$ p(Z\mid W,\alpha,\beta)=\frac{p(W,Z\mid\alpha,\beta)}{p(W\mid \alpha,\beta)}\propto p(W,Z\mid\alpha,\beta) $$由于变量$W,\alpha,\beta$已知,所以分母项都相同,可以不予考虑,接下来对联合概率分布$p(W,Z\mid\alpha,\beta)$的表达式按上一节的概率图模型进一步分解:
$$ p(W,Z\mid\alpha,\beta)=p(W\mid Z,\alpha,\beta)p(Z\mid\alpha,\beta)=p(W\mid Z,\beta)p(Z\mid\alpha) $$所以问题转化为对$p(W\mid Z,\beta),p(Z\mid\alpha)$这两个表达式的推导
对于$p(W\mid Z,\beta)$的计算,根据概率图模型,又可以拆解为两部分,首先有:
$$ p(W\mid Z,\varphi)=\prod_{k=1}^K\prod_{v=1}^V\varphi_{kv}^{n_{kv}} $$其中,$\varphi_{kv}$是第$k$个话题生成单词集合的第$v$个单词的概率,$n_{kv}$是数据中第$k$个话题生成第$v$个单词的次数。于是:
$$ p(W\mid Z,\beta)=\int p(W\mid Z,\varphi)p(\varphi\mid\beta)d\varphi\\ =\int \prod_{k=1}^K\frac{1}{B(\beta)}\prod_{v=1}^V\varphi_{kv}^{n_{kv}+\beta_v-1}d\varphi \\ =\prod_{k=1}^K\frac{1}{B(\beta)}\int \prod_{v=1}^V\varphi_{kv}^{n_{kv}+\beta_v-1}d\varphi \\ =\prod_{k=1}^K\frac{B(n_k+\beta)}{B(\beta)} $$其中,$n_k=\{n_{k1},n_{k2},...,n_{kV}\}$,$B(\beta)=\int\prod_{i=1}^V\theta_i^{\beta_i-1}d\theta$
对于$p(Z\mid\alpha)$的推导与上面类似,首先:
$$ p(Z\mid\theta)=\prod_{m=1}^M\prod_{k=1}^K\theta_{mk}^{n_{mk}} $$其中,$\theta_{mk}$是第$m$个文本生成第$k$个话题的概率,$n_{mk}$是数据中第$m$个文本生成第$k$个话题的次数,于是:
$$ p(Z\mid\alpha)=\int p(Z\mid\theta)p(\theta\mid\alpha)d\theta\\ =\int\prod_{m=1}^M\frac{1}{B(\alpha)}\prod_{k=1}^K\theta_{mk}^{n_{mk}+\alpha_k-1}d\theta\\ =\prod_{m=1}^M\frac{1}{B(\alpha)}\int\prod_{k=1}^K\theta_{mk}^{n_{mk}+\alpha_k-1}d\theta\\ =\prod_{m=1}^M\frac{B(n_m+\alpha)}{B(\alpha)} $$其中,$n_m=\{n_{m1},n_{m2},...,n_{mK}\}$,$B(\alpha)=\int\prod_{i=1}^K\theta_i^{\alpha_i-1}d\theta$
接下来,我们组合上面的两个因子:
$$ p(Z,W\mid\alpha,\beta)=\prod_{k=1}^K\frac{B(n_k+\beta)}{B(\beta)}\cdot \prod_{m=1}^M\frac{B(n_m+\alpha)}{B(\alpha)} $$所以:
$$ p(Z\mid W,\alpha,\beta)\propto\prod_{k=1}^K\frac{B(n_k+\beta)}{B(\beta)}\cdot \prod_{m=1}^M\frac{B(n_m+\alpha)}{B(\alpha)} $$分布$p(Z\mid W,\alpha,\beta)$的满条件分布可以写作:
$$ p(z_i\mid Z_{-i},W,\alpha,\beta)=\frac{1}{Z_{z_i}}p(Z\mid W,\alpha,\beta) $$这里,$w_i$表示所有文本的单词序列的第$i$个位置的单词,$z_i$表示单词$w_i$对应的话题,$i=(m,n),i=1,2,...,M,n=1,2,...,N_m$,$Z_{-i}=\{z_j:j\neq i\}$,$Z_{z_i}$表示分布$p(Z\mid W,\alpha,\beta)$对变量$z_i$的边缘化因子,所以:
$$ p(z_i\mid Z_{-i},W,\alpha,\beta)\propto \frac{n_{kv}+\beta_v}{\sum_{v=1}^V(n_{kv}+\beta_v)}\cdot \frac{n_{mk}+\alpha_k}{\sum_{k=1}^K(n_{mk}+\alpha_k)} $$其中,第$m$个文本的第$n$个位置的单词$w_i$是单词集合的第$v$个单词,其话题$z_i$是话题集合的第$k$个话题,$n_{kv}$表示第$k$个话题中第$v$个单词的计数,但减去当前单词的计数,$n_{mk}$表示第$m$个文本中第$k$个话题的计数,但减去当前单词的话题的计数。
接下来便是通过Gibbs抽样得到的分布$p(Z\mid W,\alpha,\beta)$的样本去估计变量$\theta,\varphi$了
根据LDA图模型,后验概率满足如下方程
$$ p(\theta_m\mid Z_m,\alpha)=\frac{1}{Z_{\theta_m}}\prod_{n=1}^{N_m}p(z_{mn}\mid\theta_m)p(\theta_m\mid\alpha)=Dir(\theta_m\mid n_m+\alpha) $$这里,$n_m=\{n_{m1},n_{m2},...,n_{mK}\}$是第$m$个文本的话题的计数,$Z_{\theta_m}$表示分布$p(\theta_m,Z_m\mid\alpha)$对变量$\theta_m$的边缘化因子。于是得到参数$\theta$的估计式:
$$ \theta_{mk}=\frac{n_{mk}+\alpha_k}{\sum_{k=1}^K(n_{mk}+\alpha_k)},m=1,2,...,M;k=1,2,...,K $$同样的,通过LDA图模型,可以得到下面的关系:
$$ p(\varphi_k\mid W,Z,\beta)=\frac{1}{Z_{\varphi_k}}\prod_{i=1}^Ip(w_i\mid\varphi_k)p(\varphi_k\mid\beta)=Dir(\varphi_k\mid n_k+\beta) $$这里$n_k=\{n_{k1},n_{k2},...,n_{kV}\}$是第$k$个话题的单词的计数,$Z_{\varphi_k}$表示分布$p(\varphi_k,W\mid Z,\beta)$对变量$\varphi_k$的边缘化因子,$I$是文本集合单词的序列$W$的单词总数。于是得到参数的估计:
$$ \varphi_{kv}=\frac{n_{kv}+\beta_v}{\sum_{v=1}^V(n_{kv}+\beta_v)},k=1,2,...,K;v=1,2,...,V $$通过前面的推导,我们可以梳理出LDA的Gibbs抽样算法流程了
输入:文本的单词序列$W=\{w_1,...,w_m,...,w_M\},w_m=(w_{m1},...,w_{mn},...,w_{mN_m})$;
输出:文本的话题序列$Z=\{z_1,...,z_m,...,z_M\},z_m=(z_{m1},...,z_{mn},...,z_{mN_m})$的后验概率分布$p(Z\mid W,\alpha,\beta)$的样本计数,模型的参数$\varphi$和$\theta$的估计值;
参数:超参数$\alpha$和$\beta$,话题个数$K$;
(1)设所有计数矩阵的元素$n_{mk},n_{kv}$,计数向量的元素$n_m,n_k$初始值为0;
(2)对所有文本$w_m,m=1,2,...,M$
对第$m$个文本中的所有单词$w_{mn},n=1,2,...,N_m$抽样其话题$z_{mn}=z_k\sim Mult(\frac{1}{K})$;
增加文本-话题计数$n_{mk}=n_{mk}+1$;
增加文本-话题和计数$n_m=n_m+1$;
增加话题-单词计数$n_{kv}=n_{kv}+1$;
增加话题-单词和计数$n_k=n_k+1$;
(3)循环执行以下操作,直到进入燃烧期
对所有文本$w_m,m=1,2,...,M$
对第$m$个文本中的所有单词$w_{mn},n=1,2,...,N_m$
(a)当前的单词$w_{mn}$是第$v$个单词,话题指派$z_{mn}$是第$k$个话题:减少计数$n_{mk}=n_{mk}-1,n_m=n_m-1,n_{kv}=n_{kv}-1,n_k=n_k-1$;
$$ p(z_i\mid Z_{-i},W,\alpha,\beta)\propto \frac{n_{kv}+\beta_v}{\sum_{v=1}^V(n_{kv}+\beta_v)}\cdot \frac{n_{mk}+\alpha_k}{\sum_{k=1}^K(n_{mk}+\alpha_k)} $$(b)按照如下的满条件分布进行抽样,得到新的第$k'$个话题,将其分配给$z_{mn}$
(c)增加计数$n_{mk'}=n_{mk'}+1,n_m=n_m+1,n_{k'v}=n_{k'v}+1,n_{k'}=n_{k'}+1$
(d)得到更新的两个计数矩阵$N_{K\times V}=[n_{kv}]$和$N_{M\times K}=[n_{mk}]$,表示后验概率分布$p(Z\mid W,\alpha,\beta)$的样本计数;
$$ \theta_{mk}=\frac{n_{mk}+\alpha_k}{\sum_{k=1}^K(n_{mk}+\alpha_k)},m=1,2,...,M;k=1,2,...,K\\ \varphi_{kv}=\frac{n_{kv}+\beta_v}{\sum_{v=1}^V(n_{kv}+\beta_v)},k=1,2,...,K;v=1,2,...,V $$(4)利用得到的样本计数估计模型参数:
"""
隐狄利克雷分布的代码实现,包括Gibbs采样和变分EM算法,代码封装在ml_models.latent_dirichlet_allocation
"""
import numpy as np
class LDA(object):
def __init__(self, alpha=None, beta=None, K=10, tol=1e-3, epochs=100):
"""
:param alpha: 主题分布的共轭狄利克雷分布的超参数
:param beta: 单词分布的共轭狄利克雷分布的超参数
:param K: 主题数量
:param tol:容忍度,允许tol的隐变量差异
:param epochs:最大迭代次数
"""
self.alpha = alpha
self.beta = beta
self.K = K
self.tol = tol
self.epochs = epochs
self.theta = None # 文本-主题矩阵
self.phi = None # 主题-单词矩阵
def _init_params(self, W):
"""
初始化参数
:param W:
:return:
"""
M = len(W) # 文本数
V = 0 # 词典大小
I = 0 # 单词总数
for w in W:
V = max(V, max(w))
I += len(w)
V += 1 # 包括0
# 文本话题计数
N_M_K = np.zeros(shape=(M, self.K))
N_M = np.zeros(M)
# 话题单词计数
N_K_V = np.zeros(shape=(self.K, V))
N_K = np.zeros(self.K)
# 初始化隐状态,计数矩阵
Z = [] # 隐状态,与W一一对应
p = [1 / self.K] * self.K
hidden_status = list(range(self.K))
for m, w in enumerate(W):
z = np.random.choice(hidden_status, len(w), replace=True, p=p).tolist()
Z.append(z)
for n, k in enumerate(z):
v = w[n]
N_M_K[m][k] += 1
N_M[m] += 1
N_K_V[k][v] += 1
N_K[k] += 1
# 初始化alpha和beta
if self.alpha is None:
self.alpha = np.ones(self.K)
if self.beta is None:
self.beta = np.ones(V)
return Z, N_M_K, N_M, N_K_V, N_K, M, V, I, hidden_status
def fit(self, W):
"""
:param W: 文本集合[[...],[...]]
:return:
"""
Z, N_M_K, N_M, N_K_V, N_K, M, V, I, hidden_status = self._init_params(W)
for _ in range(self.epochs):
error_num = 0
for m, w in enumerate(W):
z = Z[m]
for n, topic in enumerate(z):
word = w[n]
N_M_K[m][topic] -= 1
N_M[m] -= 1
N_K_V[topic][word] -= 1
N_K[topic] -= 1
# 采样一个新k
p = [] # 更新多项分布
for k_ in range(self.K):
p_ = (N_K_V[k_][word] + self.beta[word]) * (N_M_K[m][k_] + self.alpha[topic]) / (
(N_K[k_] + np.sum(self.beta)) * (N_M[m] + np.sum(self.alpha)))
p.append(p_)
ps = np.sum(p)
p = [p_ / ps for p_ in p]
topic_new = np.random.choice(hidden_status, 1, p=p)[0]
if topic_new != topic:
error_num += 1
Z[m][n] = topic_new
N_M_K[m][topic_new] += 1
N_M[m] += 1
N_K_V[topic_new][word] += 1
N_K[topic_new] += 1
if error_num / I < self.tol:
break
# 计算参数theta和phi
self.theta = N_M_K / np.sum(N_M_K, axis=1, keepdims=True)
self.phi = N_K_V / np.sum(N_K_V, axis=1, keepdims=True)
def transform(self, W):
rst = []
for w in W:
tmp = np.zeros(shape=self.K)
for v in w:
try:
v_ = self.phi[:, v]
except:
v_ = np.zeros(shape=self.K)
tmp += v_
if np.sum(tmp) > 0:
tmp = tmp / np.sum(tmp)
rst.append(tmp)
return np.asarray(rst)
接下来,笔者就使用https://github.com/zhulei227/Text_Representation 项目中的部分数据来做测试吧
docs=[
["有","微信","红包","的","软件"],
["微信","支付","不行","的"],
["我们","需要","稳定的","微信","支付","接口"],
["申请","公众号","认证"],
["这个","还有","几天","放","垃圾","流量"],
["可以","提供","聚合","支付","系统"]
]
word2id={}
idx=0
W=[]
for doc in docs:
tmp=[]
for word in doc:
if word in word2id:
tmp.append(word2id[word])
else:
word2id[word]=idx
idx+=1
tmp.append(word2id[word])
W.append(tmp)
W
[[0, 1, 2, 3, 4], [1, 5, 6, 3], [7, 8, 9, 1, 5, 10], [11, 12, 13], [14, 15, 16, 17, 18, 19], [20, 21, 22, 5, 23]]
lda=LDA(epochs=1000)
lda.fit(W)
trans=lda.transform(W)
#第二句和第三句应该比较近似,因为它们都含有“微信”,“支付”
trans[1].dot(trans[2])
0.2026956065114651
#而第二句和第四句的相似度显然不如第二句和第三句
trans[1].dot(trans[3])
0.037617554858934164
#当然第二句和第五句的差距也有些大
trans[1].dot(trans[4])
0.025078369905956115
#而第一句和第二句都含有“微信”,所以相似度会比第四、五句高,同时比第三句小
trans[1].dot(trans[0])
0.1355252606255012
测试结果基本符合期望...,大家可以在更大的语料库上做测试(当然训练速度会比较慢就是了,并没有做优化~~~)