【R】数字の集合をトピックモデルで分析したらなかなか良くてびっくりした話。
前回、疎行列をクラスタリングする話を書きました。
そこでふと思いついたのが、「数字の羅列を単語の集合とみなして自然言語処理のスキームに当てはめられるんじゃね？」ということです。どういうことかというと、
[[1]] [1] 246 308 149 170 161 233 291 218 260 171 155 151 [[2]] [1] 383 331 353 [[3]] [1] 310 333 [[4]] [1] 255 247 223 284
こんな感じのデータがあるとします。個々の数字はなんでも良いのですが、例えばこれがある消費者の買った商品IDの集合だとします。通常は、前回書いたとおり、これを主成分分析にかけたり、クラスタリングしたり、あるいはあらかじめ属性がわかっていれば教師あり学習させてみたり、といろいろ思いつくわけなんです。が、じゃあ、この商品IDを「単語」とみなしてトピックモデルで分析したらうまくクラスタリングできるのだろうか？というのが今回のモチベーションです。
そもそも、トピックモデルは自然言語処理のスキームの中に位置づけられるわけなんですが、トピックモデルって大きく考えれば「カテゴリカルデータの分析」なわけで、カテゴリカルデータで良いのなら単語である必要もなく、ならば上記のような消費者の購買データのようなものでも良いよね、となるわけです。
なにはなくともデータセットを作るところから
適切なデータが手元にないので、なんとなくそれっぽいデータを作ります。前回は4つに区切った一様分布を使いました。それだときれいに別れすぎて面白くなかったので、今回は分散が同じで平均値が異なる4つの正規分布から乱数を発生させてデータセットを作ります。イメージとしてはこんな感じです。
これだと重なっているところもあるし、単純だけどちょっと意地悪なデータなので、どんな結果になるか楽しみです。
# 初期設定
maxcat <- 1000
grpnum <- 4
maxdev <- maxcat / grpnum
set.seed(4756)
# 各データ内のデータ数を指数分布に従って発生させる
cnumf <- function(){
set.seed(4756)
round(rexp(maxnum, rate=0.3)+0.5)
}
cnum <- cnumf()
# グループを一様分布で割り当てる
gnumf <- function(){
set.seed(4756)
round(runif(maxnum, min=0.5, max=grpnum + 0.49999))
}
gnum <- gnumf()
# 各データに対して指定個数の乱数を割り当てる
# 割り当てる乱数はN((100*groupNum)+100, 50)の正規分布に従う
library(foreach)
library(Matrix)
# データ作成のための乱数のseedを一様乱数で設定する
seedsVecF <- function(){
set.seed(4756)
round(runif(maxnum)*1000000)
}
seedsVec <- seedsVecF()
cnormf <- function(i){
set.seed(seedsVec[i])
normNum <- rnorm(cnum[i],
mean = (gnum[i]*100)+100,
sd =50)
normNum <- unique(abs(round(normNum)))
return(normNum)
}
cnorm <- foreach(i=1:maxnum) %do% cnormf(i)
# データ確認のためヒストグラムを描く g1 <- unlist(cnorm[gnum==1]) g2 <- unlist(cnorm[gnum==2]) g3 <- unlist(cnorm[gnum==3]) g4 <- unlist(cnorm[gnum==4]) hist(g1, col="#ff634740", breaks=100, xlim=c(0,700)) hist(g2, col="#ffd70040", breaks=100, xlim=c(0,700), add=T) hist(g3, col="#228b2240", breaks=100, xlim=c(0,700), add=T) hist(g4, col="#4169e140", breaks=100, xlim=c(0,700), add=T)
これでいけます。col=で指定したカラーの後ろの2桁分が透明度になっています。重なり具合がそれでわかるようになっています。
従来方法で分析する
いきなりトピックモデルにいく前に、主成分分析したらどうなるか、試してみます。
そのためには前回と同じように各数値の出現頻度を疎行列に変換し、それを主成分分析に突っ込む、というプロセスを経ることになります。
# 重複排除してあるので各データのデータ数を再計算する
clistN <- foreach(i=1:maxnum) %do% rep(i,length(cnorm[[i]]))
# 乱数のベクトル,データ数のデータフレームを作る
cdfN <- data.frame(id = unlist(clistN),
cat = unlist(cnorm),
yn = 1)
# sparseMtrixにするときに引数i,jに0が入らないようにする
catList <- data.frame(cat=unique(cdfN$cat),
catNum=seq(1,length(unique(cdfN$cat))),
stringsAsFactors = F)
cdfN <- merge(cdfN, catList, by="cat", all=T)
# 疎行列に変換する
ctblS <- sparseMatrix(i = cdfN$id,
j = cdfN$catNum,
x = cdfN$yn)これで疎行列に変換できました。その際、疎行列の行名、列名は要素数の通し番号にしないといけないっぽいので注意が必要です。そのため、わざわざ元の要素名と疎行列にしたときの要素名の変換テーブルcatListを作っているのです。
で、頻度行列ができたところで、これを主成分分析にかけます。通常は第二主成分までの二次元プロットで十分なことが多いのですが、今回はちょっとそれだけでは物足りなかったので第三主成分までの三次元プロットにしています。三次元プロットを表示するため、パッケージ"rgl"を使っています。インストールしていない方はあらかじめインストールしておいてください。
# 主成分分析 library(rgl) pcTest <- prcomp(ctblS, scale. =T) plot(pcTest$x[,1:2]) plot3d(pcTest$x[,1:3])
第二主成分までの二次元プロットは次のようになります。
たしかに4つのグループに分けられそうですが、上方向に伸びた腕に対して下方向に伸びている腕が短い感じがします。あまり均等ではないので、第三主成分まで入れたらどうかを見るため、三次元プロットをすると、こんな感じです。
はい、これを二次元に圧縮していたから均等な腕の長さに見えなかったんですね。ちなみにrglで出したグラフはぐりぐり動かせます。今回はスクショなので動きません。ご自身のPCで試してみてください。
ここまでの手順を振り返る
ここまでやってみて、「10,000人×1,000アイテムの購買データ」みたいなやたらと要素の多いデータでも疎行列に変換してしまえば巨大なデータを扱うのがあまり得意ではないRでもそこそこなことができることがわかりました。とはいえ、主成分分析は相関行列を作ることになるので、いずれにせよ巨大な行列を作らなけらばならず、メモリが足りない、という事態に陥ることが考えられそうです。
また、上記したように、sparseMatrix関数だと行名、列名が連番で割り当てられるため、もとのデータが数値型でなければならず、加えて途中が抜けているとうまく行列に変換できなかったりすることがあります。そのあたり、ちょっと手間だなぁ、と考えてしまいます。
で、その手間を埋めてくれるのが自然言語処理のスキームなんじゃないか、と思うわけなのです。
自然言語処理のスキームへ
もう一度、元のデータを眺めます。
[[1]] [1] 246 308 149 170 161 233 291 218 260 171 155 151 [[2]] [1] 383 331 353 [[3]] [1] 310 333 [[4]] [1] 255 247 223 284
最初に作ったデータはリスト型になっていて、各リスト要素にベクトル型で数字の羅列が格納されています。これを半角スペースで結合してしまえば、ひとつのセンテンスとして扱うことができ、そのセンテンスの集合は文書集合とみなすことができる、というわけです。
文書集合への変換
というわけで、元のデータを自然言語処理のスキームに当てはめるべく、データを変換します。
# 文字列ベクトルに変換する
cVec <- sapply(seq(1,length(cnorm)),function(x){
paste(cnorm[[x]], collapse=" ")
})paste関数でリストの各要素をひとつに結合します。その結果、
[1] "246 308 149 170 161 233 291 218 260 171 155 151" [2] "383 331 353" [3] "310 333" [4] "255 247 223 284"
このような数字の羅列の集合を得ることができました。ここから先、これらの数字はもう数字ではありません。すべて「単語」です。一つ一つのベクトル要素は「センテンス」です。ということは、これを「形態素解析」して「Document Term Matrix」に変換できる、ということなのです。
Document Term Matrixに変換する
パッケージはお気に入りのtext2vecを使います。
Word Embeddingのためのパッケージではありますが、自然言語処理のための基本的な処理もできますので、いろいろ便利です。
library(text2vec) # document term matrixを作る tokenizer <- itoken(cVec, ids=seq(1,length(cVec))) vocab <- create_vocabulary(tokenizer) vectorizer <- vocab_vectorizer(vocab) dtm <- create_dtm(tokenizer, vectorizer)
今回はstop wordsは存在しませんし、punctuationも存在しません。複合語を想定する必要もないのでn_gram_windowも指定しません。生のまま、そのままじゃんじゃん処理していきます。そうしてできたdtmをみると
> dtm
10000 x 609 sparse Matrix of class "dgCMatrix"
[[ suppressing 34 column names ‘652’, ‘632’, ‘655’ ... ]]
[[ suppressing 34 column names ‘652’, ‘632’, ‘655’ ... ]]
1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......
6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ......つーかこれ、最初になんとか作ろうとしてた疎行列なんじゃね？
ああ、とうとう気付いてはいけないことに気づいてしまいましたね。。。
そうなんです、これ、もともと作ろうとしていた疎行列そのものなんです。
> colnames(dtm) [1] "652" "632" "655" "81" "47" "57" "651" "52" "625" "71" "54" "62" "627" [14] "643" "53" "608" "85" "69" "620" "34" "648" "603" "48" "644" "593" "98" [27] "635" "86" "64" "72" "636" "78" "91" "629" "127" "614" "622" "82" "50"
列名をみるとソート順こそ連番での並びにはなっていませんが、各値になっていますし、行列の各要素はその頻度になっているのです。データを加工してまでわざわざsparsMatrix関数を使う必要がなかったのです。
っても今更な話なんで、過去は過去として、潔く次に進みましょう。
いよいよトピックモデルを当てはめる
はい、余計な過去は忘れていよいよトピックモデルです。どんな感じになるかわくわくです。
ちなみに、使うトピックモデルはLDA、結果の可視化のために"LDAvis"というパッケージをつかっています。
ブラウザで表示させてぐりぐり動くので、なかなか楽しいです。
※今回はスクショだけですので動きません。ご自身のPCで試してみてください。
なお、ここではLDAの詳しい説明はしません。そもそもの原理に興味のある方は下記サイトをご参考ください。
では、実際のコードです。
library(LDAvis) # topic model LDA ldaModel <- LatentDirichletAllocation$new(n_topics=4, vocabulary=vocab) ldaRes <- ldaModel$fit(dtm, n_iter=1000, check_convergence_every_n = 1) ldaRes$plot()
まずはトピック数4でやってみます。その結果、
きれいに４つにわかれました。そのうち、topic1は
200近辺の値が多く出現すると見積もられていますから、どうやら想定した通りの別れ方になっているようです。そこで、各トピックにどんな値が出現するのか、チェックしてみましょう。そのためには$get_word_vectorを使います。
# 各トピックにおける単語の出現頻度を集計する
wordList <- ldaRes$get_word_vectors()
order <- order(as.numeric(rownames(wordList)))
wordList <- wordList[order,]
# プロット
plot("number","frequency",type="n",
xlim=c(0,max(as.numeric(rownames(wordList)))),
ylim=c(0,max(wordList)))
points(rownames(wordList),wordList[,1],col="#ff6347",pch=16)
points(rownames(wordList),wordList[,2],col="#ffd700",pch=16)
points(rownames(wordList),wordList[,3],col="#228b22",pch=16)
points(rownames(wordList),wordList[,4],col="#4169e1",pch=16)で、プロットした結果がこちら。
ほらほら、これ、すごいと思いません？ちゃんと正規分布が再現されているんですよ？しかも、結構重なりが大きいにも関わらず、重なっている部分もちゃんと分類できているっぽい。これって下手にハードクラスタリングするよりずっといいのじゃないだろうか？と考えてしまいます。
意地悪してトピック数を変えてみる
となると、LDAにおけるトピック数を変えたときにどうなるのか、というところに興味が行きます。ならばやってみよう。
トピック数=2のケース
############ topic=2で実験 # topic model LDA ldaModel <- LatentDirichletAllocation$new(n_topics=2, vocabulary=vocab) ldaRes_k2 <- ldaModel$fit(dtm, n_iter=1000, check_convergence_every_n = 1) #ldaRes_k2$plot() wordList <- ldaRes_k2$get_word_vectors() order <- order(as.numeric(rownames(wordList))) wordList <- wordList[order,]
LDAvisはトピック数が3以上でないと動かなかったため、コメントアウトしています。なので、各トピックにおける単語の出現頻度の分布だけ。
おお、きれいに二つの山になっている。左側2つの正規分布と右側2つの正規分布をそれぞれ併合した形になっているわけですな。それにしてもきれいに出るもんだなぁ。
トピック数=3のケース
############ topic=3で実験
# topic model LDA
ldaModel <- LatentDirichletAllocation$new(n_topics=3, vocabulary=vocab)
ldaRes_k3 <- ldaModel$fit(dtm, n_iter=1000, check_convergence_every_n = 1)
ldaRes_k3$plot()
wordList <- ldaRes_k3$get_word_vectors()
order <- order(as.numeric(rownames(wordList)))
wordList <- wordList[order,]
# プロット
plot("number","frequency",type="n",
xlim=c(0,max(as.numeric(rownames(wordList)))),
ylim=c(0,max(wordList)))
points(rownames(wordList),wordList[,1],col="#ff6347",pch=16)
points(rownames(wordList),wordList[,2],col="#ffd700",pch=16)
points(rownames(wordList),wordList[,3],col="#228b22",pch=16)
はい、きちんと３つに分かれました。感じ的にはどこかの二つの正規分布を無理やりひとつにまとめた、という感じでしょうか。分布で確認してみましょう。
うん、そうですね。真ん中のふたつの正規分布を一つに併合した感じ。
トピック数=7のケース
最後にトピック数を7にした場合。なぜ7にしたかというと、
もとのデータが4つの正規分布で、重なっている箇所が3箇所あるから。つまり、この重なっている部分がトピックとして独立したりしないのかなー、と思ったためです。
############ topic=7で実験
# topic model LDA
ldaModel <- LatentDirichletAllocation$new(n_topics=7, vocabulary=vocab)
ldaRes_k7 <- ldaModel$fit(dtm, n_iter=1000, check_convergence_every_n = 1)
ldaRes_k7$plot()
wordList <- ldaRes_k7$get_word_vectors()
order <- order(as.numeric(rownames(wordList)))
wordList <- wordList[order,]
# プロット
plot("number","frequency",type="n",
xlim=c(0,max(as.numeric(rownames(wordList)))),
ylim=c(0,max(wordList)))
points(rownames(wordList),wordList[,1],col="#ff6347",pch=16)
points(rownames(wordList),wordList[,2],col="#ffd700",pch=16)
points(rownames(wordList),wordList[,3],col="#228b22",pch=16)
points(rownames(wordList),wordList[,4],col="#4169e1",pch=16)
points(rownames(wordList),wordList[,5],col="#c71585",pch=16)
points(rownames(wordList),wordList[,6],col="#adff2f",pch=16)
points(rownames(wordList),wordList[,7],col="#20b2aa",pch=16)はい、結果はというと、
こんな感じで左端の正規分布のみ独立していて、それ以外は無理やり2つに分けたような感じ。
改めて考えてみると、トピックモデルって要は「ある単語ペアがあるセンテンス中に同時に現れる確率」も計算しているわけで、いくら分布が重なっているからといっても、その重なり部分が独立したトピックにはなりづらいのだなぁ、と。だから、こんな風にトピックにおける単語の出現頻度の分布を眺めてみると、適切なグルーピングになっているのだな、と。
いやぁ、なんか妙にトピックモデルの勉強になった気がする。
