明けましておめでとうございます

遅れ馳せながら、明けましておめでとうございます :)
三ヶ日も終わり、正月気分なのは僕だけでしょうか？

昨年を振り返ると、研究会、Workshopや勉強会に複数参加をする等、アクティブに動けた１年だったと思います。
今年最初の仕事は、そういった研究成果をしっかりとまとめ、業績とすること！

残り少ない学生生活…締切に追われる毎日ですが、大いに楽しみたいと思います！

追伸
卒業後は研究からは離れますが、興味のある論文やその実装、技術に関してメモをするこのblogは、続けたいなと思っています
　

Probabilistic Latent Semantic Analysis : PLSA （Rで実装）

前回のエントリからはや一ヶ月。月日が立つのは早いものです。
修論に向け、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アルゴリズムで、温度スケジューリングは簡単のため

1. 温度T=1.0からスタート
2. 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の値によってかなり推定精度にバラつきがありますが、シンプルな理論/実装で動くモデルなので良いですね。例によって、間違い、またより良い実装方法がございましたらご指摘ください。