はじめに

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

本記事は以前の衆議院記事の拡張版であり、分析対象を1950年まで遡って再分析を行いました。ただ、少し残念なことに、得られた大きな結論には目立った変化が見られませんでした。

本研究の目的は、衆議院外務委員会の議事録データをもとに、単語の意味の変化を捉えることにあります。従来からこの目的に対しては様々なモデルや手法が提案されていますが、本記事で紹介するword2vec系モデルの最大の特徴は、単語×年のペアごとにベクトルを推定するのではなく、単語の意味ベクトル自体は時間を通じて固定であり、その重要度（重み）の方が年ごとに変化するという構造にあります。

このような設計により、パラメータ数を抑えつつ、意味構造の変化を効率的かつ解釈可能な形で捉えることが可能になります。

では早速、本格的な説明に入ります！

word2vecの仕組み

急な質問で恐縮ですが、皆さんは次の2文のうち、どちらが日本語として自然だと思いますか？

戦争は外交の失敗である

戦争はチーズタッカルビの失敗である

もちろん、様々な視点から異論もあるかもしれませんが、多くの方はおそらく前者の方が日本語として妥当だと感じるのではないでしょうか。その理由は、戦争のような政治的なテーマについて語っている文脈に、突如として「チーズタッカルビ」のような韓国料理の名前が登場するのは不自然だからです。

この違和感の正体は、単語ごとの共起関係（同時出現）にあります。つまり、ある単語が登場したときに、一緒に現れると自然な単語と、現れると不自然な単語があるということです。

実際、Mikolov, Chen, Corrado and Dean（2013）が考案したword2vecをはじめとする自然言語処理のモデルでは、このような共起の自然さをデータから学習させます。たとえば、今回用いたような衆議院外務委員会の議事録データをもとに、実際に出現する単語のペアを「自然な文脈」として学習し、一方で、実際には出現しないような組み合わせ（例：「戦争」と「チーズタッカルビ」）をノイズとして挿入することで、モデルに意味の構造を学ばせることができます。

本記事で紹介するモデルは概ねword2vecの考えに基づいていますが、ベイズモデルのため、Rudolph, Ruiz, Mandt and Blei（2016）の構造も参考にしています。

モデル構造と問題意識

ここでは、モデルの構造となる尤度関数と事前分布について説明しますが、その前に、本モデルがそもそもどのような課題を解決しようとしているのかを簡単に述べます。

本モデルが取り組む第一の課題は、word2vec系のモデルにおける次元数の推定です。既存の研究では、次元数を理論的な根拠なしに固定したり、複数の次元数でモデルを個別に推定して予測精度を比較したりする方法が一般的です。しかし、ベイジアンの立場から見ると、これはあまり筋のよいアプローチとは言えません。

そこで本モデルでは、次元の「数」ではなく、それぞれの次元の「重要度」自体をパラメータとして推定するという発想を取り入れています。これにより、次元の数を事前に固定する必要がなく、自然に次元のスパース性や意味的な構造が表現されるようになります。

さらに、単語ベクトルについても、あらかじめ平均や分散を固定した単純な正規分布に従わせるのではなく、単語ベクトルの事前分布そのものをデータから柔軟に推定できる構造にしています。これにより、より多様な語義構造や語彙の分布的性質を捉えることが期待されます。

尤度関数

まず、前述の「戦争は外交の失敗である」と、（デタラメである）「戦争はチーズタッカルビの失敗である」を比較することで、以下のようなデータ構造を想定します：

正解フラグ	目標単語	文脈1	文脈2
1	外交	戦争	失敗
0	チーズタッカルビ	戦争	失敗

ここで「正解フラグ」は、その行（レコード）が実際のデータに登場したものか（= 1）、あるいは分析者が意図的に挿入したノイズか（= 0）を示します。

「目標単語」はスコアを推定したい中心語であり、「文脈」はその周囲に現れた語句です。

このデータをもとに、目標単語の無限次元の意味ベクトル（embedding）、文脈単語の無限次元のベクトル（context）、および年ごとの無限次元の次元重要度（importance）を用いて、以下のようにスコア$\theta$を定義します：

$$
\theta = \sum_{d=1}^{\infty}(importance_{year_{i}, d} \cdot embedding_{word_{i},d} \cdot \sum_{j \in c_{i}}context_{j,d})
$$

このスコア $\theta_i$ に対してロジット変換を施し、正解ラベルをベルヌーイ分布から生成するモデルです：

$$
正解フラグ_{i} \sim Bernoulli(logit(\theta_{i}))
$$

この仕組みによって、「文脈にふさわしい単語」と「文脈に不自然な単語」を区別するための意味空間と年次変化が学習されます。さらに、単語ベクトルを年ごとの次元重要度で重み付けすることで、単語の意味の変化を経時的に追うことが可能になります。

さて、数学的素養のある読者であれば、「$\theta$における無限和は本当に定義できるのか？」という疑問を持つかもしれません。確かに、無限和は常に収束するとは限りません。たとえば、数列1, -1, 1, -1, ...の和を考えると、項数が偶数のときは0、奇数のときは1になりますが、「無限大が奇数か偶数か」は定義されていないため、この無限和は収束しません。

しかし本モデルでは、後述するように、年ごとの次元重要度ベクトルに対して「棒折り過程」と呼ばれるノンパラメトリック・ベイズ構造を導入しています。これは、重要度が次元を追うごとに指数関数的にゼロへ収束していく特性を持っており、そのおかげで$\theta$の無限和も問題なく定義できるようになっています。

事前分布

ここの内容は、基本的には前回の記事の内容をそのまま利用しています。

単語埋め込みベクトル

単語ベクトルの事前分布には、正規分布を設定することが多いですが、本当に正規分布で良いのか？という疑問に回答するため、正規分布から成すディリクレ過程で設定しています。

まずは全体のハイパーパラメーター$\alpha$をガンマ分布からサンプリングします。

$$
\alpha \sim Gamma(0.001, 0.001)
$$

次に、棒折り過程でディリクレ過程を構築します。無限にあるクラスターの中のクラスターpについてこのように必要な変数をサンプリングします：

$$
\pi_{p}\sim Beta(1, \alpha)
$$

$$
p_{p} = \pi_{p} \prod\limits_{l=1}^{p - 1} (1 - \pi_{l})
$$

また、クラスターpの中心ベクトルとそのばらつきは下記のようにサンプリングされます：

$$
P_{latent,p} \sim Normal(0,1)
$$

$$
P_{\sigma,p} \sim Gamma(0.001, 0.001)
$$

単語Sのベクトルをサンプリングする際は、まず上記の棒折り過程で構築した無限次元の確率ベクトルをパラメーターとするカテゴリ分布から、インデックスをサンプリングします：

$$
\eta_{S} \sim Categorical(p)
$$

次に、サンプリングされた$\eta_{s}$を基に、ベクトルをサンプリングします：

$$
embedding_{S} \sim Normal(P_{latent, \eta_{S}}, P_{\sigma, \eta_{S}})
$$

この一連のサンプリングのプロセス分布を$G$で表します。

年ごとの次元重要度

次元数は自己回帰的な階層棒折り過程で定式化します。まず、最初の年の重要度ベクトルを棒折り過程でサンプリングします。具体的には、

$$
\beta \sim Gamma(0.001, 0.001)
$$

をサンプリングし、次に棒折り過程で$年重要度ベクトル_{1,d}$を求める

$$
\delta_{1,d}\sim Beta(1, \beta)
$$

$$
importance_{1,d} = \delta_{1,d} \prod\limits_{l=1}^{d - 1} (1 - \delta_{1,l})
$$

これをd = 1からd = $\infty$まで繰り返します。

次に、経時的変化の度合いを示すパラメータをサンプリングします：

$$
\theta\sim Gamma(0.001, 0.001)
$$

そして、Y年の重要度ベクトルを、Y-1年の重要度ベクトルを「事前分布」とする階層棒おり過程からサンプリングします。

$$
\delta_{Y,d} \sim Beta\left( \theta \delta_{Y-1,d}, \theta \left(1 - \sum_{l=1}^{d} \delta_{Y-1,l} \right) \right)
$$

$$
importance_{Y,d} = \delta_{Y,d} \prod\limits_{l=1}^{d - 1} (1 - \delta_{Y,l})
$$

これもd = 1からd = $\infty$まで繰り返します。

どうして階層化するのかというと、前の年の次元重要度のベクトルをこの年の年重要度ベクトルの事前分布にしないと、次元重要度のベクトル同士の重なりがなくなり、うまく経時的変化を捉えられません。

さて、詳細は筆者の別の論文のAppendix（Wu, 2025）を参照していただきたいですが：

年Yの次元重要度のベクトルのd番目の次元の期待値はこのように書くことができます：

$$
\mathbb{E}[importance_{Y,d}] = \frac{1}{1 + \beta} \left( \frac{\beta}{1 + \beta} \right)^{d-1}.
$$

$\beta$がガンマ分布から生成されているため、0より大きい実数で、$\frac{\beta}{1+\beta}$は0から1までの間の実数になります。したがって、年Yの次元重要度のベクトルは、次元を追うごとに指数関数的にゼロへ収束していく特性を持っていることがわかります。

このように、ノンパラメトリック・ベイズの代表的な構造である棒折り過程を用いることで、無限次元の重要度ベクトルであっても収束性とモデルの数学的健全性が保たれているだけでなく、急速に値がゼロに収束する特性は次元の変数選択の機能も果たすこともできます。

既存モデルとの比較

Rudolph and Blei（2018）は、単語の意味の変化を加味したword2vecを提案しています。具体的には、このようなモデル形式です：

$$
embedding_{Y} \sim Normal(embedding_{Y - 1}, \lambda)
$$

ただし、従来の手法では単語ベクトルを年ごとに個別に推定する必要があるため、計算資源に大きな負担がかかります。たとえば、年数を$T$、次元数を$D$、語彙サイズを$W$とした場合、Rudolph and Blei（2018）の手法では$T\times D \times W$個のパラメータを保持する必要があり、特に年数や単語数が多い場合には、メモリやストレージの限界にすぐ直面します。

それに対して、本記事で提案するモデルは、$(T + W) \times D$次元だけでデータを記述できる構造を採用しています。もちろん、理論上は$D$（次元数）を無限大とするベイズ的発想に基づいていますが、実際の推定においては、十分に大きな$D$を選ぶことで上限を設ける実装上の工夫が必要です。

この構造により、必要なパラメータ数が劇的に削減され、スケーラビリティの面でも実用的なモデルとなっています。

Stanコード

実装用のStanコードはこちらです：

dynamic_word2vec.stan

functions {
  vector stick_breaking(vector breaks){
    int length = size(breaks) + 1;
    vector[length] result;
    
    result[1] = breaks[1];
    real summed = result[1];
    for (d in 2:(length - 1)) {
      result[d] = (1 - summed) * breaks[d];
      summed += result[d];
    }
    result[length] = 1 - summed;
    
    return result;
  }
  
  real weighted_inner_product(
    vector vec_a, vector vec_b, vector weight
  ){
    vector[size(vec_a)] result_record;
    for (i in 1:size(vec_a)){
      result_record[i] = vec_a[i] * vec_b[i] * weight[i];
    }
    return sum(result_record);
  }
  
  real entropy(
    vector distribution, real threshold
  ){
    vector[size(distribution)] result;
    
    for (i in 1:size(distribution)){
      if (distribution[i] < threshold){
        result[i] = 0.0;
      } else {
        result[i] = distribution[i] * log(distribution[i]);
      }
    }
    
    return -sum(result);
  }
  
  real embedding_prior_lpmf(
    array[] int word_seq,
    int start, int end,
    
    int group_type,
    
    array[] vector word_embedding,
    
    vector group,
    array[] vector group_embedding, vector group_sigma
  ){
    vector[end - start + 1] lambda;
    int count = 1;
    for (w in start:end){
      vector[group_type] case_vector;
      for (g in 1:group_type){
        case_vector[g] = log(group[g]) + normal_lpdf(word_embedding[w] | group_embedding[g], group_sigma[g]);
      }
      lambda[count] = log_sum_exp(case_vector);
      count += 1;
    }
    return sum(lambda);
  }
  
  real partial_sum_lpmf(
    array[] int flag,
    int start, int end,
    
    array[] int year,
    array[] int word,
    array[] int word_lead_1,
    array[] int word_lag_1,
    array[] int word_lead_2,
    array[] int word_lag_2,
    array[] int word_lead_3,
    array[] int word_lag_3,
    
    array[] vector dimension_year,
    array[] vector word_embedding,
    array[] vector word_context
  ){
    vector[end - start + 1] lambda;
    int count = 1;
    for (i in start:end){
      lambda[count] = weighted_inner_product(
        word_embedding[word[i]],  
        word_context[word_lead_1[i]] + word_context[word_lag_1[i]] + 
        word_context[word_lead_2[i]] + word_context[word_lag_2[i]] + 
        word_context[word_lead_3[i]] + word_context[word_lag_3[i]], 
        dimension_year[year[i]]);
      count += 1;
    }
    return bernoulli_logit_lupmf(flag | lambda);
  }
}
data {
  int<lower=1> word_type;
  int<lower=1> year_type;
  int<lower=1> group_type;
  int<lower=1> dimension_type;
  
  array[word_type] int<lower=1,upper=word_type> word_seq;
  
  int<lower=1> N;
  array[N] int<lower=1,upper=year_type> year;
  array[N] int<lower=1,upper=word_type> word;
  array[N] int<lower=1,upper=word_type> word_lead_1;
  array[N] int<lower=1,upper=word_type> word_lag_1;
  array[N] int<lower=1,upper=word_type> word_lead_2;
  array[N] int<lower=1,upper=word_type> word_lag_2;
  array[N] int<lower=1,upper=word_type> word_lead_3;
  array[N] int<lower=1,upper=word_type> word_lag_3;
  array[N] int<lower=0,upper=1> flag;
  
  int<lower=0> val_N;
  array[val_N] int<lower=1,upper=year_type> val_year;
  array[val_N] int<lower=1,upper=word_type> val_word;
  array[val_N] int<lower=1,upper=word_type> val_word_lead_1;
  array[val_N] int<lower=1,upper=word_type> val_word_lag_1;
  array[val_N] int<lower=1,upper=word_type> val_word_lead_2;
  array[val_N] int<lower=1,upper=word_type> val_word_lag_2;
  array[val_N] int<lower=1,upper=word_type> val_word_lead_3;
  array[val_N] int<lower=1,upper=word_type> val_word_lag_3;
  array[val_N] int<lower=0,upper=1> val_flag;
}
parameters{
  real<lower=0> dimension_global_alpha;
  
  real<lower=0> dimension_across_year_alpha;
  array[year_type] vector<lower=0, upper=1>[dimension_type - 1] dimension_year_breaks;
  
  real<lower=0> group_alpha;                                       // ディリクレ過程の全体のパラメータ
  vector<lower=0, upper=1>[group_type - 1] group_breaks;  // ディリクレ過程のstick-breaking representationのためのパラメータ
  
  vector<lower=0>[group_type] group_sigma;
  array[group_type] vector[dimension_type] group_embedding;
  array[word_type] vector[dimension_type] word_embedding;
  array[word_type] vector[dimension_type] word_context;
}
transformed parameters {
  array[year_type] simplex[dimension_type] dimension_year;
  simplex[group_type] group;
  
  for (t in 1:year_type){
    dimension_year[t] = stick_breaking(dimension_year_breaks[t]);
  }
  
  group = stick_breaking(group_breaks);
}
model{
  dimension_global_alpha ~ gamma(0.001, 0.001);
  
  dimension_year_breaks[1] ~ beta(1, dimension_global_alpha);
  
  dimension_across_year_alpha ~ gamma(0.001, 0.001);
  
  for (t in 2:year_type){
    for (d in 1:(dimension_type - 1)){
      dimension_year_breaks[t, d] ~ beta(dimension_across_year_alpha * dimension_year[t - 1, d], dimension_across_year_alpha * (1 - sum(dimension_year[t - 1, 1:d])));
    }
  }
  
  group_alpha ~ gamma(0.001, 0.001);
  
  group_breaks ~ beta(1, group_alpha);
  
  group_sigma ~ gamma(0.001, 0.001);
  
  for (g in 1:group_type){
    group_embedding[g] ~ normal(0, 1);
  }
  
  int grainsize = 1;
  
  target += reduce_sum(
    embedding_prior_lupmf, word_seq, grainsize,
    
    group_type,
    
    word_embedding,
    
    group,
    group_embedding, group_sigma
  );
  
  target += reduce_sum(
    partial_sum_lupmf, flag, grainsize,
    
    year,
    word,
    word_lead_1,
    word_lag_1,
    word_lead_2,
    word_lag_2,
    word_lead_3,
    word_lag_3,
    
    dimension_year,
    word_embedding,
    word_context
  );
}
generated quantities {
  array[dimension_type] real drawn_G;
  array[word_type] vector[group_type] estimated_eta;
  vector[year_type] entropy_per_year;
  real F1_score;
  {
    int sampled_group = categorical_rng(group);
    drawn_G = normal_rng(group_embedding[sampled_group], group_sigma[sampled_group]);
    
    for (w in 1:word_type){
      vector[group_type] case_vector;
      for (g in 1:group_type){
        case_vector[g] = log(group[g]) + normal_lpdf(word_embedding[w] | group_embedding[g], group_sigma[g]);
      }
      estimated_eta[w] = softmax(case_vector);
    }
    
    for (i in 1:year_type){
      entropy_per_year[i] = entropy(dimension_year[i], 0.00001);
    }
    
    array[val_N] int predicted;
    int TP = 0;
    int FP = 0;
    int FN = 0;
    
    for (i in 1:val_N){
      predicted[i] = bernoulli_logit_rng(
        weighted_inner_product(
          word_embedding[val_word[i]],  
          word_context[val_word_lead_1[i]] + word_context[val_word_lag_1[i]] + 
          word_context[val_word_lead_2[i]] + word_context[val_word_lag_2[i]] + 
          word_context[val_word_lead_3[i]] + word_context[val_word_lag_3[i]], 
          dimension_year[val_year[i]])
      );
      if (val_flag[i] == 1 && predicted[i] == 1){
        TP += 1;
      }
      else if (val_flag[i] == 0 && predicted[i] == 1){
        FP += 1;
      }
      else if (val_flag[i] == 1 && predicted[i] == 0){
        FN += 1;
      }
    }
    F1_score = (2 * TP) * 1.0/(2 * TP + FP + FN);
  }
}

前処理とモデル推定

さて、ここからはデータをStanに渡す前の前処理を実施します：

`%>%` <- magrittr::`%>%`
mecabbing <- function(text){
  text |>
    stringr::str_replace_all("[^一-龠ぁ-んーァ-ヶーa-zA-Z]", " ") |>
    RMeCab::RMeCabC(1) |>
    unlist() %>%
    .[which(names(.) == "名詞")] |>
    stringr::str_c(collapse = " ")
}

#　データのダウンロード
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 <- quanteda::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, 1950, 2017)) |>
  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), .progress = TRUE)) |>
  dplyr::mutate(
    text_id = dplyr::row_number()
  )

sub_corp_use_tokens <- sub_corp_df |>
  dplyr::pull(text_mecabbed) |>
  quanteda::phrase() |>
  quanteda::tokens() |>
  quanteda::tokens_remove("の") |>
  quanteda::dfm() |>
  quanteda::dfm_trim(min_termfreq = 100, max_termfreq = 200000) |>
  colnames()

word_master <- tibble::tibble(
  word = sub_corp_use_tokens
) |>
  dplyr::arrange(word) |>
  dplyr::mutate(
    word_id = dplyr::row_number()
  )

year_master <- sub_corp_df |>
  dplyr::select(record_year) |>
  dplyr::distinct() |>
  dplyr::arrange(record_year) |>
  dplyr::mutate(
    year_id = dplyr::row_number()
  )

sub_corp_long_df_lag_correct <- sub_corp_df |>
  dplyr::select(text_id, record_year, text_mecabbed) |>
  dplyr::mutate(
    word = text_mecabbed |>
      quanteda::phrase() |>
      quanteda::tokens() |>
      quanteda::tokens_remove("の") |>
      as.list()
  ) |>
  tidyr::unnest_longer(word) |>
  dplyr::filter(word %in% word_master$word) |>
  dplyr::select(text_id, record_year, word) |>
  dplyr::left_join(word_master, by = "word") |>
  dplyr::left_join(year_master, by = "record_year") |>
  dplyr::group_by(text_id) |>
  dplyr::mutate(
    word_lag_1 = dplyr::lag(word_id, 1),
    word_lead_1 = dplyr::lead(word_id, 1),
    word_lag_2 = dplyr::lag(word_id, 2),
    word_lead_2 = dplyr::lead(word_id, 2),
    word_lag_3 = dplyr::lag(word_id, 3),
    word_lead_3 = dplyr::lead(word_id, 3)
  ) |>
  dplyr::ungroup() |>
  tidyr::drop_na()

set.seed(12345)

sub_corp_long_df_lag_correct_and_neg <- sub_corp_long_df_lag_correct |>
  dplyr::mutate(ans = 1) |>
  dplyr::bind_rows(
    sub_corp_long_df_lag_correct |>
      dplyr::mutate(word_id = sample(sub_corp_long_df_lag_correct$word_id, length(sub_corp_long_df_lag_correct$word_id), replace = FALSE)) |>
      dplyr::mutate(ans = 0)
  )

set.seed(12345)
val_id <- c(
  sample(which(sub_corp_long_df_lag_correct_and_neg$ans == 1), 5000, replace = FALSE),
  sample(which(sub_corp_long_df_lag_correct_and_neg$ans == 0), 5000, replace = FALSE)
)

data_list <- list(
  word_type = nrow(word_master),
  year_type = nrow(year_master),
  group_type = 10,
  dimension_type = 20,
  
  word_seq = 1:nrow(word_master),
  
  N = nrow(sub_corp_long_df_lag_correct_and_neg[-val_id,]),
  year = sub_corp_long_df_lag_correct_and_neg$year_id[-val_id],
  word = sub_corp_long_df_lag_correct_and_neg$word_id[-val_id],
  word_lead_1 = sub_corp_long_df_lag_correct_and_neg$word_lead_1[-val_id],
  word_lag_1 = sub_corp_long_df_lag_correct_and_neg$word_lag_1[-val_id],
  word_lead_2 = sub_corp_long_df_lag_correct_and_neg$word_lead_2[-val_id],
  word_lag_2 = sub_corp_long_df_lag_correct_and_neg$word_lag_2[-val_id],
  word_lead_3 = sub_corp_long_df_lag_correct_and_neg$word_lead_3[-val_id],
  word_lag_3 = sub_corp_long_df_lag_correct_and_neg$word_lag_3[-val_id],
  flag = sub_corp_long_df_lag_correct_and_neg$ans[-val_id],
  
  val_N = nrow(sub_corp_long_df_lag_correct_and_neg[val_id,]),
  val_year = sub_corp_long_df_lag_correct_and_neg$year_id[val_id],
  val_word = sub_corp_long_df_lag_correct_and_neg$word_id[val_id],
  val_word_lead_1 = sub_corp_long_df_lag_correct_and_neg$word_lead_1[val_id],
  val_word_lag_1 = sub_corp_long_df_lag_correct_and_neg$word_lag_1[val_id],
  val_word_lead_2 = sub_corp_long_df_lag_correct_and_neg$word_lead_2[val_id],
  val_word_lag_2 = sub_corp_long_df_lag_correct_and_neg$word_lag_2[val_id],
  val_word_lead_3 = sub_corp_long_df_lag_correct_and_neg$word_lead_3[val_id],
  val_word_lag_3 = sub_corp_long_df_lag_correct_and_neg$word_lag_3[val_id],
  val_flag = sub_corp_long_df_lag_correct_and_neg$ans[val_id]
)

学習データの量を確認しましょう：

> nrow(sub_corp_long_df_lag_correct_and_neg[-val_id,])
[1] 21545478

2154万件です。これは、おそらく単純なベイズモデルから一歩も踏み出したことのない方には想像もつかない、本物のビッグデータです。ここでMCMCにこだわってしまうと、おそらく計算が終わるのは私たち30代の「ひ孫」の世代になるでしょう。ですが、データ分析はウイスキーづくりではありません。百年単位の時間がかかるような分析は、ビジネスでもアカデミアでも現実的に利用できません。そのため本記事では素直に、分散処理に対応した変分推論 を用いて、効率的にモデルを推定します。

ではここでモデルのコンパイルと推定に入りましょう：

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


m_dw2v_estimate <- m_dw2v_init$variational(
  seed = 12345,
  threads = 6,
  data = data_list
)

推定のログと所要時間はこちらです：

> m_dw2v_estimate$output()

method = variational
  variational
    algorithm = meanfield (Default)
      meanfield
    iter = 10000 (Default)
    grad_samples = 1 (Default)
    elbo_samples = 100 (Default)
    eta = 1 (Default)
    adapt
      engaged = 1 (Default)
      iter = 50 (Default)
    tol_rel_obj = 0.01 (Default)
    eval_elbo = 100 (Default)
    output_samples = 1000 (Default)
id = 1 (Default)
data
  file = /var/folders/yp/7nvktzd95xxg_nbx7mdvxy1w0000gn/T/Rtmp2yo9kT/standata-23316a124e9e.json
init = 2 (Default)
random
  seed = 12345
output
  file = /var/folders/yp/7nvktzd95xxg_nbx7mdvxy1w0000gn/T/Rtmp2yo9kT/dynamic_word2vec-202507062237-1-42855f.csv
  diagnostic_file =  (Default)
  refresh = 100 (Default)
  sig_figs = -1 (Default)
  profile_file = /var/folders/yp/7nvktzd95xxg_nbx7mdvxy1w0000gn/T/Rtmp2yo9kT/dynamic_word2vec-profile-202507062237-1-500fbe.csv
  save_cmdstan_config = 0 (Default)
num_threads = 6 (Default)

------------------------------------------------------------
EXPERIMENTAL ALGORITHM:
  This procedure has not been thoroughly tested and may be unstable
  or buggy. The interface is subject to change.
------------------------------------------------------------



Gradient evaluation took 18.3348 seconds
1000 transitions using 10 leapfrog steps per transition would take 183348 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    -14614068.854             1.000            1.000
   200    -13449376.898             0.543            1.000
   300    -12694258.273             0.382            0.087
   400    -12140446.577             0.298            0.087
   500    -11805548.264             0.244            0.059
   600    -11591744.661             0.206            0.059
   700    -11442847.822             0.179            0.046
   800    -11331364.053             0.158            0.046
   900    -11246617.256             0.141            0.028
  1000    -11179069.383             0.127            0.028
  1100    -11123898.808             0.028            0.018
  1200    -11077360.960             0.020            0.013
  1300    -11038095.746             0.014            0.010   MEDIAN ELBO CONVERGED

Drawing a sample of size 1000 from the approximate posterior... 
COMPLETED

> m_dw2v_estimate$time()
$total
[1] 47343.49

47343秒、13時間かかりました。Macさん、よく頑張りました！

推定結果

年ごとの次元重要度

では、まずは年ごとの次元重要度を可視化しましょう：

m_dw2v_summary |>
  dplyr::filter(stringr::str_detect(variable, "dimension_year\\[")) |>
  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(
    year_master, by = c("id_1" = "year_id")
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = record_year, y = mean, color = as.factor(id_2))) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = record_year, ymin = q5, ymax = q95, fill = as.factor(id_2)), alpha = 0.3) +
  ggplot2::scale_x_continuous(breaks = seq(min(year_master$record_year), max(year_master$record_year), by = 2)) +
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + 
  ggplot2::labs(
    title = "Change of Dimensions' Importance Across Time",
    x = "Year",
    y = "Importance",
    color = "dimension id",
    fill = "dimension id"
  ) + 
  ggplot2::geom_vline(xintercept = 1989, linetype = "dashed", color = "blue", alpha = 0.6, linewidth = 1.5) + 
  ggplot2::annotate("text", x = 1989, y = Inf, label = "End of Cold War", vjust = 2, color = "blue")

おーーーーーーーーーーーーーーい！！！！！めっちゃ安定してるじゃん！！！！！時間と電気代返してくれ😭

1954年から1955年までのジャンプを除くと離散的な変化があまりないようです。約70年分のデータを13時間かけた推定したのでもっと急激な変化を期待していましたが、この70年の間で少なくとも衆議院外務委員会の出席者が使う単語の意味に大きな変化がなかったと言えそうです。また、最大で20次元まで使えると設定しましたが、結果として12次元から13次元しか利用されていないようです。

ちなみに、筆者が同じモデル構造で国連での投票行動を分析した際はもっと激しく上下した結果が出たので：

次元の重要度の推移がほぼ直線になっているのはこのモデルの性質ではなく、どちらかというとデータから発見した傾向として捉えた方がいいでしょう。

次に、年ごとの次元重要度ベクトルのエントロピーを計算して可視化します：

m_dw2v_summary |>
  dplyr::filter(stringr::str_detect(variable, "^entropy_per_year\\[")) |>
  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(
    year_master, by = c("id_1" = "year_id")
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = record_year, y = mean)) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = record_year, ymin = q5, ymax = q95), fill = ggplot2::alpha("blue", 0.3)) +
  ggplot2::scale_x_continuous(breaks = seq(min(year_master$record_year), max(year_master$record_year), by = 2)) +
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + 
  ggplot2::labs(
    title = "Change of Entropy Across Time",
    x = "Year",
    y = "Entropy"
  ) + 
  ggplot2::geom_vline(xintercept = 1989, linetype = "dashed", color = "blue", alpha = 0.6, linewidth = 1.5) + 
  ggplot2::annotate("text", x = 1989, y = Inf, label = "End of Cold War", vjust = 2, color = "blue")

エントロピーで見たら、1950年から1954年までは複雑性が低かったんですが、55年体制の始まりである1955年にエントロピーが急上昇して、それから長期的に下落していきます。

単語ベクトル

ここではまず、単語埋め込みベクトルが何種類のクラスターから生成されたのかを確認しましょう：

m_dw2v_summary |>
  dplyr::filter(stringr::str_detect(variable, "^group\\[")) |>
  dplyr::mutate(
    id = dplyr::row_number()
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_bar(ggplot2::aes(x = as.factor(id), y = mean), 
                    stat = "identity", fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::geom_errorbar(
    ggplot2::aes(x = as.factor(id), ymin = q5, ymax = q95),
    width = 0.2
  ) + 
  ggplot2::labs(
    title = "Posterior Distribution of Group",
    x = "Group",
    y = "distribution")

主に二つのグループですね！では実際は利用されていないグループも含めて、各グループの代表的な単語を可視化しましょう：

m_dw2v_summary |>
  dplyr::filter(stringr::str_detect(variable, "^estimated_eta\\[")) |>
  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(-mean) |>
        dplyr::slice_head(n = 20) |>
        dplyr::select(word, mean) |>
        dplyr::mutate(
          group = this_df$id_2[1],
          rid = dplyr::row_number() |>
            rev()
        )
    }
  ) |>
  dplyr::bind_rows() |>
  ggplot2::ggplot() +
  ggplot2::geom_text(ggplot2::aes(x = as.factor(group), y = rid, label = word), color = "blue", family = "HiraKakuPro-W3") +
  ggplot2::theme_minimal(base_family = "HiraKakuPro-W3") +
  ggplot2::theme(
    panel.grid.major = ggplot2::element_blank(),
    panel.grid.minor = ggplot2::element_blank(),
    axis.title.y = ggplot2::element_blank(),
    axis.text.y = ggplot2::element_blank(),
    axis.ticks.y = ggplot2::element_blank()
  ) + 
  ggplot2::labs(
    title = "Representative words for each group",
    x = "Group"
  )

グループ1は実際の外交問題に関する単語で、グループ2は人名がメインになっています。実際の外交問題の単語と登場人物の名前は意味的に大きく異なるため、別の分布が設定されるのは納得できる結果だといえるでしょう。

次に、Tsneの図で学習した単語ベクトルを可視化します。Tsneの図における距離が近ければ、単語の意味が近いという形で学習結果の妥当性を確認します；

set.seed(12345)
tsne_df <- m_dw2v_summary |>
  dplyr::filter(stringr::str_detect(variable, "^word_embedding\\[")) |>
  dplyr::pull(median) |>
  matrix(ncol = 20) |>
  rbind(m_dw2v_estimate$draws("drawn_G") |>
          tibble::as_tibble() |>
          dplyr::mutate(
            dplyr::across(
              dplyr::everything(),
              ~ as.numeric(.)
            )
          ) |>
          as.matrix()
        ) |>
  Rtsne::Rtsne() |>
  purrr::pluck("Y") |>
  tibble::as_tibble() |>
  dplyr::bind_cols(
    type = c(rep("word", nrow(word_master)), rep("sample", 1000)),
    name = c(word_master$word, rep("sample", 1000)),
  ) 

g_tsne <- tsne_df |>
  dplyr::filter(type == "word") |>
  ggplot2::ggplot() +
  ggplot2::geom_text(ggplot2::aes(x = V1, y = V2, label = name), color = ggplot2::alpha("blue", 0.5)) + 
  ggplot2::geom_density_2d(data = tsne_df |>
                             dplyr::filter(type == "sample"),
                           ggplot2::aes(x = V1, y = V2),
                           color = ggplot2::alpha("black", 0.2),
                           show.legend = FALSE
  )

plotly::ggplotly(g_tsne)

詳しく見ると、このように、新聞と曜日が近辺にあり、これは日本語の一般的な特性を捉えた結果だと思われます：

他にも、他国との領土紛争の対象になっている島が近辺にあることも確認できます：

条約系の内容も見られます：

単語の意味の経時変化

最後に、単語の経時変化を、単語ベクトルと年ごとの次元重要度の重みづけで確認します。まずは重み付けされたベクトルを保存するリストとコサイン類似度を計算するための関数をそれぞれ用意します：

word_embedding_matrix_weighted <- year_master |>
  nrow() |>
  seq_len() |>
  purrr::map(
    \(i){
      word_master |>
        nrow() |>
        seq_len() |>
        purrr::map(
          \(j){
            word_embedding_matrix[j,] * dimension_year_matrix[i,]
          }
        ) |>
        do.call(what = rbind, args = _)
    }
  )

cos_sim <- function(a, b){
  as.numeric(a %*% b)/sqrt(sum(a ^ 2) * sum(b ^ 2))
}

まずは、「安保」の意味の変化から確認しましょう。色の濃さはコサイン類似度を表しています。

year_master |>
  nrow() |>
  seq_len() |>
  purrr::map(
    \(i){
      this_matrix <- word_embedding_matrix_weighted[[i]]
      
      this_matrix |>
        nrow() |>
        seq_len() |>
        purrr::map_dbl(
          \(i){
            cos_sim(this_matrix[2665, ], this_matrix[i, ])
          }
        ) |>
        tibble::tibble(
          cos_sim = _,
          word = word_master$word
        ) |>
        dplyr::arrange(-cos_sim) |>
        dplyr::slice_head(n = 20) |>
        dplyr::mutate(
          record_year = year_master$record_year[i],
          rid = dplyr::row_number() |>
            rev()
        )
    }
  ) |>
  dplyr::bind_rows() |>
  dplyr::filter(
    record_year %in% c(1950, 1955, 1960, 1965, 1970, 1975, 1980, 1985, 1990, 1995, 2000, 2005, 2015)
  ) |>
  ggplot2::ggplot() +
  ggplot2::geom_text(ggplot2::aes(x = as.factor(record_year), y = rid, label = word, alpha = cos_sim), color = "blue", family = "HiraKakuPro-W3") +
  ggplot2::theme_minimal(base_family = "HiraKakuPro-W3") +
  ggplot2::theme(
    panel.grid.major = ggplot2::element_blank(),
    panel.grid.minor = ggplot2::element_blank(),
    axis.title.y = ggplot2::element_blank(),
    axis.text.y = ggplot2::element_blank(),
    axis.ticks.y = ggplot2::element_blank()
  ) + 
  ggplot2::labs(
    title = "安保",
    x = "年"
  )

1950年代から1970年代まで「安保」とコサイン類似度の高い単語に「核軍縮」があったが、最近になって上位の単語から消えました。これは日本外交の言説の重要な変化をモデルが捉えたものだと言えるでしょう。

次に、「軍隊」を確認しましょう：

year_master |>
  nrow() |>
  seq_len() |>
  purrr::map(
    \(i){
      this_matrix <- word_embedding_matrix_weighted[[i]]
      
      this_matrix |>
        nrow() |>
        seq_len() |>
        purrr::map_dbl(
          \(i){
            cos_sim(this_matrix[5580, ], this_matrix[i, ])
          }
        ) |>
        tibble::tibble(
          cos_sim = _,
          word = word_master$word
        ) |>
        dplyr::arrange(-cos_sim) |>
        dplyr::slice_head(n = 20) |>
        dplyr::mutate(
          record_year = year_master$record_year[i],
          rid = dplyr::row_number() |>
            rev()
        )
    }
  ) |>
  dplyr::bind_rows() |>
  dplyr::filter(
    record_year %in% c(1950, 1955, 1960, 1965, 1970, 1975, 1980, 1985, 1990, 1995, 2000, 2005, 2015)
  ) |>
  ggplot2::ggplot() +
  ggplot2::geom_text(ggplot2::aes(x = as.factor(record_year), y = rid, label = word, alpha = cos_sim), color = "blue", family = "HiraKakuPro-W3") +
  ggplot2::theme_minimal(base_family = "HiraKakuPro-W3") +
  ggplot2::theme(
    panel.grid.major = ggplot2::element_blank(),
    panel.grid.minor = ggplot2::element_blank(),
    axis.title.y = ggplot2::element_blank(),
    axis.text.y = ggplot2::element_blank(),
    axis.ticks.y = ggplot2::element_blank()
  ) + 
  ggplot2::labs(
    title = "軍隊",
    x = "年"
  )

最初の数十年はベトナム戦争時に南ベトナムで活動したゲリラ組織「南ベトナム解放民族戦線」の略称である「ベトコン」がトップに来ていたのですが、ベトナム戦争の終結に伴い「軍隊」とのコサイン類似度が下降傾向になり、2000年以降は完全に上位の単語から消えました。

最後に、「日本」を見てみましょう：

year_master |>
  nrow() |>
  seq_len() |>
  purrr::map(
    \(i){
      this_matrix <- word_embedding_matrix_weighted[[i]]
      
      this_matrix |>
        nrow() |>
        seq_len() |>
        purrr::map_dbl(
          \(i){
            cos_sim(this_matrix[3627, ], this_matrix[i, ])
          }
        ) |>
        tibble::tibble(
          cos_sim = _,
          word = word_master$word
        ) |>
        dplyr::arrange(-cos_sim) |>
        dplyr::slice_head(n = 20) |>
        dplyr::mutate(
          record_year = year_master$record_year[i],
          rid = dplyr::row_number() |>
            rev()
        )
    }
  ) |>
  dplyr::bind_rows() |>
  dplyr::filter(
    record_year %in% c(1950, 1955, 1960, 1965, 1970, 1975, 1980, 1985, 1990, 1995, 2000, 2005, 2015)
  ) |>
  ggplot2::ggplot() +
  ggplot2::geom_text(ggplot2::aes(x = as.factor(record_year), y = rid, label = word, alpha = cos_sim), color = "blue", family = "HiraKakuPro-W3") +
  ggplot2::theme_minimal(base_family = "HiraKakuPro-W3") +
  ggplot2::theme(
    panel.grid.major = ggplot2::element_blank(),
    panel.grid.minor = ggplot2::element_blank(),
    axis.title.y = ggplot2::element_blank(),
    axis.text.y = ggplot2::element_blank(),
    axis.ticks.y = ggplot2::element_blank()
  ) + 
  ggplot2::labs(
    title = "日本",
    x = "年"
  )

終戦直後、「軍隊」や（連合国軍の）「駐屯」などがあったのに対し、主権回復に伴いこうした単語が出なくなりました。

このように、本モデルは、単語の意味を時間と共に変わらない単語ベクトルと時間と共に変化する次元重要度ベクトルに分解することで、単語の経時的変化を捉えることに成功しました。

結論

いかがでしたか？このように、棒おり過程などのノンパラメトリックモデルを活用することで、データから柔軟にモデルの構造自体を学習することが可能になります。

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

参考文献

Mikolov, T., Chen, K., Corrado, G., and Dean, J. (2013). Efficient estimation of word representations in vector space. arXiv preprint arXiv:1301.3781.

Rudolph, M., & Blei, D. (2018, April). Dynamic embeddings for language evolution. In Proceedings of the 2018 world wide web conference (pp. 1003-1011).

Rudolph, M., Ruiz, F., Mandt, S., & Blei, D. (2016). Exponential family embeddings. Advances in Neural Information Processing Systems, 29.

Wu, Tung-Wen. (2025). Estimating Dynamic Dimensionality in Roll Call Data: An Application to UNGA Voting. Poster presented at the Japanese Society for Quantitative Political Science (JSQPS) Summer Meeting, Kyoto, July 2025.