pomp {pomp} | R Documentation |
Create a new pomp
object.
pomp(data, times, t0, rprocess, dprocess, rmeasure, dmeasure, initializer, ...)
data |
An array holding the data.
This array is of dimensions nobs x ntimes , where nobs is the number of observed variables and ntimes is the number of times at which observations were made.
It is also possible to specify data as a data-frame, in which case the argument times must be the name of the column of observation times.
Note that if data is a data-frame, it will be coerced to an array with storage-mode `double'.
|
times |
The times corresponding to the observations.
times must be a strictly increasing numeric vector.
If data is a data-frame, times should be the name of the column of observation times.
|
t0 |
The zero-time.
This must be no later than the time of the first observation, times[1] .
|
rprocess |
Function of prototype rprocess(xstart,times,params,...) which simulates from the unobserved process.
|
dprocess |
Function of prototype dprocess(x,times,params,log=FALSE,...) which evaluates the likelihood of a sequence of consecutive state transitions.
|
rmeasure |
Function of prototype rmeasure(x,times,params,...) which simulates from the observation process.
|
dmeasure |
Function of prototype dmeasure(y,x,times,params,log=FALSE,...) which gives the likelihood of y given x .
|
initializer |
Function of prototype initializer(params,t0,...) which yields initial conditions for the state process when given a vector, params , of parameters.
By default, i.e., if it is unspecified when pomp is called, the initializer assumes any names in params that end in ".0" correspond to initial value parameters.
These are simply copied over as initial conditions when init.state is called (see init.state-pomp ).
The names of the state variables are the same as the corresponding initial value parameters, but with the ".0" dropped.
|
... |
Any additional arguments are stored in a slot `userdata' and are passed as arguments to each of the functions rprocess , dprocess , rmeasure , dmeasure , and initializer whenever they are evaluated.
Using this mechanism, the user can store additional information necessary for the definition of the model.
|
It is the user's responsibility to ensure that the rprocess
, dprocess
, rmeasure
, dmeasure
, and initializer
functions satisfy the following conditions:
rprocess
rprocess
must have at least the following arguments:
xstart
, times
, params
, and ...
.
It can also take additional arguments.
It is guaranteed that these will be filled with the corresponding elements the user has included as additional arguments in the construction of the pomp
object.
In calls to rprocess
, xstart
will be a rank-2 array (matrix) with rows corresponding to state variables and columns corresponding to independent realizations of the process.
params
will similarly be a rank-2 array with rows corresponding to parameters and columns corresponding to independent realizations.
The columns of params
are to be matched up with those of xstart
;
in particular, they will agree in number.
Both xstart
and params
must have rownames, which are available for use by the user.
rprocess
must return a rank-3 array with rownames.
Suppose x
is the array returned.
Then dim(x)=c(nvars,nreps,ntimes)
, where nvars
is the number of state variables (=nrow(xstart)
), nreps
is the number of independent realizations simulated (=ncol(xstart)
), and ntimes
is the length of the vector times
.
x[,j,k]
is the value of the state process in the j
-th realization at time times[k]
.
In particular, x[,,1]
must be identical to xstart
.
The rownames of x
must correspond to those of xstart
.
dprocess
dprocess
must have at least the following arguments:
x
, times
, params
, log
, and ...
.
It may take additional arguments.
It is guaranteed that these will be filled with the corresponding elements the user has included as additional arguments in the construction of the pomp
object.
In calls to dprocess
, x
will be an nvars
x nreps
x ntimes
array, where these terms have the same meanings as above.
params
will be a rank-2 array with rows corresponding to individual parameters and columns corresponding to independent realizations.
The columns of params
are to be matched up with those of x
; in particular, they will agree in number.
Both x
and params
must have rownames, available for use by the user.
dprocess
must return a rank-2 array (matrix).
Suppose d
is the array returned.
Then dim(d)=c(nreps,ntimes-1)
.
d[j,k]
is the probability density of the transition from state x[,j,k-1]
at time times[k-1]
to state x[,j,k]
at time times[k]
.
If log=TRUE
, then the log of the p.d.f. is returned.
In writing this function, you may assume that the transitions are consecutive.
rmeasure
rmeasure
must have at least the arguments x
, times
, params
, and ...
.
It may take additional arguments, which will be filled with user-specified data as above.
x
must be a rank-3 array dimension c(nvars,nreps,ntimes)
, where these variables have the same meanings as above.
times
is the corresponding set of times.
params
is a rank-2 array of dimension c(npars,nreps)
as above.
rmeasure
must return a rank-3 array with rownames.
If y
is the returned array, then dim(y)=c(nobs,nreps,ntimes)
, where nobs
is the number of observable variables and nreps
, ntimes
agree with the corresponding dimensions of x
.
y[,j,k]
must be the vector of observables in the j
-th realization at time times[k]
.
dmeasure
dmeasure
must have at least the arguments y
, x
, times
, params
, log
, and ...
.
y
is a rank-2 array of observations (nobs
x ntimes
);
x
, a rank-3 array of states;
params
, a rank-2 array containing parameters;
and times
, the corresponding observation times.
It may take additional arguments which will be filled with user-specified data as above.
dmeasure
must return a rank-2 array of dimension nreps
x ntimes
.
If d
is the returned array, then d[j,k]
is the p.d.f. of y[,k]
given x[,j,k]
at time times[k]
.
If log=TRUE
, then the log of the p.d.f. is returned.
initializer
initializer
must have at least the arguments params
, t0
, and ...
.
params
is a named vector of parameters.
t0
will be the time at which initial conditions are desired.
initializer
must return a named vector of initial conditions.
An object of class pomp
.
Some error checking is done, but complete error checking is impossible. If the user-specified functions do not conform to the above specifications (see Details), then the results may be invalid.
Aaron A. King (kingaa at umich dot edu)
pomp-class, euler, vignette(intro_to_pomp)
## Simple example: a 2-D Brownian motion, observed with normal error ## first, set up the pomp object: rw2 <- pomp( rprocess = function (xstart, times, params, ...) { ## this function simulates two independent random walks with intensities s1, s2 nsims <- ncol(params) ntimes <- length(times) dt <- diff(times) x <- array(0,dim=c(2,nsims,ntimes)) rownames(x) <- rownames(xstart) noise.sds <- params[c('s1','s2'),] x[,,1] <- xstart for (j in 2:ntimes) { ## we are mimicking a continuous-time process, so the increments have SD ~ sqrt(dt) ## note that we do not have to assume that 'times' are equally spaced x[,,j] <- rnorm(n=2*nsims,mean=x[,,j-1],sd=noise.sds*dt[j-1]) } x }, dprocess = function (x, times, params, log = FALSE, ...) { ## given a sequence of consecutive states in 'x', this function computes the p.d.f. nsims <- ncol(params) ntimes <- length(times) dt <- diff(times) d <- array(0,dim=c(2,nsims,ntimes-1)) noise.sds <- params[c('s1','s2'),] for (j in 2:ntimes) d[,,j-1] <- dnorm(x[,,j]-x[,,j-1],mean=0,sd=noise.sds*dt[j-1],log=TRUE) if (log) { apply(d,c(2,3),sum) } else { exp(apply(d,c(2,3),sum)) } }, rmeasure = function (x, times, params, ...) { ## noisy observations of the two walks with common noise SD 'tau' nsims <- dim(x)[2] ntimes <- dim(x)[3] y <- array(0,dim=c(2,nsims,ntimes)) rownames(y) <- c('y1','y2') for (j in 1:nsims) { for (k in 1:ntimes) { y[,j,k] <- rnorm(2,mean=x[,j,k],sd=params['tau',j]) } } y }, dmeasure = function (y, x, times, params, log = FALSE, ...) { ## noisy observations of the two walks with common noise SD 'tau' d1 <- dnorm( x=y['y1',], mean=x['x1',,], sd=params['tau',], log=TRUE ) d2 <- dnorm( x=y['y2',], mean=x['x2',,], sd=params['tau',], log=TRUE ) if (log) { d1+d2 } else { exp(d1+d2) } }, times=1:100, data=rbind( y1=rep(0,100), y2=rep(0,100) ), t0=0 ) ## specify some parameters p <- rbind(s1=c(2,2,3),s2=c(0.1,1,2),tau=c(1,5,0),x1.0=c(0,0,5),x2.0=c(0,0,0)) ## simulate examples <- simulate(rw2,params=p) rw2 <- examples[[1]] ## by default, simulate produces a list of pomp objects plot(rw2) t <- seq(0,20) X <- simulate(rw2,params=p[,1],nsim=10,states=TRUE,obs=TRUE,times=t) ## compute the process model likelihoods f <- dprocess( rw2, x=X$states, times=t, params=matrix( p[,1], nrow=nrow(p), ncol=10, dimnames=list(rownames(p),NULL) ), log=TRUE ) apply(f,1,sum) ## compute the measurement likelihoods g <- dmeasure( rw2, y=X$obs[,7,], x=X$states, times=t, params=matrix( p[,1], nrow=nrow(p), ncol=10, dimnames=list(rownames(p),NULL) ), log=TRUE ) apply(g,1,sum) ## For more examples, see the vignettes.