[Depmix-commits] r199 - trunk/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 30 12:47:09 CEST 2008
Author: maarten
Date: 2008-06-30 12:47:09 +0200 (Mon, 30 Jun 2008)
New Revision: 199
Added:
trunk/R/depmixsim-class.R
Modified:
trunk/R/depmix-class.R
trunk/R/responseMVN.R
trunk/R/transInit.R
Log:
put simulate methods in corresponding class files
added separate depmixsim file
Modified: trunk/R/depmix-class.R
===================================================================
--- trunk/R/depmix-class.R 2008-06-30 10:21:14 UTC (rev 198)
+++ trunk/R/depmix-class.R 2008-06-30 10:47:09 UTC (rev 199)
@@ -118,7 +118,91 @@
}
)
+setMethod("is.stationary",signature(object="depmix"),
+ function(object) {
+ return(object at stationary)
+ }
+)
+setMethod("simulate",signature(object="depmix"),
+ function(object,nsim=1,seed=NULL,...) {
+ 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) {
+ 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 <- apply(object at prior@x,2,rep,nsim)
+ for(j in 1:ns) {
+ if(!is.stationary(object)) object at transition[[j]]@x <- as.matrix(apply(object at transition[[j]]@x,2,rep,nsim))
+ for(i in 1:nr) {
+ object at response[[j]][[i]]@y <- as.matrix(responses[[i]])
+ object at response[[j]][[i]]@x <- as.matrix(apply(object at response[[j]][[i]]@x,2,rep,nsim))
+ }
+ }
+ object at ntimes <- rep(object at ntimes,nsim)
+
+ # make appropriate array for transition densities
+ nt <- sum(object at ntimes)
+ if(is.stationary(object)) trDens <- array(0,c(1,ns,ns)) else trDens <- array(0,c(nt,ns,ns))
+
+ # make appropriate array for response densities
+ dns <- array(,c(nt,nr,ns))
+
+ # compute observation and transition densities
+ for(i in 1:ns) {
+ for(j in 1:nr) {
+ dns[,j,i] <- dens(object at response[[i]][[j]]) # remove this response as an argument from the call to setpars
+ }
+ trDens[,,i] <- dens(object at transition[[i]])
+ }
+
+ # compute initial state probabilties
+ object at init <- dens(object at prior)
+ object at trDens <- trDens
+ object at dens <- dns
+
+ return(object)
+ }
+)
+
+
+
#
# SUMMARY method: to do
#
Added: trunk/R/depmixsim-class.R
===================================================================
--- trunk/R/depmixsim-class.R (rev 0)
+++ trunk/R/depmixsim-class.R 2008-06-30 10:47:09 UTC (rev 199)
@@ -0,0 +1,8 @@
+# Class for simulated depmix model
+
+setClass("depmix.sim",
+ contains="depmix",
+ representation(
+ states="matrix"
+ )
+)
\ No newline at end of file
Modified: trunk/R/responseMVN.R
===================================================================
--- trunk/R/responseMVN.R 2008-06-30 10:21:14 UTC (rev 198)
+++ trunk/R/responseMVN.R 2008-06-30 10:47:09 UTC (rev 199)
@@ -16,7 +16,7 @@
}
)
-dm_dmvnorm <- function(y,mean,sigma,log=FALSE,logdet,invsigma) {
+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
@@ -41,6 +41,15 @@
if(missing(invSigma)) {
distval <- mahalanobis(x, center = mean, cov = sigma)
} else {
+ if (NCOL(x) != NCOL(invSigma)) {
+ stop("x and invSigma have non-conforming size")
+ }
+ if (NROW(invSigma) != NCOL(invSigma)) {
+ stop("invSigma must be a square matrix")
+ }
+ if (length(mean) != NROW(invSigma)) {
+ stop("mean and invSigma have non-conforming size")
+ }
distval <- mahalanobis(x, center = mean, cov = invsigma, inverted=TRUE)
}
if(missing(logdet)) logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values))
@@ -54,14 +63,14 @@
setMethod("logDens","MVNresponse",
- function(object) {
- dm_dmvnorm(x=object at y,mean=predict(object),sigma=object at parameters$Sigma,log=TRUE)
+ 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)
+ function(object,log=FALSE,...) {
+ dm_dmvnorm(x=object at y,mean=predict(object),sigma=object at parameters$Sigma,log=log,...)
}
)
Modified: trunk/R/transInit.R
===================================================================
--- trunk/R/transInit.R 2008-06-30 10:21:14 UTC (rev 198)
+++ trunk/R/transInit.R 2008-06-30 10:47:09 UTC (rev 199)
@@ -126,3 +126,27 @@
}
)
+setMethod("simulate",signature(object="transInit"),
+ function(object,nsim=1,seed=NULL,times,is.prior=FALSE,...) {
+ if(!is.null(seed)) set.seed(seed)
+ if(is.prior) {
+ pr <- dens(object)
+ sims <- array(apply(pr,1,rmultinom,n=nsim,size=1),dim=c(ncol(pr),nsim,nrow(pr)))
+ states <- t(apply(sims,c(2,3), function(x) which(x==1)))
+ return(states)
+ } else {
+ if(missing(times)) {
+ # this is likely to be a stationary model...
+ pr <- predict(object)
+ } else {
+ pr <- predict(object)[times,]
+ if(length(times)==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))
+ states <- t(apply(sims,c(2,3), function(x) which(x==1)))
+ # states <- apply(apply(pr,2,rmultinom rmultinom(nt*nsim,size=1,prob=pr),2,function(x) which(x==1))
+ return(states)
+ }
+ }
+)
\ No newline at end of file
More information about the depmix-commits
mailing list