[Depmix-commits] r135 - in trunk: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 29 12:40:53 CEST 2008
Author: maarten
Date: 2008-05-29 12:40:53 +0200 (Thu, 29 May 2008)
New Revision: 135
Added:
trunk/R/responseGLMGAMMA.r
trunk/R/responseGLMPOISSON.R
trunk/man/simulate.Rd
Modified:
trunk/R/allGenerics.R
trunk/R/responseGLMBINOM.R
trunk/R/responseGLMMULTINOM.R
trunk/R/responseMVN.R
trunk/R/responseNORM.R
trunk/R/simulate.R
Log:
- added Poisson and Gamma GLMresponse models
- reorganized simulate methods
Modified: trunk/R/allGenerics.R
===================================================================
--- trunk/R/allGenerics.R 2008-05-21 10:18:43 UTC (rev 134)
+++ trunk/R/allGenerics.R 2008-05-29 10:40:53 UTC (rev 135)
@@ -60,3 +60,5 @@
setGeneric("nresp", function(object, ...) standardGeneric("nresp"))
+setGeneric("is.stationary", function(object,...) standardGeneric("is.stationary"))
+
Modified: trunk/R/responseGLMBINOM.R
===================================================================
--- trunk/R/responseGLMBINOM.R 2008-05-21 10:18:43 UTC (rev 134)
+++ trunk/R/responseGLMBINOM.R 2008-05-29 10:40:53 UTC (rev 135)
@@ -19,3 +19,18 @@
dbinom(x=object at y,size=object at n,prob=predict(object),log=log)
}
)
+
+setMethod("simulate",signature(object="BINOMresponse"),
+ function(object,nsim=1,seed=NULL,time) {
+ if(missing(time)) {
+ # draw in one go
+ pr <- predict(object)
+ } else {
+ pr <- predict(object)[time,]
+ }
+ nt <- nrow(pr)
+ response <- rbinom(nt*nsim,size=object at n,prob=pr)
+ if(nsim > 1) response <- matrix(response,ncol=nsim)
+ return(response)
+ }
+)
Added: trunk/R/responseGLMGAMMA.r
===================================================================
--- trunk/R/responseGLMGAMMA.r (rev 0)
+++ trunk/R/responseGLMGAMMA.r 2008-05-29 10:40:53 UTC (rev 135)
@@ -0,0 +1,35 @@
+setClass("GAMMAresponse",contains="GLMresponse")
+
+# method 'fit'
+# use: in EM (M step)
+# returns: (fitted) response with (new) estimates of parameters
+
+# methods 'logDens' & dens
+# use: instead of density slot in rModel
+# returns: matrix with log(p(y|x,parameters))
+setMethod("logDens","GAMMAresponse",
+ function(object) {
+ dpois(x=object at y,shape=predict(object),log=TRUE)
+ }
+)
+
+setMethod("dens","GAMMAresponse",
+ function(object,log=FALSE) {
+ dpois(x=object at y,shape=predict(object),log=log)
+ }
+)
+
+setMethod("simulate",signature(object="GAMMAresponse"),
+ function(object,nsim=1,seed=NULL,time) {
+ if(missing(time)) {
+ # draw in one go
+ shape <- predict(object)
+ } else {
+ shape <- predict(object)[time,]
+ }
+ nt <- nrow(shape)
+ response <- rgamma(nt*nsim,shape=shape)
+ if(nsim > 1) response <- matrix(response,ncol=nsim)
+ return(response)
+ }
+)
Modified: trunk/R/responseGLMMULTINOM.R
===================================================================
--- trunk/R/responseGLMMULTINOM.R 2008-05-21 10:18:43 UTC (rev 134)
+++ trunk/R/responseGLMMULTINOM.R 2008-05-29 10:40:53 UTC (rev 135)
@@ -40,3 +40,19 @@
}
}
)
+
+setMethod("simulate",signature(object="MULTINOMresponse"),
+ function(object,nsim=1,seed=NULL,time) {
+ if(missing(time)) {
+ # draw all times in one go
+ pr <- predict(object)
+ } else {
+ pr <- predict(object)[time,]
+ if(length(time)==1) pr <- matrix(pr,ncol=length(pr))
+ }
+ nt <- nrow(pr)
+ sims <- array(apply(pr,1,rmultinom,n=nsim,size=1),dim=c(ncol(pr),nsim,nt))
+ response <- t(apply(sims,c(2,3), function(x) which(x==1)))
+ return(response)
+ }
+)
Added: trunk/R/responseGLMPOISSON.R
===================================================================
--- trunk/R/responseGLMPOISSON.R (rev 0)
+++ trunk/R/responseGLMPOISSON.R 2008-05-29 10:40:53 UTC (rev 135)
@@ -0,0 +1,35 @@
+setClass("POISSIONresponse",contains="GLMresponse")
+
+# method 'fit'
+# use: in EM (M step)
+# returns: (fitted) response with (new) estimates of parameters
+
+# methods 'logDens' & dens
+# use: instead of density slot in rModel
+# returns: matrix with log(p(y|x,parameters))
+setMethod("logDens","POISSONresponse",
+ function(object) {
+ dpois(x=object at y,lambda=predict(object),log=TRUE)
+ }
+)
+
+setMethod("dens","POISSONresponse",
+ function(object,log=FALSE) {
+ dpois(x=object at y,lambda=predict(object),log=log)
+ }
+)
+
+setMethod("simulate",signature(object="POISSONresponse"),
+ function(object,nsim=1,seed=NULL,time) {
+ if(missing(time)) {
+ # draw in one go
+ lambda <- predict(object)
+ } else {
+ lambda <- predict(object)[time,]
+ }
+ nt <- nrow(lambda)
+ response <- rpois(nt*nsim,lambda=lambda)
+ if(nsim > 1) response <- matrix(response,ncol=nsim)
+ return(response)
+ }
+)
Modified: trunk/R/responseMVN.R
===================================================================
--- trunk/R/responseMVN.R 2008-05-21 10:18:43 UTC (rev 134)
+++ trunk/R/responseMVN.R 2008-05-29 10:40:53 UTC (rev 135)
@@ -11,3 +11,68 @@
object
}
)
+
+dm_dmvnorm <- function(y,mean,sigma,log=FALSE,logdet,invsigma) {
+ # taken from mvtnorm package
+ # allows passing of logdet (sigma) and invsigma to save
+ # computation when called repeated times with same sigma
+ if (is.vector(x)) {
+ x <- matrix(x, ncol = length(x))
+ }
+ if (missing(mean)) {
+ mean <- rep(0, length = ncol(x))
+ }
+ if (missing(sigma)) {
+ sigma <- diag(ncol(x))
+ }
+ if (NCOL(x) != NCOL(sigma)) {
+ stop("x and sigma have non-conforming size")
+ }
+ if (NROW(sigma) != NCOL(sigma)) {
+ stop("sigma must be a square matrix")
+ }
+ if (length(mean) != NROW(sigma)) {
+ stop("mean and sigma have non-conforming size")
+ }
+ if(missing(invSigma)) {
+ distval <- mahalanobis(x, center = mean, cov = sigma)
+ } else {
+ distval <- mahalanobis(x, center = mean, cov = invsigma, inverted=TRUE)
+ }
+ if(missing(logdet)) logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values))
+ logretval <- -(ncol(x) * log(2 * pi) + logdet + distval)/2
+ if (log) {
+ return(logretval)
+ } else {
+ return(exp(logretval))
+ }
+}
+
+
+setMethod("logDens","MVNresponse",
+ function(object) {
+ dm_dmvnorm(x=object at y,mean=predict(object),sigma=object at parameters$Sigma,log=TRUE)
+ }
+)
+
+setMethod("dens","MVNresponse",
+ function(object,log=FALSE) {
+ dm_dmvnorm(x=object at y,mean=predict(object),sigma=object at parameters$Sigma,log=log)
+ }
+)
+
+setMethod("simulate",signature(object="MVNresponse"),
+ function(object,nsim=1,seed=NULL,time) {
+ if(missing(time)) {
+ # draw in one go
+ mu <- predict(object)
+ } else {
+ mu <- predict(object)[time,]
+ }
+ nt <- nrow(mu)
+ response <- mvrnorm(nt*nsim,mu=mu,Sigma=object at parameters$Sigma)
+ if(nsim > 1) response <- array(response,dim=c(nt,ncol(response),nsim))
+ return(response)
+ }
+)
+
Modified: trunk/R/responseNORM.R
===================================================================
--- trunk/R/responseNORM.R 2008-05-21 10:18:43 UTC (rev 134)
+++ trunk/R/responseNORM.R 2008-05-29 10:40:53 UTC (rev 135)
@@ -32,4 +32,20 @@
function(object) {
object at x%*%object at parameters$coefficients
}
-)
\ No newline at end of file
+)
+
+setMethod("simulate",signature(object="NORMresponse"),
+ function(object,nsim=1,seed=NULL,time) {
+ if(missing(time)) {
+ # draw in one go
+ mu <- predict(object)
+ } else {
+ mu <- predict(object)[time]
+ }
+ nt <- length(mu)
+ sd <- getpars(object)["sd"]
+ response <- rnorm(nt*nsim,mean=mu,sd=sd)
+ if(nsim > 1) response <- matrix(response,ncol=nsim)
+ return(response)
+ }
+)
Modified: trunk/R/simulate.R
===================================================================
--- trunk/R/simulate.R 2008-05-21 10:18:43 UTC (rev 134)
+++ trunk/R/simulate.R 2008-05-29 10:40:53 UTC (rev 135)
@@ -1,8 +1,8 @@
# simulate data from a depmix model
# TODO: move this to all generics etc...
-setGeneric("is.stationary", function(object,...) standardGeneric("is.stationary"))
- setMethod("is.stationary",signature(object="depmix"),
+#setGeneric("is.stationary", function(object,...) standardGeneric("is.stationary"))
+setMethod("is.stationary",signature(object="depmix"),
function(object) {
return(object at stationary)
}
@@ -83,43 +83,5 @@
}
)
-setMethod("simulate",signature(object="NORMresponse"),
- function(object,nsim=1,seed=NULL,time) {
- if(missing(time)) {
- # draw in one go
- mu <- predict(object)
- } else {
- mu <- predict(object)[time]
- }
- nt <- length(mu)
- sd <- getpars(object)["sd"]
- response <- rnorm(nt*nsim,mean=mu,sd=sd)
- if(nsim > 1) response <- matrix(response,ncol=nsim)
- return(response)
- }
-)
-setMethod("simulate",signature(object="MULTINOMresponse"),
- function(object,nsim=1,seed=NULL,time) {
- if(missing(time)) {
- # draw all times in one go
- pr <- predict(object)
- } else {
- pr <- predict(object)[time,]
- if(length(time)==1) pr <- matrix(pr,ncol=length(pr))
- }
- nt <- nrow(pr)
- sims <- array(apply(pr,1,rmultinom,n=nsim,size=1),dim=c(ncol(pr),nsim,nt))
- response <- t(apply(sims,c(2,3), function(x) which(x==1)))
-
-# if(nsim > 1) {
-# response <- t(apply(pr,1,function(x) apply(rmultinom(n=nsim,size=1,pr=x),2,function(y) which(y==1))))
-# } else {
-# response <- apply(pr,1,function(x) apply(rmultinom(n=nsim,size=1,pr=x),2,function(y) which(y==1)))
-# }
- return(response)
- }
-)
-
-
Added: trunk/man/simulate.Rd
===================================================================
--- trunk/man/simulate.Rd (rev 0)
+++ trunk/man/simulate.Rd 2008-05-29 10:40:53 UTC (rev 135)
@@ -0,0 +1,56 @@
+\name{simulate,depmix-method}
+
+\docType{method}
+
+\alias{simulate,GLMresponse-method}
+\alias{simulate,transInit-method}
+
+\title{Methods to simulate from depmix models}
+
+\description{
+
+Random draws from \code{depmix} objects.
+
+}
+
+\usage{
+ \S4method{simulate}{depmix}(object)
+ \S4method{simulate}{response}(object)
+ \S4method{simulate}{GLMresponse}(object)
+ \S4method{simulate}{transInit}(object)
+}
+
+\arguments{
+ \item{object}{Object to generate random draws. An object of class \code{depmix}, \code{response}, \code{transInit}}
+ \item{nsim}{The number of draws.}
+ \item{seed}{Set the seed.}
+ \item{times}{An indicator vector indicating for which times in the complete
+ series to generate the data. For internal use.}
+ \item{...}{Not used currently.}
+}
+
+\details{
+
+ For a \code{depmix} model, simulate generates \code{nsim} random state
+ sequences, each of the same length as the observation sequence in the
+ \code{depmix} model (i.e., \code{sum(ntimes(object))}. The state sequences
+ are then used to generate \code{nsim} observation sequences of thee same
+ length.
+
+ Setting the \code{times} option selects the time points in the total
+ state/observation sequence (i.e., counting is continued over ntimes).
+ Usage of the \code{times} option is not recommended.
+
+}
+
+\value{
+
+ For a \code{depmix} object, a \code{list} with two elements;
+ \code{states}, an N*nsim matrix with a state sequence in each column,
+ and \code{responses}, an N*nobs*nsim array with observation sequences, where
+ N = \code{sum(ntimes(object))}
+}
+
+\author{Maarten Speekenbrink}
+
+\keyword{methods}
More information about the depmix-commits
mailing list