[Depmix-commits] r401 - in tags/release-0.3-0: . R data inst inst/doc man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Mar 9 13:40:10 CET 2010
Author: ingmarvisser
Date: 2010-03-09 13:40:10 +0100 (Tue, 09 Mar 2010)
New Revision: 401
Added:
tags/release-0.3-0/R/
tags/release-0.3-0/R/EM.R
tags/release-0.3-0/R/allGenerics.R
tags/release-0.3-0/R/depmix-class.R
tags/release-0.3-0/R/depmix.R
tags/release-0.3-0/R/depmixAIC.R
tags/release-0.3-0/R/depmixBIC.R
tags/release-0.3-0/R/depmixfit-class.R
tags/release-0.3-0/R/depmixfit.R
tags/release-0.3-0/R/depmixsim-class.R
tags/release-0.3-0/R/fb.R
tags/release-0.3-0/R/forwardbackward.R
tags/release-0.3-0/R/freepars.R
tags/release-0.3-0/R/getpars.R
tags/release-0.3-0/R/llratio.R
tags/release-0.3-0/R/logLik.R
tags/release-0.3-0/R/lystig.R
tags/release-0.3-0/R/makeDepmix.R
tags/release-0.3-0/R/makePriorModel.R
tags/release-0.3-0/R/makeResponseModels.R
tags/release-0.3-0/R/makeTransModels.R
tags/release-0.3-0/R/mlogit.R
tags/release-0.3-0/R/multinomial.R
tags/release-0.3-0/R/nobs.R
tags/release-0.3-0/R/pa2conr.R
tags/release-0.3-0/R/response-class.R
tags/release-0.3-0/R/responseGLM.R
tags/release-0.3-0/R/responseGLMBINOM.R
tags/release-0.3-0/R/responseGLMGAMMA.R
tags/release-0.3-0/R/responseGLMMULTINOM.R
tags/release-0.3-0/R/responseGLMPOISSON.R
tags/release-0.3-0/R/responseMVN.R
tags/release-0.3-0/R/responseNORM.R
tags/release-0.3-0/R/setpars.R
tags/release-0.3-0/R/stationary.R
tags/release-0.3-0/R/transInit.R
tags/release-0.3-0/R/viterbi.R
tags/release-0.3-0/data/
tags/release-0.3-0/data/balance.rda
tags/release-0.3-0/data/speed.rda
tags/release-0.3-0/inst/
tags/release-0.3-0/inst/CITATION
tags/release-0.3-0/inst/doc/
tags/release-0.3-0/inst/doc/depmix-intro.pdf
tags/release-0.3-0/man/
tags/release-0.3-0/man/AIC.Rd
tags/release-0.3-0/man/GLMresponse.Rd
tags/release-0.3-0/man/balance.Rd
tags/release-0.3-0/man/depmix-class.Rd
tags/release-0.3-0/man/depmix-internal.Rd
tags/release-0.3-0/man/depmix-methods.Rd
tags/release-0.3-0/man/depmix.Rd
tags/release-0.3-0/man/depmix.fit.Rd
tags/release-0.3-0/man/depmix.fitted-class.Rd
tags/release-0.3-0/man/depmix.sim-class.Rd
tags/release-0.3-0/man/depmixS4-package.Rd
tags/release-0.3-0/man/forwardbackward.Rd
tags/release-0.3-0/man/llratio.Rd
tags/release-0.3-0/man/makeDepmix.Rd
tags/release-0.3-0/man/mix-class.Rd
tags/release-0.3-0/man/mix.Rd
tags/release-0.3-0/man/mix.fitted-class.Rd
tags/release-0.3-0/man/mix.sim-class.Rd
tags/release-0.3-0/man/posterior.Rd
tags/release-0.3-0/man/response-class.Rd
tags/release-0.3-0/man/response-classes.Rd
tags/release-0.3-0/man/responses.Rd
tags/release-0.3-0/man/simulate.Rd
tags/release-0.3-0/man/speed.Rd
tags/release-0.3-0/man/transInit.Rd
tags/release-0.3-0/tests/
tags/release-0.3-0/tests/test1speed.R
tags/release-0.3-0/tests/test1speed.Rout.save
tags/release-0.3-0/tests/test2getsetpars.R
tags/release-0.3-0/tests/test2getsetpars.Rout.save
tags/release-0.3-0/tests/test3responses.R
Log:
Moved pkg to tag 0.3-0 directory
Copied: tags/release-0.3-0/R/EM.R (from rev 398, pkg/depmixS4/R/EM.R)
===================================================================
--- tags/release-0.3-0/R/EM.R (rev 0)
+++ tags/release-0.3-0/R/EM.R 2010-03-09 12:40:10 UTC (rev 401)
@@ -0,0 +1,182 @@
+#
+# Maarten Speekenbrink 23-3-2008
+#
+
+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
+}
+
+# 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
+
+ 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)
+
+ 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])
+ # update dens slot of the model
+ object at dens[,k,i] <- dens(object at response[[i]][[k]])
+ }
+ }
+
+ # expectation
+ B <- apply(object at dens,c(1,3),prod)
+ gamma <- object at init*B
+ LL <- sum(log(rowSums(gamma)))
+ # normalize
+ gamma <- gamma/rowSums(gamma)
+
+ # print stuff
+ 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")
+ converge <- TRUE
+ }
+
+ LL.old <- LL
+ j <- j+1
+
+ }
+
+ class(object) <- "mix.fitted"
+
+ if(converge) object at message <- "Log likelihood converged to within tol."
+ else object at message <- "'maxit' iterations reached in EM without convergence."
+
+ # no constraints in EM
+ object at conMat <- matrix()
+ object at lin.lower <- numeric()
+ object at lin.upper <- numeric()
+
+ object
+
+}
+
+# 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'")
+
+ ns <- object at nstates
+
+ ntimes <- ntimes(object)
+ lt <- length(ntimes)
+ et <- cumsum(ntimes)
+ bt <- c(1,et[-lt]+1)
+
+ converge <- FALSE
+ j <- 0
+
+ # A <- object at trDens
+ # B <- object at dens
+ # init <- object at init
+
+ # initial expectation
+ fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
+ LL <- fbo$logLike
+ LL.old <- LL + 1
+
+ while(j <= maxit & !converge) {
+
+ # maximization
+
+ # should become object at prior <- fit(object at prior)
+ object at prior@y <- fbo$gamma[bt,,drop=FALSE]
+ object at prior <- fit(object at prior, w=NULL,ntimes=NULL)
+ object at init <- dens(object at prior)
+
+ trm <- matrix(0,ns,ns)
+ for(i in 1:ns) {
+ if(max(ntimes(object)>1)) { # skip transition parameters update in case of latent class model
+ if(!object at stationary) {
+ 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])
+ }
+ # 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)
+ }
+ # update trDens slot of the model
+ object at trDens[,,i] <- dens(object at transition[[i]])
+ }
+ }
+
+ for(i in 1:ns) {
+
+ for(k in 1:nresp(object)) {
+ object at response[[i]][[k]] <- fit(object at response[[i]][[k]],w=fbo$gamma[,i])
+ # update dens slot of the model
+ object at dens[,k,i] <- dens(object at response[[i]][[k]])
+ }
+ }
+
+ # expectation
+ fbo <- fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),stationary=object at stationary)
+ LL <- fbo$logLike
+
+ 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")
+ converge <- TRUE
+ }
+
+ LL.old <- LL
+ j <- j+1
+
+ }
+
+ #if(class(object)=="depmix") class(object) <- "depmix.fitted"
+ #if(class(object)=="mix") class(object) <- "mix.fitted"
+
+ class(object) <- "depmix.fitted"
+
+ if(converge) object at message <- "Log likelihood converged to within tol."
+ else object at message <- "'maxit' iterations reached in EM without convergence."
+
+ # no constraints in EM
+ object at conMat <- matrix()
+ object at lin.lower <- numeric()
+ object at lin.upper <- numeric()
+
+ object
+}
Copied: tags/release-0.3-0/R/allGenerics.R (from rev 398, pkg/depmixS4/R/allGenerics.R)
===================================================================
--- tags/release-0.3-0/R/allGenerics.R (rev 0)
+++ tags/release-0.3-0/R/allGenerics.R 2010-03-09 12:40:10 UTC (rev 401)
@@ -0,0 +1,76 @@
+
+#
+# Ingmar Visser, 23-3-2008
+#
+
+.First.lib <- function(lib, pkg) {
+ require(stats)
+ require(methods)
+ require(MASS)
+ require(nnet)
+ require(MCMCpack)
+}
+
+.Last.lib <- function(libpath) {}
+
+# Guess what: all generics
+
+setGeneric("depmix", function(response,data=NULL,nstates,transition=~1,family=gaussian(),prior=~1,initdata=NULL,
+ respstart=NULL,trstart=NULL,instart=NULL,ntimes=NULL, ...) standardGeneric("depmix"))
+
+setGeneric("GLMresponse", function(formula, data = NULL, family = gaussian(), pstart =
+ NULL, fixed = NULL, prob=TRUE, ...) standardGeneric("GLMresponse"))
+
+setGeneric("MVNresponse", function(formula, data = NULL,pstart=NULL,fixed=NULL,...) standardGeneric("MVNresponse"))
+
+setGeneric("transInit", function(formula, nstates, data = NULL, family = multinomial(),
+ pstart = NULL, fixed = NULL, prob=TRUE, ...) standardGeneric("transInit"))
+
+setGeneric("npar", function(object, ...) standardGeneric("npar"))
+
+setGeneric("nobs", function(object, ...) standardGeneric("nobs"))
+
+setGeneric("ntimes", function(object, ...) standardGeneric("ntimes"))
+
+setGeneric("nstates", function(object, ...) standardGeneric("nstates"))
+
+setGeneric("nresp", function(object, ...) standardGeneric("nresp"))
+
+setGeneric("freepars", function(object, ...) standardGeneric("freepars"))
+
+setGeneric("nlin", function(object, ...) standardGeneric("nlin"))
+
+# setGeneric("logLik", function(object, ...) standardGeneric("logLik"))
+
+setGeneric("fit", function(object, ...) standardGeneric("fit"))
+
+setGeneric("posterior", function(object, ...) standardGeneric("posterior"))
+
+setGeneric("forwardbackward", function(object, ...) standardGeneric("forwardbackward"))
+
+setGeneric("simulate", function(object,nsim=1,seed=NULL, ...) standardGeneric("simulate"))
+
+setGeneric("predict", function(object, ...) standardGeneric("predict"))
+
+# setGeneric("AIC", function(object, ..., k=2) standardGeneric("AIC"))
+
+setGeneric("BIC", function(object, ...) standardGeneric("BIC"))
+
+setGeneric("getdf",function(object) standardGeneric("getdf"))
+
+setGeneric("setpars", function(object,values,which="pars",...) standardGeneric("setpars"))
+
+setGeneric("getpars", function(object,which="pars",...) standardGeneric("getpars"))
+
+setGeneric("logDens",function(object,...) standardGeneric("logDens"))
+
+setGeneric("dens",function(object,...) standardGeneric("dens"))
+
+setGeneric("summary")
+
+setGeneric("ntimes", function(object, ...) standardGeneric("ntimes"))
+
+setGeneric("nresp", function(object, ...) standardGeneric("nresp"))
+
+setGeneric("is.stationary", function(object,...) standardGeneric("is.stationary"))
+
Copied: tags/release-0.3-0/R/depmix-class.R (from rev 398, pkg/depmixS4/R/depmix-class.R)
===================================================================
--- tags/release-0.3-0/R/depmix-class.R (rev 0)
+++ tags/release-0.3-0/R/depmix-class.R 2010-03-09 12:40:10 UTC (rev 401)
@@ -0,0 +1,307 @@
+
+#
+# Ingmar Visser, 11-6-2008
+#
+
+#
+# DEPMIX CLASS BELOW THE MIX CLASS
+#
+
+#
+# Class definition, accessor functions, print and summary methods
+#
+
+#
+# MIX CLASS
+#
+
+setClass("mix",
+ representation(response="list", # response models
+ prior="ANY", # the prior model (multinomial logistic)
+ dens="array", # response densities (B)
+ init="array", # usually called pi
+ nstates="numeric",
+ nresp="numeric",
+ ntimes="numeric",
+ npars="numeric" # number of parameters
+ )
+)
+
+# accessor functions
+setMethod("npar","mix",
+ function(object) return(object at npars)
+)
+
+setMethod("ntimes","mix",
+ function(object) return(object at ntimes)
+)
+
+setMethod("nstates","mix",
+ function(object) return(object at nstates)
+)
+
+setMethod("nresp","mix",
+ function(object) return(object at nresp)
+)
+
+setMethod("is.stationary",signature(object="mix"),
+ function(object) {
+ return(TRUE)
+ }
+)
+
+setMethod("simulate",signature(object="mix"),
+ function(object,nsim=1,seed=NULL,...) {
+
+ if(!is.null(seed)) set.seed(seed)
+
+ ntim <- ntimes(object)
+ nt <- sum(ntim)
+ bt <- 1:nt
+
+ nr <- nresp(object)
+ ns <- nstates(object)
+
+ # simulate state sequences first, then observations
+
+ # random generation is slow when done separately for each t, so first draw
+ # variates for all t, and then determine state sequences iteratively
+ states <- array(,dim=c(nt,nsim))
+ states[bt,] <- simulate(object at prior,n=nsim,is.prior=T)
+ sims <- array(,dim=c(nt,ns,nsim))
+
+# for(i in 1:ns) {
+# if(is.stationary(object)) {
+# # TODO: this is a temporary fix!!!
+# sims[,i,] <- simulate(object at transition[[i]],nsim=nsim,times=rep(1,nt))
+# } else {
+# sims[,i,] <- simulate(object at transition[[i]],nsim=nsim)
+# }
+# }
+# # track states
+# for(case in 1:lt) {
+# for(i in (bt[case]+1):et[case]) {
+# states[i,] <- sims[cbind(i,states[i-1,],1:nsim)]
+# }
+# }
+
+ states <- as.vector(states)
+ responses <- list(length=nr)
+ #responses <- array(,dim=c(nt,nr,nsim))
+ for(i in 1:nr) {
+ tmp <- matrix(,nrow=nt*nsim,ncol=NCOL(object at response[[1]][[i]]@y))
+ for(j in 1:ns) {
+ tmp[states==j,] <- simulate(object at response[[j]][[i]],nsim=nsim)[states==j,]
+ }
+ responses[[i]] <- tmp
+ }
+
+ # generate new mix.sim object
+ class(object) <- "mix.sim"
+ object at states <- as.matrix(states)
+
+ object at prior@x <- as.matrix(apply(object at prior@x,2,rep,nsim))
+ for(j in 1:ns) {
+# if(!is.stationary(object)) object at transition[[j]]@x <- as.matrix(apply(object at transition[[j]]@x,2,rep,nsim))
+ for(i in 1:nr) {
+ object at response[[j]][[i]]@y <- as.matrix(responses[[i]])
+ object at response[[j]][[i]]@x <- as.matrix(apply(object at response[[j]][[i]]@x,2,rep,nsim))
+ }
+ }
+ object at ntimes <- rep(object at ntimes,nsim)
+
+ # make appropriate array for transition densities
+ nt <- sum(object at ntimes)
+
+ # make appropriate array for response densities
+ dns <- array(,c(nt,nr,ns))
+
+ # compute observation and transition densities
+ for(i in 1:ns) {
+ for(j in 1:nr) {
+ dns[,j,i] <- dens(object at response[[i]][[j]]) # remove this response as an argument from the call to setpars
+ }
+ }
+
+ # compute initial state probabilties
+ object at init <- dens(object at prior)
+ object at dens <- dns
+
+ return(object)
+ }
+)
+
+
+
+#
+# PRINT method
+#
+
+setMethod("show","mix",
+ function(object) {
+ cat("Initial state probabilties model \n")
+ print(object at prior)
+ cat("\n")
+ for(i in 1:object at nstates) {
+ cat("Response model(s) for state", i,"\n\n")
+ for(j in 1:object at nresp) {
+ cat("Response model for response",j,"\n")
+ print(object at response[[i]][[j]])
+ cat("\n")
+ }
+ cat("\n")
+ }
+ }
+)
+
+#
+# SUMMARY method: to do
+#
+
+
+#
+# Ingmar Visser, 23-3-2008
+#
+
+#
+# Class definition, accessor functions, print and summary methods
+#
+
+#
+# DEPMIX CLASS
+#
+
+setClass("depmix",
+ representation(transition="list", # transition models (multinomial logistic)
+ trDens="array", # transition densities (A)
+ stationary="logical"
+ ),
+ contains="mix"
+)
+
+#
+# PRINT method
+#
+
+setMethod("show","depmix",
+ function(object) {
+ cat("Initial state probabilties model \n")
+ print(object at prior)
+ cat("\n")
+ for(i in 1:object at nstates) {
+ cat("Transition model for state (component)", i,"\n")
+ print(object at transition[[i]])
+ cat("\n")
+ }
+ cat("\n")
+ for(i in 1:object at nstates) {
+ cat("Response model(s) for state", i,"\n\n")
+ for(j in 1:object at nresp) {
+ cat("Response model for response",j,"\n")
+ print(object at response[[i]][[j]])
+ cat("\n")
+ }
+ cat("\n")
+ }
+ }
+)
+
+setMethod("is.stationary",signature(object="depmix"),
+ function(object) {
+ return(object at stationary)
+ }
+)
+
+setMethod("simulate",signature(object="depmix"),
+ function(object,nsim=1,seed=NULL,...) {
+
+ if(!is.null(seed)) set.seed(seed)
+
+ ntim <- ntimes(object)
+ nt <- sum(ntim)
+ lt <- length(ntim)
+ et <- cumsum(ntim)
+ bt <- c(1,et[-lt]+1)
+
+ nr <- nresp(object)
+ ns <- nstates(object)
+
+ # simulate state sequences first, then observations
+
+ # random generation is slow when done separately for each t, so first draw
+ # variates for all t, and then determine state sequences iteratively
+ states <- array(,dim=c(nt,nsim))
+ states[bt,] <- simulate(object at prior,n=nsim,is.prior=T)
+ sims <- array(,dim=c(nt,ns,nsim))
+ for(i in 1:ns) {
+ if(is.stationary(object)) {
+ # TODO: this is a temporary fix!!!
+ sims[,i,] <- simulate(object at transition[[i]],nsim=nsim,times=rep(1,nt))
+ } else {
+ sims[,i,] <- simulate(object at transition[[i]],nsim=nsim)
+ }
+ }
+ # track states
+ for(case in 1:lt) {
+ for(i in (bt[case]+1):et[case]) {
+ states[i,] <- sims[cbind(i,states[i-1,],1:nsim)]
+ }
+ }
+
+ states <- as.vector(states)
+ responses <- list(length=nr)
+ #responses <- array(,dim=c(nt,nr,nsim))
+ for(i in 1:nr) {
+ tmp <- matrix(,nrow=nt*nsim,ncol=NCOL(object at response[[1]][[i]]@y))
+ for(j in 1:ns) {
+ tmp[states==j,] <- simulate(object at response[[j]][[i]],nsim=nsim)[states==j,]
+ }
+ responses[[i]] <- tmp
+ }
+
+ # generate new depmix.sim object
+ class(object) <- "depmix.sim"
+ object at states <- as.matrix(states)
+
+ object at prior@x <- as.matrix(apply(object at prior@x,2,rep,nsim))
+ for(j in 1:ns) {
+ if(!is.stationary(object)) object at transition[[j]]@x <- as.matrix(apply(object at transition[[j]]@x,2,rep,nsim))
+ for(i in 1:nr) {
+ object at response[[j]][[i]]@y <- as.matrix(responses[[i]])
+ object at response[[j]][[i]]@x <- as.matrix(apply(object at response[[j]][[i]]@x,2,rep,nsim))
+ }
+ }
+ object at ntimes <- rep(object at ntimes,nsim)
+
+ # make appropriate array for transition densities
+ nt <- sum(object at ntimes)
+ if(is.stationary(object)) trDens <- array(0,c(1,ns,ns)) else trDens <- array(0,c(nt,ns,ns))
+
+ # make appropriate array for response densities
+ dns <- array(,c(nt,nr,ns))
+
+ # compute observation and transition densities
+ for(i in 1:ns) {
+ for(j in 1:nr) {
+ dns[,j,i] <- dens(object at response[[i]][[j]]) # remove this response as an argument from the call to setpars
+ }
+ trDens[,,i] <- dens(object at transition[[i]])
+ }
+
+ # compute initial state probabilties
+ object at init <- dens(object at prior)
+ object at trDens <- trDens
+ object at dens <- dns
+
+ return(object)
+ }
+)
+
+
+
+#
+# SUMMARY method: to do
+#
+
+
+
Copied: tags/release-0.3-0/R/depmix.R (from rev 398, pkg/depmixS4/R/depmix.R)
===================================================================
--- tags/release-0.3-0/R/depmix.R (rev 0)
+++ tags/release-0.3-0/R/depmix.R 2010-03-09 12:40:10 UTC (rev 401)
@@ -0,0 +1,93 @@
+#
+# Ingmar Visser, 11-6-2008
+#
+
+#
+# Main function to construct mix models
+#
+
+#
+# UNIVARIATE AND MULTIVARIATE MIXTURE OF GLM'S
+#
+
+
+setGeneric("mix", function(response,data=NULL,nstates,family=gaussian(),prior=~1,initdata=NULL,
+ respstart=NULL,instart=NULL, ...) standardGeneric("mix"))
+
+
+setMethod("mix",
+ signature(response="ANY"),
+ function(response, data=NULL, nstates, family=gaussian(), prior=~1, initdata=NULL,
+ respstart=NULL, instart=NULL, ...) {
+
+ # make response models
+ response <- makeResponseModels(response=response,data=data,nstates=nstates,family=family,values=respstart)
+
+ # FIX ME: this only works if data are actually provided ... (maybe make this obligatory ...)
+ ntimes <- rep(1,nrow(data))
+
+ # make prior model
+ prior <- makePriorModel(nstates=nstates,ncases=length(ntimes),formula=prior,data=initdata,values=instart)
+
+ # call main depmix with all these models, ntimes and stationary
+ model <- makeMix(response=response,prior=prior)
+
+ # deal with starting values here!!!!!!
+
+ return(model)
+ }
+)
+
+#
+# Ingmar Visser, 23-3-2008
+#
+
+#
+# Main function to construct depmix models
+#
+
+#
+# UNIVARIATE AND MULTIVARIATE MARKOV MIXTURE OF GLM'S
+#
+
+setMethod("depmix",
+ signature(response="ANY"),
+ function(response, data=NULL, nstates, transition=~1, family=gaussian(), prior=~1, initdata=NULL,
+ respstart=NULL, trstart=NULL, instart=NULL, ntimes=NULL, ...) {
+
+ if(is.null(data)) {
+ if(is.null(ntimes)) stop("'ntimes' must be provided if not in the data")
+ } else {
+ if(is.null(attr(data,"ntimes"))) {
+ if(is.null(ntimes)) ntimes <- nrow(data)
+ } else {
+ ntimes <- attr(data,"ntimes")
+ }
+ if(sum(ntimes)!=nrow(data)) stop("'ntimes' and data do not match")
+ }
+
+ # make response models
+ response <- makeResponseModels(response=response,data=data,nstates=nstates,family=family,values=respstart)
+
+ # make transition models
+ stationary=FALSE
+ if(transition==~1) stationary=TRUE
+ transition <- makeTransModels(nstates=nstates,formula=transition,data=data,stationary=stationary,values=trstart)
+
+ # make prior model
+ prior <- makePriorModel(nstates=nstates,ncases=length(ntimes),formula=prior,data=initdata,values=instart)
+
+ # call main depmix with all these models, ntimes and stationary
+ model <- makeDepmix(response=response,transition=transition,prior=prior,ntimes=ntimes,stationary=stationary)
+
+ # deal with starting values here!!!!!!
+
+ return(model)
+ }
+)
+
+
+
+
+
+
Copied: tags/release-0.3-0/R/depmixAIC.R (from rev 398, pkg/depmixS4/R/depmixAIC.R)
===================================================================
--- tags/release-0.3-0/R/depmixAIC.R (rev 0)
+++ tags/release-0.3-0/R/depmixAIC.R 2010-03-09 12:40:10 UTC (rev 401)
@@ -0,0 +1,13 @@
+# depends on logLik and freepars
+setMethod("AIC", signature(object="depmix"),
+ function(object, ..., k=2){
+ c(-2 * logLik(object) + freepars(object) * k)
+ }
+)
+
+# depends on logLik and freepars
+setMethod("AIC", signature(object="mix"),
+ function(object, ..., k=2){
+ c(-2 * logLik(object) + freepars(object) * k)
+ }
+)
\ No newline at end of file
Copied: tags/release-0.3-0/R/depmixBIC.R (from rev 398, pkg/depmixS4/R/depmixBIC.R)
===================================================================
--- tags/release-0.3-0/R/depmixBIC.R (rev 0)
+++ tags/release-0.3-0/R/depmixBIC.R 2010-03-09 12:40:10 UTC (rev 401)
@@ -0,0 +1,12 @@
+# depends on logLik, freepars and nobs
+setMethod("BIC", signature(object="depmix"),
+ function(object, ...){
+ c(-2 * logLik(object) + freepars(object) * log(nobs(object)))
+ }
+)
+
+setMethod("BIC", signature(object="mix"),
+ function(object, ...){
+ c(-2 * logLik(object) + freepars(object) * log(nobs(object)))
+ }
+)
Copied: tags/release-0.3-0/R/depmixfit-class.R (from rev 398, pkg/depmixS4/R/depmixfit-class.R)
===================================================================
--- tags/release-0.3-0/R/depmixfit-class.R (rev 0)
+++ tags/release-0.3-0/R/depmixfit-class.R 2010-03-09 12:40:10 UTC (rev 401)
@@ -0,0 +1,115 @@
+
+#
+# Ingmar Visser, 11-6-2008
+#
+
+# Changes
+# - added lin.upper and lin.lower slots to these objects
+
+#
+# MIX.FITTED CLASS
+#
+
+setClass("mix.fitted",
+ representation(message="character", # convergence information
+ conMat="matrix", # constraint matrix on the parameters for general linear constraints
+ lin.upper="numeric", # upper bounds for linear constraint
+ lin.lower="numeric", # lower bounds for linear constraints
+ posterior="data.frame" # posterior probabilities for the states
+ ),
+ contains="mix"
+)
+
+# accessor functions
+
+setMethod("posterior","mix.fitted",
+ function(object) {
+ return(object at posterior)
+ }
+)
+
+setMethod("show","mix.fitted",
+ function(object) {
+ cat("Convergence info:",object at message,"\n")
+ print(logLik(object))
+ cat("AIC: ", AIC(object),"\n")
+ cat("BIC: ", BIC(object),"\n")
+ }
+)
+
+setMethod("summary","mix.fitted",
+ function(object) {
+ cat("Mixture probabilities model \n")
+ print(object at prior)
+ cat("\n")
+ for(i in 1:object at nstates) {
+ cat("Response model(s) for state", i,"\n\n")
+ for(j in 1:object at nresp) {
+ cat("Response model for response",j,"\n")
+ print(object at response[[i]][[j]])
+ cat("\n")
+ }
+ cat("\n")
+ }
+ }
+)
+
+#
+# Ingmar Visser, 23-3-2008
+#
+
+#
+# DEPMIX.FITTED CLASS
+#
+
+setClass("depmix.fitted",
+ representation(message="character", # convergence information
+ conMat="matrix", # constraint matrix on the parameters for general linear constraints
+ lin.upper="numeric", # upper bounds for linear constraints
+ lin.lower="numeric", # lower bounds for linear constraints
+ posterior="data.frame" # posterior probabilities for the states
+ ),
+ contains="depmix"
+)
+
+# accessor functions
+
+setMethod("posterior","depmix.fitted",
+ function(object) {
+ return(object at posterior)
+ }
+)
+
+setMethod("show","depmix.fitted",
+ function(object) {
+ cat("Convergence info:",object at message,"\n")
+ print(logLik(object))
+ cat("AIC: ", AIC(object),"\n")
+ cat("BIC: ", BIC(object),"\n")
+ }
+)
+
+setMethod("summary","depmix.fitted",
+ function(object) {
+ cat("Initial state probabilties model \n")
+ print(object at prior)
+ cat("\n")
+ for(i in 1:object at nstates) {
+ cat("Transition model for state (component)", i,"\n")
+ print(object at transition[[i]])
+ cat("\n")
+ }
+ cat("\n")
+ for(i in 1:object at nstates) {
+ cat("Response model(s) for state", i,"\n\n")
+ for(j in 1:object at nresp) {
+ cat("Response model for response",j,"\n")
+ print(object at response[[i]][[j]])
+ cat("\n")
+ }
+ cat("\n")
+ }
+ }
+)
+
+
Copied: tags/release-0.3-0/R/depmixfit.R (from rev 398, pkg/depmixS4/R/depmixfit.R)
===================================================================
--- tags/release-0.3-0/R/depmixfit.R (rev 0)
+++ tags/release-0.3-0/R/depmixfit.R 2010-03-09 12:40:10 UTC (rev 401)
@@ -0,0 +1,131 @@
+
+setMethod("fit",
+ signature(object="mix"),
+ function(object,fixed=NULL,equal=NULL,conrows=NULL,conrows.upper=0,conrows.lower=0,method=NULL,...) {
+
+ # when there are constraints donlp should be used
+ # otherwise EM is good
+
+ # can/does EM deal with fixed constraints??? it should be possible for sure
+ if(is.null(method)) {
+ if(is.null(equal)&is.null(conrows)&is.null(fixed)) {
+ method="EM"
+ } else {
+ method="donlp"
+ }
+ }
+
+ # determine which parameters are fixed
+ if(!is.null(fixed)) {
+ if(length(fixed)!=npar(object)) stop("'fixed' does not have correct length")
+ } else {
+ if(!is.null(equal)) {
+ if(length(equal)!=npar(object)) stop("'equal' does not have correct length")
+ fixed <- !pa2conr(equal)$free
+ } else {
+ fixed <- getpars(object,"fixed")
+ }
+ }
+
+ # set those fixed parameters in the appropriate submodels
+ object <- setpars(object,fixed,which="fixed")
+
+ if(is.nan(logLik(object))) stop("Initial model infeasible, log likelihood is NaN; please provide better starting values. ")
+
+ if(method=="EM") {
+ object <- em(object,verbose=TRUE,...)
+ }
+
+ if(method=="donlp") {
+ # get the full set of parameters
+ allpars <- getpars(object)
+ # get the reduced set of parameters, ie the ones that will be optimized
+ pars <- allpars[!fixed]
+
+ # set bounds, if any
+ par.u <- rep(+Inf, length(pars))
+ par.l <- rep(-Inf, length(pars))
+
+ # make loglike function that only depends on pars
+ logl <- function(pars) {
+ allpars[!fixed] <- pars
+ object <- setpars(object,allpars)
+ -logLik(object)
+ }
+
+ if(!require(Rdonlp2)) stop("donlp optimization requires the 'Rdonlp2' package")
+
+ # make constraint matrix and its upper and lower bounds
+ lincon <- matrix(0,nr=0,nc=npar(object))
+ lin.u <- numeric(0)
+ lin.l <- numeric(0)
+
+ # incorporate equality constraints, if any
+ if(!is.null(equal)) {
+ if(length(equal)!=npar(object)) stop("'equal' does not have correct length")
+ equal <- pa2conr(equal)$conr
+ lincon <- rbind(lincon,equal)
+ lin.u <- c(lin.u,rep(0,nrow(equal)))
+ lin.l <- c(lin.l,rep(0,nrow(equal)))
+ }
+
+ # incorporate general linear constraints, if any
+ if(!is.null(conrows)) {
+ if(ncol(conrows)!=npar(object)) stop("'conrows' does not have the right dimensions")
+ lincon <- rbind(lincon,conrows)
+ if(any(conrows.upper==0)) {
+ lin.u <- c(lin.u,rep(0,nrow(conrows)))
+ } else {
+ if(length(conrows.upper)!=nrow(conrows)) stop("'conrows.upper does not have correct length")
+ lin.u <- c(lin.u,conrows.upper)
+ }
+ if(any(conrows.lower==0)) {
+ lin.l <- c(lin.l,rep(0,nrow(conrows)))
+ } else {
+ if(length(conrows.lower)!=nrow(conrows)) stop("'conrows.lower does not have correct length")
+ 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]
+
+ # set donlp2 control parameters
+ cntrl <- donlp2.control(hessian=FALSE,difftype=2,report=TRUE)
+
+ mycontrol <- function(info) {
+ return(TRUE)
+ }
+
+ # optimize the parameters
+ result <- donlp2(pars,logl,
+ par.upper=par.u,
+ par.lower=par.l,
+ A=lincon,
+ lin.upper=lin.u,
+ lin.lower=lin.l,
+ control=cntrl,
+ control.fun=mycontrol,
+ ...
+ )
+
+ if(class(object)=="depmix") class(object) <- "depmix.fitted"
+ if(class(object)=="mix") class(object) <- "mix.fitted"
+
+ object at conMat <- linconFull
+ object at message <- result$message
+ object at lin.upper <- lin.u
+ object at lin.lower <- lin.l
+
+ # put the result back into the model
+ allpars[!fixed] <- result$par
+ object <- setpars(object,allpars)
+
+ }
+
+ object at posterior <- viterbi(object)
+
+ return(object)
+ }
+)
\ No newline at end of file
Copied: tags/release-0.3-0/R/depmixsim-class.R (from rev 398, pkg/depmixS4/R/depmixsim-class.R)
===================================================================
--- tags/release-0.3-0/R/depmixsim-class.R (rev 0)
+++ tags/release-0.3-0/R/depmixsim-class.R 2010-03-09 12:40:10 UTC (rev 401)
@@ -0,0 +1,16 @@
+# Classes for simulated mix and depmix models
+
+setClass("mix.sim",
+ contains="mix",
+ representation(
+ states="matrix"
+ )
+)
+
+setClass("depmix.sim",
+ contains="depmix",
+ representation(
+ states="matrix"
+ )
+)
+
Copied: tags/release-0.3-0/R/fb.R (from rev 398, pkg/depmixS4/R/fb.R)
===================================================================
--- tags/release-0.3-0/R/fb.R (rev 0)
+++ tags/release-0.3-0/R/fb.R 2010-03-09 12:40:10 UTC (rev 401)
@@ -0,0 +1,79 @@
+#
+# Maarten Speekenbrink
+#
+# FORWARD-BACKWARD algoritme, 23-3-2008
+#
+
+fb <- function(init,A,B,ntimes=NULL,return.all=FALSE,stationary=TRUE) {
+
+ # Forward-Backward algorithm (used in Baum-Welch)
+ # Returns alpha, beta, and full data likelihood
+ # NOTE THE CHANGE IN FROM ROW TO COLUMN SUCH THAT TRANSPOSING A IS NOT NECCESSARY ANYMORE
+ # IN COMPUTING ALPHA AND BETA BUT IS NOW NECCESSARY IN COMPUTING XI
+ # A = T*K*K matrix with transition probabilities, from row to column!!!!!!!
+ # B = T*K matrix with elements ab_{ij} = P(y_i|s_j)
+ # init = K vector with initial probabilities
+
+ # NOTE: to prevent underflow, alpha and beta are scaled, using c
+
+ # NOTE: xi[t,i,j] = P(S[t] = j & S[t+1] = i)
+
+
+ B <- apply(B,c(1,3),prod)
+
+ nt <- nrow(B)
+ ns <- ncol(init)
+
+ alpha <- matrix(ncol=ns,nrow=nt)
+ beta <- matrix(ncol=ns,nrow=nt)
+ sca <- vector(length=nt)
+ xi <- array(dim=c(nt,ns,ns))
+
+ if(is.null(ntimes)) ntimes <- nt
+
+ lt <- length(ntimes)
+ et <- cumsum(ntimes)
+ bt <- c(1,et[-lt]+1)
+
+ for(case in 1:lt) {
+ alpha[bt[case],] <- init[case,]*B[bt[case],] # initialize
+ sca[bt[case]] <- 1/sum(alpha[bt[case],])
+ alpha[bt[case],] <- alpha[bt[case],]*sca[bt[case]]
+
+ if(ntimes[case]>1) {
+ for(i in bt[case]:(et[case]-1)) {
+ if(stationary) alpha[i+1,] <- (A[1,,]%*%alpha[i,])*B[i+1,]
+ else alpha[i+1,] <- (A[i,,]%*%alpha[i,])*B[i+1,]
+ sca[i+1] <- 1/sum(alpha[i+1,])
+ alpha[i+1,] <- sca[i+1]*alpha[i+1,]
+ }
+ }
+
+ beta[et[case],] <- 1*sca[et[case]] # initialize
+
+ if(ntimes[case]>1) {
+ for(i in (et[case]-1):bt[case]) {
+ if(stationary) beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[1,,]*sca[i]
+ else beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[i,,]*sca[i]
+ }
+
+ for(i in bt[case]:(et[case]-1)) {
+ if(stationary) xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[1,,])
+ else xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[i,,])
+ }
+ }
+
+ }
+
+ gamma <- alpha*beta/sca
+ like <- -sum(log(sca))
+
+ if(return.all) {
+ res <- list(alpha=alpha,beta=beta,gamma=gamma,xi=xi,sca=sca,logLike=like)
+ } else {
+ res <- list(gamma=gamma,xi=xi,logLike=like)
+ }
+
+ res
+}
+
Copied: tags/release-0.3-0/R/forwardbackward.R (from rev 398, pkg/depmixS4/R/forwardbackward.R)
===================================================================
--- tags/release-0.3-0/R/forwardbackward.R (rev 0)
+++ tags/release-0.3-0/R/forwardbackward.R 2010-03-09 12:40:10 UTC (rev 401)
@@ -0,0 +1,15 @@
+#
+# Ingmar Visser
+#
+# FORWARD-BACKWARD function, user interface, 10-06-2008
+#
+
+setMethod("forwardbackward","depmix",
+ function(object, return.all=TRUE, ...) {
+ fb(init=object at init,A=object at trDens,B=object at dens,ntimes=ntimes(object),
+ stationary=object at stationary,return.all=return.all)
+ }
+)
+
+
+
Copied: tags/release-0.3-0/R/freepars.R (from rev 398, pkg/depmixS4/R/freepars.R)
===================================================================
--- tags/release-0.3-0/R/freepars.R (rev 0)
+++ tags/release-0.3-0/R/freepars.R 2010-03-09 12:40:10 UTC (rev 401)
@@ -0,0 +1,44 @@
+# depends on getpars(object)
+setMethod("freepars","mix",
+ function(object) {
+ free <- sum(!getpars(object,which="fixed"))
+ free
+ }
+)
+
+# depends on nlin(object) and getpars(object)
+setMethod("freepars","depmix.fitted",
+ function(object) {
+ free <- sum(!getpars(object,which="fixed"))
+ free <- free-nlin(object)
+ free
+ }
+)
+
+# depends on nlin(object) and getpars(object)
+setMethod("freepars","mix.fitted",
+ function(object) {
+ free <- sum(!getpars(object,which="fixed"))
+ free <- free-nlin(object)
+ free
+ }
+)
+
+setMethod("nlin","mix.fitted",
+ function(object) {
+ conMat <- object at conMat[which(object at lin.lower==object at lin.upper),,drop=FALSE]
+ if(nrow(conMat)==0) nlin <- 0
+ else nlin <- qr(conMat)$rank
+ nlin
+ }
+)
+
+setMethod("nlin","depmix.fitted",
+ function(object) {
+ conMat <- object at conMat[which(object at lin.lower==object at lin.upper),,drop=FALSE]
+ if(nrow(conMat)==0) nlin <- 0
+ else nlin <- qr(conMat)$rank
+ nlin
+ }
+)
+
Copied: tags/release-0.3-0/R/getpars.R (from rev 398, pkg/depmixS4/R/getpars.R)
===================================================================
--- tags/release-0.3-0/R/getpars.R (rev 0)
+++ tags/release-0.3-0/R/getpars.R 2010-03-09 12:40:10 UTC (rev 401)
@@ -0,0 +1,26 @@
+setMethod("getpars","mix",
+ function(object,which="pars",...) {
+ parameters <- getpars(object at prior,which=which)
+ for(i in 1:object at nstates) {
+ for(j in 1:object at nresp) {
+ parameters <- c(parameters,getpars(object at response[[i]][[j]],which=which))
+ }
+ }
+ return(parameters)
+ }
+)
+
+setMethod("getpars","depmix",
+ function(object,which="pars",...) {
+ parameters <- getpars(object at prior,which=which)
+ for(i in 1:object at nstates) {
+ parameters <- c(parameters,getpars(object at transition[[i]],which=which))
+ }
+ for(i in 1:object at nstates) {
+ for(j in 1:object at nresp) {
+ parameters <- c(parameters,getpars(object at response[[i]][[j]],which=which))
+ }
+ }
+ return(parameters)
+ }
+)
\ No newline at end of file
Copied: tags/release-0.3-0/R/llratio.R (from rev 398, pkg/depmixS4/R/llratio.R)
===================================================================
--- tags/release-0.3-0/R/llratio.R (rev 0)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/depmix -r 401
More information about the depmix-commits
mailing list