【学習動機】
以前、テキストマイニングをRで使って遊んでみたとき, トピックモデルが出てきたが、理論を全く知らないのも悲しいので、岩田『トピックモデル』
で学ぶことにした。まずはユニグラムモデルについて学ぶ。
【学習内容】
ユニグラムモデル(unigram model)はBOW(bag of words)表現された文書集合の最も簡単な生成モデル。このモデルの理論を学んだので、Rで人工データを作成し、最尤推定とMAP推定でパラメータの推定を行う。
ライブラリの読み込みと表示桁数の指定。
library(tidyverse)
library(ggplot2)
options(digits=3)
set.seed(536329)単語分布の設定をする。語彙の種類とそれに付与される確率を決定。
df_word_dist <- data.frame(
words = c("あじ","いくら","うに","えんがわ","かれい","こはだ","さわら","さんま","しゃこ","すずき"),
prob = c(0.2,0.1,0.075,0.075,0.05,0.075,0.1,0.2,0.075,0.05)
)
numV <- nrow(df_word_dist) #単語分布のボキャブラリーの種類(number of Vocabulary)
allword <- df_word_dist$words %>% as.vector() #全単語(あとでファクタ型にしたときの並べ替えで使う)文書数と各文書の単語数を指定する。
numD <- 12 #文書の数(number of Documents)
vnumW <- c(20,70,10,80,20,80,10,80,20,80,40,60) #各文書に含まれる単語数(number of words)のベクトル
numW <- sum(vnumW) #単語の総数シミュレーションデータ作成関数が以下のようになる。引数には文書の数と各文書の単語数ベクトル、単語分布のデータフレームをとる。
sim_unigram <- function(numD, vnumW, df_word_dist){
list_word <- vector("list", numD)
for(d in 1:numD){
list_word[[d]] <- sample(df_word_dist$words, vnumW[d], replace=TRUE, prob = df_word_dist$prob)
}
return(list_word)
}これにより人工データを作成し、生成された文書1の単語を見てみる。
DW <- sim_unigram(numD=numD, vnumW = vnumW, df_word_dist = df_word_dist)
DW[[1]]## [1] あじ こはだ さんま すずき えんがわ さんま うに
## [8] さわら さんま かれい すずき しゃこ さんま あじ
## [15] さわら あじ さんま すずき あじ さんま
## 10 Levels: あじ いくら うに えんがわ かれい こはだ さわら ... すずき各語彙の頻度テーブルを作成し、生成された文書1のテーブルを見てみる。
DWtab <- DW %>% purrr::map(factor,levels = allword) %>% purrr::map(table)
DWtab[[1]]##
## あじ いくら うに えんがわ かれい こはだ さわら さんま
## 4 0 1 1 1 1 2 6
## しゃこ すずき
## 1 3ユニグラムモデルでは文書によらず単語分布が一緒なので、語彙の頻度テーブルを統一してしまって問題ない。
allWtab <- DW %>% unlist() %>% factor(levels = allword) %>% table()最尤推定(phiv = Nv/N)を行う。
MLphi <- allWtab / numW
MLphi## .
## あじ いくら うに えんがわ かれい こはだ さわら さんま
## 0.2035 0.1035 0.0842 0.0930 0.0526 0.0579 0.0860 0.1930
## しゃこ すずき
## 0.0807 0.0456事前分布にディリクレ分布を設定し、パラメータを(beta,beta,beta,…,beta)として、MAP推定(phiv = (Nv+beta-1)/(N+(beta-1)V) beta:paramofDirichlet)を行う。
beta <- 2 #観測がない単語の推定値を負にしないためにはbetaが1以上とする必要がある
MAPphi <- (allWtab + beta - 1) / (numW + (beta - 1)*numV)
MAPphi## .
## あじ いくら うに えんがわ かれい こはだ さわら さんま
## 0.2017 0.1034 0.0845 0.0931 0.0534 0.0586 0.0862 0.1914
## しゃこ すずき
## 0.0810 0.0466ちなみに、ベイズ推定では事前分布にディリクレ分布として、パラメータを(beta,beta,beta,…,beta)とすると、事後分布はDirichlet(phi|beta+N1,beta+N2,…,beta+Nv)が得られるが、ここでは求めない。
MAP推定と最尤推定(ML)と真のパラメータを棒グラフにして並べる。赤がMAP推定値、緑が最尤推定値、青が真の値を示している。
df_est <- data.frame(
voc = df_word_dist$words,
phi = df_word_dist$prob,
MLphi = MLphi %>% as.vector(),
MAPphi = MAPphi %>% as.vector()
)
df_est_wide <- gather(df_est, key = "estimator",value,-voc)
df_est_wide$voc <- df_est_wide$voc %>% factor(levels = df_word_dist$words)
ggplot(df_est_wide, aes(x=estimator, y=value, fill=estimator)) +
geom_bar(stat = "identity", position = "dodge") +
facet_wrap(~voc)
サンプルサイズを10%にすると当然ながら推定精度が落ちる。
set.seed(536329)
numD <- 12 #文書の数(number of Documents)
vnumW <- c(20,70,10,80,20,80,10,80,20,80,40,60)*0.1 #各文書に含まれる単語数(number of words)のベクトル
