[Depmix-commits] r489 - pkg/depmixS4/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Sep 2 15:56:24 CEST 2011
Author: ingmarvisser
Date: 2011-09-02 15:56:24 +0200 (Fri, 02 Sep 2011)
New Revision: 489
Modified:
pkg/depmixS4/R/allGenerics.R
pkg/depmixS4/R/depmix-class.R
pkg/depmixS4/R/depmixfit-class.R
pkg/depmixS4/R/depmixsim-class.R
pkg/depmixS4/R/fb.R
pkg/depmixS4/R/forwardbackward.R
pkg/depmixS4/R/freepars.R
pkg/depmixS4/R/getpars.R
pkg/depmixS4/R/llratio.R
pkg/depmixS4/R/logLik.R
pkg/depmixS4/R/lystig.R
pkg/depmixS4/R/nobs.R
pkg/depmixS4/R/response-class.R
pkg/depmixS4/R/responseGLM.R
pkg/depmixS4/R/responseGLMBINOM.R
pkg/depmixS4/R/responseGLMGAMMA.R
pkg/depmixS4/R/responseGLMMULTINOM.R
pkg/depmixS4/R/responseGLMPOISSON.R
pkg/depmixS4/R/responseNORM.R
pkg/depmixS4/R/setpars.R
Log:
Changed line endings to Unix consistently.
Modified: pkg/depmixS4/R/allGenerics.R
===================================================================
--- pkg/depmixS4/R/allGenerics.R 2011-09-01 15:16:58 UTC (rev 488)
+++ pkg/depmixS4/R/allGenerics.R 2011-09-02 13:56:24 UTC (rev 489)
@@ -1,85 +1,85 @@
-
-#
-# 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) {
- require(stats)
- require(methods)
- require(MASS)
- require(nnet)
- require(Rsolnp)
- require(stats4)
- 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"))
+
+#
+# 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) {
+ require(stats)
+ require(methods)
+ require(MASS)
+ require(nnet)
+ require(Rsolnp)
+ require(stats4)
+ 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"))
Modified: pkg/depmixS4/R/depmix-class.R
===================================================================
--- pkg/depmixS4/R/depmix-class.R 2011-09-01 15:16:58 UTC (rev 488)
+++ pkg/depmixS4/R/depmix-class.R 2011-09-02 13:56:24 UTC (rev 489)
@@ -1,306 +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
-#
-
-
-
+
+#
+# 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
+#
+
+
+
Modified: pkg/depmixS4/R/depmixfit-class.R
===================================================================
--- pkg/depmixS4/R/depmixfit-class.R 2011-09-01 15:16:58 UTC (rev 488)
+++ pkg/depmixS4/R/depmixfit-class.R 2011-09-02 13:56:24 UTC (rev 489)
@@ -1,138 +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")
- }
- }
- }
-)
-
-
+
+#
+# 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")
+ }
+ }
+ }
+)
+
+
Modified: pkg/depmixS4/R/depmixsim-class.R
===================================================================
--- pkg/depmixS4/R/depmixsim-class.R 2011-09-01 15:16:58 UTC (rev 488)
+++ pkg/depmixS4/R/depmixsim-class.R 2011-09-02 13:56:24 UTC (rev 489)
@@ -1,16 +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"
- )
-)
-
+# Classes for simulated mix and depmix models
+
+setClass("mix.sim",
+ contains="mix",
+ representation(
+ states="matrix"
+ )
+)
+
+setClass("depmix.sim",
+ contains="depmix",
+ representation(
+ states="matrix"
+ )
+)
+
Modified: pkg/depmixS4/R/fb.R
===================================================================
--- pkg/depmixS4/R/fb.R 2011-09-01 15:16:58 UTC (rev 488)
+++ pkg/depmixS4/R/fb.R 2011-09-02 13:56:24 UTC (rev 489)
@@ -1,124 +1,124 @@
-#
-# Maarten Speekenbrink
-#
-# FORWARD-BACKWARD algoritme, 23-3-2008
-#
-
-fb <- function(init,A,B,ntimes=NULL,return.all=FALSE,stationary=TRUE,useC=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 sca
-
- # NOTE: xi[t,i,j] = P(S[t] = j & S[t+1] = i) !!!NOTE the order of i and j!!!
-
- nt <- nrow(B)
- ns <- ncol(init)
-
-# print(head(B))
-#
-# bin <- array(0,dim=dim(B)[c(1,3)])
-# res <- .C("ddens",
-# as.double(B),
-# out=as.double(bin),
-# dim=dim(B),
-# package="depmixS4")[c("out")]
-# B <- matrix(res$out,nc=ns,byrow=TRUE)
-#
-# print(head(B))
-
- B <- apply(B,c(1,3),prod)
-
- if(is.null(ntimes)) ntimes <- nt
-
- lt <- length(ntimes)
- et <- cumsum(ntimes)
- bt <- c(1,et[-lt]+1)
-
- if(useC) {
-
- alpha <- matrix(0,ncol=ns,nrow=nt)
- beta <- matrix(0,ncol=ns,nrow=nt)
- sca <- rep(0,nt)
- xi <- array(0,dim=c(nt,ns,ns))
-
- res <- .C("forwardbackward",
- hom=as.integer(stationary),
- ns=as.integer(ns),
- lt=as.integer(lt),
- nt=as.integer(nt),
- ntimes=as.integer(ntimes),
- bt=as.integer(bt),
- et=as.integer(et),
- init=as.double(t(init)),
- A=as.double(A),
- B=as.double(t(B)),
- alpha=as.double(alpha),
- beta=as.double(beta),
- sca=as.double(sca),
- xi=as.double(xi),
- PACKAGE="depmixS4")[c("alpha","beta","sca","xi")]
-
- alpha <- matrix(res$alpha,nc=ns,byrow=TRUE)
- beta <- matrix(res$beta,nc=ns,byrow=TRUE)
- xi <- array(res$xi,dim=c(nt,ns,ns))
- xi[et,,] <- NA
- sca <- res$sca
-
- } else {
-
- alpha <- matrix(ncol=ns,nrow=nt)
- beta <- matrix(ncol=ns,nrow=nt)
- sca <- vector(length=nt)
- xi <- array(dim=c(nt,ns,ns))
-
- 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
-}
-
+#
+# Maarten Speekenbrink
+#
+# FORWARD-BACKWARD algoritme, 23-3-2008
+#
+
+fb <- function(init,A,B,ntimes=NULL,return.all=FALSE,stationary=TRUE,useC=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 sca
+
+ # NOTE: xi[t,i,j] = P(S[t] = j & S[t+1] = i) !!!NOTE the order of i and j!!!
+
+ nt <- nrow(B)
+ ns <- ncol(init)
+
+# print(head(B))
+#
+# bin <- array(0,dim=dim(B)[c(1,3)])
+# res <- .C("ddens",
+# as.double(B),
+# out=as.double(bin),
+# dim=dim(B),
+# package="depmixS4")[c("out")]
+# B <- matrix(res$out,nc=ns,byrow=TRUE)
+#
+# print(head(B))
+
+ B <- apply(B,c(1,3),prod)
+
+ if(is.null(ntimes)) ntimes <- nt
+
+ lt <- length(ntimes)
+ et <- cumsum(ntimes)
+ bt <- c(1,et[-lt]+1)
+
+ if(useC) {
+
+ alpha <- matrix(0,ncol=ns,nrow=nt)
+ beta <- matrix(0,ncol=ns,nrow=nt)
+ sca <- rep(0,nt)
+ xi <- array(0,dim=c(nt,ns,ns))
+
+ res <- .C("forwardbackward",
+ hom=as.integer(stationary),
+ ns=as.integer(ns),
+ lt=as.integer(lt),
+ nt=as.integer(nt),
+ ntimes=as.integer(ntimes),
+ bt=as.integer(bt),
+ et=as.integer(et),
+ init=as.double(t(init)),
+ A=as.double(A),
+ B=as.double(t(B)),
+ alpha=as.double(alpha),
+ beta=as.double(beta),
+ sca=as.double(sca),
+ xi=as.double(xi),
+ PACKAGE="depmixS4")[c("alpha","beta","sca","xi")]
+
+ alpha <- matrix(res$alpha,nc=ns,byrow=TRUE)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/depmix -r 489
More information about the depmix-commits
mailing list