RでPLSA(PLSI)

Probabilistic Latent Semantic AnalysisとProbabilistic Latent Semantic Indexingのどっちの名前を使えばいいのかわかりませんが，昔PLSAの実装に挑戦したときのソースが出てきたので晒しとく．
ちゃんと動いてるかどうかは不明．Tempered EMアルゴリズムではなくて普通のEMアルゴリズム．
Rのforループは遅いので，行列演算とapplyを駆使してどれだけforループを減らせるかという目的で書いてたらソースコードがカオスになって諦めた．
行列演算とapplyに頭を使うぐらいなら別の言語でforループで書いたほうが……．

#data
x <- matrix(c(1,1,1,0,0,0, 1,0,1,0,0,0, 1,1,0,0,0,0,0, 0,0,1,1,1,0, 0,0,1,0,1,0,0,0,1,1,0), 6)

#Topic Num
Z<-2
#Max iteration
MI<-50

#initialize
	D <- dim(x)[1]
	W <- dim(x)[2]

#parameters
	pz_dw <- array(0, dim=c(Z, D, W))
	pd_z <- matrix(runif(Z*D), D, Z)
	pw_z <- matrix(runif(Z*W), W, Z)
	pd_z <- t(t(pd_z) / apply(pd_z, 2, sum))
	pw_z <- t(t(pw_z) / apply(pw_z, 2, sum))
	pz <- matrix(1 / Z, Z)

################################

	Q1 <- function() {
		sum <- 0
		for(d in 1:D){
			for(w in 1:W){
				sum <- sum + x[d, w] * sum(pz * pw_z[w, ] * pd_z[d, ])
			}
		}
		sum
	}

	q1 <- Q1()

	for(t in 1:MI){
		cat("Iteration : ", t, ", Q1 : ", q1,  "\n")
#E-step 
		for(d in 1:D){
			for(w in 1:W){
				pz_dw[, d, w] <- pz * pw_z[w,] * pd_z[d, ] / sum(pz * pw_z[w,] * pd_z[d, ])
			}
		}
#M-step
		for(z in 1:Z){
			pw_z[, z] <- apply(x * pz_dw[z, , ], 2, sum) / sum(x * pz_dw[z,,])
		}
		pw_z <- t(t(pw_z) / apply(pw_z, 2, sum))
		
		for(d in 1:D){
			pd_z[d, ] <- apply(as.numeric(x[d, ]) * t(pz_dw[, d, ]), 2, sum) / sum(x * pz_dw[z,,])
		}
		pd_z <- t(t(pd_z) / apply(pd_z, 2, sum))
		
		s <- sum(x)
		for(z in 1:Z){
			pz[z] <- sum(x * pz_dw[z, , ]) / s
		}
		pz <- pz / sum(pz)

		qt1 <- Q1()
		if(abs(q1 - qt1) <= 1e-10) {
			break;
		} else {
			q1 <- qt1
		}

	}
#output
pd_z