[Depmix-commits] r98 - in trunk: . R

[noreply at r-forge.r-project.org]

Author: ingmarvisser
Date: 2008-03-22 00:45:47 +0100 (Sat, 22 Mar 2008)
New Revision: 98

Modified:
   trunk/NAMESPACE
   trunk/R/EM.R
   trunk/R/allGenerics.R
   trunk/R/depmix.fitted.R
   trunk/depmixNew-test2.R
   trunk/depmixNew-test3.R
Log:
Speeded up EM a little bit

Modified: trunk/NAMESPACE
===================================================================
--- trunk/NAMESPACE	2008-03-21 16:13:42 UTC (rev 97)
+++ trunk/NAMESPACE	2008-03-21 23:45:47 UTC (rev 98)
@@ -28,8 +28,11 @@
 	BIC,
 	fit,
 	npar,
+	freepars,
 	getdf,
 	nobs,
+	nresp,
+	ntimes,
 	depmix,
 	GLMresponse,
 	transInit,

Modified: trunk/R/EM.R
===================================================================
--- trunk/R/EM.R	2008-03-21 16:13:42 UTC (rev 97)
+++ trunk/R/EM.R	2008-03-21 23:45:47 UTC (rev 98)
@@ -8,28 +8,21 @@
 	lt <- length(ntimes)
 	et <- cumsum(ntimes)
 	bt <- c(1,et[-lt]+1)
-		
-	LL <- logLik(object)
 	
 	converge <- FALSE
 	j <- 0
 	
 	A <- object at trDens
+	B <- apply(object at dens,c(1,3),prod)
+	init <- object at init
+	
+	# initial expectation
+	fbo <- fb(init=object at init,A=object at trDens,B=apply(object at dens,c(1,3),prod),ntimes=ntimes(object))
+	LL <- fbo$logLike
+	LL.old <- LL + 1
+	
 	while(j <= maxit & !converge) {
-		
-		for(i in 1:ns) {
-			A[,,i] <- object at trDens[,,i]
-		}
-		
-		B <- apply(object at dens,c(1,3),prod)
-		init <- object at init
-		LL.old <- LL
-		j <- j+1
-		
-		# expectation
-		fbo <- fb(init=init,A=A,B=B,ntimes=ntimes(object))
-		LL <- fbo$logLike
-		
+				
 		# maximization
 				
 		# should become object at prior <- fit(object at prior)
@@ -62,13 +55,20 @@
 				object at dens[,k,i] <- dens(object at response[[i]][[k]])
 			}
 		}
+		
+		# expectation
+		fbo <- fb(init=object at init,A=object at trDens,B=apply(object at dens,c(1,3),prod),ntimes=ntimes(object))
+		LL <- fbo$logLike
 				
-		LL <- logLik(object)
 		if(verbose&((j%%5)==0)) cat("iteration",j,"logLik:",LL,"\n")
 		if( (LL >= LL.old) & (LL - LL.old < tol))  {
 			cat("iteration",j,"logLik:",LL,"\n")
 			converge <- TRUE
 		}
+		
+		LL.old <- LL
+		j <- j+1
+		
 	}
 	
 	class(object) <- "depmix.fitted"
@@ -79,7 +79,8 @@
 	object at conMat <- matrix()
 	
 	# what do we want in slot posterior?
-	object at posterior <- viterbi(object)
+	# this is moved to depmix.fit
+	# object at posterior <- viterbi(object)
 	
 	object
 }

Modified: trunk/R/allGenerics.R
===================================================================
--- trunk/R/allGenerics.R	2008-03-21 16:13:42 UTC (rev 97)
+++ trunk/R/allGenerics.R	2008-03-21 23:45:47 UTC (rev 98)
@@ -36,6 +36,8 @@
 
 setGeneric("fit", function(object, ...) standardGeneric("fit"))
 
+setGeneric("posterior", function(object, ...) standardGeneric("posterior"))
+
 setGeneric("predict", function(object, ...) standardGeneric("predict"))
 
 setGeneric("AIC", function(object, ..., k=2) standardGeneric("AIC"))

Modified: trunk/R/depmix.fitted.R
===================================================================
--- trunk/R/depmix.fitted.R	2008-03-21 16:13:42 UTC (rev 97)
+++ trunk/R/depmix.fitted.R	2008-03-21 23:45:47 UTC (rev 98)
@@ -131,10 +131,11 @@
 			# put the result back into the model
 			allpars[!fixed] <- result$par
 			object <- setpars(object,allpars)
-			object at posterior <- viterbi(object)
 			
-		}		
+		}
 		
+		object at posterior <- viterbi(object)
+		
 		return(object)
 	}
 )
@@ -190,3 +191,9 @@
 		}
 	}
 )
+
+setMethod("posterior","depmix.fitted",
+	function(object) {
+		return(object at posterior)
+	}
+)
Modified: trunk/depmixNew-test2.R
===================================================================
--- trunk/depmixNew-test2.R	2008-03-21 16:13:42 UTC (rev 97)
+++ trunk/depmixNew-test2.R	2008-03-21 23:45:47 UTC (rev 98)
@@ -10,6 +10,10 @@
 # test model with EM optimization, no covariates
 # 
 
+setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/rforge/depmix/trunk/")
+
+library(depmixS4)
+
 data(speed)
 
 trstart=c(0.899,0.101,0.084,0.916)
@@ -21,6 +25,52 @@
 
 logLik(mod)
 
+source("R/EM.R")
+
+gc()
+Rprof()
+mod1 <- em(mod,verb=TRUE)
+Rprof(NULL)
+summaryRprof()
+
+
+
+dd <- mod1 at dens[1:100,,]
+
+iter=1000
+
+gc()
+system.time(
+	for(i in 1:iter) {
+		x<-apply(dd,c(1,3),prod)
+	}
+)
+
+gc()
+system.time(
+	for(i in 1:iter) {
+		x<-exp(apply(log(dd),c(1,3),sum))
+	}
+)
+
+
+
+source("R/EM.R")
+
+gc()
+Rprof()
+mod1 <- em(mod)
+Rprof(NULL)
+summaryRprof()
+
+
+
+
+
+# 		object at parameters$coefficients <- pars$coefficients
+
+
+
 mod1 <- fit(mod)
 ll <- logLik(mod1)
 
Modified: trunk/depmixNew-test3.R
===================================================================
--- trunk/depmixNew-test3.R	2008-03-21 16:13:42 UTC (rev 97)
+++ trunk/depmixNew-test3.R	2008-03-21 23:45:47 UTC (rev 98)
@@ -9,6 +9,10 @@
 # BALANCE SCALE data example with age as covariate on class membership
 # 
 
+setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/rforge/depmix/trunk/")
+
+library(depmixS4)
+
 data(balance)
 
 # now fit some latent class models
@@ -26,10 +30,30 @@
 pars <- getpars(mod)
 fixed <- c(1,0,1,1,1,1,rep(c(1,0),8))
 
+gc()
+Rprof()
+mod1 <- em(mod)
+Rprof(NULL)
+summaryRprof()
+
+source("R/EM.R")
+source("R/responses.R")
+
+gc()
+Rprof()
 mod1 <- fit(mod,fixed=fixed)
+Rprof(NULL)
+summaryRprof()
 
+
 mod1
 
+setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/rforge/depmix/trunk/")
+
+source("R/depmix.fitted.R")
+source("R/viterbi.R")
+
+
 # 'log Lik.' -1083.036 (df=9)
 
 #