[Depmix-commits] r373 - trunk/R

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

Author: maarten
Date: 2010-03-01 16:52:44 +0100 (Mon, 01 Mar 2010)
New Revision: 373

Modified:
   trunk/R/EM.R
   trunk/R/makeTransModels.R
   trunk/R/transInit.R
Log:
- added support for identity link in transInit models
- identity link is now default for stationary models in makeTransInit

Modified: trunk/R/EM.R
===================================================================
--- trunk/R/EM.R	2010-02-25 21:55:09 UTC (rev 372)
+++ trunk/R/EM.R	2010-03-01 15:52:44 UTC (rev 373)
@@ -144,8 +144,11 @@
 					for(k in 1:ns) {
 						trm[i,k] <- sum(fbo$xi[-c(et),k,i])/sum(fbo$gamma[-c(et),i])
 					}
-					# FIX THIS; it will only work with a specific trinModel
-					object at transition[[i]]@parameters$coefficients <- object at transition[[i]]@family$linkfun(trm[i,],base=object at transition[[i]]@family$base)
+					# FIX THIS; it will only work with specific trinModels??
+					object at transition[[i]]@parameters$coefficients <- switch(object at transition[[i]]@family$link,
+					identity = object at transition[[i]]@family$linkfun(trm[i,]),
+					mlogit = object at transition[[i]]@family$linkfun(trm[i,],base=object at transition[[i]]@family$base),
+					object at transition[[i]]@family$linkfun(trm[i,]))
 				}
 				# update trDens slot of the model
 				object at trDens[,,i] <- dens(object at transition[[i]])

Modified: trunk/R/makeTransModels.R
===================================================================
--- trunk/R/makeTransModels.R	2010-02-25 21:55:09 UTC (rev 372)
+++ trunk/R/makeTransModels.R	2010-03-01 15:52:44 UTC (rev 373)
@@ -19,10 +19,10 @@
 	for(i in 1:nstates) {
 		if(tst) {
 			# in the case of stationary models we do not need data?!?!?!?!?!
-			if(stationary) models[[i]] <- transInit(formula,multinomial(base=base),data=data.frame(1),nstates=nstates,pstart=values[i,],prob=prob)
+			if(stationary) models[[i]] <- transInit(formula,multinomial(link="identity"),data=data.frame(1),nstates=nstates,pstart=values[i,],prob=prob)
 			else models[[i]] <- transInit(formula,multinomial(base=base),data=data,nstates=nstates,pstart=values[i,],prob=prob)
 		} else {
-			if(stationary) models[[i]] <- transInit(formula,multinomial(base=base),data=data.frame(1),nstates=nstates,prob=FALSE)
+			if(stationary) models[[i]] <- transInit(formula,multinomial(link="identity"),data=data.frame(1),nstates=nstates,prob=FALSE)
 			else models[[i]] <- transInit(formula,multinomial(base=base),data=data,nstates=nstates,prob=FALSE)
 		}
 	}

Modified: trunk/R/transInit.R
===================================================================
--- trunk/R/transInit.R	2010-02-25 21:55:09 UTC (rev 372)
+++ trunk/R/transInit.R	2010-03-01 15:52:44 UTC (rev 373)
@@ -106,29 +106,48 @@
 			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 # 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)) {
-			if(NCOL(y) < 3) {
-				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,weights=w,size=0,softmax=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
-			}
-		} else {
-			if(NCOL(y) < 3) {
-				fit <- nnet.default(x,y,size=0,entropy=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
-			} else {
-				fit <- nnet.default(x,y,size=0,softmax=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))
+		}
+		switch(object at family$link,
+		  mlogit = {
+    		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)) {
+    			if(NCOL(y) < 3) {
+    				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,weights=w,size=0,softmax=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+    			}
+    		} else {
+    			if(NCOL(y) < 3) {
+    				fit <- nnet.default(x,y,size=0,entropy=TRUE,skip=TRUE,mask=mask,Wts=Wts,trace=FALSE)
+    			} else {
+    				fit <- nnet.default(x,y,size=0,softmax=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))
+  		},
+  		identity = {
+  		  # object at y = fbo$xi[,,i]/fbo$gamma[,i]
+  		  # should become (see em):
+  		  #for(k in 1:ns) {
+				#		trm[i,k] <- sum(fbo$xi[-c(et),k,i])/sum(fbo$gamma[-c(et),i])
+				#	}
+				if(!is.null(w)) {
+				  sw <- sum(w)
+				  pars <- colSums(w*object at y)/sum(w)
+				} else {
+				  pars <- colMeans(object at y)
+				}
+			  object <- setpars(object,pars)
+  		},
+  		stop("link function not implemented")
+  	)
 		object
 	}
 )
@@ -156,4 +175,4 @@
 			return(states)
 		}
 	}
-)
\ No newline at end of file
+)