[Depmix-commits] r60 - trunk

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

Author: ingmarvisser
Date: 2008-03-07 10:47:57 +0100 (Fri, 07 Mar 2008)
New Revision: 60

Modified:
   trunk/EM.R
Log:
Added convergence message to return object

Modified: trunk/EM.R
===================================================================
--- trunk/EM.R	2008-03-07 09:22:06 UTC (rev 59)
+++ trunk/EM.R	2008-03-07 09:47:57 UTC (rev 60)
@@ -1,4 +1,5 @@
 em <- function(object,maxit=100,tol=1e-6,verbose=FALSE,...) {
+	
 	if(!is(object,"depmix")) stop("object is not of class 'depmix'")
 	
 	# pseudocode
@@ -24,7 +25,6 @@
 		B <- apply(object at dens,c(1,3),prod)
 		# TODO: add functionality for inModel
 		init <- object at init
-# 		print(init)
 		LL.old <- LL
 		j <- j+1
 		
@@ -34,21 +34,16 @@
 		
 		# maximization
 		
-		#object at init <- fit(object at init,ip=fbo$gamma[1,])
-		#object at init <- matrix(fbo$gamma[1,],nrow=1)
-		
-		# FIX ME for length(ntimes)>1
-		# print(fbo$gamma[1,])
 		# Here we need an average of gamma[bt[case],], which may need to be weighted ?? (see Rabiner, p283)
-		
-
-		
 		# this is without weighting
 		initprobs <- apply(fbo$gamma[bt,],2,mean)
 		
 		# should become object at prior <- fit(object at prior)
 		object at prior@y <- fbo$gamma[bt,]
-		object at prior <- fit(object at prior,w=NULL)
+		object at prior <- fit(object at prior, w=NULL)
+		
+		# init needs to be recomputed here?
+		
 		#object at initModel <- setpars(object at initModel,values=object at initModel@family$linkfun(initprobs,base=object at initModel@family$base))
 		
 		# This should become:
@@ -57,9 +52,7 @@
 		# bt <- c(1,et[-lt]+1)
 		# object at initModel@y <- fbo$gamma[bt,]
 		# object at initModel <- fit(object at initModel)
-		
-		#et <- cumsum(object at ntimes)
-		
+				
 		trm <- matrix(0,ns,ns)
 		for(i in 1:ns) {
 			
@@ -75,9 +68,6 @@
 					object at transition[[i]]@parameters$coefficients <- object at transition[[i]]@family$linkfun(trm[i,],base=object at transition[[i]]@family$base)
 				}
 				
-				#object at trModels[[i]] <- fit(object at trModels[[i]],w=NULL,ntimes=object at ntimes) # check this
-				#object at trans[,,i] <- exp(logDens(object at trModels[[i]]))
-				
 				# update trans slot of the model
 				object at trDens[,,i] <- dens(object at transition[[i]])
 			}
@@ -88,13 +78,17 @@
 				object at dens[,k,i] <- dens(object at response[[i]][[k]])
 			}
 		}
-		
-		#object <- setpars(object,getpars(object)) # set parameters and recompute (bit of a roundabout way)
-		
+				
 		LL <- logLik(object)
 		if(verbose) cat("iteration",j,"logLik:",LL,"\n")
 		if( (LL >= LL.old) & (LL - LL.old < tol))  converge <- TRUE
 	}
 	
+	if(converge) object at message <- "Log likelihood converged to within tol."
+	else object at message <- "'maxit' iterations reached in EM without convergence."
+	
+	# no constraints in EM
+	object at conMat <- NULL
+	
 	object
 }