はじめに

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

本記事は、新しくでデータサイエンスチームにジョインした社員向けの研修資料で、テキストデータモデリングの第一弾になります。

最終的には高度なトピックモデルと埋め込み表現モデルまで紹介して、メンバー全員が手元のタスクに特化したテキストモデルを構築できるようしたいですが、基礎が大事なので、古典的なLDA（Latent Dirichlet Allocation）から説明したいと思います：

LDAによる文章生成の考え方

まずは、LDA系のトピックモデルがどのように文章を捉えているかを、直感的な例を使って説明しましょう。

文章はどうやって「生まれた」のか？

例えば、以下のような3つの求人原稿があったとします。

原稿A

子供が急に熱になった場合は早めに帰れます！経験がなくても大丈夫です！
原稿B

学業と部活優先で問題ありません！初めての方でもしっかりサポートします！
原稿C

未経験者歓迎！残業多めですがやろうと思えばしっかり稼げます！

いきなりですが、少し哲学的な問いをしてみます。これらの文章は、どういう「仕組み」で生まれたと考えられるでしょうか？

LDAの世界観：文章はトピックの混合でできている

LDAでは、各文章は複数のトピック（話題・主題）の混合でできていると考えます。たとえば、上記の原稿をLDA的に解釈すると、次のようになります：

原稿A：
「親の働きやすさ」と「未経験」の2つのトピックを持っていて、各トピックから単語をサンプリングして文章を作った。
原稿B：
「学生の両立支援」と「未経験」のトピックを持っていて、各トピックから単語をサンプリングして文章を作った。
原稿C：
「未経験」と「高収入」のトピックを持っていて、各トピックから単語をサンプリングして文章を作った。

ここでいう「サンプリング」とは、トピックごとに用意された単語の分布から、単語を1つずつ選んでいくイメージです。

単語はどう選ばれる：トピック → 単語の二段階

もう少し細かく見てみましょう。

例えば、原稿Aが「親の働きやすさ」60%、「未経験」40%というトピック構成だったとします。

最初の単語を書くとき、まずどのトピックから取るかをサンプリング（例：「親の働きやすさ」）。
次に、そのトピックの中の単語分布から具体的な単語をサンプリング（例：「子供」）。
次の単語も同様に、まずトピックを選び（例：「親の働きやすさ」）、次に単語を選ぶ（例：「急に熱」）。

これを繰り返して、文章全体が構築されていきます。また、ここで強調しておきたいのは、LDAは文書の中の単語の順番を無視しています。これはかなり大胆な仮定に見えますが、トピックの抽出が目的だったらそこまで性能に影響しないことを、この後の実証分析で確認していただければと思います。

観測から構造を逆算するのがLDA

LDAは、上記のような「見えない生成過程」を前提としながら、実際に観測された文章データから次の2つを同時に推定します：

各文書がどのトピックをどの程度含んでいるか（文書ごとのトピック分布）
各トピックにどんな単語がよく使われるか（トピックごとの単語分布）

つまり、観測された文章から、背後にある「話題の混合」と「単語の出現傾向」を逆算するモデルがLDAなのです。

次は、実際にどのようにしてこれらのトピックを推定していくのかを見ていきましょう。

LDAの確率モデル

では、上の段落で説明したモデルを実際の確率モデルに書き換えてみましょう！

まず、文章dのトピック分布$\theta_{d}$をサンプリングします：

$$
\gamma_{d} \sim Gamma(0.01, 0.01)
$$

$$
\theta_{d} \sim Dirichlet(\gamma_{d})
$$

次に、トピックpの単語分布$\beta_{p}$をサンプリングします：

$$
\delta_{p} \sim Gamma(0.01, 0.01)
$$

$$
\beta_{p} \sim Dirichlet(\delta_{p})
$$

普通のLDAだったら、$\gamma$と$\delta$はハイパーパラメータとして指定しますが、せっかくStanで実装するのでパラメータとして推定します。

次に、i番目の単語がどの単語かをサンプリングします。まず、i番目の単語は文書$doc_{i}$に属していますので、$\theta_{doc_{i}}$からトピックを引きます：

$$
z_{i} \sim Categorical(\theta_{doc_{i}})
$$

次に、引いたトピックの単語分布から、単語$word_{i}$を抽出します：

$$
word_{i} \sim Categorical(\beta_{z_{i}})
$$

ただし、Stanのように自動微分に依存するパッケージでは、$z_{i}$のような離散パラメータを直接扱うことができません。しかし、$z_{i}$のような変数が取りうる値が少数である場合は、総当たり（全ての値を列挙）によって対処することが可能です。

具体的には：

$z_i$ が取りうるすべての値に対して尤度を計算する
それらの対数尤度を log_sum_exp 関数に渡すことで、$z_i$ を周辺化できる

このアプローチにより、離散パラメータを含むモデルでもStanで扱えるようになります。

LDAのStan実装

Stanでは、LDAをこのように実装できます。処理速度のため並列処理機能とよりトピックの特徴的な単語を抽出できるとされるFREXを計算するパートを入れています：

LDA.stan

functions {
  real empirical_cumulative_density(
    real x, vector distribution
  ){
    vector[size(distribution)] result;
    for (i in 1:size(distribution)){
      if (x >= distribution[i]){
        result[i] = 1.0;
      } else {
        result[i] = 0.0;
      }
    }
    return mean(result);
  }
  
  real partial_sum_lpmf(
    array[] int freq,
    
    int start, int end,
    
    int topic_type,
    
    array[] vector doc_topic, array[] vector topic_word,
    
    array[] int word, array[] int doc
  ){
    vector[end - start + 1] log_likelihood;
    int count = 1;
    
    for (i in start:end){
      vector[topic_type] case_when;
      for (j in 1:topic_type){
        case_when[j] = log(doc_topic[doc[i], j]) + log(topic_word[j, word[i]]);
      }
      log_likelihood[count] = freq[count] * log_sum_exp(case_when);
      count += 1;
    }
    
    return sum(log_likelihood);
  }
}
data {
  int word_type;
  int doc_type;
  int topic_type;
  
  int N;
  array[N] int word;
  array[N] int doc;
  array[N] int freq;
}
parameters {
  vector<lower=0>[doc_type] sigma_doc;
  vector<lower=0>[topic_type] sigma_topic;
  
  array[doc_type] simplex[topic_type] doc_topic;
  array[topic_type] simplex[word_type] topic_word;
}
model{
  sigma_doc ~ gamma(0.01, 0.01);
  sigma_topic ~ gamma(0.01, 0.01);
  
  for (i in 1:doc_type){
    doc_topic[i] ~ dirichlet(rep_vector(sigma_doc[i], topic_type));
  }
  
  for (i in 1:topic_type){
    topic_word[i] ~ dirichlet(rep_vector(sigma_topic[i], word_type));
  }
  
  target += reduce_sum(
    partial_sum_lupmf, freq, 1,
    
    topic_type,
    
    doc_topic, topic_word,
    
    word, doc
  );
}
generated quantities {
  matrix[word_type, topic_type] frex;
  {
    matrix[word_type, topic_type] B;
    
    for (w in 1:word_type){
      for (t in 1:topic_type){
        B[w, t] = topic_word[t, w];
      }
    }
    
    for (w in 1:word_type){
      for (t in 1:topic_type){
        frex[w, t] = ((0.5/empirical_cumulative_density(B[w, t], B[w,:]')) + (0.5/empirical_cumulative_density(B[w, t], B[:,t])))^(-1);
      }
    }
  }
}

R言語とQuantedaによる前処理

今回の分析では、livedoor ニュースコーパスを利用します。

まず、テキストデータの読み込みとMecabで単語分割等を実施します：

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::`%>%`

review_df <- list.files("text") %>%
  .[which(stringr::str_detect(., ".txt") == FALSE)] |>
  purrr::map(
    \(this_file){
      this_file %>%
        stringr::str_c("text/", .)  |> 
        list.files() %>%
        .[which(stringr::str_detect(., "LICENSE") == FALSE)] %>% 
        stringr::str_c("text/", this_file, "/", .) |>
        purrr::map(
          \(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()
  )

次に、単語文書行列を縦持ちにするテータフレイムと単語のIDを作成するためのマスター表を作成します：

review_dfm_df <- review_df |>
  dplyr::pull(text) |>
  quanteda::phrase() |>
  quanteda::tokens() |>
  quanteda::tokens_ngrams(1:2) |>
  quanteda::dfm() |>
  quanteda::dfm_trim(min_docfreq = 200) |>
  quanteda::convert(to = "tripletlist") |>
  tibble::as_tibble() |>
  dplyr::mutate(
    document = document |>
      stringr::str_remove_all("text") |>
      as.integer()
  )

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

最後に、単語IDを縦持ちにした単語文書行列に付与します；

review_df_with_id <- review_dfm_df |>
  dplyr::left_join(
    word_master, by = "feature"
  )

ここで、Stanにデータを渡すための整形とモデルコンパイルを実施し：

data_list <- list(
  word_type = nrow(word_master),
  doc_type = max(review_df_with_id$document),
  topic_type = 10,
  
  N = nrow(review_df_with_id),
  word = review_df_with_id$word_id,
  doc = review_df_with_id$document,
  freq = review_df_with_id$frequency
)

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

実際のモデル推定に入ります：

> m_lda_estimate <- m_lda_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.184976 seconds 
1000 transitions using 10 leapfrog steps per transition would take 1849.76 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     -5406659.205             1.000            1.000 
   200     -5322185.569             0.508            1.000 
   300     -5266398.667             0.342            0.016 
   400     -5212980.305             0.259            0.016 
   500     -5166861.354             0.209            0.011 
   600     -5129535.148             0.175            0.011 
   700     -5098921.582             0.151            0.010 
   800     -5074103.829             0.133            0.010 
   900     -5053629.621             0.119            0.009   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  415.2 seconds.

7分で終わりました！推定結果の事後分布をデータフレイムとして保存します：

m_lda_summary <- m_lda_estimate$summary()

推定結果可視化

さて、まずはトピックの代表的な単語をFREXから確認しましょう：

> m_lda_summary |>
     dplyr::filter(stringr::str_detect(variable, "^frex\\[")) |>
     dplyr::mutate(
         id = variable |>
             purrr::map(
                 \(x){
                     as.integer(stringr::str_split(x, "\\[|\\]|,")[[1]][2:3])
                 }
             )
     ) |>
     tidyr::unnest_wider(id, names_sep = "_") |>
     dplyr::left_join(
         word_master, by = c("id_1" = "word_id")
     ) |>
     split(~ id_2) |>
     purrr::map(
         \(this_df){
             this_df |>
                 dplyr::arrange(dplyr::desc(mean)) |>
                 dplyr::select(feature)
         }
     ) |>
     dplyr::bind_cols() |>
     `colnames<-`(stringr::str_c("topic_", 1:10)) |>
     dplyr::glimpse()
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`
• `feature` -> `feature...10`
Rows: 908
Columns: 10
$ topic_1  <chr> "対応", "フォン", "スマート_フォン", "スマート", "搭載", "ドコモ", "発表", "向け", "エスマックス", "夏", "開始", "製品", "通信", "モバイル", "万", "バッテリー", "端末"…
$ topic_2  <chr> "画面", "機能", "設定", "通", "カメラ", "サービス", "デジタル", "動画", "無料", "入力", "デジ", "接続", "操作", "使用", "機器", "利用", "デジ_通", "起動", "ページ", "…
$ topic_3  <chr> "名", "応募", "ソフトバンク", "会", "ユーザー", "イベント", "期", "サイト", "モデル", "ブログ", "開催", "ブランド", "連絡", "プレゼント", "ー", "木", "駅", "社", "発売…
$ topic_4  <chr> "氏", "公開", "本", "ネット", "監督", "選手", "番組", "賞", "放送", "コメント", "映像", "掲示板", "ネット_掲示板", "出演", "試合", "アメリカ", "チーム", "声", "発言", …
$ topic_5  <chr> "映画", "作品", "位", "写真", "部", "作", "特集", "シーン", "登場", "本_作", "編集", "ランキング", "劇場", "夢", "週", "最", "舞台", "アクション", "全国", "テーマ", "…
$ topic_6  <chr> "韓国", "ドラマ", "田", "子", "ちゃん", "家族", "話題", "理由", "視聴", "率", "金", "関係", "挑戦", "問題", "関連", "本人", "インタビュー", "関連_記事", "光", "意味", …
$ topic_7  <chr> "更新", "ソフトウェア", "とき", "場合", "必要", "力", "方", "お客", "点", "情報", "実施", "際", "メール", "注意", "仕事", "事", "アップデート", "企業", "失敗", "効果",…
$ topic_8  <chr> "人", "女", "会社", "代", "独_女", "歳", "肌", "友達", "転職", "心", "独", "経験", "収", "男", "笑", "系", "話", "美", "人_人", "恋", "言葉", "営業", "大切", "クリスマ…
$ topic_9  <chr> "結婚", "男性", "好き", "も", "女性", "料理", "ん_歳", "女子", "生活", "恋愛", "気", "最近", "感", "モテ", "子供", "店", "味", "体", "子ども", "印象", "家", "そう", "…
$ topic_10 <chr> "アプリ", "記事", "簡単", "万_円", "チェック", "ゲーム", "性", "データ", "検索", "ポイント", "確認", "海", "登録", "システム", "状態", "商品", "家電", "終了", "以", "…

以下は、各トピックの語彙に基づく内容の解釈です。

トピック1：スマートフォンやモバイル端末に関する技術情報や製品発表を扱う話題。新機種のスペックや通信機能、発売時期などが中心。
トピック2：カメラや画面、設定、接続などのデジタル機器の機能説明や使い方ガイドに関連する内容。ユーザー向けマニュアルやレビュー記事を想起させる。
トピック3：キャンペーンやイベント、プレゼント応募など、企業プロモーションや販促活動に関連する話題。特にソフトバンクやモデル関連の話題が多い。
トピック4：芸能人・スポーツ選手・監督などの発言や出演、ネット上での反応に関するニュース。掲示板やコメント、放送に関連する語が多い。
トピック5：映画や演劇、特集記事など、エンタメ作品に関するランキング・紹介・批評などの内容。作品の舞台やシーンに関する詳細も含まれる。
トピック6：韓国ドラマや家族、視聴率、話題性といったワードが見られ、ドラマに関するレビューや芸能人の話題、視聴者の関心を反映した記事が想定される。
トピック7：ソフトウェア更新やセキュリティ対策、トラブル対応といったサポート・お知らせ情報。企業からユーザーへの注意喚起やFAQ的な内容が中心。
トピック8：独身女性や恋愛、転職、美容に関する生活情報や自己啓発的な話題。「独女」「恋」「営業」などの語から、ライフスタイル記事が想定される。
トピック9：結婚、恋愛、家庭、料理など、男女の関係や日常生活に関する記事。女性誌などに見られる恋愛相談やライフスタイル記事の雰囲気。
トピック10：アプリやポイント、チェック、検索など、実用系の記事やアプリ紹介、買い物情報などに関連する内容。使いやすさやお得情報に焦点がある。

全体として、ニュースサイトや求人メディアによくあるカテゴリが反映されており、LDAが自然にこうした意味的なまとまりを捉えていることが確認できます。

最後に、トピック同士の相関も確認しましょう：

m_lda_summary |>
  dplyr::filter(stringr::str_detect(variable, "^doc_topic\\[")) |>
  dplyr::pull(mean) |>
  matrix(ncol = 10) |>
  cor() |>
  as.data.frame() |>
  tibble::rowid_to_column(var = "topic_a") |>
  tibble::tibble() |>
  tidyr::pivot_longer(!topic_a, names_to = "topic_b", values_to = "correlation") |>
  dplyr::mutate(
    topic_a = factor(stringr::str_c("topic_", topic_a), levels = stringr::str_c("topic_", 1:10)),
    topic_b = factor(stringr::str_c("topic_", stringr::str_remove_all(topic_b, "V")), levels = stringr::str_c("topic_", 1:10))
  ) |>
  ggplot2::ggplot() +
  ggplot2::geom_tile(ggplot2::aes(x = topic_a, y = topic_b, fill = correlation)) +
  ggplot2::scale_fill_gradient2(
    low = "red",
    mid = "white",
    high = "blue",
    midpoint = 0
  )

ここでまず、トピック1・2・10の相関が高いことがわかります。つまり、文書においてトピック1の成分が高いと、トピック2やトピック10の成分も高くなる傾向があるということです。

実際にトピック1・2・10の内容を確認すると、いずれもデジタル・IT系の話題であり、同じ文書内で言及されるのは自然だと考えられます。

次に、トピック4・5・6の相関も高く、こちらは芸能系の話題が中心のため、同じ文書で扱われる可能性が高いでしょう。

最後に、トピック7・8・9・10の相関も強いですが、こちらは生活に密着した話題が多いのが特徴です。トピック10はデジタル・IT系でありながら、アプリやポイントなど生活に密着した内容も含んでいるため、トピック1・2との相関だけでなく、トピック7・8・9との相関も強くなっています。

結論

いかがでしょうか？このように、トピックモデルは単語の順序を無視する大胆なモデルですが、しっかりと解釈できるトピックを抽出することができます。

テキストデータモデルはテキストデータだけでなく、レコメンドエンジン、政治学の投票モデル、ないしは社会学と生物学のネットワークに共通する部分しかないので、きちんと押さえていくことで、モデリングの視野が広がると思います。

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