[Depmix-commits] r157 - trunk/R

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

Author: maarten
Date: 2008-06-12 15:44:52 +0200 (Thu, 12 Jun 2008)
New Revision: 157

Modified:
   trunk/R/EM.R
Log:
- changed EM so that a simpler routine (not fb) is used for mix models

Modified: trunk/R/EM.R
===================================================================
--- trunk/R/EM.R	2008-06-11 21:06:46 UTC (rev 156)
+++ trunk/R/EM.R	2008-06-12 13:44:52 UTC (rev 157)
@@ -3,9 +3,91 @@
 # 
 
 em <- function(object,maxit=100,tol=1e-6,verbose=FALSE,...) {
+  if(!is(object,"mix")) stop("object is not of class '(dep)mix'")
+  call <- match.call()
+  if(is(object,"depmix")) {
+    call[[1]] <- as.name("em.depmix")
+  } else {
+    call[[1]] <- as.name("em.mix")
+  }
+  object <- eval(call, parent.frame())
+  object
+}
+
+
+em.mix <- function(object,maxit=100,tol=1e-6,verbose=FALSE,...) {
+  if(!is(object,"mix")) stop("object is not of class 'mix'")
+
+  ns <- object at nstates
+
+	ntimes <- ntimes(object)
+	lt <- length(ntimes)
+	et <- cumsum(ntimes)
+	bt <- c(1,et[-lt]+1)
+
+	converge <- FALSE
+	j <- 0
 	
-	if(!is(object,"mix")) stop("object is not of class '(dep)mix'")
+	# compute responsibilities
+  B <- apply(object at dens,c(1,3),prod)
+  gamma <- object at init*B
+  LL <- sum(log(rowSums(gamma)))
+  # normalize
+  gamma <- gamma/rowSums(gamma)
+
+	LL.old <- LL + 1
+
+	while(j <= maxit & !converge) {
+
+		# maximization
+
+		# should become object at prior <- fit(object at prior)
+		object at prior@y <- gamma[bt,,drop=FALSE]
+		object at prior <- fit(object at prior, w=NULL,ntimes=NULL)
+		object at init <- dens(object at prior)
+
+		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]])
+			}
+		}
+
+		# expectation
+    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(verbose&((j%%5)==0)) cat("iteration",j,"logLik:",LL,"\n")
+		if( (LL >= LL.old) & (LL - LL.old < tol))  {
+			cat("iteration",j,"logLik:",LL,"\n")
+			converge <- TRUE
+		}
+
+		LL.old <- LL
+		j <- j+1
+
+	}
+
+	class(object) <- "mix.fitted"
+
+	if(converge) object at message <- "Log likelihood converged to within tol."
+	else object at message <- "'maxit' iterations reached in EM without convergence."
+
+	# no constraints in EM
+	object at conMat <- matrix()
+
+	object
 	
+}
+
+em.depmix <- function(object,maxit=100,tol=1e-6,verbose=FALSE,...) {
+	
+	if(!is(object,"depmix")) stop("object is not of class '(dep)mix'")
+	
 	ns <- object at nstates
 	
 	ntimes <- ntimes(object)
@@ -76,9 +158,11 @@
 		
 	}
 	
-	if(class(object)=="depmix") class(object) <- "depmix.fitted"
-	if(class(object)=="mix") class(object) <- "mix.fitted"
+	#if(class(object)=="depmix") class(object) <- "depmix.fitted"
+	#if(class(object)=="mix") class(object) <- "mix.fitted"
 	
+	class(object) <- "depmix.fitted"
+	
 	if(converge) object at message <- "Log likelihood converged to within tol."
 	else object at message <- "'maxit' iterations reached in EM without convergence."