[Depmix-commits] r132 - trunk/R

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

Author: maarten
Date: 2008-05-21 02:04:29 +0200 (Wed, 21 May 2008)
New Revision: 132

Modified:
   trunk/R/EM.R
   trunk/R/fb.R
   trunk/R/mlogit.R
   trunk/R/responseGLMMULTINOM.R
   trunk/R/transInit.R
Log:
- fixed EM, now works on Test 4 in depmixNew-test1.R :-)

Modified: trunk/R/EM.R
===================================================================
--- trunk/R/EM.R	2008-05-20 16:41:14 UTC (rev 131)
+++ trunk/R/EM.R	2008-05-21 00:04:29 UTC (rev 132)
@@ -16,12 +16,12 @@
 	converge <- FALSE
 	j <- 0
 	
-	A <- object at trDens
-	B <- apply(object at dens,c(1,3),prod)
-	init <- object at init
+	# A <- object at trDens
+	# B <- object at dens
+	# init <- object at init
 	
 	# initial expectation
-	fbo <- fb(init=object at init,A=object at trDens,B=apply(object at dens,c(1,3),prod),ntimes=ntimes(object))
+	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
 	
@@ -52,6 +52,9 @@
 				# update trDens slot of the model
 				object at trDens[,,i] <- dens(object at transition[[i]])
 			}
+		}
+		
+		for(i in 1:ns) {
 			
 			for(k in 1:nresp(object)) {
 				object at response[[i]][[k]] <- fit(object at response[[i]][[k]],w=fbo$gamma[,i])
@@ -61,7 +64,7 @@
 		}
 		
 		# expectation
-		fbo <- fb(init=object at init,A=object at trDens,B=apply(object at dens,c(1,3),prod),ntimes=ntimes(object))
+		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(verbose&((j%%5)==0)) cat("iteration",j,"logLik:",LL,"\n")

Modified: trunk/R/fb.R
===================================================================
--- trunk/R/fb.R	2008-05-20 16:41:14 UTC (rev 131)
+++ trunk/R/fb.R	2008-05-21 00:04:29 UTC (rev 132)
@@ -19,7 +19,7 @@
 	# NOTE: xi[t,i,j] = P(S[t] = j & S[t+1] = i)
 	
 	
-	if(length(dim(B)) == 3) B <- apply(B,c(1,3),prod)
+	B <- apply(B,c(1,3),prod)
 	
 	nt <- nrow(B)	
 	ns <- ncol(init)
Modified: trunk/R/mlogit.R
===================================================================
--- trunk/R/mlogit.R	2008-05-20 16:41:14 UTC (rev 131)
+++ trunk/R/mlogit.R	2008-05-21 00:04:29 UTC (rev 132)
@@ -24,8 +24,9 @@
 	linkinv <- function(eta,base) {
 		linv <- function(eta,base) {
 			pp <- numeric(length(eta))
-			if(any(is.infinite(eta))) {
+			if(any(is.infinite(eta)) || any(eta > log(.Machine$double.xmax)) || any(eta < log(.Machine$double.xmin))) {
 				pp[which(is.infinite(eta))] <- 1
+				pp[which(eta > log(.Machine$double.xmax))] <- 1 # change this to something better!
 			} else {
 				expb <- exp(eta)
 				sumb <- sum(expb)

Modified: trunk/R/responseGLMMULTINOM.R
===================================================================
--- trunk/R/responseGLMMULTINOM.R	2008-05-20 16:41:14 UTC (rev 131)
+++ trunk/R/responseGLMMULTINOM.R	2008-05-21 00:04:29 UTC (rev 132)
@@ -5,12 +5,15 @@
 	function(object,w) {
 		pars <- object at parameters
 		base <- object at family$base # delete me
-		y <- object at y[,-base]
+		#y <- object at y[,-base]
 		x <- object at x
 		mask <- matrix(1,nrow=nrow(pars$coefficients),ncol=ncol(pars$coefficients))
-		mask[,base] <- 0
-		fit <- nnet.default(x,y,weights=w,size=0,entropy=TRUE,skip=TRUE,mask=mask,rang=0,trace=FALSE)
-		pars$coefficients <- matrix(fit$wts,ncol=ncol(pars$coefficients),nrow=nrow(pars$coefficients),byrow=TRUE)
+		mask[,base] <- 0 # fix base category coefficients to 0
+		mask <- rbind(0,mask) # fix "bias" nodes to 0
+		Wts <- mask
+		Wts[-1,] <- t(pars$coefficients) # set starting weights
+		fit <- nnet.default(x,y,weights=w,size=0,entropy=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+		pars$coefficients <- matrix(fit$wts,ncol=ncol(pars$coefficients),nrow=nrow(pars$coefficients)+1,byrow=TRUE)[-1,]
 		object <- setpars(object,unlist(pars))
 		object
 	}
Modified: trunk/R/transInit.R
===================================================================
--- trunk/R/transInit.R	2008-05-20 16:41:14 UTC (rev 131)
+++ trunk/R/transInit.R	2008-05-21 00:04:29 UTC (rev 132)
@@ -175,13 +175,20 @@
 			#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,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?
+
+		mask <- matrix(1,nrow=nrow(pars$coefficients),ncol=ncol(pars$coefficients))
+		mask[,base] <- 0 # fix base category coefficients to 0
+		mask <- rbind(0,mask) # fix "bias" nodes to 0
+		Wts <- mask
+		Wts[-1,] <- pars$coefficients # set starting weights
+		Wts[Wts == Inf] <- .Machine$double.max.exp # Fix this!!!!
+		Wts[Wts == -Inf] <- .Machine$double.min.exp # Fix this!!!!!
+		if(!is.null(w)) {
+      fit <- nnet.default(x,y,weights=w,size=0,entropy=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+    } else {
+      fit <- nnet.default(x,y,size=0,entropy=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+    }
+		pars$coefficients <- t(matrix(fit$wts,ncol=ncol(pars$coefficients),nrow=nrow(pars$coefficients)+1)[-1,])
 		object <- setpars(object,unlist(pars))
 		object
 	}