[Pomp-commits] r356 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Oct 4 17:43:05 CEST 2010
Author: kingaa
Date: 2010-10-04 17:43:05 +0200 (Mon, 04 Oct 2010)
New Revision: 356
Modified:
pkg/R/simulate-pomp.R
Log:
- 'simulate.internal' is now the named function for the guts of 'simulate'
Modified: pkg/R/simulate-pomp.R
===================================================================
--- pkg/R/simulate-pomp.R 2010-10-01 14:05:21 UTC (rev 355)
+++ pkg/R/simulate-pomp.R 2010-10-04 15:43:05 UTC (rev 356)
@@ -1,80 +1,78 @@
## simulate a partially-observed Markov process
-setMethod(
- "simulate",
- "pomp",
- function (object, nsim = 1, seed = NULL, params,
- states = FALSE, obs = FALSE,
- times = time(object,t0=TRUE), ...) {
- ntimes <- length(times)
- times <- as.numeric(times)
- if (ntimes<1)
- stop("if length of ",sQuote("times")," is less than 1, there is no work to do",call.=FALSE)
- if (missing(params)) {
- if (length(object at params)>0)
- params <- object at params
- else
- stop("no ",sQuote("params")," specified",call.=FALSE)
- }
- if (is.null(dim(params)))
- params <- matrix(params,ncol=1,dimnames=list(names(params),NULL))
- npars <- ncol(params)
- if (!is.null(seed)) { # set the random seed (be very careful about this)
- if (!exists('.Random.seed',envir=.GlobalEnv)) runif(1)
- save.seed <- get('.Random.seed',envir=.GlobalEnv)
- set.seed(seed)
- }
- if (!is.numeric(nsim)||(length(nsim)>1)||nsim<1)
- stop(sQuote("nsim")," must be a positive integer")
- nsim <- as.integer(nsim)
- xstart <- init.state(object,params=params,t0=times[1])
- nreps <- npars*nsim # total number of replicates
- ## we will do the process model simulations with single calls to the user functions
- if (nsim > 1) { # make nsim copies of the IC and parameter matrices
- xstart <- array(rep(xstart,nsim),dim=c(nrow(xstart),nreps),dimnames=list(rownames(xstart),NULL))
- params <- array(rep(params,nsim),dim=c(nrow(params),nreps),dimnames=list(rownames(params),NULL))
- }
- x <- rprocess(object,xstart=xstart,times=times,params=params) # simulate the process model
- ## rprocess returns a rank-3 matrix (nvars x nreps x ntimes)
- if (obs || !states) { # we need only simulate the observation process if obs=T or states=F
- y <- rmeasure(object,x=x,times=times,params=params)
- ## rmeasure returns a rank-3 matrix (nobs x nreps x ntimes)
- }
- if (!is.null(seed)) { # restore the RNG state
- assign('.Random.seed',save.seed,envir=.GlobalEnv)
- }
- if (!obs && !states) { # both obs=F and states=F, return a list of pomp objects
- nm.y <- rownames(y)
- nobs <- nrow(y)
- retval <- lapply(
- seq(length=nreps),
- function (rep) {
- po <- as(object,'pomp') # copy the pomp object
- po at t0 <- times[1]
- po at times <- times[-1]
- po at data <- array(0,dim=c(nobs,ntimes-1),dimnames=list(nm.y,NULL))
- po at data[,] <- y[,rep,-1] # replace the data
- po at params <- params[,rep] # store the parameters
- po at states <- array(dim=dim(x)[c(1,3)],dimnames=list(rownames(x),NULL))
- po at states[,] <- x[,rep,] # store the states
- po
- }
- )
- if (nreps==1) {
- retval <- retval[[1]]
- }
- } else {
- if (!obs && states) { # return states only
- retval <- x
- } else if (obs && !states) { # return only observations
- retval <- y
- } else { # return both states and observations
- retval <- list(
- states = x,
- obs = y
- )
- }
- }
- retval
- }
- )
+simulate.internal <- function (object, nsim = 1, seed = NULL, params,
+ states = FALSE, obs = FALSE,
+ times = time(object,t0=TRUE), ...) {
+ ntimes <- length(times)
+ times <- as.numeric(times)
+ if (ntimes<1)
+ stop("if length of ",sQuote("times")," is less than 1, there is no work to do",call.=FALSE)
+ if (missing(params)) {
+ if (length(object at params)>0)
+ params <- object at params
+ else
+ stop("no ",sQuote("params")," specified",call.=FALSE)
+ }
+ if (is.null(dim(params)))
+ params <- matrix(params,ncol=1,dimnames=list(names(params),NULL))
+ npars <- ncol(params)
+ if (!is.null(seed)) { # set the random seed (be very careful about this)
+ if (!exists('.Random.seed',envir=.GlobalEnv)) runif(1)
+ save.seed <- get('.Random.seed',envir=.GlobalEnv)
+ set.seed(seed)
+ }
+ if (!is.numeric(nsim)||(length(nsim)>1)||nsim<1)
+ stop(sQuote("nsim")," must be a positive integer")
+ nsim <- as.integer(nsim)
+ xstart <- init.state(object,params=params,t0=times[1])
+ nreps <- npars*nsim # total number of replicates
+ ## we will do the process model simulations with single calls to the user functions
+ if (nsim > 1) { # make nsim copies of the IC and parameter matrices
+ xstart <- array(rep(xstart,nsim),dim=c(nrow(xstart),nreps),dimnames=list(rownames(xstart),NULL))
+ params <- array(rep(params,nsim),dim=c(nrow(params),nreps),dimnames=list(rownames(params),NULL))
+ }
+ x <- rprocess(object,xstart=xstart,times=times,params=params) # simulate the process model
+ ## rprocess returns a rank-3 matrix (nvars x nreps x ntimes)
+ if (obs || !states) { # we need only simulate the observation process if obs=T or states=F
+ y <- rmeasure(object,x=x,times=times,params=params)
+ ## rmeasure returns a rank-3 matrix (nobs x nreps x ntimes)
+ }
+ if (!is.null(seed)) { # restore the RNG state
+ assign('.Random.seed',save.seed,envir=.GlobalEnv)
+ }
+ if (!obs && !states) { # both obs=F and states=F, return a list of pomp objects
+ nm.y <- rownames(y)
+ nobs <- nrow(y)
+ retval <- lapply(
+ seq(length=nreps),
+ function (rep) {
+ po <- as(object,'pomp') # copy the pomp object
+ po at t0 <- times[1]
+ po at times <- times[-1]
+ po at data <- array(0,dim=c(nobs,ntimes-1),dimnames=list(nm.y,NULL))
+ po at data[,] <- y[,rep,-1] # replace the data
+ po at params <- params[,rep] # store the parameters
+ po at states <- array(dim=dim(x)[c(1,3)],dimnames=list(rownames(x),NULL))
+ po at states[,] <- x[,rep,] # store the states
+ po
+ }
+ )
+ if (nreps==1) {
+ retval <- retval[[1]]
+ }
+ } else {
+ if (!obs && states) { # return states only
+ retval <- x
+ } else if (obs && !states) { # return only observations
+ retval <- y
+ } else { # return both states and observations
+ retval <- list(
+ states = x,
+ obs = y
+ )
+ }
+ }
+ retval
+}
+
+setMethod("simulate",signature=signature(object="pomp"),definition=simulate.internal)
More information about the pomp-commits
mailing list