[Pomp-commits] r793 - in pkg/pomp: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jan 8 04:21:33 CET 2013
Author: kingaa
Date: 2013-01-08 04:21:32 +0100 (Tue, 08 Jan 2013)
New Revision: 793
Added:
pkg/pomp/R/pomp-class.R
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/R/pfilter.R
pkg/pomp/R/pomp.R
Log:
- some changes preparatory to merging in mif2
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2013-01-07 20:15:21 UTC (rev 792)
+++ pkg/pomp/DESCRIPTION 2013-01-08 03:21:32 UTC (rev 793)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.43-5
-Date: 2012-08-22
+Version: 0.43-6
+Date: 2013-01-07
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing, Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
Maintainer: Aaron A. King <kingaa at umich.edu>
URL: http://pomp.r-forge.r-project.org
@@ -14,8 +14,9 @@
BuildVignettes: no
Collate: aaa.R authors.R version.R eulermultinom.R plugins.R
parmat.R slice-design.R profile-design.R sobol.R bsplines.R sannbox.R
- pomp-fun.R pomp.R pomp-methods.R rmeasure-pomp.R rprocess-pomp.R init-state-pomp.R
- dmeasure-pomp.R dprocess-pomp.R skeleton-pomp.R simulate-pomp.R trajectory-pomp.R plot-pomp.R
+ pomp-fun.R pomp-class.R pomp.R pomp-methods.R
+ rmeasure-pomp.R rprocess-pomp.R init-state-pomp.R dmeasure-pomp.R dprocess-pomp.R skeleton-pomp.R
+ simulate-pomp.R trajectory-pomp.R plot-pomp.R
pfilter.R pfilter-methods.R traj-match.R bsmc.R
mif-class.R particles-mif.R mif.R mif-methods.R compare-mif.R
pmcmc.R pmcmc-methods.R compare-pmcmc.R
Modified: pkg/pomp/R/pfilter.R
===================================================================
--- pkg/pomp/R/pfilter.R 2013-01-07 20:15:21 UTC (rev 792)
+++ pkg/pomp/R/pfilter.R 2013-01-08 03:21:32 UTC (rev 793)
@@ -7,6 +7,7 @@
pred.mean="array",
pred.var="array",
filter.mean="array",
+ paramMatrix="array",
eff.sample.size="numeric",
cond.loglik="numeric",
saved.states="list",
@@ -16,18 +17,38 @@
tol="numeric",
nfail="integer",
loglik="numeric"
+ ),
+ prototype=prototype(
+ pred.mean=array(data=numeric(0),dim=c(0,0)),
+ pred.var=array(data=numeric(0),dim=c(0,0)),
+ filter.mean=array(data=numeric(0),dim=c(0,0)),
+ paramMatrix=array(data=numeric(0),dim=c(0,0)),
+ eff.sample.size=numeric(0),
+ cond.loglik=numeric(0),
+ saved.states=list(),
+ saved.params=list(),
+ seed=as.integer(NA),
+ Np=as.integer(NA),
+ tol=as.double(NA),
+ nfail=as.integer(NA),
+ loglik=as.double(NA)
)
)
-## question: when pfilter.internal is called by mif, do we need to compute the prediction means and variances of the state variables each time, or only at the end?
+## question: when pfilter.internal is called by mif,
+## do we need to compute the prediction means and variances of the state variables each time,
+## or only at the end?
pfilter.internal <- function (object, params, Np,
tol, max.fail,
pred.mean, pred.var, filter.mean,
+ cooling.fraction, cooling.m, mif2 = FALSE,
.rw.sd, seed, verbose,
save.states, save.params,
- transform) {
-
+ transform,
+ .ndone = 0) {
+
+ mif2 <- as.logical(mif2)
transform <- as.logical(transform)
if (missing(seed)) seed <- NULL
@@ -38,7 +59,7 @@
save.seed <- get(".Random.seed",pos=.GlobalEnv)
set.seed(seed)
}
-
+
if (length(params)==0)
stop(sQuote("pfilter")," error: ",sQuote("params")," must be specified",call.=FALSE)
@@ -48,7 +69,7 @@
one.par <- FALSE
times <- time(object,t0=TRUE)
ntimes <- length(times)-1
-
+
if (missing(Np))
Np <- NCOL(params)
if (is.function(Np)) {
@@ -85,7 +106,7 @@
paramnames <- rownames(params)
if (is.null(paramnames))
stop(sQuote("pfilter")," error: ",sQuote("params")," must have rownames",call.=FALSE)
-
+
x <- init.state(
object,
params=if (transform) {
@@ -106,7 +127,7 @@
pparticles <- vector(mode="list",length=ntimes)
else
pparticles <- list()
-
+
random.walk <- !missing(.rw.sd)
if (random.walk) {
rw.names <- names(.rw.sd)
@@ -128,7 +149,7 @@
eff.sample.size <- numeric(ntimes)
nfail <- 0
npars <- length(rw.names)
-
+
## set up storage for prediction means, variances, etc.
if (pred.mean)
pred.m <- matrix(
@@ -138,8 +159,8 @@
dimnames=list(c(statenames,rw.names),NULL)
)
else
- pred.m <- array(dim=c(0,0))
-
+ pred.m <- array(data=numeric(0),dim=c(0,0))
+
if (pred.var)
pred.v <- matrix(
data=0,
@@ -148,7 +169,7 @@
dimnames=list(c(statenames,rw.names),NULL)
)
else
- pred.v <- array(dim=c(0,0))
+ pred.v <- array(data=numeric(0),dim=c(0,0))
if (filter.mean)
if (random.walk)
@@ -166,13 +187,31 @@
dimnames=list(statenames,NULL)
)
else
- filt.m <- array(dim=c(0,0))
+ filt.m <- array(data=numeric(0),dim=c(0,0))
+ if (mif2) {
+ if (missing(cooling.fraction))
+ stop("pfilter error: ",sQuote("cooling.fraction")," must be specified for method mif2",call.=FALSE)
+ cooling.fraction <- as.numeric(cooling.fraction)
+ }
+
for (nt in seq_len(ntimes)) {
-
+
+ if (mif2) {
+ cool.sched <- try(
+ mif2.cooling(cooling.fraction,nt,cooling.m+.ndone,ntimes),
+ silent=FALSE
+ )
+ if (inherits(cool.sched,"try-error"))
+ stop("pfilter error: cooling schedule error",call.=FALSE)
+ sigma1 <- sigma*cool.sched$alpha
+ } else {
+ sigma1 <- sigma
+ }
+
## transform the parameters if necessary
if (transform) tparams <- partrans(object,params,dir="forward")
-
+
## advance the state variables according to the process model
X <- try(
rprocess(
@@ -186,7 +225,7 @@
)
if (inherits(X,'try-error'))
stop(sQuote("pfilter")," error: process simulation error",call.=FALSE)
-
+
if (pred.var) { ## check for nonfinite state variables and parameters
problem.indices <- unique(which(!is.finite(X),arr.ind=TRUE)[,1])
if (length(problem.indices)>0) { # state variables
@@ -225,7 +264,7 @@
if (any(!is.finite(weights))) {
stop(sQuote("pfilter")," error: ",sQuote("dmeasure")," returns non-finite value",call.=FALSE)
}
-
+
## compute prediction mean, prediction variance, filtering mean,
## effective sample size, log-likelihood
## also do resampling if filtering has not failed
@@ -233,7 +272,8 @@
.Call(
pfilter_computations,
X,params,Np[nt+1],
- random.walk,sigma,
+ random.walk,
+ sigma1,
pred.mean,pred.var,
filter.mean,one.par,
weights,tol
@@ -246,17 +286,17 @@
all.fail <- xx$fail
loglik[nt] <- xx$loglik
eff.sample.size[nt] <- xx$ess
-
+
x <- xx$states
params <- xx$params
-
+
if (pred.mean)
pred.m[,nt] <- xx$pm
if (pred.var)
pred.v[,nt] <- xx$pv
if (filter.mean)
filt.m[,nt] <- xx$fm
-
+
if (all.fail) { ## all particles are lost
nfail <- nfail+1
if (verbose)
@@ -268,16 +308,16 @@
if (save.states) {
xparticles[[nt]] <- x
}
-
+
if (save.params) {
pparticles[[nt]] <- params
}
-
+
if (verbose && (nt%%5==0))
cat("pfilter timestep",nt,"of",ntimes,"finished\n")
-
+
}
-
+
if (!is.null(seed)) {
assign(".Random.seed",save.seed,pos=.GlobalEnv)
seed <- save.seed
@@ -289,6 +329,7 @@
pred.mean=pred.m,
pred.var=pred.v,
filter.mean=filt.m,
+ paramMatrix=if (mif2) params else array(data=numeric(0),dim=c(0,0)),
eff.sample.size=eff.sample.size,
cond.loglik=loglik,
saved.states=xparticles,
Added: pkg/pomp/R/pomp-class.R
===================================================================
--- pkg/pomp/R/pomp-class.R (rev 0)
+++ pkg/pomp/R/pomp-class.R 2013-01-08 03:21:32 UTC (rev 793)
@@ -0,0 +1,148 @@
+## as of version 0.37-1 'pomp' is a generic function
+setGeneric("pomp",function(data,...)standardGeneric("pomp"))
+
+## this is the initial-condition setting function that is used by default
+## it simply finds all parameters in the vector 'params' that have a name ending in '.0'
+## and returns a vector with their values with names stripped of '.0'
+default.initializer <- function (params, t0, ...) {
+ ivpnames <- grep("\\.0$",names(params),value=TRUE)
+ if (length(ivpnames)<1)
+ stop("default initializer error: no parameter names ending in ",
+ sQuote(".0")," found: see ",sQuote("pomp")," documentation")
+ x <- params[ivpnames]
+ names(x) <- sub("\\.0$","",ivpnames)
+ x
+}
+
+## define the pomp class
+setClass(
+ 'pomp',
+ representation=representation(
+ data = 'array',
+ times = 'numeric',
+ t0 = 'numeric',
+ rprocess = 'function',
+ dprocess = 'function',
+ dmeasure = 'pomp.fun',
+ rmeasure = 'pomp.fun',
+ skeleton.type = 'character',
+ skeleton = 'pomp.fun',
+ skelmap.delta.t = 'numeric',
+ initializer = 'function',
+ states = 'array',
+ params = 'numeric',
+ covar = 'matrix',
+ tcovar = 'numeric',
+ obsnames = 'character',
+ statenames = 'character',
+ paramnames = 'character',
+ covarnames = 'character',
+ zeronames = 'character',
+ has.trans = 'logical',
+ par.trans = 'pomp.fun',
+ par.untrans = 'pomp.fun',
+ PACKAGE = 'character',
+ userdata = 'list'
+ ),
+ prototype=prototype(
+ data=array(data=numeric(0),dim=c(0,0)),
+ times=numeric(0),
+ t0=numeric(0),
+ rprocess=function(xstart,times,params,...)stop(sQuote("rprocess")," not specified"),
+ dprocess=function(x,times,params,log=FALSE,...)stop(sQuote("dprocess")," not specified"),
+ dmeasure=pomp.fun(),
+ rmeasure=pomp.fun(),
+ skeleton.type="map",
+ skeleton=pomp.fun(),
+ skelmap.delta.t=as.numeric(NA),
+ initializer=default.initializer,
+ states=array(data=numeric(0),dim=c(0,0)),
+ params=numeric(0),
+ covar=array(data=numeric(0),dim=c(0,0)),
+ tcovar=numeric(0),
+ obsnames=character(0),
+ statenames=character(0),
+ paramnames=character(0),
+ covarnames=character(0),
+ zeronames=character(0),
+ has.trans=FALSE,
+ par.trans=pomp.fun(),
+ par.untrans=pomp.fun(),
+ PACKAGE="",
+ userdata=list()
+ ),
+ validity=function (object) {
+ retval <- character(0)
+ if (length(object at data)<1)
+ retval <- append(retval,paste(sQuote("data"),"is a required argument"))
+ if (length(object at times)<1)
+ retval <- append(retval,paste(sQuote("times"),"is a required argument"))
+ if (length(object at t0)<1)
+ retval <- append(retval,paste(sQuote("t0"),"is a required argument"))
+ if (!is.numeric(object at params) || (length(object at params)>0 && is.null(names(object at params))))
+ retval <- append(retval,paste(sQuote("params"),"must be a named numeric vector"))
+ if (ncol(object at data)!=length(object at times))
+ retval <- append(retval,paste("the length of",sQuote("times"),"should match the number of observations"))
+ if (length(object at t0)<1)
+ retval <- append(retval,paste(sQuote("t0"),"is a required argument"))
+ if (length(object at t0)>1)
+ retval <- append(retval,paste(sQuote("t0"),"must be a single number"))
+ if (object at t0 > object at times[1])
+ retval <- append(retval,paste("the zero-time",sQuote("t0"),"must occur no later than the first observation"))
+ if (object at skelmap.delta.t <= 0)
+ retval <- append(retval,paste(sQuote("skelmap.delta.t"),"must be positive"))
+ if (!all(c('xstart','times','params','...')%in%names(formals(object at rprocess))))
+ retval <- append(
+ retval,
+ paste(
+ sQuote("rprocess"),"must be a function of prototype",
+ sQuote("rprocess(xstart,times,params,...)")
+ )
+ )
+ if (!all(c('x','times','params','log','...')%in%names(formals(object at dprocess))))
+ retval <- append(
+ retval,
+ paste(
+ sQuote("dprocess"),"must be a function of prototype",
+ sQuote("dprocess(x,times,params,log,...)")
+ )
+ )
+ if (!all(c('params','t0','...')%in%names(formals(object at initializer))))
+ retval <- append(
+ retval,
+ paste(
+ sQuote("initializer"),"must be a function of prototype",
+ sQuote("initializer(params,t0,...)")
+ )
+ )
+ if (length(object at tcovar)!=nrow(object at covar)) {
+ retval <- append(
+ retval,
+ paste(
+ "the length of",sQuote("tcovar"),
+ "should match the number of rows of",sQuote("covar")
+ )
+ )
+ } else if (!all(object at covarnames%in%colnames(object at covar))) {
+ missing <- object at covarnames[!(object at covarnames%in%colnames(object at covar))]
+ retval <- append(
+ retval,
+ paste(
+ "covariate(s)",paste(missing,collapse=","),
+ "are not found among the columns of",sQuote("covar")
+ )
+ )
+ }
+ if (!is.numeric(object at tcovar))
+ retval <- append(
+ retval,
+ paste(
+ sQuote("tcovar"),
+ "must either be a numeric vector or must name a numeric vector in the data frame",
+ sQuote("covar")
+ )
+ )
+
+ if (length(retval)==0) TRUE else retval
+ }
+ )
Modified: pkg/pomp/R/pomp.R
===================================================================
--- pkg/pomp/R/pomp.R 2013-01-07 20:15:21 UTC (rev 792)
+++ pkg/pomp/R/pomp.R 2013-01-08 03:21:32 UTC (rev 793)
@@ -1,50 +1,3 @@
-## define the pomp class
-setClass(
- 'pomp',
- representation(
- data = 'array',
- times = 'numeric',
- t0 = 'numeric',
- rprocess = 'function',
- dprocess = 'function',
- dmeasure = 'pomp.fun',
- rmeasure = 'pomp.fun',
- skeleton.type = 'character',
- skeleton = 'pomp.fun',
- skelmap.delta.t = 'numeric',
- initializer = 'function',
- states = 'array',
- params = 'numeric',
- covar = 'matrix',
- tcovar = 'numeric',
- obsnames = 'character',
- statenames = 'character',
- paramnames = 'character',
- covarnames = 'character',
- zeronames = 'character',
- has.trans = 'logical',
- par.trans = 'pomp.fun',
- par.untrans = 'pomp.fun',
- PACKAGE = 'character',
- userdata = 'list'
- )
- )
-
-## this is the initial-condition setting function that is used by default
-## it simply finds all parameters in the vector 'params' that have a name ending in '.0'
-## and returns a vector with their values with names stripped of '.0'
-default.initializer <- function (params, t0, ...) {
- ivpnames <- grep("\\.0$",names(params),value=TRUE)
- if (length(ivpnames)<1)
- stop("default initializer error: no parameter names ending in ",sQuote(".0")," found: see ",sQuote("pomp")," documentation")
- x <- params[ivpnames]
- names(x) <- sub("\\.0$","",ivpnames)
- x
-}
-
-## as of version 0.37-1 'pomp' is a generic function
-setGeneric("pomp",function(data,...)standardGeneric("pomp"))
-
## basic constructor of the pomp class
pomp.constructor <- function (data, times, t0, ..., rprocess, dprocess,
rmeasure, dmeasure, measurement.model,
More information about the pomp-commits
mailing list