[Depmix-commits] r619 - pkg/depmixS4/R

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

Author: maarten
Date: 2014-02-05 22:29:21 +0100 (Wed, 05 Feb 2014)
New Revision: 619

Modified:
   pkg/depmixS4/R/EM.R
Log:
- fixed bug in classification = "hard" option; now uses proper classification likelihood criterion to determine convergence
- speed improvement in emviterbi help function in em 

Modified: pkg/depmixS4/R/EM.R
===================================================================
--- pkg/depmixS4/R/EM.R	2014-02-04 13:50:43 UTC (rev 618)
+++ pkg/depmixS4/R/EM.R	2014-02-05 21:29:21 UTC (rev 619)
@@ -24,6 +24,7 @@
     out
 }
 
+
 emviterbi <- function(A,B,init,ntimes,nstates,homogeneous,na.allow=TRUE) {
     # used for EM with hard classification, so that we don't need to change the object...
     # returns the most likely state sequence
@@ -36,6 +37,7 @@
     
     delta <- psi <- matrix(nrow=nt,ncol=ns)
     state <- vector(length=nt)
+    pstate <- vector(length=nt)
     
     prior <- init
     
@@ -67,10 +69,8 @@
           
         }
       }
-      
       # trace maximum likely state
       state[et[case]] <- which.max(delta[et[case],])
-      
       # this doesn't need a for loop does it???? FIX ME
       if(ntimes[case]>1) {
         for(i in (et[case]-1):bt[case]) {
@@ -78,8 +78,17 @@
         }
       }
     }
-    
-    delta <- data.frame(state,delta) 	
+    # compute the unconditional probability of the state sequence
+    pstate[bt] <- prior[cbind(1:lt,state[bt])]
+    btt <- bt[ntimes>1]
+    ett <- et[ntimes>1]
+    idx <- unlist(mapply(seq,btt+1,ett))
+    if(!homogeneous) {
+      pstate[idx] <- A[cbind(idx,state[idx],state[idx-1])]
+    } else {
+      pstate[idx] <- A[cbind(1,state[idx],state[idx-1])]
+    }
+    delta <- data.frame(state,pstate,delta)
     return(delta)
 }
 
@@ -138,7 +147,7 @@
 		  vstate <- apply(gamma,1,which.max)
 		  B <- dens
 		  if(na.allow) B[is.na(B)] <- 1
-		  fbo$logLike <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
+		  fbo$logLike <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes),vstate)])) + sum(log(init[cbind(1:lt,vstate)]))
 		} else {
 		  fbo <- fb(init=init,matrix(0,1,1),B=dens,ntimes=ntimes(object))
 		}
@@ -154,7 +163,7 @@
 		  vstate <- apply(fbo$gamma,1,which.max)
 		  B <- object at dens
 		  if(na.allow) B[is.na(B)] <- 1
-		  fbo$logLike <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
+		  fbo$logLike <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes),vstate)])) + sum(log(init[cbind(1:lt,vstate)]))
           
 		}
 		LL <- fbo$logLike
@@ -185,7 +194,7 @@
 		  vstate <- apply(fbo$gamma,1,which.max)
 		  B <- dens
 		  if(na.allow) B[is.na(B)] <- 1
-		  fbo$logLike <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
+		  fbo$logLike <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes),vstate)])) + sum(log(init[cbind(1:lt,vstate)]))
 		}
 		LL <- fbo$logLike
 		# print stuff
@@ -288,13 +297,15 @@
 	# initial expectation
   if(clsf == "hard") {
     fbo <- list()
-	  vstate <- emviterbi(A=trDens,B=dens,init=init,ntimes=object at ntimes,nstates=ns,homogeneous=object at homogeneous,na.allow=na.allow)[,1]
+    vit <- emviterbi(A=trDens,B=dens,init=init,ntimes=object at ntimes,nstates=ns,homogeneous=object at homogeneous,na.allow=na.allow)
+	  vstate <- vit[,1]
+    pstate <- vit[,2]
 	  fbo$gamma <- as.matrix(model.matrix(~ factor(vstate,levels=1:ns) - 1))
 	  fbo$xi <- array(0,dim=c(sum(ntimes),ns,ns))
 	  fbo$xi[cbind(1:(sum(ntimes)- 1),vstate[-1],vstate[-length(vstate)])] <- 1
 	  B <- dens
 	  if(na.allow) B[is.na(B)] <- 1
-	  fbo$logLike <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
+	  fbo$logLike <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes),vstate)])) + sum(log(pstate))
 	} else {
 	  fbo <- fb(init=init,A=trDens,B=dens,ntimes=ntimes(object),homogeneous=object at homogeneous)
   }
@@ -333,22 +344,26 @@
 		
 		for(i in 1:ns) {
 			for(k in 1:nresp(object)) {
-				response[[i]][[k]] <- fit(response[[i]][[k]],w=fbo$gamma[,i])
-				# update dens slot of the model
-				dens[,k,i] <- dens(response[[i]][[k]])
+				if(sum(fbo$gamma[,i])>0) {
+          response[[i]][[k]] <- fit(response[[i]][[k]],w=fbo$gamma[,i])
+  				# update dens slot of the model
+  				dens[,k,i] <- dens(response[[i]][[k]])
+				}
 			}
 		}
 		
 		if(clsf == "hard") {
       fbo <- list()
-      vstate <- emviterbi(A=trDens,B=dens,init=init,ntimes=object at ntimes,nstates=ns,homogeneous=object at homogeneous,na.allow=na.allow)[,1]
+      vit <- emviterbi(A=trDens,B=dens,init=init,ntimes=object at ntimes,nstates=ns,homogeneous=object at homogeneous,na.allow=na.allow)
+      vstate <- vit[,1]
+      pstate <- vit[,2]
 		  #vstate <- viterbi(object)[,1]
 		  fbo$gamma <- as.matrix(model.matrix(~ factor(vstate,levels=1:ns) - 1))
 		  fbo$xi <- array(0,dim=c(sum(ntimes),ns,ns))
 		  fbo$xi[cbind(1:(sum(ntimes)- 1),vstate[-1],vstate[-length(vstate)])] <- 1
 		  B <- dens
 		  if(na.allow) B[is.na(B)] <- 1
-		  fbo$logLike <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes),vstate)]))
+		  fbo$logLike <- sum(log((apply(B,c(1,3),prod))[cbind(1:sum(ntimes),vstate)])) + sum(log(pstate))
 		} else {
 		  # expectation
 		  fbo <- fb(init=init,A=trDens,B=dens,ntimes=ntimes(object),homogeneous=object at homogeneous)