StanでLDAをまわしてみた　（自動微分変分推論）

# はじめに
LDAはトピックモデルと呼ばれる自然言語処理の分野で広く使われている数理モデルです。モデルの詳細は後ほど説明しますが、ざっくり言うと書かれた単語からその文章のトピックを複数推定するモデルです。

LDAのパラメータ推定には一般的にはRのtopicmodelsパッケージやPythonのgensimモジュールが使われることが多いですが、今回はモデルの拡張や仕組みの理解を目的として、stanでモデルをゼロから構築して推定することを試みます。

# モデルの定義
コードは[松浦さんのブログ](http://statmodeling.hatenablog.com/entry/topic-model-4)から借りてきました。ただし、stanのバージョンアップに合わせて、一部コードを修正しています。

```stan
data {
  int<lower=1> K;
  int<lower=1> M;
  int<lower=1> V;
  int<lower=1> N;
  int<lower=1,upper=V> W[N];
  int<lower=1,upper=N> Offset[M,2]; // range of word index per doc
  
  vector<lower=0>[K] Alpha; // topic prior
  vector<lower=0>[V] Beta; // word prior
}
parameters {
  simplex[K] theta[M];
  simplex[V] phi[K];
}
model {
  // prior
  for (m in 1:M)
    theta[m] ~ dirichlet(Alpha);
  for (k in 1:K)
    phi[k] ~ dirichlet(Beta);
  
  // likelihood
  for (m in 1:M) {
    for (n in Offset[m,1]:Offset[m,2]) {
      real gamma[K];
      for (k in 1:K)
        gamma[k] = log(theta[m,k]) + log(phi[k,W[n]]);
      target += log_sum_exp(gamma);
    }
  }
}
```

# データと下処理
今回は20_newsコーパスを使いました。これは20個のジャンルからなる20,000本のニュース記事の集合で、トピックモデルの評価にしばしば使われるコーパスです。

予め単語にインデックスを付与し、文章は`1 10 8 9`などの整数のリストに変換しています。したがって別ファイル`vocab.txt`にインデックスと単語の対応関係`1 --> "news"`などを記録しています。

```r
library(tidyverse)
library(tidytext)

corpus = read_csv("./corpus_20news.txt", col_names = "words")
corpus = corpus %>% sample_n(1000)
corpus = corpus %>% mutate(doc = seq(1:nrow(corpus)))
corpus.tidy = corpus %>% unnest_tokens(word, words)
corpus.tidy = corpus.tidy %>% mutate(word = as.integer(word))
```

上のコードでコーパスをtidytext形式に変換しています。

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> doc </th>
   <th style="text-align:right;"> word </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 178 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 878 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 866 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 3579 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 2045 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 522 </td>
  </tr>
</tbody>
</table>

ここから、頻出語の削除とそれに伴うインデックスの振り直しを行っていきます。

```r
stop_words = corpus.tidy %>% count(word, sort=T) %>% filter(n > 1000 | n < 15)
corpus.tidy = corpus.tidy %>% anti_join(stop_words)

appeared_word = unique(corpus.tidy$word)

vocab = read_csv("./vocab_20news.txt", col_names = "word")
vocab = vocab %>% mutate(index = seq(nrow(vocab)))
vocab.new = vocab[appeared_word,]

vocab.new = vocab.new %>% mutate(new_ind = seq(nrow(vocab.new)))

corpus.lda = corpus.tidy %>% left_join(vocab.new, by = c("word" = "index"))
```
あとは、文書iに属する単語がtidy形式でインデックス何番目から何番目にあるのかを表す情報をoffset変数に格納していきます。

```r
M = nrow(corpus)
N = nrow(corpus.lda)
V = nrow(vocab.new)
K = 20

offset = matrix(nrow = M, ncol = 2)

last_ind = 0
for (m in 1:M) {
  n_words = as.integer(corpus.lda %>% filter(doc == m) %>% count())
  offset[m,1] = last_ind + 1
  offset[m,2] = offset[m,1] + (n_words - 1)
  last_ind = offset[m,2]
}
```

# パラメータ推定
モデルとデータがそろったのでパラメータ推定を行っていきましょう。アルゴリズムにはver.2.9からstanに搭載されたADVI（自動微分変分推論）を使いました。stanといえばNUTSというアルゴリズムが有名ですが、MCMCだとLDAのようなパラメータの多いアルゴリズムでは推定にとてつもない時間が必要になってしまいます。ADVIであれば、MCMCより精度は劣るものの、高速に推論を行うことができます。

```r
library(rstan)

load("./20news.RData")

data = list(
  K = K,
  M = M,
  V = V,
  N = N,
  W = corpus.lda$new_ind,
  Offset = offset,
  Alpha = rep(1, K),
  Beta = rep(0.5, V)
)

sm = stan_model(file = "./lda.stan")

fit.vb = vb(
  sm,
  data = data,
  output_samples = 2000,
  adapt_engaged = FALSE,
  eta = .1
)
```

etaは自動的に調整してくれる機能もありますが、時間がかかるのでここでは天下り的に適当な値を指定しています。

# 結果
結果は以下のようになりました。

|         |           |            |            |         |        |        |           |         |       |            |
|:--------|:----------|:-----------|:-----------|:--------|:-------|:-------|:----------|:--------|:------|:-----------|
|topic 1  |university |ca          |time        |apr      |problem |made    |people     |fact     |law    |israel      |
|topic 2  |image      |years       |time        |make     |things  |people  |set        |god      |file   |color       |
|topic 3  |ftp        |information |windows     |system   |file    |version |files      |jpeg     |scsi   |gif         |
|topic 4  |image      |list        |system      |graphics |bit     |format  |line       |run      |images |jpeg        |
|topic 5  |time       |apr         |make        |message  |problem |car     |back       |point    |god    |drug        |
|topic 6  |man        |day         |time        |apr      |work    |people  |car        |called   |ve     |drive       |
|topic 7  |image      |mail        |pc          |people   |find    |price   |software   |file     |files  |openwindows |
|topic 8  |good       |ftp         |state       |space    |apr     |made    |system     |point    |year   |team        |
|topic 9  |good       |ca          |time        |support  |make    |people  |jews       |religion |ve     |god         |
|topic 10 |good       |image       |university  |key      |ca      |work    |people     |system   |bit    |jpeg        |
|topic 11 |cs         |good        |state       |time     |make    |ve      |government |god      |part   |give        |
|topic 12 |good       |time        |bill        |made     |people  |person  |place      |life     |book   |government  |
|topic 13 |good       |ca          |time        |apr      |make    |work    |people     |thing    |file   |jpeg        |
|topic 14 |man        |years       |information |uk       |problem |things  |people     |car      |thing  |god         |
|topic 15 |good       |university  |years       |ed       |time    |support |make       |work     |ve     |god         |
|topic 16 |good       |makes       |apr         |make     |people  |back    |ve         |ll       |point  |god         |
|topic 17 |good       |question    |time        |make     |windows |work    |people     |god      |run    |file        |
|topic 18 |cs         |good        |question    |time     |make    |work    |read       |system   |ve     |data        |
|topic 19 |good       |information |mail        |time     |apr     |make    |work       |problem  |ve     |game        |
|topic 20 |question   |apr         |david       |make     |windows |problem |put        |people   |back   |god       
なんとなくトピックらしきものが抽出されてる気もしますが、重なるものが多かったりあんまり品質はよくないようです。

# 感想
ある程度予想はしていたのですが、やはり既存のパッケージに比べると推定時間・精度ともに悪いです。ただ、LDAには無数の派生系があり、そうしたモデルの中にはどうしても自分で実装するしかないものも多いので、その後の拡張のしやすさという意味ではstanを使う意義もあるのではないでしょうか。僕はとりあえずgaussian-LDAをstanでやってみようと思っています。
コードおよびデータは[ここ](https://github.com/Hiroyuki1993/lda_stan)にも保存しています。