Latent Dirichlet Allocations の Python 実装

LDA とは "Latent Dirichlet Allocation"。文書中の単語の「トピック」を確率的に求める言語モデル
潜在的ディリクレ配分法」と訳されていることもあるが、その名前だと「それってなんだっけ?」という人のほうが多そうw。


各単語が「隠れトピック」(話題、カテゴリー)から生成されている、と想定して、そのトピックを文書集合から教師無しで推定することができる。特徴は、果物の apple と音楽の apple とコンピュータ関連の apple を区別することが出来る(ことが期待される)という点。そのために、どのトピックを生成しやすいかという分布を各文章も持つ。細かい話は略。


結果の見方としては、定量的にはパープレキシティを見るし(一般に小さいほどいい)、定性的には各トピックがどのような単語を生成するか、その確率上位のものを見てふむふむする。この「各トピックが生成する単語」で音楽トピックか、スポーツトピックか、というのを判断できる。
ただし判断するのは人間。そこらへんは教師無しのつらいところ。


で、勉強を兼ねて LDA まわりを以前からあれこれ実装してみたりしているので、その第一弾として素の LDA を Python で実装したものを MIT License で公開。

必要なモジュールは numpy と nltk。
nltk はブラウンコーパスと lemmatize(原形への変換) に使用している。


実行例は次の通り。

$ ./lda.py -c 0:20 -s 0 -k 10 --alpha=0.5 --beta=0.5 -i 50

これで「ブラウンコーパスの 0 版から 19 版までを文書集合とし、stop words は除外して、トピック数は 10、トピックの事前ディリクレ分布のパラメータは α=0.5、単語の事前ディリクレ分布のパラメータは β=0.5 で、推論を50回イテレーション」して、単語の生成確率上位 20 個をトピックごとに出力する。
オプションの意味は -? とか付けてもらえば出力するので、それ見てね。
推論は collapsed Gibbs サンプリング。イテレーションごとにパープレキシティを出力するが、確率的な推論なので必ずしも毎回下がるとは限らない。


実装のキモの部分は以下の通り。70 行もないので読むのも簡単?
numpy をできるだけ使って、繰り返し部分 (inference メソッド) は特に極力高速化している。が、そうは言っても C++ とかではないので、実用的な速度は期待しないでね。いろいろ改造して実験するためのもの、という感じ。

class LDA:
    """LDA : collapsed Gibbs で推論"""
    def __init__(self, K, alpha, beta):
        self.K = K
        self.alpha = alpha # parameter of topics prior
        self.beta = beta   # parameter of words prior

    def set_corpus(self, corpus, stopwords):
        """コーパスの読み込みとカウンタの初期化"""
        voca = vocabulary.Vocabulary(stopwords==0)
        self.docs = [voca.doc_to_ids(doc) for doc in corpus]

        M = len(self.docs)
        self.V = voca.size()

        self.z_m_n = [] # topics of words of documents
        self.n_m_z = numpy.zeros((M, self.K)) + self.alpha     # word count of each document and topic
        self.n_z_t = numpy.zeros((self.K, self.V)) + self.beta # word count of each topic and vocabulary
        self.n_z = numpy.zeros(self.K) + self.V * self.beta    # word count of each topic

        self.N = 0
        for m, doc in enumerate(self.docs):
            self.N += len(doc)
            z_n = numpy.random.randint(0, self.K, len(doc))
            self.z_m_n.append(z_n)
            for t, z in zip(doc, z_n):
                self.n_m_z[m, z] += 1
                self.n_z_t[z, t] += 1
                self.n_z[z] += 1
        return voca

    def inference(self):
        """イテレーション1回分の推論を行う"""
        for m, doc in enumerate(self.docs):
            z_n = self.z_m_n[m]
            for n, t in enumerate(doc):
                # n 番目の単語 t (トピック z)についてカウンタを減算
                z = z_n[n]
                self.n_m_z[m, z] -= 1
                self.n_z_t[z, t] -= 1
                self.n_z[z] -= 1

                # トピックの事後分布からサンプリング
                p_z = self.n_z_t[:, t] * self.n_m_z[m] / self.n_z
                new_z = numpy.random.multinomial(1, p_z / p_z.sum()).argmax()

                # サンプリングされた新トピックを設定&カウンタを増加
                z_n[n] = new_z
                self.n_m_z[m, new_z] += 1
                self.n_z_t[new_z, t] += 1
                self.n_z[new_z] += 1

    def worddist(self):
        """トピック-単語分布を返す"""
        return self.n_z_t / self.n_z[:, numpy.newaxis]

    def perplexity(self):
        """パープレキシティを計算"""
        phi = self.worddist()
        log_per = 0
        Kalpha = self.K * self.alpha
        for m, doc in enumerate(self.docs):
            theta = self.n_m_z[m,:] / (len(doc) + Kalpha)
            for w in doc:
                log_per -= numpy.log(numpy.inner(phi[:,w], theta))
        return numpy.exp(log_per / self.N)

カウンタにあらかじめパラメータαやβなどが加算されていて、inference() ではベクトルの掛け算割り算だけで事後分布が計算できているところがミソ。
そうしない方がアルゴリズムとしては見通しがずっといいのだけど、結構速度に差があったので速い方を採用。


LDA に限らず stop words をどのように扱うか、というのがよくポイントになるので、そこらへんを制御できるようにいくつか細工している。
というわけで条件をいくつか変えてみての実行結果とかは次回。