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