[Depmix-commits] r269 - / trunk/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed May 13 15:07:39 CEST 2009
Author: ingmarvisser
Date: 2009-05-13 15:07:39 +0200 (Wed, 13 May 2009)
New Revision: 269
Modified:
todo
trunk/R/EM.R
trunk/R/depmixfit.R
trunk/R/response-class.R
Log:
Updated todo (and minor comments in some R files)
Modified: todo
===================================================================
--- todo 2009-04-21 10:24:35 UTC (rev 268)
+++ todo 2009-05-13 13:07:39 UTC (rev 269)
@@ -6,8 +6,8 @@
done in the test files)
- similarly for actual markov models and mixture models
-2) add example of adding new response model to JSS paper (still looking
-for a good candidate model ...)
+2) add example of adding new response model to JSS paper: exgaus models
+for speed data
3) check models with boundary values in the transition and initial
state probabilities (possibly work on identity link for multinomial
@@ -23,15 +23,24 @@
6) Tests:
- make output files for test files
+7) Bugs:
+ - check error sent by Rita Gaio
+ - check error by Maartje
+
TODO Medium term
1) Discrimination learning example, need covariates with these data!
-
-2) Factor depmix model as example of user specified model
-3) Discrimination learning example: find data set with interesting covariates
+2) BB or other optimization routines??
+3) Speed issues:
+ - add gradients
+ - run Rprofile tests to determine speed limits
+ - possibly move lystig and forward/backward to C routines
+
+4) Get Hessian for standard errors (computed by e.g. Louis' method)
+
TODO long term
@@ -40,25 +49,23 @@
1) Multi group possibilities: use group factor in call to depmix??
2) Missing data options?
+
+3) Mixture of depmix models?
Other capabilities
1) Look at log and exp of matrices to have continuous time observations
-(also see mlm)
+(also see mlm, see Bockenholt 2005
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
+4) Check glmc pack for constrained glm models
-5) Check glmc pack for constrained glm models
+5) stochmod pack for stochastic models
-6) stochmod pack for stochastic models
-
-7) use relevel function to recode factors into different base categories
-
-8) Standard errors (Hessian) for EM (computed by e.g. Louis' method)
+6) use relevel function to recode factors into different base categories
Modified: trunk/R/EM.R
===================================================================
--- trunk/R/EM.R 2009-04-21 10:24:35 UTC (rev 268)
+++ trunk/R/EM.R 2009-05-13 13:07:39 UTC (rev 269)
@@ -3,49 +3,49 @@
#
em <- function(object,maxit=100,tol=1e-6,verbose=FALSE,...) {
- if(!is(object,"mix")) stop("object is not of class '(dep)mix'")
- call <- match.call()
- if(is(object,"depmix")) {
- call[[1]] <- as.name("em.depmix")
- } else {
- call[[1]] <- as.name("em.mix")
- }
- object <- eval(call, parent.frame())
- object
+ if(!is(object,"mix")) stop("object is not of class '(dep)mix'")
+ call <- match.call()
+ if(is(object,"depmix")) {
+ call[[1]] <- as.name("em.depmix")
+ } else {
+ call[[1]] <- as.name("em.mix")
+ }
+ object <- eval(call, parent.frame())
+ object
}
-
+# em for lca and mixture models
em.mix <- function(object,maxit=100,tol=1e-6,verbose=FALSE,...) {
- if(!is(object,"mix")) stop("object is not of class 'mix'")
-
- ns <- object at nstates
-
+ if(!is(object,"mix")) stop("object is not of class 'mix'")
+
+ ns <- object at nstates
+
ntimes <- ntimes(object)
lt <- length(ntimes)
et <- cumsum(ntimes)
bt <- c(1,et[-lt]+1)
-
+
converge <- FALSE
j <- 0
# compute responsibilities
- B <- apply(object at dens,c(1,3),prod)
- gamma <- object at init*B
- LL <- sum(log(rowSums(gamma)))
- # normalize
- gamma <- gamma/rowSums(gamma)
-
+ B <- apply(object at dens,c(1,3),prod)
+ gamma <- object at init*B
+ LL <- sum(log(rowSums(gamma)))
+ # normalize
+ gamma <- gamma/rowSums(gamma)
+
LL.old <- LL + 1
-
+
while(j <= maxit & !converge) {
-
+
# maximization
-
+
# should become object at prior <- fit(object at prior)
object at prior@y <- gamma[bt,,drop=FALSE]
object at prior <- fit(object at prior, w=NULL,ntimes=NULL)
object at init <- dens(object at prior)
-
+
for(i in 1:ns) {
for(k in 1:nresp(object)) {
object at response[[i]][[k]] <- fit(object at response[[i]][[k]],w=gamma[,i])
@@ -62,7 +62,9 @@
gamma <- gamma/rowSums(gamma)
# print stuff
- if(verbose&((j%%5)==0)) cat("iteration",j,"logLik:",LL,"\n")
+ if(verbose&((j%%5)==0)) {
+ cat("iteration",j,"logLik:",LL,"\n")
+ }
if( (LL >= LL.old) & (LL - LL.old < tol)) {
cat("iteration",j,"logLik:",LL,"\n")
@@ -88,6 +90,7 @@
}
+# em for hidden markov models
em.depmix <- function(object,maxit=100,tol=1e-6,verbose=FALSE,...) {
if(!is(object,"depmix")) stop("object is not of class '(dep)mix'")
Modified: trunk/R/depmixfit.R
===================================================================
--- trunk/R/depmixfit.R 2009-04-21 10:24:35 UTC (rev 268)
+++ trunk/R/depmixfit.R 2009-05-13 13:07:39 UTC (rev 269)
@@ -29,7 +29,7 @@
# set those fixed parameters in the appropriate submodels
object <- setpars(object,fixed,which="fixed")
-
+
if(method=="EM") {
object <- em(object,verbose=TRUE,...)
}
@@ -84,7 +84,7 @@
lin.l <- c(lin.l,conrows.lower)
}
}
-
+
# select only those columns of the constraint matrix that correspond to non-fixed parameters
linconFull <- lincon
lincon <- lincon[,!fixed,drop=FALSE]
@@ -95,7 +95,7 @@
mycontrol <- function(info) {
return(TRUE)
}
-
+
# optimize the parameters
result <- donlp2(pars,logl,
par.upper=par.u,
Modified: trunk/R/response-class.R
===================================================================
--- trunk/R/response-class.R 2009-04-21 10:24:35 UTC (rev 268)
+++ trunk/R/response-class.R 2009-05-13 13:07:39 UTC (rev 269)
@@ -33,3 +33,6 @@
return(sum(!object at fixed))
}
)
+
+
+
More information about the depmix-commits
mailing list