はじめに

こんにちは、事業会社で働くデータサイエンティストです。

今回の記事では、多項分散回帰（Taddy, 2015）をベイズ機械学習用のStan言語で実装します。

多項分散回帰については、こちらの記事をご確認ください：

なぜやるのか

多項分散回帰を提案した著者は、すでに実装用のR言語のパッケージまで用意してくれました。

では、どうして車輪の再発明をしてStan言語で実装しますか？どうして読者の皆さんこの記事を読まないといけないですか？

答えは簡単です。モデル改造の土台作りのためです。

分散多項回帰もそうですが、他の人が用意したパッケージの中の、特にC、C++、pytorchなどで実装したアルゴリズムを、ドメイン知識を駆使して目の前の分析タスクに特化させるのは難しいです。

わかりやすい実例でいうと、distromパッケージには単語の傾きに時間的な変動を加える機能がないです。ただし、長いタイムスパンのデータの場合、どうしても単語の意味の変化を確認したいニーズが出てきます。この時、ベースとなる多項分散回帰をStanで実装し、そこから事前分布の設計などを通じて改造すれば、実現する可能性が高まります。

では、早速多項分散回帰の仕組みの説明に入ります。

モデル構造

今回の分析では、livedoor ニュースコーパスを利用して、各カテゴリーの特徴的な単語を識別します。

多項分散回帰は、LDA（Blei, Ng and Jordan, 2003）をはじめとするトピックモデルとword2vec（Mikolov et al., 2013）をはじめとする多くのニューラルネットワーク系・深層学習系のモデルと同じように、テキストデータの生成モデルです。

ここでひとつ強調しておきたいのは、本記事の目的は（Taddy, 2015）の完全な再現ではなく、あくまでもモデルのエッセンスの部分だけを取り出して同じ機能を再現できるかを検証します。元の論文に載っているモデルに興味ある方はこちらをご確認ください：

多項分布との関係は省略しますが、数式でいうと

全ての単語wと全てのカテゴリーcについて

$$\sigma_{w, c} \sim Gamma(0.01, 0.01)$$

$$\beta_{w, c} \sim Laplace(0, \sigma_{w, c})$$

次に、全てのドキュメントdと単語wについて

$$出現回数_{d,w} \sim Poisson(exp(latent_{文書_{d}} + latent_{単語_{w}} + \beta_{w, カテゴリ_{d}}))$$

になります。

要するに、文書の中のある単語の出現回数は、

文書の「隠れ長さ」（$latent_{文書_{d}}$）

この文書は長いから、どんな単語でもそもそも出てきやすいし、出現回数も多い（極端な事例：今日の日経新聞の全文を一つの文書にするとき）

単語の全体的な「人気度」（$latent_{単語_{w}}$）

この単語は使われやすいから、どんな文書でも登場しやすい（極端な事例：を、は、に、が）

カテゴリ内の単語の「人気度」（$\beta_{w, カテゴリ_{d}}$）

文書の「隠れ長さ」と単語の全体的な「人気度」だけで説明できない要素。文書が特定のカテゴリに該当すると利用頻度が上がる単語

を変数とするポワソン回帰でモデリングできます。

特に、文書が特定のカテゴリに該当すると単語の利用頻度が上がる度合いを示す$\beta$が高い単語をそのカテゴリの特徴的な単語にして、主な分析対象とします。

Taddy（2015）は$latent_{単語_{w}}$をパラメータとして推定するのではなく、$log(文書の長さ_{d})$という定数でモデルに渡しても、モデルの性能に悪影響がないことを実際のデータ分析で示しました。

これによって多項分布を独立したポワソン分布に置き換えて分散処理できると著者が説明しますが、個人的には文書の量とともに増加するパラメータを推定対象から外せるのもありがたいポイントだと思います。

ビジネスと歴史研究になると、数十万から数百万・数千万レベルの文書を分析しないといけない場面も多々あり、その際に$latent_{単語_{w}}$を全て推定しないといけなくなると、メモリなどで問題が起き、高度なインフラエンジニアリングのスキルで対処する必要が出てしまいます。なので、$log(文書の長さ_{d})$を一つの説明変数としてモデルに渡して$latent_{単語_{w}}$の代用品として使うことでモデルのスケーラビリティが一気に上がります。

モデル実装

Stanの実装コードはこちらになります。分散多項回帰には離散確率変数がないため、直感的にモデルの数式からコードに変換できます。

dmr.stan

functions {
  real partial_sum_lpmf(
    array[] int frequency,
    int start, int end,
    
    array[] int word,
    array[] int category,
    array[] real log_m,
    
    vector intercept_word,
    matrix beta_word
  ){
    vector[end - start + 1] lambda;
    int count = 1;
    for (i in start:end){
      lambda[count] = poisson_log_lpmf(frequency[count] | log_m[i] + intercept_word[word[i]] + beta_word[category[i], word[i]]);
      count += 1;
    }
    return sum(lambda);
  }
}
data {
  int category_type;
  int word_type;
  
  int N;
  array[N] int word;
  array[N] int category;
  array[N] real log_m;
  array[N] int frequency;
}
parameters {
  vector[word_type] intercept_word;
  matrix<lower=0>[category_type, word_type] sigma_word;
  matrix[category_type, word_type] beta_word;
}
model {
  to_vector(sigma_word) ~ gamma(0.01, 0.01);
  to_vector(beta_word) ~ double_exponential(0, to_vector(sigma_word));
  
  int grainsize = 1;
  
  target += reduce_sum(
    partial_sum_lupmf, frequency, grainsize,
    
    word,
    category,
    log_m,
    
    intercept_word,
    beta_word
  );
}

前処理

まずは、日本語の分かち書き用の関数とmagrittrパッケージのパイプを用意します。

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

`%>%` <- magrittr::`%>%`

magrittrのパイプ（%>%）はbase Rのパイプ（|>）より性能が豊富ですが、使い過ぎるとかえって可読性が下がることもありますので、base Rのパイプを基本にし、magrittrのパイプの利用を最小限にします。

次に、こちらからlivedoor ニュースコーパスのデータをダウンロードして、前処理をはじめます。

最初は、.txtファイルに格納されたテキストを読み込み、分かち書き処理を実施してtibbleデータフレイムに変換します。

review_df <- list.files("text") %>%
  # CHANGES.txtとREADME.txtなどを削除
  .[which(stringr::str_detect(., ".txt") == FALSE)] |>
  purrr::map(
    # フォルダごとの処理
    \(this_file){
      this_file %>%
        # フォルダのパスを作成
        stringr::str_c("text/", .)  |> 
        list.files() %>%
        # フォルダ内のLICENSE.txtファイルを削除
        .[which(stringr::str_detect(., "LICENSE") == FALSE)] %>% 
        # テキストファイルのパスを作成
        stringr::str_c("text/", this_file, "/", .) |>
        purrr::map(
          # .txtファイルごとの処理
          \(this_text){
            this_text |>
              readr::read_lines() |>
              stringr::str_c(collapse = " ") |>
              mecabbing() |>
              stringr::str_remove_all(
                quanteda::stopwords("ja", source = "marimo") |> 
                  c("的", "の", "こと", "さ", "ら", "今", "いま") |>
                  stringr::str_c(collapse = "|")
              ) |>
              tibble::tibble(
                category = this_file,
                text = _
              )
          }
        )
    },
    .progress = TRUE
  ) |>
  dplyr::bind_rows() |>
  dplyr::mutate(
  　# 文書のカテゴリなどのメタ情報のマスターテーブルの役割を演じるため、キーを付与
    doc_id = dplyr::row_number()
  )

次に、単語文書行列を作成して、縦持ちにします。

review_dfm_df <- review_df |>
  dplyr::pull(text) |>
  quanteda::phrase() |>
  quanteda::tokens() |>
  quanteda::tokens_ngrams(1:2) |>
  quanteda::dfm() |>
  quanteda::dfm_trim(min_docfreq = 50) |>
  quanteda::convert(to = "tripletlist") |>
  tibble::as_tibble() |>
  dplyr::mutate(
    # documentを元の"text1", "text2"から1, 2に変換し、
    # review_dfとdoc_idでleft joinできるようにする
    document = document |>
      stringr::str_remove_all("text") |>
      as.integer()
  )

単語とカテゴリをidに変換するためのマスター表も作成します：

word_master <- review_dfm_df |>
  dplyr::select(feature) |>
  dplyr::distinct() |>
  dplyr::arrange(feature) |>
  dplyr::mutate(
    word_id = dplyr::row_number()
  )

category_master <- review_df |>
  dplyr::select(category) |>
  dplyr::distinct() |>
  dplyr::arrange(category) |>
  dplyr::mutate(
    category_id = dplyr::row_number()
  )

ここでは、review_dfm_dfの負例をサンプリングします。

単語文書行列を確認すればわかるように、行列内の要素はほとんどゼロです。

> review_df |>
     dplyr::pull(text) |>
     quanteda::phrase() |>
     quanteda::tokens() |>
     quanteda::tokens_ngrams(1:2) |>
     quanteda::dfm() |>
     quanteda::dfm_trim(min_docfreq = 50)
     
Document-feature matrix of: 7,367 documents, 3,960 features (97.49% sparse) and 0 docvars.
       features
docs    友人 代表 独 女 式 状態 人 出席 回数 お願い
  text1    6    3  5  5  2    2 10    1    1      5
  text2    1    0  2  2  0    0  2    0    0      0
  text3    0    0  5  6  0    0  2    0    0      0
  text4    2    0  0  0  0    0  7    0    0      0
  text5    0    0  2  2  0    0  9    0    0      0
  text6    1    0  0  0  0    0  4    0    0      0
[ reached max_ndoc ... 7,361 more documents, reached max_nfeat ... 3,950 more features ]

これはこのコンパスの特性ではなく、テキストデータ、ないしはより抽象的にいうと遺伝子データや購買データのような高次元のカテゴリデータに共通する現象です。このようにデータがゼロまみれの状態でゼロの情報を格納するのは極めて効率が悪いので、quanteda::convert(to = "tripletlist")はゼロのレコードを全て省略します。

ただ、ゼロが全くないとモデルに入れられないので、ここでは、ゼロを一部だけ「復元」します。具体的には、全てのゼロを補完するのではなく、文書の長さと同じ量のゼロを補完するやり方を紹介します。

これは私の理論的な根拠のないアイデアなので、興味ある方はぜひ大規模データか統計理論で妥当性を検証してください。

ついでに、idと$log(文書の長さ)$も計算します。

set.seed(12345)
review_dfm_neg_df <- review_dfm_df |>
  split(~ document) |>
  purrr::map(
    \(this_df){
      words_not_used <- word_master$feature[which(word_master$feature %in% this_df$feature == FALSE)]
      
      this_df |>
        dplyr::bind_rows(
          tibble::tibble(
            document = this_df$document[1],
            feature = sample(words_not_used, sum(this_df$frequency), replace = FALSE),
            frequency = 0
          )
        )
    }
  ) |>
  dplyr::bind_rows() |>
  dplyr::group_by(document) |>
  dplyr::mutate(
    log_m = log(sum(frequency))
  ) |>
  dplyr::ungroup() |>
  dplyr::left_join(
    review_df |>
      dplyr::select(doc_id, category),
    by = c("document" = "doc_id")
  ) |>
  dplyr::left_join(
    word_master, by = "feature"
  ) |>
  dplyr::left_join(
    category_master, by = "category"
  )

モデル推定

ここでは、Stanに渡すためモデル推定用のデータをlistに入れます。

data_list <- list(
  category_type = nrow(category_master),
  word_type = nrow(word_master),
  
  N = nrow(review_dfm_neg_df),
  word = review_dfm_neg_df$word_id,
  category = review_dfm_neg_df$category_id,
  log_m = review_dfm_neg_df$log_m,
  frequency = review_dfm_neg_df$frequency
)

次に、モデルをコンパイルします。

m_dmr_init <- cmdstanr::cmdstan_model("dmr.stan",
                                  cpp_options = list(
                                    stan_threads = TRUE
                                  )
)

最後に、実際の推定に入ります：

> m_dmr_estimate <- m_dmr_init$variational(
     seed = 12345,
     threads = 20,
     iter = 50000,
     data = data_list
 )
------------------------------------------------------------ 
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.045152 seconds 
1000 transitions using 10 leapfrog steps per transition would take 451.52 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   -281352741.767             1.000            1.000 
   200    -82829929.261             1.698            2.397 
   300    -30658642.199             1.699            1.702 
   400    -13274304.150             1.602            1.702 
   500     -6742206.004             1.475            1.310 
   600     -4114604.893             1.336            1.310 
   700     -3014265.021             1.197            1.000 
   800     -2535836.041             1.071            1.000 
   900     -2321733.373             0.962            0.969 
  1000     -2221423.233             0.871            0.969 
  1100     -2171287.969             0.794            0.639   MAY BE DIVERGING... INSPECT ELBO 
  1200     -2144325.219             0.729            0.639   MAY BE DIVERGING... INSPECT ELBO 
  1300     -2128373.529             0.673            0.365   MAY BE DIVERGING... INSPECT ELBO 
  1400     -2118453.814             0.625            0.365   MAY BE DIVERGING... INSPECT ELBO 
  1500     -2111813.533             0.584            0.189   MAY BE DIVERGING... INSPECT ELBO 
  1600     -2107187.370             0.547            0.189   MAY BE DIVERGING... INSPECT ELBO 
  1700     -2103740.303             0.515            0.092   MAY BE DIVERGING... INSPECT ELBO 
  1800     -2101053.222             0.487            0.092 
  1900     -2099140.874             0.461            0.045 
  2000     -2097381.564             0.438            0.045 
  2100     -2096064.506             0.417            0.023 
  2200     -2094808.775             0.398            0.023 
  2300     -2093823.669             0.381            0.013 
  2400     -2092890.451             0.365            0.013 
  2500     -2092187.455             0.351            0.007   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  236.7 seconds.

4分くらいで終わりました。

ここで可視化に必要な$\beta$を格納します：

beta_word <- m_dmr_estimate$summary("beta_word")

beta_word_with_meta <- beta_word |>
  dplyr::mutate(
    id = variable |>
      purrr::map(
        \(x){
          stringr::str_split(x, "\\[|\\]|,")[[1]][2:3]
        }
      )
  ) |>
  tidyr::unnest_wider(id, names_sep = "_") |>
  dplyr::mutate(
    dplyr::across(dplyr::starts_with("id"),
                  ~ as.integer(.))
  ) |>
  dplyr::left_join(word_master, by = c("id_2" = "word_id")) |>
  dplyr::left_join(category_master, by = c("id_1" = "category_id"))

推定結果の可視化

ここでは、各カテゴリの代表的な単語を上位30件抽出します。

> purrr::map(
     1:nrow(category_master),
     \(i){
         beta_word_with_meta |>
             dplyr::filter(category == category_master$category[i]) |>
             dplyr::arrange(dplyr::desc(mean)) |>
             dplyr::select(feature)
     }
 ) |>
     dplyr::bind_cols() |>
     `colnames<-`(category_master$category) |>
     dplyr::glimpse(width = 210)
New names:
• `feature` -> `feature...1`
• `feature` -> `feature...2`
• `feature` -> `feature...3`
• `feature` -> `feature...4`
• `feature` -> `feature...5`
• `feature` -> `feature...6`
• `feature` -> `feature...7`
• `feature` -> `feature...8`
• `feature` -> `feature...9`
Rows: 3,960
Columns: 9
$ `dokujo-tsushin` <chr> "アラサー", "広沢", "女子_力", "歳_会社", "ん_仮名", "エッセイ", "栗", "せつこ", "友人_たち", "柔軟", "人_結婚", "女子_会", "専業_主婦", "嘘", "覚悟", "カワイイ…
$ `it-life-hack`   <chr> "キーボード", "ブログ_虎巻", "岡本_奈知", "インク_キャノンインクタンク", "情報_皆ん", "満載_虎巻", "ハードディスク", "トップページ", "ロゴ", "大島_克彦", "バイ…
$ `kaden-channel`  <chr> "光学", "プロジェクター", "食品", "階段", "自動車", "天国_階段", "手入れ", "高画質", "発売_売れ筋", "機種_発売", "連載_天国", "グーグル", "ビデオカメラ", "フェ…
$ `livedoor-homme` <chr> "シャツ", "特集_関連", "ドライバー", "辛口_説教", "お答え", "トレーニング", "転職_将来", "自動車", "求人_転職", "ゴルフ", "ゴルファー", "仕事_若手", "求人", "若…
$ `movie-enter`    <chr> "映画_まとめ", "ジョン", "ジャパン_プレミア", "ピカ_デリー", "銃", "本_公開", "監督_最新", "回_公開", "オブ", "戦争", "話題_作", "劇場_話", "士", "ダーク_ナイト…
$ peachy           <chr> "鍋", "週_ランキング", "職_兼", "バレンタイン", "ワイン", "賞品_発送", "野菜", "脚", "花", "プラン", "ハリ", "ボタン_応募", "ホテル", "チョコレート", "了承_当選…
$ smax             <chr> "サムスン_電子", "画素_カメラ", "更新", "ソフトウェア_更新", "取扱_説明", "画面_部", "ディスプレイ_インチ", "電池_残", "以_注意", "キーボード", "連続_通話", "執…
$ `sports-watch`   <chr> "本田", "解説_者", "野村", "なでしこ", "サポーター", "岡田", "香川", "く", "金メダル", "問い", "格闘技", "サッカー_解説", "逆転", "松井", "球団", "斎藤", "陸", …
$ `topic-news`     <chr> "香川", "掲示板_物議", "高橋", "擁護", "逮捕", "市民", "ニュース_関連", "北朝鮮", "フジテレビ", "紗", "罪", "生徒", "週刊_誌", "島", "韓国_人", "疑惑", "自殺", …

独女通信？に特徴的な単語に「結婚」、「女子力」、「イケメン」などが入っています。movie-enterに特徴的な単語は「映画」、「戦争」、「話題作」など、映画関連の単語が入っています。topic newsには、「北朝鮮」、「韓国」などの単語が上位に来ています。

抽出結果から分かるように、このStanで再現した分散多項回帰モデルは、カテゴリに特徴的な単語を抽出することに成功しています。

結論

いかがでしたか？このように、Stanで実装した分散多項回帰があれば、簡単な改造だけで手元のデータに特化したモデルを構築することができます。やり方については、また別の記事で説明します。

ちなみに、2024年9月上旬に弊社のデータサイエンティスト新卒採用関連イベントに登壇します。興味ある学生の方はぜひご応募ください：

私たちと一緒に働きたい方はぜひ下記のリンクもご確認ください：

参考文献

Blei, David M., Andrew Y. Ng, and Michael I. Jordan. "Latent dirichlet allocation." Journal of Machine Learning Research 3.Jan (2003): 993-1022.

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

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