はじめに

　『トピックモデル』(MLPシリーズ)の勉強会資料のまとめです。各種モデルやアルゴリズムを「数式」と「プログラム」を用いて解説します。
　本の補助として読んでください。

　この記事では、トピックモデルの生成モデルをR言語でスクラッチ実装します。

【前節の内容】

【他の節の内容】

【この節の内容】

4.1 トピックモデルの生成モデルの実装

　トピックモデル(topic models)の生成モデル(generative model)・生成過程(generative process)を実装して、簡易的な文書データ(トイデータ)を作成する。
　トピックモデルについては「4.1：トピックモデルの生成モデルの導出【青トピックモデルのノート】 - からっぽのしょこ」を参照のこと。

　利用するパッケージを読み込む。

# 利用パッケージ
library(tidyverse)
library(MCMCpack)

　この記事では、基本的に パッケージ名::関数名() の記法を使うので、パッケージの読み込みは不要である。ただし、作図コードについてはパッケージ名を省略するので、ggplot2 を読み込む必要がある。
　また、ネイティブパイプ演算子 |> を使う。magrittrパッケージのパイプ演算子 %>% に置き換えられるが、その場合は magrittr を読み込む必要がある。

パラメータの設定

　文書集合に関する値を設定する。

# 文書数を指定
D <- 10

# 語彙数を指定
V <- 6

# トピック数を指定
true_K <- 4

　文書数 $D$ 、語彙(単語の種類)数 $V$ 、トピック数 $K$ を指定する。

　トピック分布のハイパーパラメータを設定する。

# トピック分布のハイパーパラメータを指定
true_alpha   <- 1
true_alpha_k <- rep(true_alpha, times = true_K) # 一様な値の場合
#true_alpha_k <- runif(n = true_K, min = 1, max = 2) # 多様な値の場合
true_alpha_k

　 $K$ 個の値(ベクトル) $\boldsymbol{\alpha} = (\alpha_1, \cdots, \alpha_K)$ を指定する。各値は0より大きい値 $\alpha_k \gt 0$ を満たす必要がある。全て同じ値 $\alpha_1 = \cdots = \alpha_K = \alpha$ とする場合は、1個の値(スカラ) $\alpha$ で表す。

　単語分布のパラメータを設定する。

# トピック分布のパラメータを生成
#true_theta_dk <- matrix(1/K, nrow = D, ncol = true_K) # 一様な値の場合
true_theta_dk <- MCMCpack::rdirichlet(n = D, alpha = true_alpha_k) # 多様な値の場合
true_theta_dk; dim(true_theta_dk); rowSums(true_theta_dk)

　 $D$ 個のトピック分布のパラメータ $\boldsymbol{\Theta} = \{\boldsymbol{\theta}_1, \cdots, \boldsymbol{\theta}_D\}$ を作成する。パラメータは $K$ 個の値(ベクトル) $\boldsymbol{\theta}_d = (\theta_{d1}, \cdots, \theta_{dK})$ であり、各値は0から1の値 $0 \leq \theta_{dk} \leq 1$ で、総和が1になる値 $\sum_{k=1}^K \theta_{dk} = 1$ を満たす必要がある。
　ディリクレ分布の乱数は、rdirichlet() で生成できる。サンプルサイズの引数 n に文書数 D、パラメータの引数 alpha にハイパーパラメータ $\boldsymbol{\alpha}$ を指定する。

　単語分布のハイパーパラメータを設定する。

# 単語分布のハイパーパラメータを指定
true_beta   <- 1
true_beta_v <- rep(true_beta, times = V) # 一様な値の場合
#true_beta_v <- runif(n = V, min = 1, max = 2) # 多様な値の場合
true_beta_v

　 $V$ 個の値(ベクトル) $\boldsymbol{\beta} = (\beta_1, \cdots, \beta_V)$ を指定する。各値は0より大きい値 $\beta_v \gt 0$ を満たす必要がある。全て同じ値 $\beta_1 = \cdots = \beta_V = \beta$ とする場合は、1個の値(スカラ) $\beta$ で表す。

　単語分布のパラメータを設定する。

# 単語分布のパラメータを生成
#true_phi_kv <- matrix(1/V, nrow = true_K, ncol = V) # 一様な値の場合
true_phi_kv <- MCMCpack::rdirichlet(n = true_K, alpha = true_beta_v) # 多様な値の場合
true_phi_kv; dim(true_phi_kv); rowSums(true_phi_kv)

　 $K$ 個の単語分布(語彙分布)のパラメータ $\boldsymbol{\Phi} = \{\boldsymbol{\phi}_1, \cdots, \boldsymbol{\phi}_K\}$ を作成する。パラメータは $V$ 個の値(ベクトル) $\boldsymbol{\phi}_k = (\phi_{k1}, \cdots, \phi_{kV})$ であり、各値は0から1の値 $0 \leq \phi_{kv} \leq 1$ で、総和が1になる値 $\sum_{v=1}^V \phi_{kv} = 1$ を満たす必要がある。
　rdirichlet() の n 引数にトピック数 $K$ 、alpha 引数にハイパーパラメータ $\boldsymbol{\beta}$ を指定する。
　true_K 行 V 列のマトリクスとして扱う。

頻度データの生成

　トピックモデルでは単語が出現した順番(語順)は考慮しないので、語順情報の有無によって2通り方法で頻度情報を作成する。どちらの方法でも同じ形式のオブジェクトを作成できる。

語順を記録する場合

　語彙の出現順と共に語彙ごとの出現頻度を作成する。こちらの方法では、混合ユニグラムモデルの生成過程をより再現するように処理している。

　設定したパラメータに従う文書データを作成する。

# 文書集合を初期化
w_lt <- list()

# トピック集合を初期化
true_z_lt <- list()

# 各文書の単語数を初期化
N_d <- rep(NA, times = D)

# 文書ごとの各語彙の単語数を初期化
N_dv <- matrix(0, nrow = D, ncol = V)

# 文書データを生成
for(d in 1:D) { # 文書ごと
  
  # 単語数を生成
  N_d[d] <- sample(x = 10:20, size = 1) # 下限:上限を指定
  
  # 各単語のトピックを初期化
  true_z_n <- rep(NA, times = N_d[d]) # トピック集合
  
  # 各単語の語彙を初期化
  w_n <- rep(NA, times = N_d[d])
  
  for(n in 1:N_d[d]) { # 単語ごと
    
    # トピックを生成
    onehot_k <- rmultinom(n = 1, size = 1, prob = true_theta_dk[d, ]) |> # (カテゴリ乱数)
      as.vector() # one-hot符号
    k <- which(onehot_k == 1) # トピック番号
    
    # トピックを割当
    true_z_n[n] <- k
    
    # 語彙を生成
    onehot_v <- rmultinom(n = 1, size = 1, prob = true_phi_kv[k, ]) |> # (カテゴリ乱数)
      as.vector() # one-hot符号
    v <- which(onehot_v == 1) # 語彙番号
    
    # 語彙を割当
    w_n[n] <- v # 単語集合
    
    # 頻度をカウント
    N_dv[d, v] <- N_dv[d, v] + 1
  }
  
  # トピック集合を格納
  true_z_lt[[d]] <- true_z_n
  
  # 単語集合を格納
  w_lt[[d]] <- w_n # 文書集合
  
  # 途中経過を表示
  print(
    paste0(
      "document: ", d, ", ", 
      "words: ", N_d[d], ", ", 
      "topics: ", paste0(sort(unique(true_z_n)), collapse = ", ")
    )
  )
}

個の文書について、順番に単語集合を生成していく。
- 文書 $d$ の単語数 $N_d$ をランダムに決める。この例では、下限数・上限数を指定して一様乱数を生成する。
- 個の単語について、順番に語彙を生成していく。
  - 文書 $d$ のトピック分布(カテゴリ分布) $p(k | \boldsymbol{\theta}_d)$ に従い、単語のトピック $z_{dn}$ に $K$ 個のトピック $k = 1, \dots, K$ からトピック番号を割り当てる。
  - 割り当てられたトピック $k$ の単語分布(カテゴリ分布) $p(v | \boldsymbol{\phi}_k)$ に従い、単語 $w_{dn}$ に $V$ 個の語彙 $v = 1, \dots, V$ から語彙番号を割り当てる。
  - 生成した語彙 $v$ に応じて出現回数 $N_{dv}$ をカウントする。

　文書間で単語数が異なる(マトリクスとして扱えない)ので、文書集合(単語集合の集合) $\mathbf{W} = \{\mathbf{w}_1, \cdots, \mathbf{w}_D\}$ をリストとして、単語集合 $\mathbf{w}_d = \{w_{d1}, \cdots, w_{dN_d}\}$ のベクトルを格納する。

　カテゴリ分布の乱数は、rmultinom() で生成できる。サンプルサイズの引数 n に 1、試行回数の引数 size に 1、パラメータの引数 prob に各文書のトピック分布のパラメータ $\boldsymbol{\theta}_d$ または各トピックの単語分布の $\boldsymbol{\phi}_k$ を指定する。
　one-hotベクトルが生成されるので、which() で値が 1 の要素番号を得る。

　作成したデータを確認する。

# 観測データを確認
head(w_lt)
head(N_dv)

　w_lt[[d]] が $\mathbf{w}_d$ 、w_lt[[d]][n] が $w_{dn}$ 、N_dv[d, v] が $N_{dv}$ に対応する。

語順を記録しない場合

　語順情報(文書集合)を作成せずに語彙ごとの出現頻度を作成する。こちらの方法では、より簡易的に処理している。

　設定したパラメータに従う文書データを作成する。

# トピック集合を初期化
true_z_lt <- list()

# 各文書の単語数を生成
N_d <- sample(x = 10:20, size = D, replace = TRUE) # 下限:上限を指定

# 文書ごとの各語彙の単語数を初期化
N_dv <- matrix(0, nrow = D, ncol = V)

# 文書データを生成
for(d in 1:D) { # 文書ごと
  
  # 各単語のトピックを生成
  true_z_n <- sample(x = 1:true_K, size = N_d[d], replace = TRUE, prob = true_theta_dk[d, ]) # (カテゴリ乱数)
  
  for(n in 1:N_d[d]) { # 単語ごと
    
    # 語彙を生成
    v <- sample(x = 1:V, size = 1, prob = true_phi_kv[true_z_n[n], ]) # (カテゴリ乱数)
    
    # 頻度をカウント
    N_dv[d, v] <- N_dv[d, v] + 1
  }
  
  # トピック集合を格納
  true_z_lt[[d]] <- true_z_n
  
  # 途中経過を表示
  print(
    paste0(
      "document: ", d, ", ", 
      "words: ", N_d[d], ", ", 
      "topics: ", paste0(sort(unique(true_z_n)), collapse = ", ")
    )
  )
}

各文書の単語数 $(N_1, \cdots, N_D)$ をランダムに決める。この例では、下限数・上限数を指定して一様乱数を生成する。
個の文書について、順番に個分の語彙頻度を生成していく。
- 文書 $d$ のトピック分布(カテゴリ分布) $p(k | \boldsymbol{\theta}_d)$ に従い、各単語のトピック $\mathbf{z}_d$ を生成する。
- 割り当てられたトピック $z_d$ の単語分布(カテゴリ分布) $p(v | \boldsymbol{\phi}_{z_d})$ に従い、 $V$ 個の語彙 $v = 1, \dots, V$ ごとの出現回数 $(N_{d1}, \cdots, N_{dV})$ を生成する。

　カテゴリ分布の乱数のインデックスは、sample() で生成できる。
　トピック番号を生成する場合は、確率変数の引数 x にトピック番号 $\{1, \dots, K\}$ 、試行回数の引数 size に単語数 $N_d$ 、生成確率の引数 prob に各文書のトピック分布のパラメータ $\boldsymbol{\theta}_d$ を指定する。
　語彙番号を生成する場合は、x 引数に語彙番号 $\{1, \dots, V\}$ 、size 引数に 1、prob 引数に各トピックの単語分布のパラメータ $\boldsymbol{\phi}_k$ を指定する。

　生成したデータを確認する。

# 観測データを確認
head(N_dv)

　作図などで語順を使う場合は、ランダムに語順を決める。

# 文書集合を初期化
w_lt <- list()

# 文書ごとに単語集合を作成
for(d in 1:D) {
  
  # 語彙番号を作成
  v_vec <- mapply(
    FUN = \(v, n) {rep(v, times = n)}, # 語彙番号を複製
    1:V,      # 語彙番号
    N_dv[d, ] # 出現回数
  ) |> 
    unlist()
  
  # 語彙を割当
  w_n <- sample(x = v_vec, size = N_d[d], replace = FALSE) # 単語集合
  
  # 単語集合を格納
  w_lt[[d]] <- w_n
}
head(w_lt)

　mapply() を使って 1 から V の語彙番号を出現回数 N_dv[d, ] に応じてそれぞれ複製する。リストが出力されるので unlist() でベクトルに変換する。
　N_d[d] 個の語彙番号を sample() でランダムに入れ替えて単語集合のベクトルとする。
　文書ごとに単語集合のベクトルを作成して、文書集合のリストに格納していく。

　各種推論アルゴリズムでは、文書集合 $\mathbf{W}$ から得られる語彙頻度データ $(N_{11}, \cdots, N_{DV})$ のみを用いてパラメータ類を推定する。

データの可視化

　次は、生成した真の分布や文書データのグラフを作成する。

データの集計

　文書や語彙に関する単語数を集計する。

# 各文書の単語数を取得
N_d <- rowSums(N_dv)
N_d

　各文書の単語数 $(N_1, \cdots, N_D)$ 、 $N_d = \sum_{v=1}^V N_{dv}$ は N_dv の行ごとの和に対応する。

　文書集合に関する値を集計する。

# 文書数を取得
D <- nrow(N_dv)
D <- length(N_d)

# 語彙数を取得
V <- ncol(N_dv)

# 全文書の単語数を取得
N <- sum(N_dv)
N <- sum(N_d)
D; V; N

　文書数 $D$ は N_dv の行数または N_d の要素数、語彙数 $V$ は N_dv の列数に対応する。
　全文書の単語数(総単語数) $N = \sum_{d=1}^D \sum_{v=1}^V N_{dv}$ は N_dv, N_d それぞれの全ての要素の和に対応する。

真の分布の図

　真の分布(パラメータ)をグラフで確認する。

続きは追々書きます…

文書データの図

　文書集合(頻度データ)とトピック集合をグラフで確認する。

　この記事では、ユニグラムモデルの生成モデルを実装して、各種推論アルゴリズムに用いる人工データを作成した。次の記事からは、3つの推論アルゴリズムを確認していく。

参考書籍

作者:岩田具治
講談社

おわりに

　またブログが自転車操業になってる…

　2024年5月29日は、えびちゅうこと私立恵比寿中学の仲村悠菜さんの17歳のお誕生日です！

　第一印象はお人形みたいに可愛いなんですが、ちょっと見てると面白い人だの方が印象が強くなっていきます。そして歌うと格好良いとなります。
　スタプラってそういう人が各グループにいませんか？

【次節の内容】

スクラッチ実装編

　トピックモデルに対する最尤推定をプログラムで確認します。

www.anarchive-beta.com

数式読解編

　トピックモデルに対する最尤推定を数式で確認します。