[Depmix-commits] r47 - trunk

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

Author: ingmarvisser
Date: 2008-03-06 11:38:48 +0100 (Thu, 06 Mar 2008)
New Revision: 47

Modified:
   trunk/fb.R
Log:
Added support for dealing with timeseries of length 1, ie for latent class models

Modified: trunk/fb.R
===================================================================
--- trunk/fb.R	2008-03-05 16:15:20 UTC (rev 46)
+++ trunk/fb.R	2008-03-06 10:38:48 UTC (rev 47)
@@ -18,10 +18,7 @@
 
 	nt <- nrow(B)	
 	ns <- ncol(init)
-		
 	
-# 	print(init)
-#  	stop("till here")
 	
 	alpha <- matrix(ncol=ns,nrow=nt)
 	beta <- matrix(ncol=ns,nrow=nt)
@@ -38,25 +35,32 @@
 		alpha[bt[case],] <- init[case,]*B[bt[case],] # initialize
 		sca[bt[case]] <- 1/sum(alpha[bt[case],])
 		alpha[bt[case],] <- alpha[bt[case],]*sca[bt[case]]
-		for(i in bt[case]:(et[case]-1)) {
-			if(stationary) alpha[i+1,] <- (A[1,,]%*%alpha[i,])*B[i+1,]
-			else alpha[i+1,] <- (A[i,,]%*%alpha[i,])*B[i+1,]
-			sca[i+1] <- 1/sum(alpha[i+1,])
-			alpha[i+1,] <- sca[i+1]*alpha[i+1,]
+		
+		if(ntimes[case]>1) {
+			for(i in bt[case]:(et[case]-1)) {
+				if(stationary) alpha[i+1,] <- (A[1,,]%*%alpha[i,])*B[i+1,]
+				else alpha[i+1,] <- (A[i,,]%*%alpha[i,])*B[i+1,]
+				sca[i+1] <- 1/sum(alpha[i+1,])
+				alpha[i+1,] <- sca[i+1]*alpha[i+1,]
+			}
 		}
-
+		
 		beta[et[case],] <- 1*sca[et[case]] # initialize
-		for(i in (et[case]-1):bt[case]) {
-			#beta[i,] <- (A[i,,]%*%(B[i+1,]*beta[i+1,]))*sca[i]
-			if(stationary) beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[1,,]*sca[i]
-			else beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[i,,]*sca[i]
+		
+		if(ntimes[case]>1) {
+			for(i in (et[case]-1):bt[case]) {
+				#beta[i,] <- (A[i,,]%*%(B[i+1,]*beta[i+1,]))*sca[i]
+				if(stationary) beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[1,,]*sca[i]
+				else beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[i,,]*sca[i]
+			}
+			
+			for(i in bt[case]:(et[case]-1)) {
+				#xi[i+1,,] <- (alpha[i,]%*%t(B[i+1,]*beta[i+1,]))*t(A[i,,])
+				if(stationary) xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[1,,])
+				else xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[i,,])
+			}
 		}
 		
-		for(i in bt[case]:(et[case]-1)) {
-			#xi[i+1,,] <- (alpha[i,]%*%t(B[i+1,]*beta[i+1,]))*t(A[i,,])
-			if(stationary) xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[1,,])
-			else xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[i,,])
-		}
 	}
 
 	gamma <- alpha*beta/sca