[Depmix-commits] r451 - in tags/release-1.0-1: . R data inst inst/doc man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jan 25 15:47:41 CET 2011
Author: ingmarvisser
Date: 2011-01-25 15:47:41 +0100 (Tue, 25 Jan 2011)
New Revision: 451
Added:
tags/release-1.0-1/R/
tags/release-1.0-1/R/EM.R
tags/release-1.0-1/R/allGenerics.R
tags/release-1.0-1/R/depmix-class.R
tags/release-1.0-1/R/depmix.R
tags/release-1.0-1/R/depmixAIC.R
tags/release-1.0-1/R/depmixBIC.R
tags/release-1.0-1/R/depmixfit-class.R
tags/release-1.0-1/R/depmixfit.R
tags/release-1.0-1/R/depmixsim-class.R
tags/release-1.0-1/R/em.control.R
tags/release-1.0-1/R/fb.R
tags/release-1.0-1/R/forwardbackward.R
tags/release-1.0-1/R/freepars.R
tags/release-1.0-1/R/getpars.R
tags/release-1.0-1/R/llratio.R
tags/release-1.0-1/R/logLik.R
tags/release-1.0-1/R/lystig.R
tags/release-1.0-1/R/makeDepmix.R
tags/release-1.0-1/R/makePriorModel.R
tags/release-1.0-1/R/makeResponseModels.R
tags/release-1.0-1/R/makeTransModels.R
tags/release-1.0-1/R/mlogit.R
tags/release-1.0-1/R/multinomial.R
tags/release-1.0-1/R/nobs.R
tags/release-1.0-1/R/pa2conr.R
tags/release-1.0-1/R/response-class.R
tags/release-1.0-1/R/responseGLM.R
tags/release-1.0-1/R/responseGLMBINOM.R
tags/release-1.0-1/R/responseGLMGAMMA.R
tags/release-1.0-1/R/responseGLMMULTINOM.R
tags/release-1.0-1/R/responseGLMPOISSON.R
tags/release-1.0-1/R/responseMVN.R
tags/release-1.0-1/R/responseNORM.R
tags/release-1.0-1/R/setpars.R
tags/release-1.0-1/R/stationary.R
tags/release-1.0-1/R/transInit.R
tags/release-1.0-1/R/viterbi.R
tags/release-1.0-1/data/
tags/release-1.0-1/data/balance.rda
tags/release-1.0-1/data/speed.rda
tags/release-1.0-1/inst/
tags/release-1.0-1/inst/CITATION
tags/release-1.0-1/inst/doc/
tags/release-1.0-1/inst/doc/baldist.pdf
tags/release-1.0-1/inst/doc/depmixS4.Rnw
tags/release-1.0-1/inst/doc/depmixS4.bib
tags/release-1.0-1/inst/doc/depmixS4.pdf
tags/release-1.0-1/man/
tags/release-1.0-1/man/AIC.Rd
tags/release-1.0-1/man/GLMresponse.Rd
tags/release-1.0-1/man/balance.Rd
tags/release-1.0-1/man/depmix-class.Rd
tags/release-1.0-1/man/depmix-internal.Rd
tags/release-1.0-1/man/depmix-methods.Rd
tags/release-1.0-1/man/depmix.Rd
tags/release-1.0-1/man/depmix.fit.Rd
tags/release-1.0-1/man/depmix.fitted-class.Rd
tags/release-1.0-1/man/depmix.sim-class.Rd
tags/release-1.0-1/man/depmixS4-package.Rd
tags/release-1.0-1/man/em.control.Rd
tags/release-1.0-1/man/forwardbackward.Rd
tags/release-1.0-1/man/llratio.Rd
tags/release-1.0-1/man/makeDepmix.Rd
tags/release-1.0-1/man/mix-class.Rd
tags/release-1.0-1/man/mix.Rd
tags/release-1.0-1/man/mix.fitted-class.Rd
tags/release-1.0-1/man/mix.sim-class.Rd
tags/release-1.0-1/man/posterior.Rd
tags/release-1.0-1/man/response-class.Rd
tags/release-1.0-1/man/response-classes.Rd
tags/release-1.0-1/man/responses.Rd
tags/release-1.0-1/man/simulate.Rd
tags/release-1.0-1/man/speed.Rd
tags/release-1.0-1/man/transInit.Rd
tags/release-1.0-1/tests/
tags/release-1.0-1/tests/test1speed.R
tags/release-1.0-1/tests/test1speed.Rout.save
tags/release-1.0-1/tests/test2getsetpars.R
tags/release-1.0-1/tests/test2getsetpars.Rout.save
tags/release-1.0-1/tests/test3responses.R
Log:
Added files and subdirectories to release 1.0-1 tag.
Added: tags/release-1.0-1/R/EM.R
===================================================================
--- tags/release-1.0-1/R/EM.R (rev 0)
+++ tags/release-1.0-1/R/EM.R 2011-01-25 14:47:41 UTC (rev 451)
@@ -0,0 +1,233 @@
+#
+# Maarten Speekenbrink 23-3-2008
+#
+
+em <- function(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-8,crit="relative",random.start=TRUE,verbose=FALSE,...) {
+
+ if(!is(object,"mix")) stop("object is not of class 'mix'")
+
+ ns <- nstates(object)
+ ntimes <- ntimes(object)
+ lt <- length(ntimes)
+ et <- cumsum(ntimes)
+ bt <- c(1,et[-lt]+1)
+
+ converge <- FALSE
+ j <- 0
+
+ # compute response probabilities
+ B <- apply(object at dens,c(1,3),prod)
+ gamma <- object at init*B
+ LL <- sum(log(rowSums(gamma)))
+ # normalize
+ gamma <- gamma/rowSums(gamma)
+
+ if(random.start) {
+ nr <- sum(ntimes(object))
+ gamma <- matrix(runif(nr*ns,min=.0001,max=.9999),nr=nr,nc=ns)
+ gamma <- gamma/rowSums(gamma)
+ # add stuff to reestimate the model here and compute LL again
+ # based on these starting values
+ }
+
+ 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) {
+ if((crit == "absolute" && LL - LL.old < tol) || (crit == "relative" && (LL.old - LL)/LL.old < tol)) {
+ cat("iteration",j,"logLik:",LL,"\n")
+ converge <- TRUE
+ }
+ } else {
+ # this should not really happen...
+ if(j > 0) warning("likelihood decreased on iteration",j)
+ }
+
+ LL.old <- LL
+ j <- j+1
+
+ }
+
+ class(object) <- "mix.fitted"
+
+ if(converge) {
+ object at message <- switch(crit,
+ relative = "Log likelihood converged to within tol. (relative change)",
+ absolute = "Log likelihood converged to within tol. (absolute change)"
+ )
+ } else object at message <- "'maxit' iterations reached in EM without convergence."
+
+ # no constraints in EM, except for the standard constraints ...
+ # which are produced by the following (only necessary for getting df right in logLik and such)
+ constraints <- getConstraints(object)
+ object at conMat <- constraints$lincon
+ object at lin.lower <- constraints$lin.l
+ object at lin.upper <- constraints$lin.u
+
+ object
+
+}
+
+# em for hidden markov models
+em.depmix <- function(object,maxit=100,tol=1e-8,crit="relative",random.start=TRUE,verbose=FALSE,...) {
+
+ if(!is(object,"depmix")) stop("object is not of class 'depmix'")
+
+ ns <- nstates(object)
+
+ ntimes <- ntimes(object)
+ lt <- length(ntimes)
+ et <- cumsum(ntimes)
+ bt <- c(1,et[-lt]+1)
+
+ converge <- FALSE
+ j <- 0
+
+
+ if(random.start) {
+
+ nr <- sum(ntimes(object))
+ gamma <- matrix(runif(nr*ns,min=.0001,max=.9999),nr=nr,nc=ns)
+ gamma <- gamma/rowSums(gamma)
+ LL <- -1e10
+
+ 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]])
+ }
+ }
+
+ # 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
+
+ if(is.nan(LL)) stop("Cannot find suitable starting values; please provide them.")
+
+ } else {
+ # 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, gamma)
+ 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(!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 specific trinModels
+ # should become object at transition = fit(object at transition, xi, gamma)
+ object at transition[[i]]@parameters$coefficients <- switch(object at transition[[i]]@family$link,
+ identity = object at transition[[i]]@family$linkfun(trm[i,]),
+ mlogit = object at transition[[i]]@family$linkfun(trm[i,],base=object at transition[[i]]@family$base),
+ object at transition[[i]]@family$linkfun(trm[i,])
+ )
+ }
+ # 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)) {
+ if((crit == "absolute" && LL - LL.old < tol) || (crit == "relative" && (LL.old - LL)/LL.old < tol)) {
+ cat("iteration",j,"logLik:",LL,"\n")
+ converge <- TRUE
+ }
+ } else {
+ # this should not really happen...
+ if(j > 0) warning("likelihood decreased on iteration",j)
+ }
+
+ LL.old <- LL
+ j <- j+1
+
+ }
+
+ class(object) <- "depmix.fitted"
+
+ if(converge) {
+ object at message <- switch(crit,
+ relative = "Log likelihood converged to within tol. (relative change)",
+ absolute = "Log likelihood converged to within tol. (absolute change)"
+ )
+ } else object at message <- "'maxit' iterations reached in EM without convergence."
+
+ # no constraints in EM, except for the standard constraints ...
+ # which are produced by the following (only necessary for getting df right in logLik and such)
+ constraints <- getConstraints(object)
+ object at conMat <- constraints$lincon
+ object at lin.lower <- constraints$lin.l
+ object at lin.upper <- constraints$lin.u
+
+ object
+}
Property changes on: tags/release-1.0-1/R/EM.R
___________________________________________________________________
Added: svn:eol-style
+ native
Added: tags/release-1.0-1/R/allGenerics.R
===================================================================
--- tags/release-1.0-1/R/allGenerics.R (rev 0)
+++ tags/release-1.0-1/R/allGenerics.R 2011-01-25 14:47:41 UTC (rev 451)
@@ -0,0 +1,80 @@
+
+#
+# Ingmar Visser, 23-3-2008
+#
+
+.First.lib <- function(lib, pkg) {
+ require(stats)
+ require(methods)
+ require(MASS)
+ require(nnet)
+ require(Rsolnp)
+}
+
+.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("getModel", function(object, ...) standardGeneric("getModel"))
+
+# setGeneric("logLik", function(object, ...) standardGeneric("logLik"))
+
+setGeneric("fit", function(object, ...) standardGeneric("fit"))
+
+setGeneric("getConstraints", function(object, ...) standardGeneric("getConstraints"))
+
+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"))
+
Property changes on: tags/release-1.0-1/R/allGenerics.R
___________________________________________________________________
Added: svn:eol-style
+ native
Added: tags/release-1.0-1/R/depmix-class.R
===================================================================
--- tags/release-1.0-1/R/depmix-class.R (rev 0)
+++ tags/release-1.0-1/R/depmix-class.R 2011-01-25 14:47:41 UTC (rev 451)
@@ -0,0 +1,306 @@
+
+#
+# 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)
+ 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))
+
+ 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) {
+ 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)
+ }
+)
+
+# setMethod("getModel",signature(object="mix"),
+# function(object,which="response",...) {
+# res <- switch(which,
+# "prior"=object at prior,
+# "response"=object at response)
+# res
+# }
+# )
+
+#
+# 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)
+ }
+)
+
+# setMethod("getModel",signature(object="depmix"),
+# function(object,which="response",...) {
+# res <- switch(which,
+# "prior"=object at prior,
+# "response"=object at response,
+# "transition"=object at transition)
+# res
+# }
+# )
+
+#
+# SUMMARY method: to do
+#
+
+
+
Property changes on: tags/release-1.0-1/R/depmix-class.R
___________________________________________________________________
Added: svn:eol-style
+ native
Added: tags/release-1.0-1/R/depmix.R
===================================================================
--- tags/release-1.0-1/R/depmix.R (rev 0)
+++ tags/release-1.0-1/R/depmix.R 2011-01-25 14:47:41 UTC (rev 451)
@@ -0,0 +1,94 @@
+#
+# 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)
+
+ 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)
+})
Property changes on: tags/release-1.0-1/R/depmix.R
___________________________________________________________________
Added: svn:eol-style
+ native
Added: tags/release-1.0-1/R/depmixAIC.R
===================================================================
--- tags/release-1.0-1/R/depmixAIC.R (rev 0)
+++ tags/release-1.0-1/R/depmixAIC.R 2011-01-25 14:47:41 UTC (rev 451)
@@ -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
Property changes on: tags/release-1.0-1/R/depmixAIC.R
___________________________________________________________________
Added: svn:eol-style
+ native
Added: tags/release-1.0-1/R/depmixBIC.R
===================================================================
--- tags/release-1.0-1/R/depmixBIC.R (rev 0)
+++ tags/release-1.0-1/R/depmixBIC.R 2011-01-25 14:47:41 UTC (rev 451)
@@ -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)))
+ }
+)
Property changes on: tags/release-1.0-1/R/depmixBIC.R
___________________________________________________________________
Added: svn:eol-style
+ native
Added: tags/release-1.0-1/R/depmixfit-class.R
===================================================================
--- tags/release-1.0-1/R/depmixfit-class.R (rev 0)
+++ tags/release-1.0-1/R/depmixfit-class.R 2011-01-25 14:47:41 UTC (rev 451)
@@ -0,0 +1,138 @@
+
+#
+# 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,which="all") {
+ ans=switch(which,
+ "all" = 1,
+ "response" = 2,
+ "prior" = 3,
+ stop("Invalid 'which' argument in summary of fitted mix model")
+ )
+ if(ans==1|ans==3) {
+ cat("Mixture probabilities model \n")
+ print(object at prior)
+ cat("\n")
+ }
+ if(ans==1|ans==2) {
+ 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,which="all") {
+ ans=switch(which,
+ "all" = 1,
+ "response" = 2,
+ "prior" = 3,
+ "transition" = 4,
+ stop("Invalid 'which' argument in summary of fitted depmix model")
+ )
+ if(ans==1|ans==3) {
+ cat("Initial state probabilties model \n")
+ print(object at prior)
+ cat("\n")
+ }
+ if(ans==1|ans==4) {
+ for(i in 1:object at nstates) {
+ cat("Transition model for state (component)", i,"\n")
+ print(object at transition[[i]])
+ cat("\n")
+ }
+ cat("\n")
+ }
+ if(ans==1|ans==2) {
+ 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")
+ }
+ }
+ }
+)
+
+
Property changes on: tags/release-1.0-1/R/depmixfit-class.R
___________________________________________________________________
Added: svn:eol-style
+ native
Added: tags/release-1.0-1/R/depmixfit.R
===================================================================
--- tags/release-1.0-1/R/depmixfit.R (rev 0)
+++ tags/release-1.0-1/R/depmixfit.R 2011-01-25 14:47:41 UTC (rev 451)
@@ -0,0 +1,298 @@
+
+
+setMethod("fit",
+ signature(object="mix"),
+ function(object,fixed=NULL,equal=NULL,conrows=NULL,conrows.upper=0,conrows.lower=0,method=NULL,emcontrol=em.control(),verbose=TRUE,...) {
+
+ fi <- !is.null(fixed)
+ cr <- !is.null(conrows)
+ eq <- !is.null(equal)
+
+ constr <- any(c(fi,cr,eq))
+
+ # when there are constraints donlp/solnp should be used
+ # otherwise EM is good
+ if(is.null(method)) {
+ if(constr) {
+ method="rsolnp"
+ } else {
+ method="EM"
+ }
+ } else {
+ if(method=="EM") {
+ if(constr) {
+ warning("EM not applicable for constrained models; optimization method changed to 'rsolnp'")
+ method="rsolnp"
+ }
+ }
+ }
+
+ if(method=="EM") {
+ object <- em(object,maxit=emcontrol$maxit,tol=emcontrol$tol,crit=emcontrol$crit,random.start=emcontrol$random.start,verbose=verbose,...)
+ }
+
+ if(method=="donlp"||method=="rsolnp") {
+
+ # check feasibility of starting values
+ if(is.nan(logLik(object))) stop("Initial model infeasible, log likelihood is NaN; please provide better starting values. ")
+
+ # determine which parameters are fixed
+ if(fi) {
+ if(length(fixed)!=npar(object)) stop("'fixed' does not have correct length")
+ } else {
+ if(eq) {
+ 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")
+
+ # 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]
+
+ constraints <- getConstraints(object)
+
+ lincon=constraints$lincon
+ lin.u=constraints$lin.u
+ lin.l=constraints$lin.l
+ par.u=constraints$par.u
+ par.l=constraints$par.l
+
+ # incorporate equality constraints provided with the fit function, if any
+ if(eq) {
+ 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(cr) {
+ 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]
+
+ # remove redundant rows in lincon (all zeroes)
+ allzero <- which(apply(lincon,1,function(y) all(y==0)))
+ if(length(allzero)>0) {
+ lincon <- lincon[-allzero,,drop=FALSE]
+ lin.u <- lin.u[-allzero]
+ lin.l <- lin.l[-allzero]
+ }
+
+ # make loglike function that only depends on pars
+ logl <- function(pars) {
+ allpars[!fixed] <- pars
+ object <- setpars(object,allpars)
+ ans = -as.numeric(logLik(object))
+ if(is.na(ans)) ans = 100000
+ ans
+ }
+
+ if(method=="donlp") {
+
+ reqdon <- require(Rdonlp2,quietly=TRUE)
+
+ if(!reqdon) stop("Rdonlp2 not available.")
+
+ # 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[!fixed],
+ par.lower=par.l[!fixed],
+ 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"
+
+ # convergence info
+ object at message <- result$message
+
+ # put the result back into the model
+ allpars[!fixed] <- result$par
+ object <- setpars(object,allpars)
+ }
+
+ if(method=="rsolnp") {
+
+ if(!(require(Rsolnp,quietly=TRUE))) stop("Method 'rsolnp' requires package 'Rsolnp'")
+
+ # separate equality and inequality constraints
+ ineq <- which(lin.u!=lin.l)
+ if(length(ineq)>0) lineq <- lincon[-ineq, ,drop=FALSE]
+ else lineq <- lincon
+
+ # returns the evaluated equality constraints
+ if(nrow(lineq)>0) {
+ eqfun <- function(pp) {
+ ans = as.vector(lineq%*%pp)
+ ans
+ }
+ # select the boundary values for the equality constraints
+ if(length(ineq)>0) lineq.bound = lin.l[-ineq]
+ else lineq.bound = lin.l
+ } else {
+ eqfun=NULL
+ lineq.bound=NULL
+ }
+
+ # select the inequality constraints
+ if(length(ineq)>0) {
+ linineq <- lincon[ineq, ,drop=FALSE]
+ ineqLB <- lin.l[ineq]
+ ineqUB <- lin.u[ineq]
+ } else {
+ ineqfun = NULL
+ ineqLB=NULL
+ ineqUB=NULL
+ }
+
+ # call to solnp
+ res <- solnp(pars,
+ logl,
+ eqfun = eqfun,
+ eqB = lineq.bound,
+ ineqfun = ineqfun,
+ ineqLB = ineqLB,
+ ineqUB = ineqUB,
+ LB = par.l[!fixed],
+ UB = par.u[!fixed],
+ control = list(delta = 1e-5, tol = 1e-6, trace = 1),
+ ...
+ )
+
+ if(class(object)=="depmix") class(object) <- "depmix.fitted"
+ if(class(object)=="mix") class(object) <- "mix.fitted"
+
+ object at message <- c(res$convergence," (0 is good in Rsolnp, check manual for other values)")
+
+ # put the result back into the model
+ allpars[!fixed] <- res$pars
+ object <- setpars(object,allpars)
+
+ }
+
+ object at conMat <- linconFull
+ object at lin.upper <- lin.u
+ object at lin.lower <- lin.l
+
+ }
+
+ object at posterior <- viterbi(object)
+
+ return(object)
+ }
+)
+
+
+setMethod("getConstraints",
+ signature(object="mix"),
+ function(object) {
+
+ # set bounds, if any (should add bounds for eg sd parameters at some point ...)
+ par.u <- rep(+Inf, npar(object))
+ par.l <- rep(-Inf, npar(object))
+
+ # 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)
+
+ ns <- nstates(object)
+ nrsp <- nresp(object)
+
+ # get bounds from submodels
+ # get internal linear constraints from submodels
+
+ # first for the prior model
+ bp <- 1
+ ep <- npar(object at prior)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/depmix -r 451
More information about the depmix-commits
mailing list