修論に向け、bag-of-featuresの実装をもくろんでおりますが、その一環としてPLSAを試してみました。
参考文献はこちら（リンク先pdf）。
T. Hofmann. Probabilistic latent semantic analysis. In Proceedings of the 15th Conference on Uncertainty in AI, 1999.
ちょうど10年前に提案されたモデルですが、LDAの元となったり、現在でも多くの論文が発表されたりと、良い言語モデルのようです。
これをRで素直に実装したのがこちら。
plsi <- function(x, K=10, eps=0.9, max_itr=200,...){
#logsumexp
logsumexp<-function(x,y,flg){
if(flg){
return(y)
}
if(x == y){
return(x + 0.69314718055)
}
vmin = min(x,y)
vmax = max(x,y)
if(vmax > (vmin + 50.0)){
return(vmax)
}else{
return(vmax + log(exp(vmin - vmax) + 1.0))
}
}
#random variable from Dirichlet distribution
rdirichlet<-function(alpha,prec=0,...){
if(prec != 0){
x<-rgamma(length(alpha),alpha*prec,1.0)
}else{
x<-rgamma(length(alpha),alpha,1.0)
}
return (x/sum(x))
}
#initialize
D <- dim(x)[1]
W <- dim(x)[2]
T <- 1 #temperature
pllik <- 0
#parameters
pz_dw <- list()
for(k in 1:K){
pz_dw[[k]] <- matrix(0, D, W)
}
pz <- rep(1/K, length=K)
pd_z <- matrix(0, K, D)
pw_z <- matrix(0, K, W)
unigram <- 1/W * apply(x,2,sum)
for(k in 1:K){
pd_z[k,] <- rdirichlet(rep(1/D, length=D))
pw_z[k,] <- rdirichlet(unigram)
}
##############################
for(t in 1:max_itr){
cat("Iteration : ",t," Temperature : ",T," ")
#E-step
#initialize expected count
cz <- rep(0, length=K)
cd_z <- matrix(0, K, D)
cw_z <- matrix(0, K, W)
for(d in 1:D){
for(w in 1:W){
z <- 0
for(k in 1:K){
pz_dw[[k]][d,w] <- T * ( log(pz[k]) + log(pd_z[k,d]) + log(pw_z[k,w]) )
z <- logsumexp(z, pz_dw[[k]][d,w], (k==1))
}
for(k in 1:K){
pz_dw[[k]][d,w] <- exp(pz_dw[[k]][d,w] - z)
if(x[d,w] != 0){
cz[k] <- cz[k] + x[d,w] * pz_dw[[k]][d,w]
cd_z[k,d] <- cd_z[k,d] + x[d,w] * pz_dw[[k]][d,w]
cw_z[k,w] <- cw_z[k,w] + x[d,w] * pz_dw[[k]][d,w]
}
}
}
}
#M-step
pz <- cz / sum(cz)
for(k in 1:K){
pd_z[k, ] <- cd_z[k, ] / sum(cd_z[k, ])
pw_z[k, ] <- cw_z[k, ] / sum(cw_z[k, ])
}
#converged?
llik <- 0
for(d in 1:D){
for(w in 1:W){
if(x[d,w] != 0){
l <- 0
pdw <- 0
for(k in 1:K){
l <- log(pz[k]) + log(pd_z[k,d]) + log(pw_z[k,w])
pdw <- logsumexp(pdw, l, (k==1))
}
llik <- llik + x[d,w] * pdw
}
}
}
cat(" log-likelihood : ",llik, " PPL : ",- 1/W * llik,"\n")
if(t > 1){
if( abs((llik - pllik) / llik) < 1e-5 || pllik > llik){
cat("Converged.\n")
break
}
}
pllik <- llik
T <- eps * T
}
return(list("pz_dw"=pz_dw, "pz"=pz, "pd_z"=pd_z, "pw_z"=pw_z))
}
素直に参考文献の式(3)~(6)を実装したものになります。推定アルゴリズムはTemperature EMアルゴリズムで、温度スケジューリングは簡単のため
- 温度T=1.0からスタート
- EMアルゴリズムの各反復において, T = eps × T (eps < 1.0)で更新
用いたデータは以下のような12×9の行列データです。
1 2 2 0 0 0 0 0 0 2 1 3 0 0 0 0 0 0 2 0 2 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 1 1 2 0 0 0 0 0 0 2 1 1 0 0 0 0 0 0 1 3 1 0 0 0 0 0 0 1 2 2 0 0 0 0 0 0 0 0 0 2 3 3 0 0 0 0 0 0 3 2 2 0 0 0 0 0 0 2 1 3 0 0 0 0 0 0 1 1 4これがテキストデータであるとすると、データ数が12、単語数が9、行列要素はデータ内で観測された頻度に相当します。
目視で確認できるように、左側・中央・右側に出現頻度が偏った、いわば３クラスが含まれるデータになります。
これを、以下のようにしてplsiを実行します。
> opt <- plsi(x, K=3, eps=0.95)オプションKはクラス数、epsはTEMの温度に関するパラメータです。これにより、以下のような出力がされます。
Iteration : 1 Temperature : 1 log-likelihood : -277.8714 Iteration : 2 Temperature : 0.9 log-likelihood : -275.7771 Iteration : 3 Temperature : 0.81 log-likelihood : -274.7316 Iteration : 4 Temperature : 0.729 log-likelihood : -272.8197 Iteration : 5 Temperature : 0.6561 log-likelihood : -268.9356 Iteration : 6 Temperature : 0.59049 log-likelihood : -265.5543 Iteration : 7 Temperature : 0.531441 log-likelihood : -264.1783 Iteration : 8 Temperature : 0.4782969 log-likelihood : -264.3499 Converged.最後の反復で若干対数尤度が悪化しているのは玉に傷です。
各クラスからの文章の出現確率、すなわちp(d|z)は、以下によって取得できます（見やすさのため小数点以下二桁で丸めています）。
> round(opt$pd_z,2)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.3 0.26 0.22 0.22
[2,] 0.28 0.33 0.22 0.17 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00
[3,] 0.00 0.01 0.00 0.00 0.22 0.22 0.27 0.27 0.0 0.00 0.00 0.00
行がクラス、列がデータを表しています。概ね正しく３つのクラスに分類できていることが伺えます。同様にして、各クラスからの単語の出現確率p(w|z)は、以下によって取得できます。
> round(opt$pw_z,2)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.00 0.00 0.00 0.00 0.00 0.00 0.3 0.26 0.44
[2,] 0.33 0.22 0.44 0.00 0.00 0.00 0.0 0.00 0.00
[3,] 0.01 0.00 0.01 0.27 0.38 0.33 0.0 0.00 0.00
行がクラス、列が単語を表しています。こちらも概ね正しく各クラスからの単語の出現確率を推定できていることが伺えます。初期値やK/epsの値によってかなり推定精度にバラつきがありますが、シンプルな理論/実装で動くモデルなので良いですね。例によって、間違い、またより良い実装方法がございましたらご指摘ください。
から、n個の確率
が生成されたとします。
のベータ分布
は以下のようになります。
であることを表しています。概ねよさげですが、ちょっと精度が悪い？
