[Depmix-commits] r57 - trunk
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 7 00:39:47 CET 2008
Author: maarten
Date: 2008-03-07 00:39:46 +0100 (Fri, 07 Mar 2008)
New Revision: 57
Modified:
trunk/EM.R
trunk/depmix-test3EM.R
trunk/responses.R
Log:
- em function in EM.R now works with depmix objects
- depmix-test3EM.R is now example with depmix models, one with intercept only (works ok) and one with a covariate (doesn't work)
Modified: trunk/EM.R
===================================================================
--- trunk/EM.R 2008-03-06 17:34:46 UTC (rev 56)
+++ trunk/EM.R 2008-03-06 23:39:46 UTC (rev 57)
@@ -1,30 +1,35 @@
em <- function(object,maxit=100,tol=1e-6,verbose=FALSE,...) {
- if(!is(object,"hmModel")) stop("object must be 'hmModel'")
+ if(!is(object,"depmix")) stop("object is not of class 'depmix'")
# pseudocode
ns <- object at nstates
+ ntimes <- ntimes(object)
+ lt <- length(ntimes)
+ et <- cumsum(ntimes)
+ bt <- c(1,et[-lt]+1)
+
LL <- logLik(object)
converge <- FALSE
j <- 0
- A <- object at trans
+ A <- object at trDens
while(j <= maxit & !converge) {
for(i in 1:ns) {
- A[,,i] <- predict(object at trModels[[i]])
+ A[,,i] <- object at trDens[,,i]
}
- B <- exp(apply(object at logdens,c(1,3),sum))
+ B <- apply(object at dens,c(1,3),prod)
# TODO: add functionality for inModel
- init <- exp(logDens(object at initModel))
+ init <- object at init
# print(init)
LL.old <- LL
j <- j+1
# expectation
- fbo <- fb(init=init,A=A,B=B,ntimes=object at ntimes)
+ fbo <- fb(init=init,A=A,B=B,ntimes=ntimes(object))
LL <- fbo$logLike
# maximization
@@ -36,15 +41,15 @@
# print(fbo$gamma[1,])
# Here we need an average of gamma[bt[case],], which may need to be weighted ?? (see Rabiner, p283)
- ntimes <- object at ntimes
- lt <- length(ntimes)
- et <- cumsum(ntimes)
- bt <- c(1,et[-lt]+1)
+
# this is without weighting
initprobs <- apply(fbo$gamma[bt,],2,mean)
- object at initModel <- setpars(object at initModel,values=object at initModel@family$linkfun(initprobs,base=object at initModel@family$base))
+ # should become object at prior <- fit(object at prior)
+ object at prior@y <- fbo$gamma[bt,]
+ object at prior <- fit(object at prior,w=NULL)
+ #object at initModel <- setpars(object at initModel,values=object at initModel@family$linkfun(initprobs,base=object at initModel@family$base))
# This should become:
# lt <- length(object at ntimes)
@@ -53,33 +58,34 @@
# object at initModel@y <- fbo$gamma[bt,]
# object at initModel <- fit(object at initModel)
- et <- cumsum(object at ntimes)
+ #et <- cumsum(object at ntimes)
trm <- matrix(0,ns,ns)
for(i in 1:ns) {
- if(max(object at ntimes)>1) { #skip transition parameters update in case of latent class model
+ if(max(ntimes(object)>1)) { #skip transition parameters update in case of latent class model
if(!object at stationary) {
- object at trModels[[i]]@y <- fbo$xi[,,i]/fbo$gamma[,i]
- object at trModels[[i]] <- fit(object at trModels[[i]],w=as.matrix(fbo$gamma[,i]),ntimes=object at ntimes) # check this
+ object at transition[[i]]@y <- fbo$xi[,,i]/fbo$gamma[,i]
+ object at transition[[i]] <- fit(object at transition[[i]],w=as.matrix(fbo$gamma[,i]),ntimes=ntimes(object)) # check this
} else {
for(k in 1:ns) {
trm[i,k] <- sum(fbo$xi[-c(et),k,i])/sum(fbo$gamma[-c(et),i])
}
- object at trModels[[i]]@parameters$coefficients <- object at trModels[[i]]@family$linkfun(trm[i,],base=object at initModel@family$base)
+ # 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)
}
#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]]))
# update trans slot of the model
- object at trans[,,i] <- predict(object at trModels[[i]])
+ object at trDens[,,i] <- dens(object at transition[[i]])
}
- for(k in 1:object at nresp) {
- object at rModels[[i]][[k]] <- fit(object at rModels[[i]][[k]],w=fbo$gamma[,i])
+ for(k in 1:nresp(object)) {
+ object at response[[i]][[k]] <- fit(object at response[[i]][[k]],w=fbo$gamma[,i])
# update logdens slot of the model
- object at logdens[,k,i] <- logDens(object at rModels[[i]][[k]])
+ object at dens[,k,i] <- dens(object at response[[i]][[k]])
}
}
Modified: trunk/depmix-test3EM.R
===================================================================
--- trunk/depmix-test3EM.R 2008-03-06 17:34:46 UTC (rev 56)
+++ trunk/depmix-test3EM.R 2008-03-06 23:39:46 UTC (rev 57)
@@ -1,38 +1,26 @@
setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/rforge/depmix/trunk/")
-source("depmixS4.R")
-source("classes.R")
-source("hmModel.R")
-source("trGLM.r")
+
+source("responses.R")
source("lystig.R")
+source("depmix.R")
source("fb.r")
source("EM.R")
load("data/speed.Rda")
-rModels <- list(
- list(
- rModel(formula=rt~1,data=speed,family=gaussian(),pstart=c(5.52,.2)),
- rModel(formula=corr~1,data=speed,family=multinomial())),
- list(
- rModel(formula=rt~1,data=speed,family=gaussian(),pstart=c(6.39,.24)),
- rModel(formula=corr~1,data=speed,family=multinomial(),pstart=c(.098,.902)))
-)
-
trstart=c(0.896,0.104,0.084,0.916)
-instart=c(.5,.5)
+trstart=c(trstart[1:2],0,0.01,trstart[3:4],0,0.01)
+instart=c(0,1)
+resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
-trstart=c(trstart[1:2],0,0.01,trstart[3:4],0,0)
+# intercept only (NOTE: we should fix the transInit etc function so that
+# transition=~1 returns an x matrix of correct dimension!
+mod <- depmix(list(rt~1,corr~1),data=speed,family=list(gaussian(),multinomial()),transition=~rep(1,439)-1,trstart=c(.9,.1,.1,.9),instart=instart,respst=resp,nst=2)
+fmod.int <- em(mod)
+logLik(fmod.int)
-mod <- depmix(rModels=rModels,data=speed,transition=~Pacc,trstart=trstart,instart=instart)
-
-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)
-
-maxit=100
-tol=1e-5
-
-fmod <- em(mod,verbose=T)
\ No newline at end of file
+# now with a covariate
+mod <- depmix(list(rt~1,corr~1),data=speed,family=list(gaussian(),multinomial()),transition=~Pacc,trstart=trstart,instart=instart,respst=resp,nst=2)
+fmod <- em(mod,verbose=T)
+logLik(fmod)
Modified: trunk/responses.R
===================================================================
--- trunk/responses.R 2008-03-06 17:34:46 UTC (rev 56)
+++ trunk/responses.R 2008-03-06 23:39:46 UTC (rev 57)
@@ -149,7 +149,7 @@
setClass("transInit",contains="GLMresponse")
-setGeneric("transInit", function(formula, ...) standardGeneric("transInit"))
+setGeneric("transInit", function(formula, ... ) standardGeneric("transInit"))
# FIX ME: data is a necessary argument to determine the dimension of x, even when there
# are no covariates (and there are by definition no responses ...)
@@ -165,6 +165,10 @@
mf <- eval(mf, parent.frame())
x <- model.matrix(attr(mf, "terms"),mf)
y <- matrix(1,ncol=1) # y is not needed in the transition and init models
+
+ # use ntimes
+ #if(!is.null(ntimes)) if(nrow(x) < sum(ntimes)) x <- apply(x,2,function(y) rep(y,length=sum(ntimes)))
+
parameters <- list()
if(is.null(nstates)) stop("'nstates' must be provided in call to trinModel")
if(family$family=="multinomial") {
@@ -543,88 +547,31 @@
setMethod("fit","transInit",
function(object,w,ntimes) {
- pars <- object at parameters
- oldfit <- function() {
- #fit.trMultinom(object,w,ntimes)
- tol <- 1e-5 # TODO: check global options
- 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),]
-
- if(!is.matrix(x)) x <- matrix(x,ncol=ncol(object at 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)
-
- mt <- as.numeric(t(mu[,-base]))
- Dl <- Sigmal <- Wl <- list()
-
- converge <- FALSE
- while(!converge) {
- b.old <- b
- for(i in 1:nt) {
- Dl[[i]] <- object at family$mu.eta(mu[i,-base])
- Sigmal[[i]] <- object at family$variance(mu[i,-base])
- Wl[[i]] <- Dl[[i]]%*%solve(Sigmal[[i]])%*%t(Dl[[i]]) # TODO:
- }
- Sigma <- bdiag(Sigmal)
- D <- bdiag(Dl)
- W <- bdiag(Wl)
-
- b[,-base] <- as.numeric(b[,-base]) + solve(t(Z)%*%W%*%Z)%*%(t(Z)%*%D%*%solve(Sigma)%*%(y-mt))
- if(abs(sum(b-b.old)) < tol) converge <- TRUE
- mu <- object at family$linkinv(x%*%b,base=base)
- mt <- as.numeric(t(mu[,-base]))
- }
- pars$coefficients <- t(b) # TODO: setpars gets matrix in wrong order!!! Fix this in setpars.
- pars
- }
-
- vglmfit <- function() {
- base <- object at family$base
- w <- cbind(w[,-base],w[,base])
- x <- slot(object,"x")
- fam <- slot(object,"family")
- fit <- vglm(w~x,fam)
- pars$coefficients[,-base] <- t(slot(fit,coefficients)) # TODO: setpars gets matrix in wrong order!!! Fix this in setpars.
- pars
- }
-
- 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)))
- na <- c(na,which(is.na(w)))
- y <- as.matrix(y)
- x <- as.matrix(x)
- na <- unique(na)
- mask <- matrix(1,nrow=nrow(pars$coefficients),ncol=ncol(pars$coefficients))
- mask[,base] <- 0
- fit <- nnet.default(x=x[-na,],y=y[-na,],weights=w[-na],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)
- #object <- setpars(object,unlist(pars))
- #object
- pars
- }
- pars <- nnetfit()
+ 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]) # why do we need to transpose?
object <- setpars(object,unlist(pars))
object
+
}
)
\ No newline at end of file
More information about the depmix-commits
mailing list