[Depmix-commits] r562 - in tags: . release-1.2-1 release-1.2-1/R release-1.2-1/data release-1.2-1/inst release-1.2-1/inst/doc release-1.2-1/man release-1.2-1/src release-1.2-1/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Aug 9 11:21:21 CEST 2012
Author: maarten
Date: 2012-08-09 11:21:21 +0200 (Thu, 09 Aug 2012)
New Revision: 562
Added:
tags/release-1.2-1/
tags/release-1.2-1/DESCRIPTION
tags/release-1.2-1/NAMESPACE
tags/release-1.2-1/NEWS
tags/release-1.2-1/R/
tags/release-1.2-1/R/EM.R
tags/release-1.2-1/R/allGenerics.R
tags/release-1.2-1/R/depmix-class.R
tags/release-1.2-1/R/depmix.R
tags/release-1.2-1/R/depmixfit-class.R
tags/release-1.2-1/R/depmixfit.R
tags/release-1.2-1/R/depmixsim-class.R
tags/release-1.2-1/R/em.control.R
tags/release-1.2-1/R/fb.R
tags/release-1.2-1/R/forwardbackward.R
tags/release-1.2-1/R/freepars.R
tags/release-1.2-1/R/getpars.R
tags/release-1.2-1/R/llratio.R
tags/release-1.2-1/R/logLik.R
tags/release-1.2-1/R/lystig.R
tags/release-1.2-1/R/makeDepmix.R
tags/release-1.2-1/R/makePriorModel.R
tags/release-1.2-1/R/makeResponseModels.R
tags/release-1.2-1/R/makeTransModels.R
tags/release-1.2-1/R/mlogit.R
tags/release-1.2-1/R/multinomial.R
tags/release-1.2-1/R/nobs.R
tags/release-1.2-1/R/pa2conr.R
tags/release-1.2-1/R/response-class.R
tags/release-1.2-1/R/responseGLM.R
tags/release-1.2-1/R/responseGLMBINOM.R
tags/release-1.2-1/R/responseGLMGAMMA.R
tags/release-1.2-1/R/responseGLMMULTINOM.R
tags/release-1.2-1/R/responseGLMPOISSON.R
tags/release-1.2-1/R/responseMVN.R
tags/release-1.2-1/R/responseNORM.R
tags/release-1.2-1/R/setpars.R
tags/release-1.2-1/R/stationary.R
tags/release-1.2-1/R/transInit.R
tags/release-1.2-1/R/viterbi.R
tags/release-1.2-1/README
tags/release-1.2-1/data/
tags/release-1.2-1/data/balance.rda
tags/release-1.2-1/data/speed.rda
tags/release-1.2-1/inst/
tags/release-1.2-1/inst/CITATION
tags/release-1.2-1/inst/doc/
tags/release-1.2-1/inst/doc/baldist.pdf
tags/release-1.2-1/inst/doc/depmixS4.Rnw
tags/release-1.2-1/inst/doc/depmixS4.bib
tags/release-1.2-1/inst/doc/depmixS4.pdf
tags/release-1.2-1/man/
tags/release-1.2-1/man/GLMresponse.Rd
tags/release-1.2-1/man/balance.Rd
tags/release-1.2-1/man/depmix-class.Rd
tags/release-1.2-1/man/depmix-internal.Rd
tags/release-1.2-1/man/depmix-methods.Rd
tags/release-1.2-1/man/depmix.Rd
tags/release-1.2-1/man/depmix.fit.Rd
tags/release-1.2-1/man/depmix.fitted-class.Rd
tags/release-1.2-1/man/depmix.sim-class.Rd
tags/release-1.2-1/man/depmixS4-package.Rd
tags/release-1.2-1/man/em.control.Rd
tags/release-1.2-1/man/forwardbackward.Rd
tags/release-1.2-1/man/llratio.Rd
tags/release-1.2-1/man/makeDepmix.Rd
tags/release-1.2-1/man/mix-class.Rd
tags/release-1.2-1/man/mix.Rd
tags/release-1.2-1/man/mix.fitted-class.Rd
tags/release-1.2-1/man/mix.sim-class.Rd
tags/release-1.2-1/man/posterior.Rd
tags/release-1.2-1/man/response-class.Rd
tags/release-1.2-1/man/response-classes.Rd
tags/release-1.2-1/man/responses.Rd
tags/release-1.2-1/man/simulate.Rd
tags/release-1.2-1/man/speed.Rd
tags/release-1.2-1/man/transInit.Rd
tags/release-1.2-1/src/
tags/release-1.2-1/src/fb.cc
tags/release-1.2-1/src/fb.h
tags/release-1.2-1/tests/
tags/release-1.2-1/tests/test1speed.R
tags/release-1.2-1/tests/test1speed.Rout.save
tags/release-1.2-1/tests/test2getsetpars.R
tags/release-1.2-1/tests/test2getsetpars.Rout.save
tags/release-1.2-1/tests/test3responses.R
Log:
fixed bugs in handling of missing data, to be released as version 1.2-1
Added: tags/release-1.2-1/DESCRIPTION
===================================================================
--- tags/release-1.2-1/DESCRIPTION (rev 0)
+++ tags/release-1.2-1/DESCRIPTION 2012-08-09 09:21:21 UTC (rev 562)
@@ -0,0 +1,12 @@
+Package: depmixS4
+Version: 1.2-1
+Date: 2012-08-09
+Title: Dependent Mixture Models - Hidden Markov Models of GLMs and Other Distributions in S4
+Author: Ingmar Visser <i.visser at uva.nl>, Maarten Speekenbrink <m.speekenbrink at ucl.ac.uk>
+Maintainer: Ingmar Visser <i.visser at uva.nl>
+Depends: R (>= 2.15.0), stats, nnet, methods, MASS, Rsolnp, stats4
+Suggests: gamlss, gamlss.dist, TTR
+Description: Fit latent (hidden) Markov models on mixed categorical and continuous (timeseries)
+ data, otherwise known as dependent mixture models
+License: GPL (>=2)
+URL: http://depmix.r-forge.r-project.org/
Property changes on: tags/release-1.2-1/DESCRIPTION
___________________________________________________________________
Added: svn:executable
+ *
Added: tags/release-1.2-1/NAMESPACE
===================================================================
--- tags/release-1.2-1/NAMESPACE (rev 0)
+++ tags/release-1.2-1/NAMESPACE 2012-08-09 09:21:21 UTC (rev 562)
@@ -0,0 +1,61 @@
+import(methods)
+
+importFrom(stats, predict, simulate)
+
+importFrom(stats4, AIC, BIC, logLik, nobs, summary)
+
+export(
+ makeDepmix,
+ makeMix,
+ lystig,
+ fb,
+ forwardbackward,
+ MVNresponse,
+ llratio,
+ multinomial,
+ em,
+ em.control,
+ viterbi,
+ mlogit,
+ logLik
+)
+
+exportClasses(
+ depmix,
+ depmix.sim,
+ mix,
+ mix.sim,
+ depmix.fitted,
+ mix.fitted,
+ response,
+ GLMresponse,
+ MVNresponse,
+ transInit
+)
+
+exportMethods(
+ fit,
+ getConstraints,
+ npar,
+ freepars,
+ nlin,
+ getdf,
+ nobs,
+ nresp,
+ ntimes,
+ nstates,
+ depmix,
+ mix,
+ posterior,
+ GLMresponse,
+ MVNresponse,
+ transInit,
+ setpars,
+ getpars,
+ predict,
+ dens,
+ show,
+ simulate,
+ summary,
+ logLik
+)
Added: tags/release-1.2-1/NEWS
===================================================================
--- tags/release-1.2-1/NEWS (rev 0)
+++ tags/release-1.2-1/NEWS 2012-08-09 09:21:21 UTC (rev 562)
@@ -0,0 +1,229 @@
+Changes in depmixS4 version 1.2-1
+
+ o Fixed a bug in handling of missing values for mix models
+
+Changes in depmixS4 version 1.2-0
+
+Major changes
+
+ o Missing values for responses are now allowed. Note that missing values
+ in covariates will cause errors.
+
+Changes in depmixS4 version 1.1-0
+
+ Major changes
+
+ o The main loop computing the forward and backward variables for use in
+ the EM algorithm is now implemented in C. Depending on model specifics
+ this results in a 2-4 fold speed increase when fitting models.
+
+ o The Changes for each release (in the NEWS file) is now split into two
+ sections: Major and Minor changes.
+
+ o Added several examples on the ?responses page (Poisson change point
+ model, similar model with a single multinomial response) and the
+ ?depmix page (model for S&P 500 returns; thanks to Chen Haibo for
+ sending this).
+
+ Minor changes
+
+ o Corrected a typo in the vignette in Equation 1; the first occurrence
+ of S read S_t instead of S_1 (thanks to Peng Yu for reporting this).
+
+ o Added a sensible error message when the data contains missings (depmixS4
+ can not handle missing data yet).
+
+ o Fixed a bug in the relative stopping criterion for EM (which resulted
+ in immediate indication of convergence for positive log likelihoods;
+ thanks to Chen Haibo for sending the S&P 500 example which brought out
+ this problem).
+
+ o Function forwardbackward now has a useC argument to determine whether
+ C code is used, the default, or not (the R code is mostly left in place
+ for easy debugging).
+
+ o Added a fix for models without covariates/intercepts. In responseGLM
+ and responseMVN the function setpars now exits when length(value) == 0.
+ In setpars.depmix, a check is added whether npar > 0.
+
+Changes in depmixS4 version 1.0-4
+
+ o Added examples of the use of ntimes argument on ?depmix and ?fit
+ help pages using the ?speed data (which now has the full reference
+ to the accompanying publication).
+
+ o Using nobs generic from stats/stats4 rather than defining them
+ anew (which gave clashes with other packages that did the same).
+
+ o Fixed a bug in simulation of gaussian response model, which was
+ returning NaNs due to an error in assignment of the sd parameter
+ (introduced in version x). Thanks to Jeffrey Arnold for reporting
+ this (bug #1365).
+
+Changes in depmixS4 version 1.0-3
+
+ o Using AIC/BIC/logLik generics from stats/stats4 rather than
+ defining them anew (which gave clashes with other packages that did
+ the same).
+
+Changes in depmixS4 version 1.0-2
+
+ o fixed a bug in simulation of binomial response model data (the response
+ consists of the number of successes, and the number of failures; in
+ simulation, the number of failures was an exact copy of the number of
+ successes).
+
+ o added a meaningful error message in the EM algorithm for lca/mixture
+ models in case the initial log likelihood is NA (thanks to Matthias
+ Ihrke for pointing this out).
+
+Changes in depmixS4 version 1.0-1
+
+ o minor changes in documentation to conform to R 2.12.0 standards.
+
+ o fixed a bug concerning random start values (the argument to specify
+ this was not passed to the EM algorithm and hence was completely
+ ineffective ...).
+
+ o changed the emcontrol argument to the fit function; it now calls
+ a function em.control which returns the list of control parameters, which
+ now also includes maxit, the max number of iterations of the EM algorithm.
+ This makes future additions to EM control parameters easier.
+
+ o random parameter initialization is now the default when using EM
+ to fit models.
+
+ o fixed a bug in multinomial models with n>1; the parameters are now
+ normalized such that they sum to unity (this bug was introduced in
+ version 0.9-0 in multinomial models with identity link).
+
+ o added an error message for multinomial response models with n>1 and
+ link='mlogit' as this case is not handled; n>1 multinomial can use the
+ 'identity' link function.
+
+Changes in depmixS4 version 1.0-0
+
+ o added a vignette to the package and upped the version number 1.0-0 to
+ celebrate publication in the Journal of Statistical Software.
+
+Changes in depmixS4 version 0.9-1
+
+ o fixed a bug in setting the lower and upper bounds for GLMresponse
+ models (the number of bounds was wrong for models with covariates/
+ predictors; these bounds are only used in constrained optimization in
+ which case they produced an error immediately; in EM optimization these
+ bounds are not used).
+
+Changes in depmixS4 version 0.9-0
+
+ o added optimization using Rsolnp, which can be invoked by using
+ method="rsolnp" in calling fit on (dep-)mix objects. Note that
+ this is meant for fitting models with additional constraints.
+ Method="rsolnp" is now the default when fitting constrained
+ models, method="donlp" is still supported.
+
+ o added documentation for control arguments that can be passed to
+ em algorithm, particularly for controlling the tolerance in
+ optimization.
+
+ o added multinomial models with identity link for transition and prior
+ probabilities. These are now the default when no covariates are
+ present.
+
+ o added bounds and constraints for multinomial identity models such
+ that these constraints are satisfied when fitting models with
+ method="rsolnp" or "donlp". Also, variance and sd parameters in
+ gaussian and multivariate normal models are given bounds to prevent
+ warnings and errors in optimization of such models using rsolnp or
+ donlp.
+
+ o added option to generate starting values as part of the EM
+ algorithm.
+
+ o fixed a bug in multinomial response models with n>1; the response for
+ these models can be specified as a k-column matrix with the number of
+ observed responses for each category; the loglikelihood for these
+ models in which there was more than 1 observation per row was
+ incorrect; note that such models may lead to some numerical
+ instabilities when n is large.
+
+Changes in depmixS4 version 0.3-0
+
+ o added multinomial response function with identity link (no covariates
+ allowed in such a model); useful when (many) boundary values occur;
+ currently no constraints are used for such models, and hence only EM
+ can be used for optimization, or alternatively, if and when Rdonlp2
+ is used, sum constraints need to be added when fitting the model.
+ See ?GLMresponse for details.
+
+ o added an example of how to specify a model with multivariate normal
+ responses (and fixed a bug in MVNresponse that prevented such models
+ from being specified in the first place). See ?makeDepmix for an
+ example.
+
+Changes in depmixS4 version 0.2-2
+
+ o fixed a warning produced when specifying conrows.upper and .lower in
+ the fit function
+
+ o added error message in case the initial log likelihood is infeasible
+
+ o fixed a bug in the fit function for multinomial response models with
+ covariates (thanks to Gilles Dutilh for spotting this)
+
+Changes in depmixS4 version 0.2-1
+
+ o fixed a bug in the Viterbi algorithm used to compute posterior states
+ (this bug was introduced in version 0.2-0)
+
+ o restructured test files somewhat
+
+ o fixed a bug in the use of the conrows argument in the fit function (a
+ missing drop=FALSE statement)
+
+ o updated help files for mix classes
+
+ o fixed a bug in setting the starting values of regression coefficients in
+ prior and transInit models with covariates (thanks to Verena Schmittmann
+ for reporting this)
+
+ o added newx argument to predict function of transInit objects, to be used
+ for predicting probabilities depending on covariates (useful in eg plotting
+ transition probabilities as function of a covariate)
+
+ o added example of the use of conrows argument in fitting functions and other
+ minor updates in documentation
+
+Changes in depmixS4 version 0.2-0
+
+ o restructured R and Rd (help) files; added depmixS4 help with a short
+ overview of the package and links to appropriate help files
+
+ o added function 'simulate' to generate new data from a (fitted) model
+
+ o added function 'forwardbackward' to access the forward and backward
+ variables as well as the smoothed transition and state variables
+
+ o added new glm distributions: gamma, poisson
+
+ o added multivariate normal distribution
+
+ o freepars now works correctly on both depmix and depmix.fitted objects
+
+ o added function 'nlin' to compute the number of linear constraints in
+ a fitted model object
+
+ o added mix class for mixture and latent class models; the depmix class
+ extends this mix class and adds a transition model to it
+
+ o added help file for makeDepmix to provide full control in specifying
+ models
+
+ o minor changes to make depmixS4 compatible with R 2.7.1
+
+
+Changes in depmixS4 version 0.1-1
+
+ o adjusted for R 2.7.0
+
+First version released on CRAN: 0.1-0
Added: tags/release-1.2-1/R/EM.R
===================================================================
--- tags/release-1.2-1/R/EM.R (rev 0)
+++ tags/release-1.2-1/R/EM.R 2012-08-09 09:21:21 UTC (rev 562)
@@ -0,0 +1,295 @@
+#
+# Maarten Speekenbrink 23-3-2008
+#
+
+rdirichlet <- function(n, alpha) {
+ # taken from gtools...
+ l <- length(alpha)
+ x <- matrix(rgamma(l * n, alpha), ncol = l, byrow = TRUE)
+ sm <- x %*% rep(1, l)
+ x/as.vector(sm)
+}
+
+which.is.max <- function(x) {
+ # taken from MASS
+ y <- seq_along(x)[x == max(x)]
+ if(length(y) > 1L) sample(y, 1L) else y
+}
+
+
+ind.max <- function(x) {
+ out <- rep(0,length(x))
+ out[which.is.max(x)] <- 1
+ out
+}
+
+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,classification=c("soft","hard"),...) {
+
+ clsf <- match.arg(classification)
+
+ 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
+
+ if(random.start) {
+
+ nr <- sum(ntimes(object))
+ gamma <- rdirichlet(nr,alpha=rep(.01,ns))
+
+ if(clsf == "hard") {
+ gamma <- t(apply(gamma,1,ind.max))
+ }
+
+ 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,matrix(0,1,1),B=object at dens,ntimes=ntimes(object))
+ 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=matrix(0,1,1),B=object at dens,ntimes=ntimes(object))
+ LL <- fbo$logLike
+ if(is.nan(LL)) stop("Starting values not feasible; please provide them.")
+
+ #B <- apply(object at dens,c(1,3),prod)
+ #gamma <- object at init*B
+ #LL <- sum(log(rowSums(gamma)))
+ #if(is.nan(LL)) stop("Starting values not feasible; please provide them.")
+ #gamma <- gamma/rowSums(gamma)
+ }
+
+ LL.old <- LL + 1
+
+ while(j <= maxit & !converge) {
+
+ # maximization
+
+ # should become object at prior <- fit(object at prior)
+
+ if(clsf == "hard") {
+ fbo$gamma <- t(apply(gamma,1,ind.max))
+ }
+ 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)
+
+ 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=matrix(0,1,1),B=object at dens,ntimes=ntimes(object))
+ LL <- fbo$logLike
+
+ #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 - LL.old)/abs(LL.old) < tol)) {
+ cat("iteration",j,"logLik:",LL,"\n")
+ converge <- TRUE
+ }
+ } else {
+ # this should not really happen...
+ if(j > 0 && (LL.old - LL) > tol) 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, classification=c("soft","hard"),...) {
+
+ if(!is(object,"depmix")) stop("object is not of class 'depmix'")
+
+ clsf <- match.arg(classification)
+
+ clsf="soft"
+
+ 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),nrow=nr,ncol=ns)
+ #gamma <- gamma/rowSums(gamma)
+ gamma <- rdirichlet(nr,alpha=rep(.01,ns))
+ if(clsf == "hard") {
+ gamma <- t(apply(gamma,1,ind.max))
+ }
+ 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
+ if(is.nan(LL)) stop("Starting values not feasible; please provide them.")
+ }
+
+ 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 - LL.old)/abs(LL.old) < tol)) {
+ cat("iteration",j,"logLik:",LL,"\n")
+ converge <- TRUE
+ }
+ } else {
+ # this should not really happen...
+ if(j > 0 && (LL.old - LL) > tol) 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
+}
Added: tags/release-1.2-1/R/allGenerics.R
===================================================================
--- tags/release-1.2-1/R/allGenerics.R (rev 0)
+++ tags/release-1.2-1/R/allGenerics.R 2012-08-09 09:21:21 UTC (rev 562)
@@ -0,0 +1,79 @@
+
+#
+# Ingmar Visser, 23-3-2008
+#
+
+# 17-6-2011: added dynamic lib statement to include the C code
+# version of forward backward routine
+
+.onLoad <- function(lib, pkg) {
+ library.dynam("depmixS4", pkg, lib)
+}
+
+.onUnLoad <- function(libpath) {
+ library.dynam.unload("depmixS4",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"))
+
+# extractors/set functions
+
+setGeneric("npar", function(object, ...) standardGeneric("npar"))
+
+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("getdf",function(object) standardGeneric("getdf"))
+
+setGeneric("nlin", function(object, ...) standardGeneric("nlin"))
+
+setGeneric("getConstraints", function(object, ...) standardGeneric("getConstraints"))
+
+setGeneric("is.stationary", function(object,...) standardGeneric("is.stationary"))
+
+setGeneric("setpars", function(object,values,which="pars",...) standardGeneric("setpars"))
+
+setGeneric("getpars", function(object,which="pars",...) standardGeneric("getpars"))
+
+
+# functions
+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("logDens",function(object,...) standardGeneric("logDens"))
+
+setGeneric("dens",function(object,...) standardGeneric("dens"))
+
+setGeneric("predict", function(object, ...) standardGeneric("predict"))
+
+
+# redundant??
+
+# setGeneric("getModel", function(object, ...) standardGeneric("getModel"))
+
+# these are imported from stats4
+# setGeneric("nobs", function(object, ...) standardGeneric("nobs"))
+# setGeneric("logLik", function(object, ...) standardGeneric("logLik"))
+# setGeneric("AIC", function(object, ..., k=2) standardGeneric("AIC"))
+# setGeneric("BIC", function(object, ...) standardGeneric("BIC"))
Added: tags/release-1.2-1/R/depmix-class.R
===================================================================
--- tags/release-1.2-1/R/depmix-class.R (rev 0)
+++ tags/release-1.2-1/R/depmix-class.R 2012-08-09 09:21:21 UTC (rev 562)
@@ -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) <- c("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) <- c("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
+#
+
+
+
Added: tags/release-1.2-1/R/depmix.R
===================================================================
--- tags/release-1.2-1/R/depmix.R (rev 0)
+++ tags/release-1.2-1/R/depmix.R 2012-08-09 09:21:21 UTC (rev 562)
@@ -0,0 +1,97 @@
+#
+# 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, ...) {
+
+ #if(!is.null(data) & any(is.na(data))) stop("'depmixS4' does not currently handle missing data.")
+
+ # 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
+ if(is.null(initdata)) initdata=data
+ 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(any(is.na(data))) stop("'depmixS4' does not currently handle missing data.")
+ 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")
+ }
+
+ if(nstates==1&transition!=~1) {stop("1-state model can not have transition covariate")}
+
+ # 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
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/depmix -r 562
More information about the depmix-commits
mailing list