Stanでちゃんと動くトピックモデルを考案してみた

# はじめに

皆さんの中には、私のように、StanでLDAをはじめとするトピックモデルを実装しようとしたがうまくいかなかった経験をした方もいるかもしれない。

本記事では、LDAやSTMがStanでうまくいかない理由を考察しながら、私が新しく考案したトピックモデルをStanで変分推論して性能を紹介する。


# Stanの変分推論がうまくLDA系の手法を扱えない理由と解決策

検証時のデータは残っていないため、定性的な紹介にはなってしまうが、Stanでトピックモデル系の手法を回すとき、大体確率的勾配上昇法でiterationが300回のところで止まり、あまり解釈性のないトピックが事後分布としてサンプリングされる。

Srivastava and Sutton（2017）の論文ではADVIというStanが利用する自動微分変分推論は質の良いトピックを抽出できない現象を報告し、ディリクレ分布がlocation scale族の分布ではないため、再パラメータ化の際に問題が生じると説明する。

Srivastava and Sutton（2017）は、LDAの問題を解消すべく、ProdLDAという手法を提案する。

本記事では、Srivastava and Sutton（2017）のモデルをベースに独自の手法の提案を試みる。

# 提案モデル

まずは、本記事が提案するモデルが仮定するデータ生成過程を説明する。

全てのドキュメントについて

```math
doc\_latent\sim Normal(0,1)
```

doc_latentはK次元のベクトルで、0から1までに収まるように標準化はしていないが、ドキュメントのトピックとして解釈できる。

全ての単語について

```math
word\_offset\sim Normal(0,1)
```

```math
word\_latent\sim Normal(0,1)
```

word_offsetは一次元の実数で、単語の全体的な人気度を示す。

word_latentはdoc_latentのようにK次元のベクトルで、標準化はしていないが、単語のトピック要素の強さとして解釈できる。

最後に、全ての文章内の単語iについて

```math
freq_{i} \sim Poisson(exp(word\_offset_{word\space id_{i}} + word\_latent_{word\space id_{i}} '* doc\_latent_{doc\space id_{i}}))
```

言葉で説明すると、これは単語が所属するドキュメントのdoc_latentとその単語のword_latentで内積を取って、単語のword_offsetを足した後に指数関数に入れる処理である。

この処理で得られた数字はポワソン分布のパラメータとなり、単語の出現回数freqを予測する。

ProdLDAは多項分布を尤度関数にするのに対し、本記事ではTaddy（2015）に倣ってポワソン分布を利用するところを強調したい。

Stanのコードは下記の通りである：

```c++: prodlda.stan
functions {
  real partial_sum_lpmf(
    array[] int freq,
    int start, int end,
    
    array[] int doc, array[] int word,
    
    array[] vector doc_latent, 
    array[] real word_offset, array[] vector word_latent
  ){
    vector[end - start + 1] lambda;
    int count = 1;
    for (i in start:end){
      lambda[count] = word_offset[word[i]] + doc_latent[doc[i]] '* word_latent[word[i]];
      count += 1;
    }
    return(
      poisson_log_lupmf(freq | lambda)
    );
  }
}
data {
  int<lower=1> N;
  int<lower=1> K;
  int<lower=1> doc_type;
  int<lower=1> word_type;
  
  array[N] int<lower=1,upper=doc_type> doc;
  array[N] int<lower=1,upper=word_type> word;
  array[N] int<lower=0> freq;
}
parameters{
  array[doc_type] vector[K] doc_latent;
  array[word_type] real word_offset;
  array[word_type] vector[K] word_latent;
}
model{
  for (d in 1:doc_type){
    doc_latent[d] ~ normal(0, 1);
  }
  for (w in 1:word_type){
    word_offset[w] ~ normal(0, 1);
    word_latent[w] ~ normal(0, 1);
  }
  
  int grainsize = 1;
  
  target += reduce_sum(
    partial_sum_lupmf, freq, grainsize,
    
    doc, word,
    
    doc_latent, 
    word_offset, word_latent
  );
}
```

# データと前処理

本記事の前処理は以前の別の記事に準拠しているため、詳細はここで確認していただきたい：

https://qiita.com/Gotoubun_taiwan/items/0d58402eb84f5508c96a

本記事では、モデルの推定時間を節約するため、データ量をかなり減らして分析を行う。

まずは単語分割用の関数を定義して、文書行列を作成してからこれを長持ちにするsub_corp_use_dfm_dfを作成する。

```r
mecabbing <- function(text){
  this_review <- stringr::str_replace_all(text, "[^一-龠ぁ-んーァ-ヶーa-zA-Z]", " ")
  mecab_output <- unlist(RMeCab::RMeCabC(this_review, 1))
  mecab_combined <- stringr::str_c(mecab_output[which(names(mecab_output) == "名詞")], collapse = " ")
  return(mecab_combined)
}

#　データのダウンロード
full_corp <- quanteda.corpora::download("data_corpus_foreignaffairscommittee") |>
  quanteda::corpus_subset(speaker != "")

# 肩書きを含む人名の抽出
capacity_full <- full_corp |>
  stringr::str_sub(1, 20) |>
  stringr::str_replace_all("\\s+.+|\n", "") |> # 冒頭の名前部分の取り出し
  stringr::str_replace( "^.+[参事|政府特別補佐人|内閣官房|会計検査院|最高裁判所長官代理者|主査|議員|副?大臣|副?議長|委員|参考人|分科員|公述人|]君(（.+）)?$", "\\1")

# 肩書きの抽出
capacity <- full_corp |>
  stringr::str_sub(1, 20) |>
  stringr::str_replace_all("\\s+.+|\n", "") |> # 冒頭の名前部分の取り出し
  stringr::str_replace( "^.+?(参事|政府特別補佐人|内閣官房|会計検査院|最高裁判所長官代理者|主査|議員|副?大臣|副?議長|委員|参考人|分科員|公述人|君(（.+）)?$)", "\\1") |>  # 冒頭の○から，名前部分までを消去
  stringr::str_replace("（.+）", "") |>
  stringr::str_replace("^○.+", "Other")

# 年の抽出
record_year <- docvars(full_corp, "date") |>
  lubridate::year() |>
  as.numeric()

# 並列処理の準備
future::plan(future::multisession, workers = 16)

sub_corp_df <- tibble::tibble(
  capacity_full = capacity_full,
  capacity = capacity,
  record_year = record_year,
  # テキストデータ形式をquantedaからbase Rのcharacter型に変換する
  text = full_corp |>
    as.list() |>
    unlist()
) |>
  dplyr::filter(dplyr::between(record_year, 2013, 2015)) |>
  dplyr::filter(capacity %in% c("委員", "大臣", "副大臣")) |>
  # 発言の冒頭の名前と肩書きを削除
  dplyr::mutate(
    text = stringr::str_remove_all(text, stringr::fixed(capacity_full))
  ) |>
  dplyr::mutate(text_mecabbed = furrr::future_map_chr(text, ~ mecabbing(.x))) |>
  dplyr::mutate(
    text_id = dplyr::row_number()
  )


sub_corp_use_dfm_df <- sub_corp_df |>
  dplyr::pull(text_mecabbed) |>
  quanteda::phrase() |>
  quanteda::tokens() |>
  quanteda::tokens_remove("の") |>
  quanteda::dfm() |>
  quanteda::dfm_trim(min_termfreq = 50, max_termfreq = 10000) |>
  quanteda::convert(to = "tripletlist") |>
  tibble::as_tibble()
```

中身はこうなっている

```r
> sub_corp_use_dfm_df
# A tibble: 250,824 × 3
   document feature frequency
   <chr>    <chr>       <dbl>
 1 text1    発言            2
 2 text50   発言            2
 3 text132  発言            1
 4 text158  発言            1
 5 text166  発言            1
 6 text197  発言            1
 7 text204  発言            1
 8 text267  発言            1
 9 text359  発言            3
10 text363  発言            1
# ℹ 250,814 more rows
# ℹ Use `print(n = ...)` to see more rows
```

次に、文書マスターと単語マスターという二枚のマスター表を作成し、sub_corp_use_dfm_dfにJOINすることで単語と文書名をindex化する

```r
doc_master <- sub_corp_use_dfm_df |>
  dplyr::pull(document) |>
  unique() |>
  tibble::tibble(
    document = _
  ) |>
  dplyr::mutate(
    document_id = dplyr::row_number()
  )


word_master <- sub_corp_use_dfm_df |>
  dplyr::pull(feature) |>
  unique() |>
  tibble::tibble(
    feature = _
  ) |>
  dplyr::mutate(
    feature_id = dplyr::row_number()
  )


sub_corp_use_dfm_df_id <- sub_corp_use_dfm_df |>
  dplyr::left_join(doc_master, by = "document") |>
  dplyr::left_join(word_master, "feature") |>
  dplyr::select(document_id, feature_id, frequency)
```

ここでは、word2vec（Mikolov et al. 2013）を参考に、ネガティブサンプル（文書行列で0になっている要素）を抽出する。

ここで気を付けていただきたいのは、本記事のネガティブサンプル抽出方法は少し乱暴で、本当にその単語が該当のドキュメントに現れていないのかを確認していない。

理由としては、文書行列というのは基本的にはかなりスパースなものなので、デタラメにつけたドキュメントと単語の組み合わせは多分本当のデータ内で現れていないのであろうと思われる。より厳密なネガティブサンプル抽出方法でトピックモデルの性能が上がるかの確認はまた今後別の記事でまとめる。

```r
sub_corp_use_dfm_df_id_neg <- sub_corp_use_dfm_df_id |>
  dplyr::bind_rows(
    sub_corp_use_dfm_df_id |>
      dplyr::mutate(
        feature_id = sample(word_master$feature_id, dplyr::n(), replace = TRUE),
        frequency = 0
      )
  )
```

では最後に、データをlist型にまとめて、Stanに渡すための準備をする

```r
data_list <- list(
  N = nrow(sub_corp_use_dfm_df_id_neg),
  K = 20,
  doc_type = nrow(doc_master),
  word_type = nrow(word_master),
  
  doc = sub_corp_use_dfm_df_id_neg$document_id,
  word = sub_corp_use_dfm_df_id_neg$feature_id,
  freq = sub_corp_use_dfm_df_id_neg$frequency
)
```

# モデル推定


ここでは実際にモデルを推定する。


```r
m_prodlda_init <- cmdstanr::cmdstan_model("prodlda.stan",
                                          cpp_options = list(
                                            stan_threads = TRUE
                                          )
                                          )
m_prodlda_estimate <- m_prodlda_init$variational(
  seed = 12345,
  threads = 10,
  iter = 50000,
  data = data_list
)
```

変分推論のプロセスは下記の通りで

```r
------------------------------------------------------------ 
EXPERIMENTAL ALGORITHM: 
  This procedure has not been thoroughly tested and may be unstable 
  or buggy. The interface is subject to change. 
------------------------------------------------------------ 
Gradient evaluation took 0.444313 seconds 
1000 transitions using 10 leapfrog steps per transition would take 4443.13 seconds. 
Adjust your expectations accordingly! 
Begin eta adaptation. 
Iteration:   1 / 250 [  0%]  (Adaptation) 
Iteration:  50 / 250 [ 20%]  (Adaptation) 
Iteration: 100 / 250 [ 40%]  (Adaptation) 
Iteration: 150 / 250 [ 60%]  (Adaptation) 
Iteration: 200 / 250 [ 80%]  (Adaptation) 
Iteration: 250 / 250 [100%]  (Adaptation) 
Success! Found best value [eta = 0.1]. 
Begin stochastic gradient ascent. 
  iter             ELBO   delta_ELBO_mean   delta_ELBO_med   notes  
   100  -137957903483599549628416.000             1.000            1.000 
   200  -3058087880879991947264.000            22.556           44.112 
   300  -32209002491675979776.000            46.353           44.112 
   400  -79290342117543690240.000            34.913           44.112 
   500  -4797321287981455360.000            31.036           15.528 
   600  -535696928667365184.000            27.189           15.528 
   700  -110690755012574736.000            23.853            7.955 
   800  -2223296200649030144.000            20.991            7.955 
   900  -197284427548197504.000            19.799            7.955 
  1000  -9983370092215800.000            19.696           10.269 
  1100  -152838372577714048.000            17.990            7.955   MAY BE DIVERGING... INSPECT ELBO 
  1200  -4369644052809880.500            19.322           10.269   MAY BE DIVERGING... INSPECT ELBO 
  1300  -87165778264286.391            21.615           10.269   MAY BE DIVERGING... INSPECT ELBO 
  1400  -562085497218604.500            20.132           10.269   MAY BE DIVERGING... INSPECT ELBO 
  1500  -43412809067628.031            19.586           10.269   MAY BE DIVERGING... INSPECT ELBO 
  1600  -8591406884017.072            18.615           10.269   MAY BE DIVERGING... INSPECT ELBO 
  1700  -7831583667540.095            17.526            7.955   MAY BE DIVERGING... INSPECT ELBO 
  1800  -1083930033048.644            16.898            7.955   MAY BE DIVERGING... INSPECT ELBO 
  1900  -2424396703666.407            16.038            6.225   MAY BE DIVERGING... INSPECT ELBO 
  2000  -224322499061.449            15.726            7.955   MAY BE DIVERGING... INSPECT ELBO 
  2100  -3031861796065.904            15.022            6.225   MAY BE DIVERGING... INSPECT ELBO 
  2200  -13909428568684.396            14.374            6.225   MAY BE DIVERGING... INSPECT ELBO 
  2300  -40596569467.562            28.603            6.225   MAY BE DIVERGING... INSPECT ELBO 
  2400  -1610468792905.306            27.451            6.225   MAY BE DIVERGING... INSPECT ELBO 
  2500  -7477819349.315            34.928            6.225   MAY BE DIVERGING... INSPECT ELBO 
  2600  -39826017915.808            33.616            6.225   MAY BE DIVERGING... INSPECT ELBO 
  2700  -1657563451.102            33.224            6.225   MAY BE DIVERGING... INSPECT ELBO 
  2800  -9085607331.989            32.066            6.225   MAY BE DIVERGING... INSPECT ELBO 
  2900  -20392808265.381            30.980            4.053   MAY BE DIVERGING... INSPECT ELBO 
  3000  -3737378671.543            30.096            4.456   MAY BE DIVERGING... INSPECT ELBO 
  3100   -112196006.633            30.167            4.456   MAY BE DIVERGING... INSPECT ELBO 
  3200    -70928085.622            29.243            4.456   MAY BE DIVERGING... INSPECT ELBO 
  3300    -29755741.921            28.398            4.053   MAY BE DIVERGING... INSPECT ELBO 
  3400    -19221889.772            27.579            4.053   MAY BE DIVERGING... INSPECT ELBO 
  3500    -35179154.386            26.804            3.840   MAY BE DIVERGING... INSPECT ELBO 
  3600     -8593471.468            26.146            3.840   MAY BE DIVERGING... INSPECT ELBO 
  3700     -6669308.185            25.447            3.094   MAY BE DIVERGING... INSPECT ELBO 
  3800     -4093079.203            24.794            3.094   MAY BE DIVERGING... INSPECT ELBO 
  3900     -1832432.475            24.190            1.384   MAY BE DIVERGING... INSPECT ELBO 
  4000     -1573847.888            23.589            1.384   MAY BE DIVERGING... INSPECT ELBO 
  4100      -958555.014            23.029            1.234   MAY BE DIVERGING... INSPECT ELBO 
  4200     -1008332.468            22.482            1.234   MAY BE DIVERGING... INSPECT ELBO 
  4300      -816182.287            21.965            1.000   MAY BE DIVERGING... INSPECT ELBO 
  4400      -774975.834            21.467            1.000   MAY BE DIVERGING... INSPECT ELBO 
  4500      -748160.862            20.990            0.975   MAY BE DIVERGING... INSPECT ELBO 
  4600      -759949.116            20.535            0.975   MAY BE DIVERGING... INSPECT ELBO 
  4700      -703336.775            20.099            0.950   MAY BE DIVERGING... INSPECT ELBO 
  4800      -697667.822            19.681            0.950   MAY BE DIVERGING... INSPECT ELBO 
  4900      -671021.832            19.280            0.935   MAY BE DIVERGING... INSPECT ELBO 
  5000      -661322.128            18.895            0.935   MAY BE DIVERGING... INSPECT ELBO 
  5100      -652177.203            18.875            0.926   MAY BE DIVERGING... INSPECT ELBO 
  5200      -644466.137            17.993            0.845   MAY BE DIVERGING... INSPECT ELBO 
  5300      -637961.217            16.114            0.818   MAY BE DIVERGING... INSPECT ELBO 
  5400      -631553.157            16.102            0.818   MAY BE DIVERGING... INSPECT ELBO 
  5500      -626735.462            15.792            0.812   MAY BE DIVERGING... INSPECT ELBO 
  5600      -621690.661            15.633            0.782   MAY BE DIVERGING... INSPECT ELBO 
  5700      -617344.445            15.556            0.642   MAY BE DIVERGING... INSPECT ELBO 
  5800      -613563.867            15.538            0.629   MAY BE DIVERGING... INSPECT ELBO 
  5900      -610186.666            15.332            0.582   MAY BE DIVERGING... INSPECT ELBO 
  6000      -607015.977            14.957            0.554   MAY BE DIVERGING... INSPECT ELBO 
  6100      -603992.092            14.939            0.553   MAY BE DIVERGING... INSPECT ELBO 
  6200      -601539.971            14.259            0.548   MAY BE DIVERGING... INSPECT ELBO 
  6300      -599236.463            13.277            0.454   MAY BE DIVERGING... INSPECT ELBO 
  6400      -597178.248            13.260            0.289   MAY BE DIVERGING... INSPECT ELBO 
  6500      -595296.593            13.021            0.235   MAY BE DIVERGING... INSPECT ELBO 
  6600      -593513.828            12.940            0.164   MAY BE DIVERGING... INSPECT ELBO 
  6700      -591940.697            12.938            0.164   MAY BE DIVERGING... INSPECT ELBO 
  6800      -590506.703            12.814            0.080   MAY BE DIVERGING... INSPECT ELBO 
  6900      -589193.432            12.803            0.053   MAY BE DIVERGING... INSPECT ELBO 
  7000      -587990.095            12.606            0.049   MAY BE DIVERGING... INSPECT ELBO 
  7100      -586693.297            12.588            0.040   MAY BE DIVERGING... INSPECT ELBO 
  7200      -585823.773            12.572            0.036   MAY BE DIVERGING... INSPECT ELBO 
  7300      -584778.416             5.740            0.016   MAY BE DIVERGING... INSPECT ELBO 
  7400      -583868.067             5.720            0.015   MAY BE DIVERGING... INSPECT ELBO 
  7500      -583034.319             1.433            0.014   MAY BE DIVERGING... INSPECT ELBO 
  7600      -582430.143             1.417            0.012   MAY BE DIVERGING... INSPECT ELBO 
  7700      -581745.874             0.956            0.010   MAY BE DIVERGING... INSPECT ELBO 
  7800      -581093.631             0.940            0.010   MAY BE DIVERGING... INSPECT ELBO 
  7900      -580435.844             0.929            0.008   MEDIAN ELBO CONVERGED   MAY BE DIVERGING... INSPECT ELBO 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  4050.8 seconds.
```

うまく収束したように見える。


# トピック内容の確認


まずはdoc_latentとword_latentを抽出する

```r
m_prodlda_summary <- m_prodlda_estimate$summary()

word_latent <- m_prodlda_summary |>
  dplyr::filter(stringr::str_detect(variable, "word_latent")) |>
  dplyr::pull(mean) |>
  matrix(ncol = 20)

doc_latent <- m_prodlda_summary |>
  dplyr::filter(stringr::str_detect(variable, "doc_latent")) |>
  dplyr::pull(mean) |>
  matrix(ncol = 20)
```

次に、各トピックのtop wordsを確認する

```r
> purrr::map(
     1:20,
     \(x) word_master |>
         dplyr::bind_cols(e = word_latent[,x]) |>
         dplyr::arrange(-e) |>
         dplyr::select(feature)
 ) |>
     purrr::reduce(dplyr::bind_cols)
New names:
• `feature` -> `feature...1`
• `feature` -> `feature...2`
New names:
• `feature` -> `feature...3`
New names:
• `feature` -> `feature...4`
New names:
• `feature` -> `feature...5`
New names:
• `feature` -> `feature...6`
New names:
• `feature` -> `feature...7`
New names:
• `feature` -> `feature...8`
New names:
• `feature` -> `feature...9`
New names:
• `feature` -> `feature...10`
New names:
• `feature` -> `feature...11`
New names:
• `feature` -> `feature...12`
New names:
• `feature` -> `feature...13`
New names:
• `feature` -> `feature...14`
New names:
• `feature` -> `feature...15`
New names:
• `feature` -> `feature...16`
New names:
• `feature` -> `feature...17`
New names:
• `feature` -> `feature...18`
New names:
• `feature` -> `feature...19`
New names:
• `feature` -> `feature...20`
# A tibble: 1,272 × 20
   feature...1 feature...2 feature...3 feature...4 feature...5 feature...6 feature...7 feature...8 feature...9 feature...10 feature...11 feature...12 feature...13 feature...14
   <chr>       <chr>       <chr>       <chr>       <chr>       <chr>       <chr>       <chr>       <chr>       <chr>        <chr>        <chr>        <chr>        <chr>       
 1 パッケージ  シリア      海          農業        旅券        品          公務員      捕鯨        試算        官民         炉           監護         ジャパン     濃縮        
 2 嘉手納      海域        空域        水          水域        郵便        和平        司法        輸送        人事         取り決め     子           ベトナム     市民        
 3 市          慰安        北          仕組み      中南米      特許        カンボジア  参拝        牛肉        所得         燃料         親権         改革         処理        
 4 トン        トン        海兵        緊急        操業        知的        め          だめ        ページ      研修         漁業         規則         選挙         トルコ      
 5 部隊        電話        配置        メンバー    テロ        処分        プロセス    葉          関税        金額         連邦         保存         多国         障害        
 6 海兵        沿岸        移転        率          包括        出願        パワー      割          農林        事務所       原子力       不法         登録         損害        
 7 移設        遺憾        中央        親権        困難        構想        スピーチ    骨子        隻          インフラ     北           締約         インド       密接        
 8 グアム      ドイツ      費用        人材        平壌        批准        総会        海上保安庁  生産        委託         意向         保全         豪           司法        
 9 以南        大使        自己        英語        唯一        討論        制約        文化        空域        アフリカ     物質         ハーグ       条項         保険        
10 ヨルダン    婦          隊          ガバナンス  南シナ海    条項        構造        秩序        線          勤務         濃縮         郵便         ベース       艦          
# ℹ 1,262 more rows
# ℹ 6 more variables: feature...15 <chr>, feature...16 <chr>, feature...17 <chr>, feature...18 <chr>, feature...19 <chr>, feature...20 <chr>
# ℹ Use `print(n = ...)` to see more rows
```

見てわかるように、トピック1は当時の記事移設問題トピックで、トピック2は慰安婦問題とシリア情勢が混ざったトピックで、トピック8は捕鯨関連のトピックになっている。

このように、完璧ではないものの、LDAをStanで回すときより大幅に改善したトピックの質だと言っても過言ではない。


最後に、それぞれのトピックの最も代表的なドキュメントを確認する。


```r
> doc_master |>
     dplyr::bind_cols(
         e = doc_latent[,1]
     ) |>
     dplyr::mutate(
         text_id = as.numeric(stringr::str_remove_all(document, "text"))
     ) |>
     dplyr::mutate(
         doc_content = purrr::map_chr(text_id,
                                      \(x) sub_corp_df$text[x]
         )
     ) |>
     dplyr::arrange(-e)
# A tibble: 6,253 × 5
   document document_id     e text_id doc_content                                                                                                                                         
   <chr>          <int> <dbl>   <dbl> <chr>                                                                                                                                               
 1 text5489        2020  1.69    5489 "　そうして前向きにやっていくということであれば、何回も言いますけれども、五月九日のものは政府としてしっかり対応しなきゃいけないと思いますよ。行った…
 2 text1018         976  1.58    1018 "　特に嘉手納町におきましては、町域の約八三％に及ぶ十二・四〇平方キロメートルが米軍基地に占用されているという実態があり、私の手元には、平成二十四年…
 3 text4820        1843  1.54    4820 "　基地の固定化を避けるという目的に沿って現在の宜野湾市から普天間飛行場の移設先に予定されているアメリカ海兵隊キャンプ・シュワブがある名護市では、本…
 4 text1388        1064  1.40    1388 "　生活の党の玉城デニーでございます。\n　国際的な子の奪取の民事上の側面に関する条約、ハーグ条約について質問させていただきます。\n　初めに、大臣、こ…
 5 text5916        4310  1.37    5916 "　そこで、次の質問なのでありますが、先ほど、全体としてという表現についても同じことを言いましたが、既存の政府見解の基本的な論理という表現を使って今…
 6 text3727        1590  1.33    3727 "　パッケージを切り離すことになった経緯について御質問をいただきましたが、かつての在日米軍のパッケージは、沖縄の負担軽減を早期に実現しなければいけな…
 7 text2002          63  1.33    2002 "　この条約の承認に関する件は以上とさせていただきますが、やはりまだ明確になっていないといいますか、明確にすべきであるというふうなことが、今現在、福…
 8 text1841        1173  1.26    1841 "　ですから、大臣がおっしゃったように、嘉手納以南の土地の返還計画を早期に進めますという、これをとっても物すごく時間のかかるスケジュールなんですね。…
 9 text4562        2531  1.24    4562 "　生活の党の玉城デニーでございます。\n　早速質問に入らせていただきます。\n　きょうは、米軍関連の事件、事故について幾つかただしたいと思います。\n　…
10 text1821        2699  1.21    1821 "　米軍が駐留していた岐阜県各務原では、米軍による傍若無人な行為が繰り返されて、住民の米軍基地反対闘争が起こりました。農民の入会地だった山梨の演習場…
# ℹ 6,243 more rows
# ℹ Use `print(n = ...)` to see more rows
```

このように、米軍関連と基地関連が多い。

トピック20を確認すると

```r
> doc_master |>
     dplyr::bind_cols(
         e = doc_latent[,20]
     ) |>
     dplyr::mutate(
         text_id = as.numeric(stringr::str_remove_all(document, "text"))
     ) |>
     dplyr::mutate(
         doc_content = purrr::map_chr(text_id,
                                      \(x) sub_corp_df$text[x]
         )
     ) |>
     dplyr::arrange(-e)
# A tibble: 6,253 × 5
   document document_id     e text_id doc_content                                                                                                                                         
   <chr>          <int> <dbl>   <dbl> <chr>                                                                                                                                               
 1 text1214         653  1.95    1214 "　ただいま議題となりました脱税の防止のための情報の交換及び個人の所得に対する租税に関する二重課税の回避のための日本国政府とジャージー政府との間の協…
 2 text662          645  1.94     662 "　ただいま議題となりました投資の促進及び保護に関する日本国政府とパプアニューギニア独立国政府との間の協定の締結について承認を求めるの件につきまして…
 3 text5603        3710  1.80    5603 "　これまで外務委員会でも、ＴＰＰの内容についていろいろ質問があったと思うんですけれども、その際、これは外交交渉であるということと、ＴＰＰ参加国に課…
 4 text840          648  1.80     840 "　ただいま議題となりました投資の促進及び保護に関する日本国政府とパプアニューギニア独立国政府との間の協定の締結について承認を求めるの件につきまして…
 5 text3183         450  1.71    3183 "　日本維新の会の石関貴史です。質問の機会をいただきまして、ありがとうございます。\n　対ロシア、今、制裁中でございますが、このことの関係で幾つか御質…
 6 text848          934  1.63     848 "　これは、実際やってみればよくわかるんです。私も仲裁を何回もやったことがあるんですけれども、仲裁は、一人はこっちがいい人、一人は向こうがいい人を選…
 7 text6261        2779  1.62    6261 "　沖縄県が今回求めてきたのは、臨時制限区域内での立入調査であります。ボーリング調査に関して、岩礁破砕の許可ではなく協議の手続がとられたのは、ブイや…
 8 text2727        1388  1.60    2727 "　ただいま議題となりました所得に対する租税に関する二重課税の回避及び脱税の防止のための日本国とアラブ首長国連邦との間の条約の締結について承認を求め…
 9 text5304         261  1.57    5304 "　まず、ワシントンＤＣでの西村内閣府副大臣の発言に関しましては、制度上の違いから米国と同一の対応は困難であるが、そのような前提や制約の中で、我が国…
10 text1356        2269  1.57    1356 "　十五件というのは決して少なくない数だと思うんです。国際司法裁判所というのは、双方がそれによって解決することについて同意をしない限りは、この手続に…
# ℹ 6,243 more rows
# ℹ Use `print(n = ...)` to see more rows
```

関税関連のトピックになっていることがわかる。

# 最後に

このように、ディリクレ分布など、値が0から1までしか取らない分布ではなく、正規分布を設定してword2vecのように隠れ変数の内積を取ることで単語の出現回数を予測する形でモデルを定式化すれば、Stanでも問題なくトピックモデルを回すことができることを確認した。

# 参考文献

Mikolov, Tomas, et al. "Distributed representations of words and phrases and their compositionality." Advances in neural information processing systems 26 (2013).

Srivastava, Akash, and Charles Sutton. "Autoencoding variational inference for topic models." arXiv preprint arXiv:1703.01488 (2017).

Taddy, Matt. "Distributed multinomial regression." The Annals of Applied Statistics 9.3 (2015): 1394-1414.