[Depmix-commits] r44 - trunk

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

Author: maarten
Date: 2008-03-05 01:32:03 +0100 (Wed, 05 Mar 2008)
New Revision: 44

Added:
   trunk/viterbi.r
Modified:
   trunk/EM.R
   trunk/classes.R
   trunk/depmix-test3EM.R
Log:
- fixed fit(trinModel), so EM works, but likelihood still goes down!
- added viterbi algorithm (will still need ntimes etc)

Modified: trunk/EM.R
===================================================================
--- trunk/EM.R	2008-03-05 00:31:15 UTC (rev 43)
+++ trunk/EM.R	2008-03-05 00:32:03 UTC (rev 44)
@@ -44,6 +44,7 @@
       #object at trModels[[i]] <- fit(object at trModels[[i]],w=NULL,ntimes=object at ntimes) # check this
       #object at trans[,,i] <- exp(logDens(object at trModels[[i]]))
       object at trans[,,i] <- predict(object at trModels[[i]])
+      
       for(k in 1:object at nresp) {
         object at rModels[[i]][[k]] <- fit(object at rModels[[i]][[k]],w=fbo$gamma[,i])
         object at logdens[,k,i] <- logDens(object at rModels[[i]][[k]])

Modified: trunk/classes.R
===================================================================
--- trunk/classes.R	2008-03-05 00:31:15 UTC (rev 43)
+++ trunk/classes.R	2008-03-05 00:32:03 UTC (rev 44)
@@ -380,27 +380,18 @@
 			pars <- object at parameters
 			b <- pars$coefficients
 			base <- object at family$base
-			
 			if(is.matrix(w)) nan <- which(is.na(rowSums(w))) else nan <- which(is.na(w))
-			
 			#vgam(cbind(w[,-base],w[,base]) ~ ) # what is this?
-			
 			y <- as.vector(t(object at family$linkinv(w[-c(nan,ntimes),-base],base=object at family$base)))
-			
-			x <- object at x[-c(nan,ntimes),]
-			
+			x <- object at x[-c(nan,ntimes),]			
 			if(!is.matrix(x)) x <- matrix(x,ncol=ncol(object at x))
-			nt <- nrow(x)
-			
+			nt <- nrow(x)			
 			Z <- matrix(ncol=length(b))
 			Z <- vector()
-			for(i in 1:nt) Z <- rbind(Z,t(bdiag(rep(list(x[i,]),ncol(w)-1))))
-			
-			mu <- object at family$linkinv(x%*%b,base=base)
-			
+			for(i in 1:nt) Z <- rbind(Z,t(bdiag(rep(list(x[i,]),ncol(w)-1))))			
+			mu <- object at family$linkinv(x%*%b,base=base)			
 			mt <- as.numeric(t(mu[,-base]))
-			Dl <- Sigmal <- Wl <- list()
-			
+			Dl <- Sigmal <- Wl <- list()			
 			converge <- FALSE
 			while(!converge) {
 				b.old <- b
@@ -433,6 +424,33 @@
     }
     
     nnetfit <- function() {
+  		require(nnet)
+  		pars <- object at parameters
+  		base <- object at family$base # delete me
+  		#y <- object at y[,-base]
+  		y <- object at y
+  		x <- object at x
+  		if(is.matrix(y)) na <- unlist(apply(y,2,function(x) which(is.na(x)))) else na <- which(is.na(y))
+  		if(is.matrix(x)) na <- c(na,unlist(apply(x,2,function(x) which(is.na(x))))) else na <- c(na,which(is.na(x)))
+  		if(!is.null(w)) na <- c(na,which(is.na(w)))
+  		y <- as.matrix(y)
+  		x <- as.matrix(x)
+  		na <- unique(na)
+  		x <- x[-na,]
+  		y <- y[-na,]
+  		y <- round(y) # delete me
+  		if(!is.null(w)) w <- w[-na]
+  		#mask <- matrix(1,nrow=nrow(pars$coefficients),ncol=ncol(pars$coefficients))
+  		#mask[,base] <- 0
+  		if(!is.null(w)) fit <- multinom(y~x-1,weights=w,trace=FALSE) else fit <- multinom(y~x-1,weights=w,trace=FALSE)
+  		ids <- vector(,length=ncol(y))
+  		ids[base] <- 1
+  		ids[-base] <- 2:ncol(y)
+  		pars$coefficients <- t(matrix(fit$wts,ncol=ncol(y))[-1,ids])
+  		object <- setpars(object,unlist(pars))
+  		#object
+  		pars
+    }
 		require(nnet)
 		pars <- object at parameters
 		base <- object at family$base # delete me
@@ -447,16 +465,15 @@
 		na <- unique(na)
 		x <- x[-na,]
 		y <- y[-na,]
+		#y <- round(y) # delete me
 		if(!is.null(w)) w <- w[-na]
 		#mask <- matrix(1,nrow=nrow(pars$coefficients),ncol=ncol(pars$coefficients))
 		#mask[,base] <- 0
-		if(!is.null(w)) fit <- multinom(y~x-1,weights=w) else fit <- multinom(y~x-1,weights=w)
-		pars$coefficients <- matrix(fit$wts,ncol=ncol(y))[-1,-1]
-		#object <- setpars(object,unlist(pars))
-		#object
-		pars
-    }
-    pars <- nnetfit()
+		if(!is.null(w)) fit <- multinom(y~x-1,weights=w,trace=FALSE) else fit <- multinom(y~x-1,weights=w,trace=FALSE)
+		ids <- vector(,length=ncol(y))
+		ids[base] <- 1
+		ids[-base] <- 2:ncol(y)
+		pars$coefficients <- t(matrix(fit$wts,ncol=ncol(y))[-1,ids]) # why do we need to transpose?
 		object <- setpars(object,unlist(pars))
 		object
 	}

Modified: trunk/depmix-test3EM.R
===================================================================
--- trunk/depmix-test3EM.R	2008-03-05 00:31:15 UTC (rev 43)
+++ trunk/depmix-test3EM.R	2008-03-05 00:32:03 UTC (rev 44)
@@ -28,11 +28,10 @@
 
 logLik(mod)
 
-mod at trModels[[1]] <- mod at trModels[[2]] <- trGLM(~Pacc,data=speed,nstate=2)
-mod at trModels[[1]]@parameters$coefficients[,2] <- c(-2.153550,.01)
-mod at trModels[[2]]@parameters$coefficients[,2] <- c(2.389200,0)
+#mod at trModels[[1]] <- mod at trModels[[2]] <- trGLM(~Pacc,data=speed,nstate=2)
+#mod at trModels[[1]]@parameters$coefficients[,2] <- c(-2.153550,.01)
+#mod at trModels[[2]]@parameters$coefficients[,2] <- c(2.389200,0)
 
-object <- mod
 maxit=100
 tol=1e-5
 

Added: trunk/viterbi.r
===================================================================
--- trunk/viterbi.r	                        (rev 0)
+++ trunk/viterbi.r	2008-03-05 00:32:03 UTC (rev 44)
@@ -0,0 +1,25 @@
+viterbi <- function(A,B,prior) {
+    # returns the most likely state sequence
+    nt <- nrow(B)
+    ns <- ncol(A)
+    delta <- psi <- matrix(nrow=nt,ncol=ns)
+    state <- vector(length=nt)
+    # initialization
+    delta[1,] <- - (log(prior) + log(B[1,]))
+    psi[1,] <- 0
+    # recursion
+    for(i in 2:nt) {
+        for(j in 1:ns) {
+            delta[i,j] <- min(delta[i-1,] - log(A[,j])) - log(B[i,j])
+            k <- which.min(delta[i-1,] - log(A[,j]))
+            if(length(k) == 0) k <- 0
+            psi[i,j] <- k
+        }
+    }
+    #trace maximum likely state
+    state[nt] <- which.min(delta[nt,])
+    for(i in (nt-1):1) {
+        state[i] <- psi[i+1,state[i+1]]
+    }
+    return(state)
+}
\ No newline at end of file