LDA(主题模型）

gibbs采样法解决LDA问题, 其中超参数$\alpha$， $\beta$各自的每个元素都相同, 便于计算

公式推导, 推荐阅读:https://www.cnblogs.com/pinard/p/6867828.html

input：主题数$T$, 词汇表$tokens$, 语料库$texts$,
output：词的主题分布$\varphi$, 文档的主体分布$\theta$
process：

1. 初始化超参数 $alpha$, $\beta$;
2. 构建一张语料库单词和所在文档的词汇l链表$word-document$, 该表的第$i$个索引为第$i$个单词所在的文档;
3. gibbs采样直到收敛:
   1. 建立两个计数矩阵:
      a) 词汇表$tokens$中的每个单词$v$所属主题的计数$C-wt$
      b) 语料库文档$texts$的每篇文档$d$所属主题计数$C-dt$
   2. 随机为语料库的每个单词$w$赋予主题编号$t$, 查询该词所属文档$d$, 并更新两个计数矩阵
      C_wt[w:t]+=1， C_dt[d:t]+=1
   3. 迭代， 遍历语料库每一个单词$w$, 查询该词所属文档$d$, 计算刨除该词之后的该文档该单词所在的主题分布, 并随机取样, 赋予该单词主题编号, 更新两个计数矩阵, 记录在每一次迭代中预料库每个单词的主题编号
   4. 计算每个单词的主题分布$\varphi$， 和每一个文档的主题分布$\theta$
      ${\varphi }_{w,t}=\frac{C-wt[w,t]+{\beta }_{wt}}{{\sum }_{v=1}^{V}C-wt[v,t]+{\beta }_{vt}}$
      ${\theta }_{d,t}=\frac{C-dt[d,t]+{\alpha }_t}{{\sum }_{s=1}^{T}C-dt[d,s]+{\alpha }_s}$

如何计算某个文d档的某个单词w的主题分布

$p(d_w=t|C-wt,C-dt)=\frac{C-dt[d,t]+{\alpha }_t}{{\sum }_{s=1}^T(C-dt[d,s]+{\alpha }_s)}\ast \frac{C-wt[w,t]+{\beta }_{wt}}{{\sum }_{v=1}^V(C-wt[v,t]+{\beta }_{vt})}$

其中要注意我们在gibbs采样时要计算刨除该词时,该文档该词的主题分布, 刨除该词体现在相应的两个计数矩阵减1

---

import numpy as np


class SmoothedLDA(object):

    def __init__(self, T, **kwargs):

        '''
        T: 主题数
        V: 词汇表数目
        N: 所有文档词的总数
        D： 文档的总数
        alpha：每一篇文档的主题先验dirichlet分布的参数 (1, T)
        beta： 每一个词的主题先验dirichlet分布参数 (V, T)
        phi： 主题的词分布
        theta： 文档的主题分布 (D, T)
        初始化两个分布内先验参数取值均是相同的，便于计算
        '''

        self.T = T  # 主题数

        self.alpha = (50.0 / self.T) * np.ones(self.T)
        if "alpha" in kwargs.keys():
            self.alpha = (kwargs["alpha"] * np.ones(self.T))

        self.beta = 0.01
        if "beta" in kwargs.keys():
            self.beta = kwargs["beta"]

    def _init_params(self, texts, tokens):
        # texts：语料库
        # tokens：词汇表
        self.tokens = tokens # 词汇表 (V,1)
        self.D = len(texts)  # 文档总数
        self.V = len(np.unique(self.tokens))  # 词汇表数目
        self.N = np.sum(np.array([len(doc) for doc in texts]))  # 所有文档总的单词数
        self.beta = self.beta * np.ones(self.V)

        self.word_document = np.zeros(self.N)  # 语料库词的文档索引

        count = 0
        for doc_idx, doc in enumerate(texts):
            for word_idx, word in enumerate(doc):
                word_idx = word_idx + count
                self.word_document[word_idx] = doc_idx
            count = count + len(doc)

    def train(self, texts, tokens, n_gibbs=2000):
        self._init_params(texts, tokens)
        C_wt, C_dt, assignments = self._gibbs_sampler(n_gibbs, texts)
        self.fit_params(C_wt, C_dt)
        return C_wt, C_dt, assignments

    def what_did_you_learn(self, top_n=10):
        for tt in range(selt.T):
            top_idx = np.argsort(self.phi[:, tt])[::-1][:top_n]
            top_tokens = self.tokens[top_idx]
            print("\nTop Words for Topic %s:\n" %str(tt))
            for token in top_tokens:
                print("\t%s\n" % str(token))

    def fit_params(self, C_wt, C_dt):
        # 更新每个主题的词分布
        # 更新每个文档的主题分布
        self.phi = np.zeros([self.V, self.T])
        self.theta = np.zeros([self.D, self.T])
        b, a = self.beta[0], self.alpha[0]

        for ii in range(self.V):
            for jj in range(self.T):
                self.phi[ii, jj] = (C_wt[ii, jj] + b) / (np.sum(C_wt[:, jj]) + self.V * b)

        for dd in range(self.D):
            for jj in range(self.T):
                self.theta[dd, jj] = (C_dt[dd, jj] + a) / (np.sum(C_dt[dd, :]) + self.T * a)
        return self.phi, self.theta

    def _estimate_topic_prob(self, ii, d, C_wt, C_dt):

        p_vec = np.zeros(self.T)
        b, a = self.beta[0], self.alpha[0]
        for jj in range(self.T):
            frac1 = (C_wt[ii, jj] + b) / (np.sum(C_wt[:, jj]) + self.V * b)
            frac2 = (C_dt[d, jj] + a) / (np.sum(C_dt[d, :]) + self.T * a)
            p_vec[jj] = frac1 * frac2
        return p_vec / np.sum(p_vec)

    def _gibbs_sampler(self, n_gibbs, texts):
        # 计数矩阵
        C_wt = np.zeros([self.V, self.T])  # 词v的第t个主题个数
        C_dt = np.zeros([self.D, self.T])  # 文档d的第t个主题个数
        assignments = np.zeros([self.N, n_gibbs + 1])  # 第k轮迭代 语料库第i个词的主题编号

        # 为语料库的每个词随机设置主题编号
        for ii in range(self.N):
            token_idx = np.concatenate(texts)[ii]  # 语料库第i个词的索引
            assignments[ii, 0] = np.random.randint(0, self.T)  # 给第i个词随机赋予主题编号

            doc = self.word_document[ii]  # 获取第i个词的文档编号
            # 如果该文档的第i个词属于第t个主题 则该文档的第t个主题个数加1
            # 如果该词属于第t个主题， 则该词的第t个主题个数加1
            C_dt[doc, assignments[ii, 0]] += 1
            C_wt[token_idx, assignments[ii, 0]] += 1

        for gg in range(n_gibbs):
            print("Gibbs iteration {} of {}".format(gg + 1, n_gibbs))

            for jj in range(self.N):
                token_idx = np.concatenate(texts)[jj]  # 语料库第j个词的索引
                doc = self.word_document[jj]   # 语料库第j个词所属文档

                C_wt[token_idx, assignments[jj, gg]] -= 1
                C_dt[doc, assignments[jj, gg]] -= 1

                # 计算文档doc某一个词token_idx的主题分布
                p_topics = self._estimate_topic_prob(token_idx, doc, C_wt, C_dt)
                # 根据分布采样一个主题, 之后更新计数矩阵
                sampled_topic = np.nonzero(np.random.multinomial(1, p_topics))[0][0]

                C_wt[token_idx, sampled_topic] += 1
                C_dt[doc, sampled_topic] += 1
                assignments[jj, gg + 1] = sampled_topic
        return C_wt, C_dt, assignments