参考了Daichi Mochihashi 的c版本的lda代码
数学原理还是半懂不懂，要继续看
数据格式
train.txt (一行一个文档，单词下标:出现的次数)
1:5 2:5 2:10
9:2 5:23 8:9
vocabulary.txt
a //index 为1
happy //index 为2
hello
处理beta结果，提取题目的代码
ProcessBeta.scala
<pre name="code" class="java">package lda_scala
import scala.collection.mutable.ArrayBuffer
import scala.io.Source
object ProcessBeta {
def main(args: Array[String]): Unit = {
val topWordArray = processBetaDocuments("result_beta.txt", 50, 10474)
mapVocabulary("vocab.txt", topWordArray)
}
def mapVocabulary(vocabularyPath: String, topWordArray: ArrayBuffer[ArrayBuffer[Int]]) {
var vocaArray = ArrayBuffer[String]()
for(line <- Source.fromFile(vocabularyPath).getLines()){
vocaArray += line
}
//println(vocaArray.take(10))
topWordArray.foreach{
x => {
print("topic: ")
x.map(x=> print(vocaArray(x) + " "))
println("")
}
}
}
def processBetaDocuments(filepath: String, nClass: Int, nLex: Int, topNumber: Int = 10): ArrayBuffer[ArrayBuffer[Int]] ={
//topicVocab(prob, Wordindex)
var topicVocab = ArrayBuffer.fill(nClass, nLex)(0.0, 0)
//println(topicVocab.size)
var wordIndex = 0
for(line <- Source.fromFile(filepath).getLines()){
// println("line: " + line)
val wArrProb = line.split("\\s+").map(x=>x.toDouble)
for(topIndex <- 0 to nClass - 1){
//for(topIndex <- 0 to ){
//println("wArrProb: " + wArrProb(topIndex))
topicVocab(topIndex)(wordIndex) = (wArrProb(topIndex), wordIndex)
// println()
}
wordIndex += 1
}
//println(topicVocab mkString " ")
var topWordArray = ArrayBuffer.fill(nClass, topNumber)(0)
//println(topicVocab.size)
var tIndex = 0
topicVocab.foreach{
x => {
// println(x.sortWith(_ > _).take(10))
println(x.sortBy(_._1).reverse.take(10))
topWordArray(tIndex) = x.sortBy(_._1).reverse.take(10).map(x => x._2)
tIndex += 1
}
}
println(topWordArray.size)
topWordArray
}
}lda的代码
lda2.scala
<pre name="code" class="java">package lda_scala
import scala.io.Source
import scala.collection.mutable.ArrayBuffer
import scala.util.Random
import java.io.FileWriter
import scala.util.control.Breaks._
//import Document
object lda2 {
val zeroBegin = true
def main(args: Array[String]): Unit = {
val s = "a b c d e"
println(s.split(" ").size)
//readDocuments("testTrain.txt")
val aptha = Array[Double]() // with k member
lda("ap.dat")
}
def lda(datapath: String, nClass: Int = 50, emmax: Int =100 , demmax: Int = 20, epsilon: Double = 1.0e-4): Unit = {
val fileAP = "result_alpha.txt"
val fileBP = "result_beta.txt"
val(dataV, maxIdV, maxLenV) = readDocuments2(datapath)
//val(dataV, maxIdV, maxLenV) = readDocuments(datapath)
var data = dataV
var nLex = maxIdV
var dLenMax = maxLenV
println("nLex: " + nLex + " " + "dLenMax: " + dLenMax)
//var alpha = ArrayBuffer.fill(nClass)(Random.nextDouble)
//var beta = ArrayBuffer.fill(nLex, nClass)(1.0 / nLex)
val(alpha_result, beta_result) = lda_learn(data,nClass, nLex, dLenMax, emmax, demmax, epsilon)
lda_write(fileAP, fileBP, alpha_result, beta_result, nClass, nLex)
}
def lda_learn(data: ArrayBuffer[Document], nClass: Int, nLex: Int, dLenMax: Int,
emmax: Int, demmax: Int, epsilon: Double): (ArrayBuffer[Double], ArrayBuffer[ArrayBuffer[Double]]) = {
//val gamma: ArrayBuffer[Double] = ArrayBuffer[Double]()
val numberOfDocument = data.size
//初始化alpha
var alpha = ArrayBuffer.fill(nClass)(Random.nextInt(Int.MaxValue).toDouble / Int.MaxValue)
val z = alpha.reduce(_+_)
alpha = alpha.map(x => x / z)
//排序alpha
alpha = alpha.sortWith(_ > _)
//初始化beta
//println("nlex: "+ nLex)
//var beta = ArrayBuffer.fill(nLex, nClass)(1.0 / nLex)
var beta = ArrayBuffer.fill(nLex, nClass)(Random.nextInt(Int.MaxValue).toDouble / Int.MaxValue *10)
var tot = 0.0
for(j <- 0 to nClass - 1){
for(i <- 0 to nLex -1){
beta(i)(j) = Random.nextInt(Int.MaxValue).toDouble / Int.MaxValue *10
tot += beta(i)(j)
}
for(i <- 0 to nLex -1)
beta(i)(j) = beta(i)(j) / tot
tot = 0.0
}
//初始化后验dirichlet的参数
var gammas = ArrayBuffer.fill(numberOfDocument, nClass)(0.0)
var betas = ArrayBuffer.fill(nLex, nClass)(0.0)
//初始化要用到的缓冲区
var q = ArrayBuffer.fill(dLenMax, nClass)(0.0)
var gamma = ArrayBuffer.fill(nClass)(0.0)
var ap = ArrayBuffer.fill(nClass)(0.0)
var nt = ArrayBuffer.fill(nClass)(0.0)
var pnt = ArrayBuffer.fill(nClass)(0.0)
println("Number of documents = " + numberOfDocument)
println("Number of words = " + nLex)
println("Number of latent Classes = " + nClass)
println("Number of outer EM iteration = " + emmax)
println("Number of inner EM iteration = " + demmax)
println("Convergence threshold = " + epsilon)
var pppl = 0.0
for(t <- 0 to emmax-1){
printf("iteration %d/%d..\n", t+1, emmax)
//variational inference
for(i <- 0 to numberOfDocument-1){
val (gammaVal, qVal) = vbem(data, i, gamma, q, nt, pnt, ap, alpha, beta, data(i).len, nClass, demmax)
gamma = gammaVal
q = qVal
val gammasVal = accum_gammas(gammas, gamma, i, nClass)
gammas = gammasVal
val betasVal = accum_betas(betas, q, nClass, data, i)
betas = betasVal
}
//VB M-step
//Newton-Raphson for alpha
val(alphaV, gammasV) = newton_alpha(alpha, gammas, numberOfDocument, nClass, 0)
alpha = alphaV
gammas = gammasV
val(betaV, betasV) = normalize_matrix_col(beta, betas, nLex, nClass)
beta = betaV
betas = betasV
//clean Buffer
betas = ArrayBuffer.fill(nLex, nClass)(0.0)
//converge?
val ppl = lda_ppl(data, beta, gammas, numberOfDocument, nClass)
if((t > 1) && (Math.abs((ppl - pppl)/pppl) < epsilon)){
if(t < 5){
printf("\nearly convergence. restarting..\n")
val (alpha_re, beta_re) = lda_learn(data, nClass, nLex, dLenMax, emmax, demmax, epsilon)
return (alpha_re, beta_re)
}else{
printf("converged")
return (alpha, beta)
}
}
pppl = ppl
}
//println(alpha mkString ",")
return (alpha,beta)
}
//correct
def lda_ppl(data: ArrayBuffer[Document], beta: ArrayBuffer[ArrayBuffer[Double]],
gammas: ArrayBuffer[ArrayBuffer[Double]], m: Int, nClass: Int): Double = {
var n = 0.0
val ds = data.size
for(i <- 0 to ds -1)
for(j <- 0 to data(i).len-1)
n += data(i).cnt(j)
return Math.exp(- lda_lik(data, beta, gammas, m, nClass) / n)
}
//correct
def lda_lik(data: ArrayBuffer[Document], beta: ArrayBuffer[ArrayBuffer[Double]],
gammasVal: ArrayBuffer[ArrayBuffer[Double]], m: Int, nClass: Int): Double = {
val numberOfDocument = data.size
var gammas = gammasVal
var egammas = ArrayBuffer.fill(m, nClass)(0.0)
val(egammasV, gammasV) = normalize_matrix_row(egammas, gammas, m, nClass)
egammas = egammasV
gammas = gammasV
var lik = 0.0
for(i <- 0 to numberOfDocument-1){
for(j <- 0 to data(i).len -1){
var z = 0.0
for(k <- 0 to nClass -1)
if(zeroBegin)
z += beta(data(i).id(j))(k) * egammas(i)(k)
else
z += beta(data(i).id(j)-1)(k) * egammas(i)(k)
lik += data(i).cnt(j) * Math.log(z)
}
}
return lik
}
//correct
def normalize_matrix_row(dst: ArrayBuffer[ArrayBuffer[Double]], src: ArrayBuffer[ArrayBuffer[Double]],
rows: Int, cols: Int): (ArrayBuffer[ArrayBuffer[Double]], ArrayBuffer[ArrayBuffer[Double]]) = {
for(i <- 0 to rows-1){
var z = 0.0
for(j <- 0 to cols -1)
z += src(i)(j)
for(j <- 0 to cols -1)
dst(i)(j) = src(i)(j) / z
}
return (dst, src)
}
//correct
def normalize_matrix_col(dst: ArrayBuffer[ArrayBuffer[Double]], src: ArrayBuffer[ArrayBuffer[Double]],
rows: Int, cols: Int): (ArrayBuffer[ArrayBuffer[Double]], ArrayBuffer[ArrayBuffer[Double]]) = {
for(j <- 0 to cols-1){
var z = 0.0
for(i <- 0 to rows -1)
z += src(i)(j)
for(i <- 0 to rows -1)
dst(i)(j) = src(i)(j) / z
}
return (dst, src)
}
//M: number of Document, K: number of Topic
def newton_alpha(alpha: ArrayBuffer[Double], gammas: ArrayBuffer[ArrayBuffer[Double]],
M: Int, K: Int, level: Int): (ArrayBuffer[Double], ArrayBuffer[ArrayBuffer[Double]]) = {
var g = ArrayBuffer.fill(K)(0.0)
var h = ArrayBuffer.fill(K)(0.0)
var pg = ArrayBuffer.fill(K)(0.0)
var palpha = ArrayBuffer.fill(K)(0.0)
if(level == 0){
for( i <- 0 to K-1){
var z = 0.0
for(j <- 0 to M-1 )
z += gammas(j)(i)
alpha(i) = z / (M * K)
}
}else{
for(i <- 0 to K-1){
var z = 0.0
for(j <- 0 to M-1)
z += gammas(j)(i)
alpha(i) = z / (M*K*Math.pow(10, level))
}
}
var psg = 0.0
for(i <- 0 to M-1){
var gs = 0.0
for(j <- 0 to K-1)
gs += gammas(i)(j)
psg += digamma(gs)
}
for(i <- 0 to K-1){
var spg = 0.0
for(j <- 0 to M-1)
spg += digamma(gammas(j)(i))
pg(i) = spg - psg
}
val MAX_NEWTON_ITERATION = 20
for(t <- 0 to MAX_NEWTON_ITERATION){
var alpha0 = 0.0
for(i <- 0 to K-1)
alpha0 += alpha(i)
var palpha0 = digamma(alpha0)
for(i <- 0 to K-1)
g(i) = M*(palpha0 - digamma(alpha(i))) + pg(i)
for(i <- 0 to K-1)
h(i) = -1 / trigamma(alpha(i))
var sh = 0.0
for(i <- 0 to K-1)
sh += h(i)
var hgz = 0.0
for(i <- 0 to K-1)
hgz += g(i) * h(i)
hgz /= (1 / trigamma(alpha0) + sh)
for(i <- 0 to K-1)
alpha(i) = alpha(i) - h(i) * (g(i) - hgz) / M
val MAX_RECURSION_LIMIT = 10
for(i <- 0 to K-1)
if(alpha(i) < 0){
if(level >= MAX_RECURSION_LIMIT){
println("error")
}else{
return newton_alpha(alpha, gammas, M, K, 1+level)
}
}
if((t>0) && converged(alpha, palpha, K, 1.0e-4)){
return (alpha, gammas)
}else{
for(i <- 0 to K-1)
palpha(i) = alpha(i)
}
}
return (alpha, gammas)
}
//tested correct
def accum_gammas(gammas: ArrayBuffer[ArrayBuffer[Double]], gamma: ArrayBuffer[Double],
n: Int, K: Int): ArrayBuffer[ArrayBuffer[Double]] = {
for(k <- 0 to K-1)
gammas(n)(k) = gamma(k)
return gammas
}
// (beta, betas) tested correct
def accum_betas(betas: ArrayBuffer[ArrayBuffer[Double]], q: ArrayBuffer[ArrayBuffer[Double]], K: Int,
data: ArrayBuffer[Document], indexDocument: Int): ArrayBuffer[ArrayBuffer[Double]] = {
val n = data(indexDocument).len
//println("document len: " + n)
for(i <- 0 to n-1)
for(k <- 0 to K-1){
//println("i: "+ i + " k: " + k)
//println("id: " + data(indexDocument).id(i))
//需要减1，因为id是以1开始
if(zeroBegin)
betas(data(indexDocument).id(i))(k) += q(i)(k) * data(indexDocument).cnt(i)
else
betas(data(indexDocument).id(i)-1)(k) += q(i)(k) * data(indexDocument).cnt(i)
}
return betas
}
def vbem(data: ArrayBuffer[Document], indexDocument: Int, gamma: ArrayBuffer[Double], q: ArrayBuffer[ArrayBuffer[Double]],
nt: ArrayBuffer[Double], pnt: ArrayBuffer[Double], ap: ArrayBuffer[Double],
alpha: ArrayBuffer[Double], beta: ArrayBuffer[ArrayBuffer[Double]], L: Int, K: Int, emmax: Int)
: (ArrayBuffer[Double], ArrayBuffer[ArrayBuffer[Double]]) = {
//var ap = apVal
//每次调用都是新的
for(k <- 0 to K-1){
nt(k) = L.toDouble / K
}
var isConverged = false
//var j = 0
breakable{
for(j <- 0 to (emmax-1)){
//每次调用都是重新赋值， 与alpha有关
for(k <- 0 to K-1){
ap(k) = Math.exp(digamma(alpha(k) + nt(k)))
}
//accumulate q
//每次都重新赋值，每个文档一个q矩阵
for(l <- 0 to L-1)
for(k <- 0 to K-1)
if(zeroBegin)
q(l)(k) = beta(data(indexDocument).id(l))(k) * ap(k)
else
q(l)(k) = beta(data(indexDocument).id(l) - 1)(k) * ap(k)
//normalize
for(l <- 0 to L-1){
var z = 0.0
for(k <- 0 to K-1)
z += q(l)(k)
for(k <- 0 to K-1)
q(l)(k) = q(l)(k) / z
}
//vb-mstep
for(k <- 0 to K -1){
var z =0.0
for(l <- 0 to L - 1)
z += q(l)(k) * data(indexDocument).cnt(l)
nt(k) = z
}
//converge
if(j > 0 && converged(nt, pnt, K, 1.0e-2))
break//isConverged = true
for(k <- 0 to K-1)
pnt(k) = nt(k)
}
}
for(k <- 0 to K-1)
gamma(k) = alpha(k) + nt(k)
return (gamma, q)
}
//len 9:2 1:8 以文本长度开头这种格式 readDocuments tested correct
def readDocuments2(filepath: String): (ArrayBuffer[Document], Int, Int) ={
var documents = ArrayBuffer[Document]()
var d = 0
var maxid = -1
var maxlen = 0
for(line <- Source.fromFile(filepath).getLines()){
println(line)
val len = line.split(" ").size - 1
val idArray = line.split(" ").tail.map(x=>x.split(":")(0).toInt)
val cntArray = line.split(" ").tail.map(x=>x.split(":")(1).toDouble)
val document = new Document(len, idArray, cntArray)
val temp_maxid = idArray.reduce((x,y) => if(x>y) x else y)
if(len > maxlen)
maxlen = len
if(temp_maxid > maxid)
maxid = temp_maxid
documents += document
}
if(zeroBegin)
maxid += 1
println(documents mkString "\n")
//Array(new Document(1,2,3))
(documents, maxid, maxlen)
}
//readDocuments tested correct
def readDocuments(filepath: String): (ArrayBuffer[Document], Int, Int) ={
var documents = ArrayBuffer[Document]()
var d = 0
var maxid = -1
var maxlen = 0
for(line <- Source.fromFile(filepath).getLines()){
println(line)
val len = line.split(" ").size
val idArray = line.split(" ").map(x=>x.split(":")(0).toInt)
val cntArray = line.split(" ").map(x=>x.split(":")(1).toDouble)
val document = new Document(len, idArray, cntArray)
val temp_maxid = idArray.reduce((x,y) => if(x>y) x else y)
if(len > maxlen)
maxlen = len
if(temp_maxid > maxid)
maxid = temp_maxid
documents += document
}
println(documents mkString "\n")
//Array(new Document(1,2,3))
(documents, maxid, maxlen)
}
def converged(u: ArrayBuffer[Double], v: ArrayBuffer[Double], n: Int, threshold: Double): Boolean = {
var us = 0.0
var ds = 0.0
var d = 0.0
for(i <- 0 to n-1)
us += u(i) * u(i)
for(i <- 0 to n-1){
d = u(i) - v(i)
ds += d*d
}
if(Math.sqrt(ds / us) < threshold)
return true
else
return false
}
//tested correct
def digamma(xval: Double): Double = {
var x = xval
var result = 0.0
val neginf = -1.0 / 0.0
val c = 12
val s = 1e-6
val d1 = -0.57721566490153286
val d2 = 1.6449340668482264365
val s3 = 1.0/12
val s4 = 1.0/120
val s5 = 1.0/252
val s6 = 1.0/240
val s7 = 1.0/132
val s8 = 691/32760
val s9 = 1/12
val s10 = 3617/8160
/* Illegal arguments */
/*
if((x == neginf) || isnan(x)) {
return 0.0/0.0;
}
*
*/
/* Singularities */
if((x <= 0) && (Math.floor(x) == x)) {
return neginf;
}
/* Negative values */
/* Use the reflection formula (Jeffrey 11.1.6):
* digamma(-x) = digamma(x+1) + pi*cot(pi*x)
*
* This is related to the identity
* digamma(-x) = digamma(x+1) - digamma(z) + digamma(1-z)
* where z is the fractional part of x
* For example:
* digamma(-3.1) = 1/3.1 + 1/2.1 + 1/1.1 + 1/0.1 + digamma(1-0.1)
* = digamma(4.1) - digamma(0.1) + digamma(1-0.1)
* Then we use
* digamma(1-z) - digamma(z) = pi*cot(pi*z)
*/
if(x < 0) {
return digamma(1-x) + Math.PI/Math.tan(-Math.PI*x);
}
/* Use Taylor series if argument <= S */
if(x <= s) return d1 - 1/x + d2*x;
/* Reduce to digamma(X + N) where (X + N) >= C */
result = 0;
while(x < c) {
result -= 1/x;
x = x + 1;
}
/* Use de Moivre's expansion if argument >= C */
/* This expansion can be computed in Maple via asympt(Psi(x),x) */
if(x >= c) {
var r = 1/x;
result += Math.log(x) - 0.5*r;
r *= r;
result -= r * (s3 - r * (s4 - r * (s5 - r * (s6 - r * s7))));
}
return result;
}
//tested correct
def trigamma(xVal: Double): Double = {
var x = xVal
var result = 0.0
val neginf = -1.0/0.0
val small = 1e-4
val large = 8
val c = 1.6449340668482264365 /* pi^2/6 = Zeta(2) */
val c1 = -2.404113806319188570799476 /* -2 Zeta(3) */
val b2 = 1./6
val b4 = -1./30
val b6 = 1./42
val b8 = -1./30
val b10 = 5./66
/* Illegal arguments */
/*
if((x == neginf) || isnan(x)) {
return 0.0/0.0;
}
* */
/* Singularities */
if((x <= 0) && (Math.floor(x) == x)) {
return -neginf;
}
/* Negative values */
/* Use the derivative of the digamma reflection formula:
* -trigamma(-x) = trigamma(x+1) - (pi*csc(pi*x))^2
*/
if(x < 0) {
result = Math.PI/Math.sin(-Math.PI*x);
return -trigamma(1-x) + result*result;
}
/* Use Taylor series if argument <= small */
if(x <= small) {
return 1/(x*x) + c + c1*x;
}
result = 0;
/* Reduce to trigamma(x+n) where ( X + N ) >= B */
while(x < large) {
result += 1/(x*x);
x = x+1;
}
/* Apply asymptotic formula when X >= B */
/* This expansion can be computed in Maple via asympt(Psi(1,x),x) */
if(x >= large) {
var r = 1/(x*x);
result += 0.5*r + (1 + r*(b2 + r*(b4 + r*(b6 + r*(b8 + r*b10)))))/x;
}
return result;
}
//tested correct
def lda_write(fileAP: String, fileBP: String, alpha: ArrayBuffer[Double],
beta: ArrayBuffer[ArrayBuffer[Double]], nClass: Int, nLex: Int){
printf("writing model..\n")
write_vector(fileAP, alpha, nClass)
write_matrix(fileBP, beta, nLex, nClass)
printf("done.\n")
}
//tested correct
def write_vector(filePath: String, vector: ArrayBuffer[Double], n: Int){
println("vector: " + vector mkString " ")
val fw = new FileWriter(filePath, false)
fw.write(vector mkString " ")
fw.close()
}
//test correct
def write_matrix(filePath: String, matrix: ArrayBuffer[ArrayBuffer[Double]],
rows: Int, cols: Int){
val fw = new FileWriter(filePath, false)
for(i <- 0 to rows -1){
for(j <- 0 to cols -1){
//fw.write("This line appended to file!") ;
fw.append(matrix(i)(j).toString + " ")
}
fw.append("\n")
}
fw.close()
}
}
topic: attorney case investigation trial government judge department justice documents office
topic: trade united states japan japanese countries world foreign international economic
topic: mrs i heart hospital family years husband children medical daughter
topic: church pope catholic roman john people vatican religious women paul
topic: abortion senate bill vote rights committee house humphrey i law
topic: space shuttle nasa launch mission earth two venus test nuclear
topic: eastern pilots district court judge union two sinhalese mexican immigration
topic: dukakis bush campaign jackson democratic republican presidential president i new
topic: stock market index points stocks trading exchange million prices shares
topic: bill house rep senate committee legislation sen federal measure i
topic: army military navy british killed soldiers ireland force died two
topic: police people killed government reported city two today injured news
topic: fire water people miles officials river ship area two three
topic: prison court convicted police charges drug two judge guilty trial
topic: bush budget house president congress reagan defense billion administration white
topic: million oil company years officials hospitals state two year new
topic: prices economy price new percent stores industry sales business market
topic: dollar cents yen late lower cent futures prices higher oil
topic: west german east germany bush soviet summit united treaty europe
topic: workers union strike contract labor plant unions pay company two
topic: state new waste states federal epa environmental officials creek county
topic: percent year billion rate sales last increase rose million months
topic: south africa government african rebels black president mandela peace war
topic: school students student university schools education college board teachers teacher
topic: air plane flight aircraft planes airport pilot force two airlines
topic: china chinese government prisoners rights people hong last kong human
topic: aid million president new service administration officials government department policy
topic: farmers million tons farm agriculture corn drought year food agricultural
topic: computer company business new years att news animals technology year
topic: i people think going solidarity walesa american thats president says
topic: aids health drug disease children patients virus study percent blood
topic: company inc corp new stock billion co offer million share
topic: iraq kuwait iraqi iran gulf united war military hostages saddam
topic: police hospital yearold death two found i man mother shot
topic: panama noriega mecham government president defense military panamanian state senate
topic: population people time spain sweden cuba first game spanish year
topic: saudi troops arabia pentagon united soldiers states american military honduran
topic: i dont like just get think im people going says
topic: court case state ruling supreme appeals law federal courts drug
topic: soviet gorbachev union moscow soviets gorbachevs news president mikhail republics
topic: south korea north korean two president roh new nations games
topic: new york year years john music american film world first
topic: party government political elections opposition minister election communist parties national
topic: million bank billion money year banks federal estate last first
topic: people smoking i years last front square government azerbaijan cigarettes
topic: earthquake damage quake rights people human philippines new aquino military
topic: barry white people i reed mayor new black blacks state
topic: cbs nbc fair northern southern new snow texas inches temperatures
topic: tax housing people homeless taxes income new cities irs report
245
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
