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

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Nov 11 16:12:13 CET 2010


Author: ingmarvisser
Date: 2010-11-11 16:12:13 +0100 (Thu, 11 Nov 2010)
New Revision: 445

Modified:
   pkg/depmixS4/R/EM.R
Log:
Added a check for bad likelihoods when generating starting values.

Modified: pkg/depmixS4/R/EM.R
===================================================================
--- pkg/depmixS4/R/EM.R	2010-10-20 11:01:22 UTC (rev 444)
+++ pkg/depmixS4/R/EM.R	2010-11-11 15:12:13 UTC (rev 445)
@@ -39,6 +39,8 @@
 		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.old <- LL + 1
@@ -122,17 +124,36 @@
 	converge <- FALSE
 	j <- 0
 	
-	# 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
-	LL.old <- LL + 1
 	
 	if(random.start) {
+				
 		nr <- sum(ntimes(object))
-		fbo$gamma <- matrix(runif(nr*ns,min=.0001,max=.9999),nr=nr,nc=ns)
-		fbo$gamma <- fbo$gamma/rowSums(fbo$gamma)
+		gamma <- matrix(runif(nr*ns,min=.0001,max=.9999),nr=nr,nc=ns)
+		gamma <- gamma/rowSums(gamma)
+		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,A=object at trDens,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,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
+		LL <- fbo$logLike
 	}
 	
+	LL.old <- LL + 1
+	
 	while(j <= maxit & !converge) {
 		
 		# maximization



More information about the depmix-commits mailing list