[Depmix-commits] r213 - pkg/R trunk/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 2 12:00:57 CEST 2008
Author: maarten
Date: 2008-07-02 12:00:57 +0200 (Wed, 02 Jul 2008)
New Revision: 213
Modified:
pkg/R/depmix-class.R
trunk/R/depmix-class.R
Log:
simulate(depmix) now forces init at x to be a matrix
Modified: pkg/R/depmix-class.R
===================================================================
--- pkg/R/depmix-class.R 2008-07-01 22:07:39 UTC (rev 212)
+++ pkg/R/depmix-class.R 2008-07-02 10:00:57 UTC (rev 213)
@@ -125,81 +125,81 @@
)
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)
- }
+ 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) {
+ 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)
+ }
)
Modified: trunk/R/depmix-class.R
===================================================================
--- trunk/R/depmix-class.R 2008-07-01 22:07:39 UTC (rev 212)
+++ trunk/R/depmix-class.R 2008-07-02 10:00:57 UTC (rev 213)
@@ -126,6 +126,7 @@
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)
@@ -167,7 +168,7 @@
class(object) <- "depmix.sim"
object at states <- as.matrix(states)
- object at prior@x <- apply(object at prior@x,2,rep,nsim)
+ 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) {
More information about the depmix-commits
mailing list