[Depmix-commits] r45 - trunk

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

Author: ingmarvisser
Date: 2008-03-05 17:14:16 +0100 (Wed, 05 Mar 2008)
New Revision: 45

Modified:
   trunk/EM.R
Log:
Added option for stationary models in EM update of transition matrix, see example file: depmix-test3EM2.R

Modified: trunk/EM.R
===================================================================
--- trunk/EM.R	2008-03-05 00:32:03 UTC (rev 44)
+++ trunk/EM.R	2008-03-05 16:14:16 UTC (rev 45)
@@ -1,59 +1,76 @@
-em <- function(object,maxit=100,tol=1e-5,verbose=FALSE,...) {
-  if(!is(object,"hmModel")) stop("object must be 'hmModel'")
-
-  # pseudocode
-  ns <- object at nstates
-
-  LL <- logLik(object)
-  
-  converge <- FALSE
-  j <- 0
-
-  A <- object at trans
-  while(j <= maxit & !converge) {
-    for(i in 1:ns) {
-        A[,,i] <- predict(object at trModels[[i]])
-    }
-    B <- exp(apply(object at logdens,c(1,3),sum))
-    # TODO: add functionality for inModel
-    init <- predict(object at initModel)
-    LL.old <- LL
-    j <- j+1
-
-    # expectation
-    fbo <- fb(init=init,A=A,B=B,ntimes=object at ntimes)
-    LL <- fbo$logLike
-    
-    # maximization
-    
-    #object at init <- fit(object at init,ip=fbo$gamma[1,])
-    #object at init <- matrix(fbo$gamma[1,],nrow=1)
+em <- function(object,maxit=100,tol=1e-9,verbose=FALSE,...) {
+	if(!is(object,"hmModel")) stop("object must be 'hmModel'")
+	
+	# pseudocode
+	ns <- object at nstates
+	
+	LL <- logLik(object)
+	
+	converge <- FALSE
+	j <- 0
+	
+	A <- object at trans
+	while(j <= maxit & !converge) {
 		
-
-    object at initModel <- setpars(object at initModel,values=object at initModel@family$linkfun(fbo$gamma[1,],base=object at initModel@family$base))
-    # This should become:
-    # lt <- length(object at ntimes)
+		for(i in 1:ns) {
+			A[,,i] <- predict(object at trModels[[i]])
+		}
+		
+		B <- exp(apply(object at logdens,c(1,3),sum))
+		# TODO: add functionality for inModel
+		init <- predict(object at initModel)
+		LL.old <- LL
+		j <- j+1
+		
+		# expectation
+		fbo <- fb(init=init,A=A,B=B,ntimes=object at ntimes)
+		LL <- fbo$logLike
+		
+		# maximization
+		
+		#object at init <- fit(object at init,ip=fbo$gamma[1,])
+		#object at init <- matrix(fbo$gamma[1,],nrow=1)
+		
+		object at initModel <- setpars(object at initModel,values=object at initModel@family$linkfun(fbo$gamma[1,],base=object at initModel@family$base))
+		# This should become:
+		# lt <- length(object at ntimes)
 		# et <- cumsum(object at ntimes)
 		# bt <- c(1,et[-lt]+1)
-    # object at initModel@y <- fbo$gamma[bt,]
-    # object at initModel <- fit(object at initModel)
-    
-    for(i in 1:ns) {
-      object at trModels[[i]]@y <- fbo$xi[,,i]/fbo$gamma[,i]
-      object at trModels[[i]] <- fit(object at trModels[[i]],w=as.matrix(fbo$gamma[,i]),ntimes=object at ntimes) # check this
-      #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]]))
-      object at trans[,,i] <- predict(object at trModels[[i]])
-      
-      for(k in 1:object at nresp) {
-        object at rModels[[i]][[k]] <- fit(object at rModels[[i]][[k]],w=fbo$gamma[,i])
-        object at logdens[,k,i] <- logDens(object at rModels[[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
-  }
-  object
+		# 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) {
+			
+			if(!object at stationary) {
+				object at trModels[[i]]@y <- fbo$xi[,,i]/fbo$gamma[,i]
+				object at trModels[[i]] <- fit(object at trModels[[i]],w=as.matrix(fbo$gamma[,i]),ntimes=object at ntimes) # check this
+			} else {
+				for(k in 1:ns) {
+					trm[i,k] <- sum(fbo$xi[-c(et),k,i])/sum(fbo$gamma[-c(et),i])
+				}
+				object at trModels[[i]]@parameters$coefficients <- object at trModels[[i]]@family$linkfun(trm[i,],base=object at initModel@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]]))
+			
+			object at trans[,,i] <- predict(object at trModels[[i]])
+			
+			for(k in 1:object at nresp) {
+				object at rModels[[i]][[k]] <- fit(object at rModels[[i]][[k]],w=fbo$gamma[,i])
+				object at logdens[,k,i] <- logDens(object at rModels[[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
+	}
+	
+	object
 }