#!/usr/bin/env python # coding: utf-8 # ## 收缩的Gibbs抽样 # 这一节我们的目的是利用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)$的表达式推导 # 首先： # # $$ # 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,\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\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$个话题的计数，但减去当前单词的话题的计数。 # ### 三.参数$\theta,\varphi$估计 # 接下来便是通过Gibbs抽样得到的分布$p(Z\mid W,\alpha,\beta)$的样本去估计变量$\theta,\varphi$了 # #### 对参数$\theta$的估计 # 根据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 # $$ # # #### 对参数$\varphi$的估计 # 同样的，通过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$； # # >>> （b）按照如下的满条件分布进行抽样，得到新的第$k'$个话题，将其分配给$z_{mn}$ # $$ # 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)} # $$ # # >>> （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)$的样本计数； # # > （4）利用得到的样本计数估计模型参数： # # $$ # \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 # $$ # ### 五.代码实现 # In[1]: """ 隐狄利克雷分布的代码实现，包括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 项目中的部分数据来做测试吧 # In[2]: docs=[ ["有","微信","红包","的","软件"], ["微信","支付","不行","的"], ["我们","需要","稳定的","微信","支付","接口"], ["申请","公众号","认证"], ["这个","还有","几天","放","垃圾","流量"], ["可以","提供","聚合","支付","系统"] ] # In[3]: 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 # In[4]: lda=LDA(epochs=1000) lda.fit(W) trans=lda.transform(W) # In[5]: #第二句和第三句应该比较近似，因为它们都含有“微信”，“支付” trans[1].dot(trans[2]) # In[6]: #而第二句和第四句的相似度显然不如第二句和第三句 trans[1].dot(trans[3]) # In[7]: #当然第二句和第五句的差距也有些大 trans[1].dot(trans[4]) # In[8]: #而第一句和第二句都含有“微信”，所以相似度会比第四、五句高，同时比第三句小 trans[1].dot(trans[0]) # 测试结果基本符合期望...，大家可以在更大的语料库上做测试（当然训练速度会比较慢就是了，并没有做优化~~~） # In[ ]: