[Depmix-commits] r449 - in pkg/depmixS4: . R

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

Author: ingmarvisser
Date: 2011-01-25 15:37:21 +0100 (Tue, 25 Jan 2011)
New Revision: 449

Modified:
   pkg/depmixS4/NEWS
   pkg/depmixS4/R/EM.R
   pkg/depmixS4/R/depmixfit.R
   pkg/depmixS4/R/fb.R
   pkg/depmixS4/R/lystig.R
Log:
Added a meaningful error message in the EM algorithm for lca/mixture 
    models in case the initial log likelihood is NA (thanks to Matthias
    Ihrke for pointing this out).

Modified: pkg/depmixS4/NEWS
===================================================================
--- pkg/depmixS4/NEWS	2010-11-23 13:13:55 UTC (rev 448)
+++ pkg/depmixS4/NEWS	2011-01-25 14:37:21 UTC (rev 449)
@@ -1,4 +1,9 @@
+Changes in depmixS4 version ???
 
+  o added a meaningful error message in the EM algorithm for lca/mixture 
+    models in case the initial log likelihood is NA (thanks to Matthias
+    Ihrke for pointing this out). 
+
 Changes in depmixS4 version 1.0-1
 
   o minor changes in documentation to conform to R 2.12.0 standards. 

Modified: pkg/depmixS4/R/EM.R
===================================================================
--- pkg/depmixS4/R/EM.R	2010-11-23 13:13:55 UTC (rev 448)
+++ pkg/depmixS4/R/EM.R	2011-01-25 14:37:21 UTC (rev 449)
@@ -28,20 +28,36 @@
 	converge <- FALSE
 	j <- 0
 	
-	# compute response probabilities
-	B <- apply(object at dens,c(1,3),prod)
-	gamma <- object at init*B
-	LL <- sum(log(rowSums(gamma)))
-	# normalize
-	gamma <- gamma/rowSums(gamma)	
-	
 	if(random.start) {
+				
 		nr <- sum(ntimes(object))
 		gamma <- matrix(runif(nr*ns,min=.0001,max=.9999),nr=nr,nc=ns)
 		gamma <- gamma/rowSums(gamma)
-		# add stuff to reestimate the model here and compute LL again
-		# based on these starting values
-	} 
+		LL <- -1e10
+		
+		for(i in 1:ns) {
+			for(k in 1:nresp(object)) {
+				object at response[[i]][[k]] <- fit(object at response[[i]][[k]],w=gamma[,i])
+				# update dens slot of the model
+				object at dens[,k,i] <- dens(object at response[[i]][[k]])
+			}
+		}
+		
+		# initial expectation
+		fbo <- fb(init=object at init,matrix(0,1,1),B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
+		LL <- fbo$logLike
+		
+		if(is.nan(LL)) stop("Cannot find suitable starting values; please provide them.")
+		
+	} else {
+		# initial expectation
+# 		fbo <- fb(init=object at init,matrix(0,1,1),B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
+		B <- apply(object at dens,c(1,3),prod)
+		gamma <- object at init*B
+		LL <- sum(log(rowSums(gamma)))
+		if(is.nan(LL)) stop("Starting values not feasible; please provide them.")
+		gamma <- gamma/rowSums(gamma)
+	}
 	
 	LL.old <- LL + 1
 	
@@ -66,6 +82,7 @@
 		B <- apply(object at dens,c(1,3),prod)
 		gamma <- object at init*B
 		LL <- sum(log(rowSums(gamma)))
+
 		# normalize
 		gamma <- gamma/rowSums(gamma)
 		
@@ -124,7 +141,6 @@
 	converge <- FALSE
 	j <- 0
 	
-	
 	if(random.start) {
 				
 		nr <- sum(ntimes(object))
@@ -150,6 +166,7 @@
 		# initial expectation
 		fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
 		LL <- fbo$logLike
+		if(is.nan(LL)) stop("Starting values not feasible; please provide them.")
 	}
 	
 	LL.old <- LL + 1

Modified: pkg/depmixS4/R/depmixfit.R
===================================================================
--- pkg/depmixS4/R/depmixfit.R	2010-11-23 13:13:55 UTC (rev 448)
+++ pkg/depmixS4/R/depmixfit.R	2011-01-25 14:37:21 UTC (rev 449)
@@ -109,7 +109,7 @@
 				allpars[!fixed] <- pars
 				object <- setpars(object,allpars)
 				ans = -as.numeric(logLik(object))
-				if(is.na(ans)) ans = 100000
+				if(is.na(ans)) ans = 100000 # remove magic number here
 				ans
 			}
 			

Modified: pkg/depmixS4/R/fb.R
===================================================================
--- pkg/depmixS4/R/fb.R	2010-11-23 13:13:55 UTC (rev 448)
+++ pkg/depmixS4/R/fb.R	2011-01-25 14:37:21 UTC (rev 449)
@@ -8,17 +8,17 @@
 
 	# Forward-Backward algorithm (used in Baum-Welch)
 	# Returns alpha, beta, and full data likelihood
+	
 	# NOTE THE CHANGE IN FROM ROW TO COLUMN SUCH THAT TRANSPOSING A IS NOT NECCESSARY ANYMORE
 	# IN COMPUTING ALPHA AND BETA BUT IS NOW NECCESSARY IN COMPUTING XI
 	# A = T*K*K matrix with transition probabilities, from row to column!!!!!!!
 	# B = T*K matrix with elements ab_{ij} = P(y_i|s_j)
 	# init = K vector with initial probabilities
 
-	# NOTE: to prevent underflow, alpha and beta are scaled, using c
+	# NOTE: to prevent underflow, alpha and beta are scaled, using sca
 	
-	# NOTE: xi[t,i,j] = P(S[t] = j & S[t+1] = i)
+	# NOTE: xi[t,i,j] = P(S[t] = j & S[t+1] = i) !!!NOTE the order of i and j!!!
 	
-	
 	B <- apply(B,c(1,3),prod)
 	
 	nt <- nrow(B)	
Modified: pkg/depmixS4/R/lystig.R
===================================================================
--- pkg/depmixS4/R/lystig.R	2010-11-23 13:13:55 UTC (rev 448)
+++ pkg/depmixS4/R/lystig.R	2011-01-25 14:37:21 UTC (rev 449)
@@ -21,6 +21,7 @@
 
 	# B = T*K*nresp matrix with elements ab_{tij} = P(y_t_i|s_j)
 	# init = K vector with initial probabilities !!!OR!!! K*length(ntimes) matrix with initial probs per case
+	
 	# Returns: 'sca'le factors, recurrent variables 'phi', loglikelihood
 		
 	nt <- nrow(B)