[Depmix-commits] r471 - in pkg/depmixS4: R src

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

Author: ingmarvisser
Date: 2011-06-21 22:50:33 +0200 (Tue, 21 Jun 2011)
New Revision: 471

Modified:
   pkg/depmixS4/R/fb.R
   pkg/depmixS4/src/fb.cc
Log:
C-code for backward variables now correct.

Modified: pkg/depmixS4/R/fb.R
===================================================================
--- pkg/depmixS4/R/fb.R	2011-06-21 13:29:11 UTC (rev 470)
+++ pkg/depmixS4/R/fb.R	2011-06-21 20:50:33 UTC (rev 471)
@@ -82,10 +82,18 @@
 			
 			beta[et[case],] <- 1*sca[et[case]] # initialize
 			
+# 			print(et[case])
+# 			print(round(beta[et[case],],6))
+			
 			if(ntimes[case]>1) {
 				for(i in (et[case]-1):bt[case]) {
+# 					cat("t: ",i,"\n")
 					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]
+# 					print(round(sca[i],6))
+# 					print(round(beta[i,],6))
+#  					print(round(B[i+1,],6))
+#  					print(round(A[i,,],6))
 				}
 				
 				for(i in bt[case]:(et[case]-1)) {
Modified: pkg/depmixS4/src/fb.cc
===================================================================
--- pkg/depmixS4/src/fb.cc	2011-06-21 13:29:11 UTC (rev 470)
+++ pkg/depmixS4/src/fb.cc	2011-06-21 20:50:33 UTC (rev 471)
@@ -92,12 +92,47 @@
 		// compute backward variables and xi
 		matrix betatp1(ns[0]);
 		matrix betat(ns[0]);
+// 		Rprintf("et[cas]: %d \n",et[cas]);		
+// 		Rprintf("sca et[cas]: %f \n",sca[et[cas]-1]);
 		// compute initial beta, ie for t=T (for each case)
 		for(int i=0; i<ns[0]; i++) {
-			betatp1[et[cas],i] <- sca[et[cas]]
+			betatp1(i+1) = sca[et[cas]-1];
 		}
-		betatp1.print();
+// 		betatp1.print();
+		for(int i=0; i<ns[0]; i++) {
+			beta[(et[cas]-1)*ns[0]+i] = betatp1(i+1);
+		}
+ 		if(ntimes[cas]>1) {
+			// loop from T-1 to 1 for each case
+ 			for(int t=(et[cas]-1); t>=bt[cas]; t--) {
+// 				Rprintf("t: %d \n", t);
+				for(int i=0; i<ns[0]; i++) {
+					for(int j=0; j<ns[0]; j++) {
+						trans(i+1,j+1) = trdens[(i*ns[0]+j)*nt[0]+t-1]; // A_t
+					}
+  					denst(i+1) = dens[t*ns[0]+i]; //B_t+1
+//  				if(stationary) beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[1,,]*sca[i]
+				}
+//  				transpose(trans).print();
+//  				denst.print();
+// 				Rprintf("sca t: %f \n",sca[t-1]);
+				betat = trans*had(denst,betatp1)*sca[t-1];
+// 				betat.print();
+				// store betat somewhere
+				for(int i=0; i<ns[0]; i++) {
+					beta[(t-1)*ns[0]+i] = betat(i+1);
+				}
+				betatp1 = betat;
+ 			}
+ 			
+ 			
+//  			for(i in bt[case]:(et[case]-1)) {
+//  				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,,])
+// 			}
+ 		}
 		
+		
 // 		R-code for the backwards/xi loop
 // 		if(ntimes[case]>1) {
 // 			for(i in (et[case]-1):bt[case]) {