[Depmix-commits] r114 - / trunk trunk/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 27 23:25:06 CET 2008
Author: ingmarvisser
Date: 2008-03-27 23:25:06 +0100 (Thu, 27 Mar 2008)
New Revision: 114
Modified:
todo
trunk/R/EM.R
trunk/R/depmix.R
trunk/R/fb.R
trunk/R/lystig.R
trunk/depmixNew-test2.R
trunk/depmixNew-test3.R
trunk/depmixNew-test4.R
Log:
Moved apply(dens etc) into fb and lystig, updated todo
Modified: todo
===================================================================
--- todo 2008-03-26 08:49:01 UTC (rev 113)
+++ todo 2008-03-27 22:25:06 UTC (rev 114)
@@ -1,76 +1,41 @@
-TODO list for depmix package release 0.1-0
+TODO list for depmix package release 0.2-0
(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: done!!!!!
- 5) transInit.Rd: detailed description on constructing transition models: done!!!!!
- 6) others: lystig, fb: done!!!!
+1) restructure R directory and help file directory by use of package.skeleton
+
+2) make C versions of fb and lystig to speed things up a little bit
+(have a look at the profile code in depmixNew-test4.R to see where
+the speed problems are)
+
+3) add other distributions:
+ - glm options
+ - ???
-- rework em function and add a working example for it: Maarten
+
+TODO Medium term
+
+1) look at initial value generation?
-- balance scale data example with covariate on initial probabilities
-(this involves a latent class model rather than a hidden Markov model)
-Ingmar: done!!!!!
+2) Discrimination learning example, need covariates with these data!
-- review summary/print functions for depmix and depmix.fitted: done!!!!!
+3) Factor depmix model as example of user specified model
-- viterbi algorithm to get maximum a posteriori states for depmix model: Maarten, done !!!!!
-- viterbi algorithm to get posterior probabilities as part of depmix.fitted: skip!!
-
-- general documentation: depmix-introduction.pdf needs updating:
- 1) check credentials
- 2) add help for posterior
- 3)
-
-
-
-
-
-
-
-
-
-
-
TODO long term
-
-- generate initial values?!?!?
-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
+2) Analytical gradients in lystig + standard errors of parameters
-3) Speed optimization: lystig in C?
+3) Multi group possibilities: use group factor in call to depmix
-4) Multi group possibilities: use group factor in call to depmix
+4) Missing data options?
-5) Missing data options?
-6) Standard errors of parameters
-
-
Other capabilities
1) Look at log and exp of matrices to have continuous time observations
Modified: trunk/R/EM.R
===================================================================
--- trunk/R/EM.R 2008-03-26 08:49:01 UTC (rev 113)
+++ trunk/R/EM.R 2008-03-27 22:25:06 UTC (rev 114)
@@ -17,11 +17,11 @@
j <- 0
A <- object at trDens
- B <- apply(object at dens,c(1,3),prod)
+# B <- apply(object at dens,c(1,3),prod)
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))
LL <- fbo$logLike
LL.old <- LL + 1
@@ -61,7 +61,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))
LL <- fbo$logLike
if(verbose&((j%%5)==0)) cat("iteration",j,"logLik:",LL,"\n")
Modified: trunk/R/depmix.R
===================================================================
--- trunk/R/depmix.R 2008-03-26 08:49:01 UTC (rev 113)
+++ trunk/R/depmix.R 2008-03-27 22:25:06 UTC (rev 114)
@@ -291,8 +291,8 @@
# depends on getpars and nobs
setMethod("logLik",signature(object="depmix"),
function(object,method="lystig") {
- if(method=="fb") ll <- fb(object at init,object at trDens,apply(object at dens,c(1,3),prod),object at ntimes,object at stationary)$logLike
- if(method=="lystig") ll <- lystig(object at init,object at trDens,apply(object at dens,c(1,3),prod),object at ntimes,object at stationary)$logLike
+ if(method=="fb") ll <- fb(object at init,object at trDens,object at dens,object at ntimes,object at stationary)$logLike
+ if(method=="lystig") ll <- lystig(object at init,object at trDens,object at dens,object at ntimes,object at stationary)$logLike
attr(ll, "df") <- freepars(object)
attr(ll, "nobs") <- nobs(object)
class(ll) <- "logLik"
Modified: trunk/R/fb.R
===================================================================
--- trunk/R/fb.R 2008-03-26 08:49:01 UTC (rev 113)
+++ trunk/R/fb.R 2008-03-27 22:25:06 UTC (rev 114)
@@ -17,7 +17,10 @@
# NOTE: to prevent underflow, alpha and beta are scaled, using c
# NOTE: xi[i,j] = P(S[t] = j & S[t+1] = i)
-
+
+
+ B <- apply(B,c(1,3),prod)
+
nt <- nrow(B)
ns <- ncol(init)
Modified: trunk/R/lystig.R
===================================================================
--- trunk/R/lystig.R 2008-03-26 08:49:01 UTC (rev 113)
+++ trunk/R/lystig.R 2008-03-27 22:25:06 UTC (rev 114)
@@ -16,12 +16,14 @@
# IN COMPUTING ALPHA AND BETA BUT IS NOW NECCESSARY IN COMPUTING XI
# A = K*K matrix with transition probabilities, from column to row !!!!!!!
# change to T*K*K
+
+ B <- apply(B,c(1,3),prod)
- # B = T*K matrix with elements ab_{ij} = P(y_i|s_j)
+ # B = T*K*nresp matrix with elements ab_{tij} = P(y_t_i|s_j)
# init = K vector with initial probabilities !!!OR!!! K*length(ntimes) matrix with initial probs per case
# Returns: 'sca'le factors, recurrent variables 'phi', loglikelihood
-
- nt <- nrow(B)
+
+ nt <- nrow(B)
ns <- ncol(init)
if(!is.null(ntimes)) {
Modified: trunk/depmixNew-test2.R
===================================================================
--- trunk/depmixNew-test2.R 2008-03-26 08:49:01 UTC (rev 113)
+++ trunk/depmixNew-test2.R 2008-03-27 22:25:06 UTC (rev 114)
@@ -25,66 +25,9 @@
logLik(mod)
-source("R/EM.R")
-
-gc()
-Rprof()
-mod1 <- em(mod,verb=TRUE)
-Rprof(NULL)
-summaryRprof()
-
-
-
-dd <- mod1 at dens[1:100,,]
-
-iter=1000
-
-gc()
-system.time(
- for(i in 1:iter) {
- x<-apply(dd,c(1,3),prod)
- }
-)
-
-gc()
-system.time(
- for(i in 1:iter) {
- x<-exp(apply(log(dd),c(1,3),sum))
- }
-)
-
-
-
-source("R/EM.R")
-
-gc()
-Rprof()
-mod1 <- em(mod)
-Rprof(NULL)
-summaryRprof()
-
-
-
-
-
-# object at parameters$coefficients <- pars$coefficients
-
-
-
mod1 <- fit(mod)
ll <- logLik(mod1)
-source("R/depmix.fitted.R")
-source("R/viterbi.R")
-
-x1=viterbi(mod1)
-x2=viterbi2(mod1)
-round(head(x1),3)
-round(head(x2),3)
-cor(x1[,1],x2[,1])
-
-
-
#
# Test optimization using Rdonlp2
#
Modified: trunk/depmixNew-test3.R
===================================================================
--- trunk/depmixNew-test3.R 2008-03-26 08:49:01 UTC (rev 113)
+++ trunk/depmixNew-test3.R 2008-03-27 22:25:06 UTC (rev 114)
@@ -48,26 +48,7 @@
# mod4 <- fit(mod2,fixed=fixed,method="donlp")
-# profiling code
-gc()
-Rprof()
-mod1 <- em(mod)
-Rprof(NULL)
-summaryRprof()
-
-source("R/EM.R")
-source("R/responses.R")
-
-gc()
-Rprof()
-mod1 <- fit(mod,fixed=fixed)
-Rprof(NULL)
-summaryRprof()
-
-
-
-
Modified: trunk/depmixNew-test4.R
===================================================================
--- trunk/depmixNew-test4.R 2008-03-26 08:49:01 UTC (rev 113)
+++ trunk/depmixNew-test4.R 2008-03-27 22:25:06 UTC (rev 114)
@@ -2,26 +2,54 @@
#
# Started by Ingmar Visser 26-2-2008
#
-# Usage: go to trunk directory and source this file in R
+# Usage: go to trunk directory and source("depmixNew-test4.R")
#
#
# BALANCE SCALE data example with age as covariate on class membership
#
-library(depmixS4)
+# library(depmixS4)
setwd("/Users/ivisser/Documents/projects/depmixProject/depmixNew/rforge/depmix/trunk/")
+#
+# optimization speed profile: case 1: latent class data
+#
+
+# source("depmixNew-test4.R")
+
data(balance)
-# 'log Lik.' -1083.036 (df=9)
+# now fit some latent class models
+trstart=c(1,0,0,1) # as this is a latent class model, the transition are not optimized
+instart=c(0.5,0.5)
+set.seed(1)
+respstart=runif(16)
+# note that ntimes argument is used to make this a mixture model
+mod1 <- depmix(list(d1~1,d2~1,d3~1,d4~1), data=balance, nstates=2,
+ family=list(multinomial(),multinomial(),multinomial(),multinomial()),
+ respstart=respstart,trstart=trstart,instart=instart,
+ ntimes=rep(1,nrow(balance)))
+logLik(mod1)
+
+# mod1 <- fit(mod)
+
+gc()
+Rprof(file="lca1")
+mod1 <- fit(mod1)
+Rprof(NULL)
+summaryRprof("lca1")
+
+
#
-# Add age as covariate on class membership
+# optimization speed profile: case 1: latent class data with cov on prior
#
+data(balance)
+
instart=c(0.5,0.5,0,0)
respstart=c(rep(c(0.1,0.9),4),rep(c(0.9,0.1),4))
trstart=c(1,0,0,1)
@@ -30,8 +58,50 @@
trstart=trstart, instart=instart, respstart=respstart,
ntimes=rep(1,nrow(balance)), prior=~age, initdata=balance)
-fixed <- c(1,0,1,0,1,1,1,1,rep(c(1,0),8))
-mod3 <- fit(mod2,fixed=fixed)
+gc()
+Rprof(file="lca2")
+mod2 <- fit(mod2)
+Rprof(NULL)
+summaryRprof("lca2")
+#
+# speed data, no covariates
+#
+data(speed)
+
+set.seed(1)
+trstart=runif(4)
+
+mod <- depmix(list(rt~1,corr~1),data=speed,nstates=2,family=list(gaussian(),multinomial()), trstart=trstart)
+
+logLik(mod)
+
+gc()
+Rprof("speed1")
+mod1 <- fit(mod)
+Rprof(NULL)
+summaryRprof("speed1")
+
+
+trstart=c(0.899,0.101,0,0.01,0.084,0.916,0,0)
+instart=c(0.5,0.5)
+resp <- c(5.52,0.202,0.472,0.528,6.39,0.24,0.098,0.902)
+
+mod <- depmix(list(rt~1,corr~1),data=speed,transition=~Pacc,nstates=2,family=list(gaussian(),multinomial()),
+ respstart=resp,trstart=trstart,instart=instart)
+
+logLik(mod)
+
+gc()
+Rprof("speed2")
+mod1 <- fit(mod)
+Rprof(NULL)
+summaryRprof("speed2")
+
+summaryRprof("lca1")
+summaryRprof("lca2")
+summaryRprof("speed1")
+summaryRprof("speed2")
+
More information about the depmix-commits
mailing list