[Depmix-commits] r55 - / trunk
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 6 18:32:42 CET 2008
Author: maarten
Date: 2008-03-06 18:32:42 +0100 (Thu, 06 Mar 2008)
New Revision: 55
Added:
trunk/EM3.R
trunk/depmix-test3EM3.R
Modified:
todo
trunk/viterbi.r
Log:
- viterbi now handles hmmModel
- added EM3.R and depmix-test3EM3.R as example of EM with viterbi states for maximisation of trinModel...doesn't work either :(
Modified: todo
===================================================================
--- todo 2008-03-06 15:44:05 UTC (rev 54)
+++ todo 2008-03-06 17:32:42 UTC (rev 55)
@@ -1,85 +1,86 @@
-
-TODO list for depmix package release 0.0-1
-
-(priority to be determined)
-
-- add help files for all functions:
- 1) depmix.Rd: help on constructing models: done!!!!!
- 2) depmix-class.Rd: detailed description of depmix class: Ingmar, done!!!!!
- 3) depmix.fit.Rd: help on fitting models: Ingmar
- Other methods for depmix.fitted objects include:
- a) posterior: get posterior state sequence: we still need the code for this!
- b) logLik: done!!!
- c) AIC, BIC: done!!!
- d) getpars, setpars: done!!!!
- e) llratio: done!!!!!
-
- 4) responses.Rd: detailed description on constructing response models
- 5) transInit.Rd: detailed description on constructing transition models
- 6) others: lystig, fb
-
-- rework em function and add a working example for it: Maarten
-
-- balance scale data example with covariate on initial probabilities
-(this involves a latent class model rather than a hidden Markov model)
-Ingmar
-
-- review summary/print functions for depmix and depmix.fitted
-
-
-- viterbi algorithm to get posterior probabilities as part of depmix.fitted
-
-
-- general documentation: depmix-introduction.pdf needs updating
-
-
-- generate initial values?!?!?
-
-
-
-
-
-
-
-
-
-TODO long term
-
-Examples
-
-1) Discrimination learning example, need covariates with these data!
-
-2) Factor depmix model as example of user specified model
-
-
-Design issues
-
-1) Data simulation according to specified model
-
-2) Analytical gradients in lystig
-
-3) Speed optimization: lystig in C?
-
-4) Multi group possibilities: use group factor in call to depmix
-
-5) Missing data options?
-
-6) Standard errors of parameters
-
-
-Other capabilities
-
-1) Look at log and exp of matrices to have continuous time observations
-(also see mlm)
-
-2) Look at HMISC for mapply -> multivariate apply for mv densities
-
-3) Look at package ks for mv normal mixture densities
-
-4) Check flexclust package for an example of extending models
-
-5) Check glmc pack for constrained glm models!!!!!
-
-6) stochmod pack for stochastic models
-
-7) use relevel function to recode factors into different base
+
+TODO list for depmix package release 0.0-1
+
+(priority to be determined)
+
+- add help files for all functions:
+ 1) depmix.Rd: help on constructing models: done!!!!!
+ 2) depmix-class.Rd: detailed description of depmix class: Ingmar, done!!!!!
+ 3) depmix.fit.Rd: help on fitting models: Ingmar
+ Other methods for depmix.fitted objects include:
+ a) posterior: get posterior state sequence: we still need the code for this!
+ b) logLik: done!!!
+ c) AIC, BIC: done!!!
+ d) getpars, setpars: done!!!!
+ e) llratio: done!!!!!
+
+ 4) responses.Rd: detailed description on constructing response models
+ 5) transInit.Rd: detailed description on constructing transition models
+ 6) others: lystig, fb
+
+- rework em function and add a working example for it: Maarten
+
+- balance scale data example with covariate on initial probabilities
+(this involves a latent class model rather than a hidden Markov model)
+Ingmar
+
+- review summary/print functions for depmix and depmix.fitted
+
+- viterbi algorithm to get maximum a posteriori states for hmmModel: Maarten, done !!!!!
+
+- viterbi algorithm to get posterior probabilities as part of depmix.fitted
+
+
+- general documentation: depmix-introduction.pdf needs updating
+
+
+- generate initial values?!?!?
+
+
+
+
+
+
+
+
+
+TODO long term
+
+Examples
+
+1) Discrimination learning example, need covariates with these data!
+
+2) Factor depmix model as example of user specified model
+
+
+Design issues
+
+1) Data simulation according to specified model
+
+2) Analytical gradients in lystig
+
+3) Speed optimization: lystig in C?
+
+4) Multi group possibilities: use group factor in call to depmix
+
+5) Missing data options?
+
+6) Standard errors of parameters
+
+
+Other capabilities
+
+1) Look at log and exp of matrices to have continuous time observations
+(also see mlm)
+
+2) Look at HMISC for mapply -> multivariate apply for mv densities
+
+3) Look at package ks for mv normal mixture densities
+
+4) Check flexclust package for an example of extending models
+
+5) Check glmc pack for constrained glm models!!!!!
+
+6) stochmod pack for stochastic models
+
+7) use relevel function to recode factors into different base
Added: trunk/EM3.R
===================================================================
--- trunk/EM3.R (rev 0)
+++ trunk/EM3.R 2008-03-06 17:32:42 UTC (rev 55)
@@ -0,0 +1,105 @@
+em <- function(object,maxit=100,tol=1e-6,verbose=FALSE,...) {
+ if(!is(object,"hmModel")) stop("object must be 'hmModel'")
+
+ # pseudocode
+ ns <- object at nstates
+
+ LL <- logLik(object)
+
+ converge <- FALSE
+ j <- 0
+
+ A <- object at trans
+ while(j <= maxit & !converge) {
+
+ for(i in 1:ns) {
+ A[,,i] <- predict(object at trModels[[i]])
+ }
+
+ B <- exp(apply(object at logdens,c(1,3),sum))
+ # TODO: add functionality for inModel
+ init <- exp(logDens(object at initModel))
+# print(init)
+ LL.old <- LL
+ j <- j+1
+
+ # expectation
+ fbo <- fb(init=init,A=A,B=B,ntimes=object at ntimes)
+ LL <- fbo$logLike
+
+ # maximization
+
+ #object at init <- fit(object at init,ip=fbo$gamma[1,])
+ #object at init <- matrix(fbo$gamma[1,],nrow=1)
+
+ # FIX ME for length(ntimes)>1
+ # 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))
+
+ # This should become:
+ # lt <- length(object at ntimes)
+ # et <- cumsum(object at ntimes)
+ # bt <- c(1,et[-lt]+1)
+ # object at initModel@y <- fbo$gamma[bt,]
+ # object at initModel <- fit(object at initModel)
+
+ et <- cumsum(object at ntimes)
+
+ trm <- matrix(0,ns,ns)
+
+ vit <- viterbi(object)
+ emp_trans <- y <- matrix(NA,nrow=sum(ntimes),ncol=2)
+ for(case in 1:lt) {
+ emp_trans[bt[case]:(et[case]-1),] <- cbind(vit[bt[case]:(et[case]-1)],vit[(bt[case]+1):(et[case])])
+ }
+
+ for(i in 1:ns) {
+
+ if(max(object at ntimes)>1) { #skip transition parameters update in case of latent class model
+ if(!object at stationary) {
+ # y is now indicator matrix with 1 in column j if transition to state j (from any state) and 0 otherwise
+ object at trModels[[i]]@y <- matrix(as.numeric(unlist(lapply(as.list(1:ns),function(x) emp_trans[,2] == x))),ncol=ns)
+ # w below is indicator vector with 1 if transition was from state i
+ object at trModels[[i]] <- fit(object at trModels[[i]],w=as.numeric(emp_trans[,1]==i),ntimes=object at ntimes) # check this
+ #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
+ } 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)
+ }
+
+ #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]])
+ }
+
+ for(k in 1:object at nresp) {
+ object at rModels[[i]][[k]] <- fit(object at rModels[[i]][[k]],w=fbo$gamma[,i])
+ # update logdens slot of the model
+ object at logdens[,k,i] <- logDens(object at rModels[[i]][[k]])
+ }
+ }
+
+ #object <- setpars(object,getpars(object)) # set parameters and recompute (bit of a roundabout way)
+
+ LL <- logLik(object)
+ if(verbose) cat("iteration",j,"logLik:",LL,"\n")
+ if( (LL >= LL.old) & (LL - LL.old < tol)) converge <- TRUE
+ }
+
+ object
+}
Added: trunk/depmix-test3EM3.R
===================================================================
--- trunk/depmix-test3EM3.R (rev 0)
+++ trunk/depmix-test3EM3.R 2008-03-06 17:32:42 UTC (rev 55)
@@ -0,0 +1,39 @@
+setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/rforge/depmix/trunk/")
+
+source("depmixS4.R")
+source("classes.R")
+source("hmModel.R")
+source("trGLM.r")
+source("lystig.R")
+source("fb.r")
+source("viterbi.r")
+source("EM3.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)
+
+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
Modified: trunk/viterbi.r
===================================================================
--- trunk/viterbi.r 2008-03-06 15:44:05 UTC (rev 54)
+++ trunk/viterbi.r 2008-03-06 17:32:42 UTC (rev 55)
@@ -1,5 +1,44 @@
-viterbi <- function(A,B,prior) {
+viterbi <- function(hmm) {
# returns the most likely state sequence
+ nt <- sum(hmm at ntimes)
+ lt <- length(hmm at ntimes)
+ et <- cumsum(hmm at ntimes)
+ bt <- c(1,et[-lt]+1)
+
+ ns <- hmm at nstates
+
+ delta <- psi <- matrix(nrow=nt,ncol=ns)
+ state <- vector(length=nt)
+
+ prior <- exp(logDens(hmm at initModel))
+
+ A <- hmm at trans
+ B <- apply(hmm at logdens,c(1,3),sum)
+
+ for(case in 1:lt) {
+ # initialization
+ delta[bt[case],] <- - (log(prior[case,]) + B[bt[case],])
+ psi[bt[case],] <- 0
+ # recursion
+ for(i in ((bt[case]+1):et[case])) {
+ for(j in 1:ns) {
+ delta[i,j] <- min(delta[i-1,] - log(A[i,,j])) - B[i,j]
+ k <- which.min(delta[i-1,] - log(A[i,,j]))
+ if(length(k) == 0) k <- 0
+ psi[i,j] <- k
+ }
+ }
+ #trace maximum likely state
+ state[et[case]] <- which.min(delta[et[case],])
+ for(i in (et[case]-1):bt[case]) {
+ state[i] <- psi[i+1,state[i+1]]
+ }
+ }
+ return(state)
+}
+
+viterbi.fb <- function(A,B,prior) {
+ # returns the most likely state sequence
nt <- nrow(B)
ns <- ncol(A)
delta <- psi <- matrix(nrow=nt,ncol=ns)
More information about the depmix-commits
mailing list