[Pomp-commits] r891 - in pkg/pomp: . R demo inst inst/include man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Mar 16 22:57:33 CET 2014
Author: kingaa
Date: 2014-03-16 22:57:33 +0100 (Sun, 16 Mar 2014)
New Revision: 891
Added:
pkg/pomp/R/dprior-pomp.R
pkg/pomp/R/rprior-pomp.R
pkg/pomp/man/prior-pomp.Rd
pkg/pomp/src/dprior.c
pkg/pomp/src/rprior.c
pkg/pomp/tests/prior.R
pkg/pomp/tests/prior.Rout.save
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/NAMESPACE
pkg/pomp/R/aaa.R
pkg/pomp/R/abc-methods.R
pkg/pomp/R/abc.R
pkg/pomp/R/bsmc.R
pkg/pomp/R/mif-methods.R
pkg/pomp/R/mif.R
pkg/pomp/R/pmcmc.R
pkg/pomp/R/pomp-class.R
pkg/pomp/R/pomp-methods.R
pkg/pomp/R/pomp.R
pkg/pomp/R/rmeasure-pomp.R
pkg/pomp/demo/sir.R
pkg/pomp/inst/NEWS
pkg/pomp/inst/include/pomp.h
pkg/pomp/man/abc.Rd
pkg/pomp/man/bsmc.Rd
pkg/pomp/man/pmcmc-methods.Rd
pkg/pomp/man/pmcmc.Rd
pkg/pomp/man/pomp-class.Rd
pkg/pomp/man/pomp.Rd
pkg/pomp/src/pomp.h
pkg/pomp/tests/abc.R
pkg/pomp/tests/abc.Rout.save
pkg/pomp/tests/bbs-trajmatch.Rout.save
pkg/pomp/tests/bbs.R
pkg/pomp/tests/bbs.Rout.save
pkg/pomp/tests/blowflies.Rout.save
pkg/pomp/tests/dacca.R
pkg/pomp/tests/dacca.Rout.save
pkg/pomp/tests/dimchecks.Rout.save
pkg/pomp/tests/fhn.Rout.save
pkg/pomp/tests/filtfail.Rout.save
pkg/pomp/tests/gillespie.Rout.save
pkg/pomp/tests/gompertz.Rout.save
pkg/pomp/tests/logistic.Rout.save
pkg/pomp/tests/ou2-bsmc.R
pkg/pomp/tests/ou2-bsmc.Rout.save
pkg/pomp/tests/ou2-forecast.Rout.save
pkg/pomp/tests/ou2-icfit.R
pkg/pomp/tests/ou2-icfit.Rout.save
pkg/pomp/tests/ou2-kalman.Rout.save
pkg/pomp/tests/ou2-mif-fp.R
pkg/pomp/tests/ou2-mif-fp.Rout.save
pkg/pomp/tests/ou2-mif.R
pkg/pomp/tests/ou2-mif.Rout.save
pkg/pomp/tests/ou2-mif2.R
pkg/pomp/tests/ou2-mif2.Rout.save
pkg/pomp/tests/ou2-nlf.Rout.save
pkg/pomp/tests/ou2-pmcmc.R
pkg/pomp/tests/ou2-pmcmc.Rout.save
pkg/pomp/tests/ou2-probe.Rout.save
pkg/pomp/tests/ou2-procmeas.Rout.save
pkg/pomp/tests/ou2-simulate.Rout.save
pkg/pomp/tests/ou2-trajmatch.Rout.save
pkg/pomp/tests/partrans.Rout.save
pkg/pomp/tests/pfilter.Rout.save
pkg/pomp/tests/pomppomp.Rout.save
pkg/pomp/tests/ricker-bsmc.R
pkg/pomp/tests/ricker-bsmc.Rout.save
pkg/pomp/tests/ricker-probe.R
pkg/pomp/tests/ricker-probe.Rout.save
pkg/pomp/tests/ricker-spect.Rout.save
pkg/pomp/tests/ricker.Rout.save
pkg/pomp/tests/rw2.Rout.save
pkg/pomp/tests/sir.Rout.save
pkg/pomp/tests/skeleton.Rout.save
pkg/pomp/tests/steps.Rout.save
pkg/pomp/tests/synlik.Rout.save
pkg/pomp/tests/verhulst.Rout.save
Log:
- new 'rprior' and 'dprior' slots in 'pomp' objects
- use environment variable POMP_FULL_TESTS to suppress some test-script execution
- change in 'bsmc' interface: now 'rprior' is used
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2014-03-11 13:23:39 UTC (rev 890)
+++ pkg/pomp/DESCRIPTION 2014-03-16 21:57:33 UTC (rev 891)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.47-5
-Date: 2014-03-10
+Version: 0.48-1
+Date: 2014-03-16
Authors at R: c(person(given=c("Aaron","A."),family="King",
role=c("aut","cre"),email="kingaa at umich.edu"),
person(given=c("Edward","L."),family="Ionides",role=c("aut")),
@@ -29,6 +29,7 @@
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
+ dprior-pomp.R rprior-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
Modified: pkg/pomp/NAMESPACE
===================================================================
--- pkg/pomp/NAMESPACE 2014-03-11 13:23:39 UTC (rev 890)
+++ pkg/pomp/NAMESPACE 2014-03-16 21:57:33 UTC (rev 891)
@@ -25,6 +25,8 @@
do_dprocess,
do_rmeasure,
do_dmeasure,
+ do_rprior,
+ do_dprior,
do_skeleton,
do_init_state
)
@@ -51,15 +53,15 @@
pomp,
plot,show,print,coerce,summary,logLik,window,"$",
dprocess,rprocess,rmeasure,dmeasure,init.state,skeleton,
- data.array,obs,partrans,coef,"coef<-",time,"time<-",timezero,"timezero<-",
+ dprior,rprior,
+ data.array,obs,partrans,coef,"coef<-",
+ time,"time<-",timezero,"timezero<-",
simulate,pfilter,
eff.sample.size,cond.logLik,
particles,mif,continue,states,trajectory,
pred.mean,pred.var,filter.mean,conv.rec,
- bsmc,
- pmcmc,dprior,
- spect,probe,
- probe.match,abc,
+ bsmc,pmcmc,abc,
+ spect,probe,probe.match,
traj.match
)
Modified: pkg/pomp/R/aaa.R
===================================================================
--- pkg/pomp/R/aaa.R 2014-03-11 13:23:39 UTC (rev 890)
+++ pkg/pomp/R/aaa.R 2014-03-16 21:57:33 UTC (rev 891)
@@ -19,6 +19,7 @@
setGeneric("filter.mean",function(object,...)standardGeneric("filter.mean"))
setGeneric("cond.logLik",function(object,...)standardGeneric("cond.logLik"))
setGeneric("eff.sample.size",function(object,...)standardGeneric("eff.sample.size"))
+setGeneric("conv.rec",function(object,...)standardGeneric("conv.rec"))
if (!exists("paste0",where="package:base")) {
paste0 <- function(...) paste(...,sep="")
Modified: pkg/pomp/R/abc-methods.R
===================================================================
--- pkg/pomp/R/abc-methods.R 2014-03-11 13:23:39 UTC (rev 890)
+++ pkg/pomp/R/abc-methods.R 2014-03-16 21:57:33 UTC (rev 891)
@@ -23,18 +23,3 @@
}
}
)
-
-setMethod(
- "dprior",
- signature=signature(object="abc"),
- function (object, params, log = FALSE, ...) {
- do.call(
- object at dprior,
- list(
- params=params,
- hyperparams=object at hyperparams,
- log=log
- )
- )
- }
- )
Modified: pkg/pomp/R/abc.R
===================================================================
--- pkg/pomp/R/abc.R 2014-03-11 13:23:39 UTC (rev 890)
+++ pkg/pomp/R/abc.R 2014-03-16 21:57:33 UTC (rev 891)
@@ -6,11 +6,9 @@
pars = 'character',
transform = 'logical',
Nabc = 'integer',
- dprior = 'function',
probes='list',
scale = 'numeric',
epsilon = 'numeric',
- hyperparams = 'list',
random.walk.sd = 'numeric',
conv.rec = 'matrix'
)
@@ -20,19 +18,20 @@
setGeneric('abc',function(object,...)standardGeneric("abc"))
abc.internal <- function (object, Nabc,
- start, pars, dprior.fun,
+ start, pars,
rw.sd, probes,
- hyperparams,
epsilon, scale,
verbose, transform,
- .ndone, .getnativesymbolinfo = TRUE,
+ .ndone = 0L,
+ .getnativesymbolinfo = TRUE,
...) {
+ object <- as(object,'pomp')
gnsi <- as.logical(.getnativesymbolinfo)
-
transform <- as.logical(transform)
-
Nabc <- as.integer(Nabc)
+ epsilon <- as.numeric(epsilon)
+ epssq <- epsilon*epsilon
if (length(start)==0)
stop(
@@ -64,7 +63,7 @@
stop(
"abc error: ",
sQuote("pars"),
- " must be a mutually disjoint subset of ",
+ " must be a subset of ",
sQuote("names(start)"),
" and must correspond to positive random-walk SDs specified in ",
sQuote("rw.sd"),
@@ -101,7 +100,7 @@
}
theta <- start
- log.prior <- dprior.fun(params=theta,hyperparams=hyperparams,log=TRUE)
+ log.prior <- dprior(object,params=theta,log=TRUE)
## we suppose that theta is a "match", which does the right thing for continue() and
## should have negligible effect unless doing many short calls to continue()
@@ -161,8 +160,8 @@
## ABC update rule
distance <- sum(((datval-simval)/scale)^2)
- if( (is.finite(distance)) && (distance<epsilon) ){
- log.prior.prop <- dprior.fun(params=theta.prop,hyperparams=hyperparams,log=TRUE)
+ if( (is.finite(distance)) && (distance<epssq) ){
+ log.prior.prop <- dprior(object,params=theta.prop,log=TRUE)
if (runif(1) < exp(log.prior.prop-log.prior)) {
theta <- theta.prop
log.prior <- log.prior.prop
@@ -180,15 +179,13 @@
'abc',
po,
params=theta,
- pars = pars,
- transform = transform,
- Nabc = Nabc,
- dprior = dprior.fun,
+ pars=pars,
+ transform=transform,
+ Nabc=Nabc,
probes=probes,
- scale = scale,
- epsilon = epsilon,
- hyperparams = hyperparams,
- random.walk.sd = rw.sd,
+ scale=scale,
+ epsilon=epsilon,
+ random.walk.sd=rw.sd,
conv.rec=conv.rec
)
@@ -199,21 +196,19 @@
signature=signature(object="pomp"),
function (object, Nabc = 1,
start, pars, rw.sd,
- dprior, probes, scale, epsilon, hyperparams,
+ probes, scale, epsilon,
verbose = getOption("verbose"),
transform = FALSE,
...) {
- transform <- as.logical(transform)
+ if (missing(start))
+ start <- coef(object)
- if (missing(start)) start <- coef(object,transform=transform)
-
if (missing(rw.sd))
stop("abc error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
- if (missing(pars)) {
+ if (missing(pars))
pars <- names(rw.sd)[rw.sd>0]
- }
if (missing(probes))
stop("abc error: ",sQuote("probes")," must be specified",call.=FALSE)
@@ -222,41 +217,19 @@
stop("abc error: ",sQuote("scale")," must be specified",call.=FALSE)
if (missing(epsilon))
- stop("abc error: ",sQuote("abc match criterion, epsilon,")," must be specified",call.=FALSE)
+ stop("abc error: abc match criterion, ",sQuote("epsilon"),", must be specified",call.=FALSE)
- if (missing(hyperparams))
- stop("abc error: ",sQuote("hyperparams")," must be specified",call.=FALSE)
-
- if (missing(dprior)) { # use default flat improper prior
- dprior <- function (params, hyperparams, log) {
- if (log) 0 else 1
- }
- } else {
- dprior <- match.fun(dprior)
- if (!all(c('params','hyperparams','log')%in%names(formals(dprior))))
- stop(
- "abc error: ",
- sQuote("dprior"),
- " must be a function of prototype ",
- sQuote("dprior(params,hyperparams,log)"),
- call.=FALSE
- )
- }
-
abc.internal(
object=object,
Nabc=Nabc,
start=start,
pars=pars,
- dprior.fun=dprior,
probes=probes,
scale=scale,
epsilon=epsilon,
rw.sd=rw.sd,
- hyperparams=hyperparams,
verbose=verbose,
- transform=transform,
- .ndone=0
+ transform=transform
)
}
)
@@ -264,62 +237,19 @@
setMethod(
"abc",
signature=signature(object="probed.pomp"),
- function (object, Nabc = 1,
- start, pars, rw.sd,
- dprior, probes, scale, epsilon, hyperparams,
+ function (object, probes,
verbose = getOption("verbose"),
transform = FALSE,
...) {
- transform <- as.logical(transform)
-
- if (missing(start)) start <- coef(object,transform=transform)
-
- if (missing(rw.sd))
- stop("abc error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
-
- if (missing(pars)) {
- pars <- names(rw.sd)[rw.sd>0]
- }
-
if (missing(probes)) probes <- object at probes
- if (missing(scale)) probes <- object at scale
- if (missing(epsilon)) probes <- object at epsilon
- if (missing(hyperparams))
- stop("abc error: ",sQuote("hyperparams")," must be specified",call.=FALSE)
-
- if (missing(dprior)) { # use default flat improper prior
- dprior <- function (params, hyperparams, log) {
- if (log) 0 else 1
- }
- } else {
- dprior <- match.fun(dprior)
- if (!all(c('params','hyperparams','log')%in%names(formals(dprior))))
- stop(
- "abc error: ",
- sQuote("dprior"),
- " must be a function of prototype ",
- sQuote("dprior(params,hyperparams,log)"),
- call.=FALSE
- )
- }
-
- abc.internal(
- object=as(object,"pomp"),
- Nabc=Nabc,
- start=start,
- pars=pars,
- dprior.fun=dprior,
- probes=probes,
- scale=scale,
- epsilon=epsilon,
- rw.sd=rw.sd,
- hyperparams=hyperparams,
- verbose=verbose,
- transform=transform,
- .ndone=0
- )
+ abc(
+ object=as(object,"pomp"),
+ probes=probes,
+ transform=transform,
+ ...
+ )
}
)
@@ -328,7 +258,7 @@
signature=signature(object="abc"),
function (object, Nabc,
start, pars, rw.sd,
- dprior, probes, scale, epsilon, hyperparams,
+ probes, scale, epsilon,
verbose = getOption("verbose"),
transform,
...) {
@@ -336,30 +266,25 @@
if (missing(Nabc)) Nabc <- object at Nabc
if (missing(start)) start <- coef(object)
if (missing(pars)) pars <- object at pars
- if (missing(rw.sd)) pars <- object at random.walk.sd
- if (missing(dprior)) dprior <- object at dprior
+ if (missing(rw.sd)) rw.sd <- object at random.walk.sd
if (missing(probes)) probes <- object at probes
- if (missing(scale)) probes <- object at scale
- if (missing(epsilon)) probes <- object at epsilon
- if (missing(hyperparams)) hyperparams <- object at hyperparams
+ if (missing(scale)) scale <- object at scale
+ if (missing(epsilon)) epsilon <- object at epsilon
if (missing(transform)) transform <- object at transform
- transform <- as.logical(transform)
- abc.internal(
- object=as(object,"pomp"),
- Nabc=Nabc,
- start=start,
- pars=pars,
- dprior.fun=dprior,
- rw.sd=rw.sd,
- probes=probes,
- scale=scale,
- epsilon=epsilon,
- hyperparams=hyperparams,
- verbose=verbose,
- transform=transform,
- .ndone=0
- )
+ abc(
+ object=as(object,"pomp"),
+ Nabc=Nabc,
+ start=start,
+ pars=pars,
+ rw.sd=rw.sd,
+ probes=probes,
+ scale=scale,
+ epsilon=epsilon,
+ verbose=verbose,
+ transform=transform,
+ ...
+ )
}
)
@@ -382,6 +307,7 @@
obj at conv.rec[-1,]
)
obj at Nabc <- as.integer(ndone+Nabc)
+
obj
}
)
Modified: pkg/pomp/R/bsmc.R
===================================================================
--- pkg/pomp/R/bsmc.R 2014-03-11 13:23:39 UTC (rev 890)
+++ pkg/pomp/R/bsmc.R 2014-03-16 21:57:33 UTC (rev 891)
@@ -68,6 +68,9 @@
if (missing(Np)) Np <- NCOL(params)
else if (is.matrix(params)&&(Np!=ncol(params)))
warning(sQuote("Np")," is ignored when ",sQuote("params")," is a matrix")
+
+ if ((!is.matrix(params)) && (Np > 1))
+ params <- rprior(object,params=parmat(params,Np))
if (transform)
params <- partrans(object,params,dir="inverse",
Added: pkg/pomp/R/dprior-pomp.R
===================================================================
--- pkg/pomp/R/dprior-pomp.R (rev 0)
+++ pkg/pomp/R/dprior-pomp.R 2014-03-16 21:57:33 UTC (rev 891)
@@ -0,0 +1,12 @@
+## evaluate the prior probability density
+setGeneric("dprior",function(object,...)standardGeneric("dprior"))
+
+dprior.internal <- function (object, params, log = FALSE,
+ .getnativesymbolinfo = TRUE, ...) {
+ .Call(do_dprior,object,params,log,.getnativesymbolinfo)
+}
+
+setMethod("dprior","pomp",
+ function (object, params, log = FALSE, ...)
+ dprior.internal(object=object,params=params,log=log,...)
+ )
Modified: pkg/pomp/R/mif-methods.R
===================================================================
--- pkg/pomp/R/mif-methods.R 2014-03-11 13:23:39 UTC (rev 890)
+++ pkg/pomp/R/mif-methods.R 2014-03-16 21:57:33 UTC (rev 891)
@@ -38,8 +38,6 @@
}
}
-setGeneric("conv.rec",function(object,...)standardGeneric("conv.rec"))
-
setMethod('conv.rec','mif',
function (object, pars, transform = FALSE, ...) {
conv.rec.internal(object=object,pars=pars,transform=transform,...)
Modified: pkg/pomp/R/mif.R
===================================================================
--- pkg/pomp/R/mif.R 2014-03-11 13:23:39 UTC (rev 890)
+++ pkg/pomp/R/mif.R 2014-03-16 21:57:33 UTC (rev 891)
@@ -404,7 +404,6 @@
transform = FALSE,
...) {
- transform <- as.logical(transform)
method <- match.arg(method)
if (missing(start)) start <- coef(object)
@@ -515,7 +514,6 @@
if (missing(cooling.fraction)) cooling.fraction <- object at cooling.fraction
if (missing(method)) method <- object at method
if (missing(transform)) transform <- object at transform
- transform <- as.logical(transform)
if (missing(Np)) Np <- object at Np
if (missing(tol)) tol <- object at tol
Modified: pkg/pomp/R/pmcmc.R
===================================================================
--- pkg/pomp/R/pmcmc.R 2014-03-11 13:23:39 UTC (rev 890)
+++ pkg/pomp/R/pmcmc.R 2014-03-16 21:57:33 UTC (rev 891)
@@ -6,8 +6,6 @@
pars = 'character',
transform = 'logical',
Nmcmc = 'integer',
- dprior = 'function',
- hyperparams = 'list',
random.walk.sd = 'numeric',
conv.rec = 'matrix',
log.prior = 'numeric'
@@ -17,28 +15,19 @@
## PMCMC algorithm functions
setGeneric('pmcmc',function(object,...)standardGeneric("pmcmc"))
-setGeneric('dprior', function(object,...)standardGeneric("dprior"))
-
-setMethod(
- "dprior",
- signature=signature(object="pmcmc"),
- function (object, params, log = FALSE, ...) {
- do.call(object at dprior,list(params=params,hyperparams=object at hyperparams,log=log))
- }
- )
-
pmcmc.internal <- function (object, Nmcmc,
- start, pars, dprior.fun,
+ start, pars,
rw.sd, Np,
- hyperparams,
tol, max.fail,
verbose, transform,
.ndone = 0L,
.prev.pfp = NULL, .prev.log.prior = NULL,
.getnativesymbolinfo = TRUE) {
+ object <- as(object,"pomp")
gnsi <- as.logical(.getnativesymbolinfo)
transform <- as.logical(transform)
+ .ndone <- as.integer(.ndone)
if (missing(start))
stop(sQuote("start")," must be specified",call.=FALSE)
@@ -96,9 +85,6 @@
rw.sd <- rw.sd[pars]
rw.names <- names(rw.sd)
- if (missing(dprior.fun))
- stop("pmcmc error: ",sQuote("dprior")," must be specified",call.=FALSE)
-
ntimes <- length(time(object))
if (missing(Np))
stop("pmcmc error: ",sQuote("Np")," must be specified",call.=FALSE)
@@ -120,9 +106,6 @@
stop(sQuote("Np")," must be a number, a vector of numbers, or a function")
Np <- as.integer(Np)
- if (missing(hyperparams))
- stop("pmcmc error: ",sQuote("hyperparams")," must be specified",call.=FALSE)
-
if (missing(Nmcmc))
stop("pmcmc error: ",sQuote("Nmcmc")," must be specified",call.=FALSE)
Nmcmc <- as.integer(Nmcmc)
@@ -161,23 +144,11 @@
)
}
- obj <- as(object,"pomp")
-
- if (Nmcmc>0)
- tmp.pmcmc <- new("pmcmc",obj,dprior=dprior.fun,hyperparams=hyperparams)
- else
- pfp <- obj
-
- if (.ndone==0) { ## compute prior and likelihood on initial parameter vector
+ if (.ndone==0L) { ## compute prior and likelihood on initial parameter vector
pfp <- try(
pfilter.internal(
- object=obj,
- params=if (transform) {
- partrans(obj,theta,dir="forward",
- .getnativesymbolinfo=gnsi)
- } else {
- theta
- },
+ object=object,
+ params=theta,
Np=Np,
tol=tol,
max.fail=max.fail,
@@ -194,7 +165,8 @@
)
if (inherits(pfp,'try-error'))
stop("pmcmc error: error in ",sQuote("pfilter"),call.=FALSE)
- log.prior <- dprior(object=tmp.pmcmc,params=theta,log=TRUE)
+ log.prior <- dprior(object,params=theta,log=TRUE,.getnativesymbolinfo=gnsi)
+ gnsi <- FALSE
} else { ## has been computed previously
pfp <- .prev.pfp
log.prior <- .prev.log.prior
@@ -205,18 +177,17 @@
for (n in seq_len(Nmcmc)) { # main loop
theta.prop <- theta
+ if (transform)
+ theta <- partrans(object,theta.prop,dir='inverse',.getnativesymbolinfo=gnsi)
theta.prop[pars] <- rnorm(n=length(pars),mean=theta.prop[pars],sd=rw.sd)
+ if (transform)
+ theta <- partrans(object,theta.prop,dir='forward',.getnativesymbolinfo=gnsi)
## run the particle filter on the proposed new parameter values
pfp.prop <- try(
pfilter.internal(
object=pfp,
- params=if (transform) {
- partrans(obj,theta.prop,dir="forward",
- .getnativesymbolinfo=gnsi)
- } else {
- theta.prop
- },
+ params=theta.prop,
Np=Np,
tol=tol,
max.fail=max.fail,
@@ -233,7 +204,7 @@
)
if (inherits(pfp.prop,'try-error'))
stop("pmcmc error: error in ",sQuote("pfilter"),call.=FALSE)
- log.prior.prop <- dprior(object=tmp.pmcmc,params=theta.prop,log=TRUE)
+ log.prior.prop <- dprior(object,params=theta.prop,log=TRUE,.getnativesymbolinfo=gnsi)
gnsi <- FALSE
## PMCMC update rule (OK because proposal is symmetric)
@@ -258,10 +229,8 @@
transform=transform,
Nmcmc=Nmcmc,
pars=pars,
- dprior=dprior.fun,
random.walk.sd=rw.sd,
Np=Np,
- hyperparams=hyperparams,
tol=tol,
conv.rec=conv.rec,
log.prior=log.prior
@@ -272,47 +241,26 @@
"pmcmc",
signature=signature(object="pomp"),
function (object, Nmcmc = 1,
- start, pars, rw.sd,
- dprior, Np, hyperparams,
+ start, pars, rw.sd, Np,
tol = 1e-17, max.fail = 0,
verbose = getOption("verbose"),
transform = FALSE,
...) {
- transform <- as.logical(transform)
if (missing(start)) start <- coef(object,transform=transform)
if (missing(rw.sd))
stop("pmcmc error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
if (missing(pars)) pars <- names(rw.sd)[rw.sd>0]
if (missing(Np))
stop("pmcmc error: ",sQuote("Np")," must be specified",call.=FALSE)
- if (missing(hyperparams))
- stop("pmcmc error: ",sQuote("hyperparams")," must be specified",call.=FALSE)
- if (missing(dprior)) { # use default flat improper prior
- dprior <- function (params, hyperparams, log) {
- if (log) 0 else 1
- }
- } else {
- dprior <- match.fun(dprior)
- if (!all(c('params','hyperparams','log')%in%names(formals(dprior))))
- stop(
- "pmcmc error: ",
- sQuote("dprior"),
- " must be a function of prototype ",
- sQuote("dprior(params,hyperparams,log)"),
- call.=FALSE
- )
- }
pmcmc.internal(
object=object,
Nmcmc=Nmcmc,
start=start,
pars=pars,
- dprior.fun=dprior,
rw.sd=rw.sd,
Np=Np,
- hyperparams=hyperparams,
tol=tol,
max.fail=max.fail,
verbose=verbose,
@@ -325,8 +273,7 @@
setMethod(
"pmcmc",
signature=signature(object="pfilterd.pomp"),
- function (object, Nmcmc = 1,
- Np, tol, ...) {
+ function (object, Nmcmc = 1, Np, tol, ...) {
if (missing(Np)) Np <- object at Np
if (missing(tol)) tol <- object at tol
@@ -346,8 +293,7 @@
signature=signature(object="pmcmc"),
function (object, Nmcmc,
start, pars, rw.sd,
- dprior, Np, hyperparams,
- tol, max.fail = 0,
+ Np, tol, max.fail = 0,
verbose = getOption("verbose"),
transform,
...) {
@@ -356,22 +302,17 @@
if (missing(start)) start <- coef(object)
if (missing(pars)) pars <- object at pars
if (missing(rw.sd)) rw.sd <- object at random.walk.sd
- if (missing(dprior)) dprior <- object at dprior
if (missing(Np)) Np <- object at Np
- if (missing(hyperparams)) hyperparams <- object at hyperparams
if (missing(tol)) tol <- object at tol
if (missing(transform)) transform <- object at transform
- transform <- as.logical(transform)
pmcmc(
object=as(object,"pomp"),
Nmcmc=Nmcmc,
start=start,
pars=pars,
- dprior=dprior,
rw.sd=rw.sd,
Np=Np,
- hyperparams=hyperparams,
tol=tol,
max.fail=max.fail,
verbose=verbose,
@@ -384,8 +325,7 @@
setMethod(
'continue',
signature=signature(object='pmcmc'),
- function (object, Nmcmc = 1,
- ...) {
+ function (object, Nmcmc = 1, ...) {
ndone <- object at Nmcmc
Modified: pkg/pomp/R/pomp-class.R
===================================================================
--- pkg/pomp/R/pomp-class.R 2014-03-11 13:23:39 UTC (rev 890)
+++ pkg/pomp/R/pomp-class.R 2014-03-16 21:57:33 UTC (rev 891)
@@ -25,6 +25,8 @@
dprocess = 'function',
dmeasure = 'pomp.fun',
rmeasure = 'pomp.fun',
+ dprior = 'pomp.fun',
+ rprior = 'pomp.fun',
skeleton.type = 'character',
skeleton = 'pomp.fun',
skelmap.delta.t = 'numeric',
@@ -52,6 +54,8 @@
dprocess=function(x,times,params,log=FALSE,...)stop(sQuote("dprocess")," not specified"),
dmeasure=pomp.fun(),
rmeasure=pomp.fun(),
+ dprior=pomp.fun(),
+ rprior=pomp.fun(),
skeleton.type="map",
skeleton=pomp.fun(),
skelmap.delta.t=as.numeric(NA),
Modified: pkg/pomp/R/pomp-methods.R
===================================================================
--- pkg/pomp/R/pomp-methods.R 2014-03-11 13:23:39 UTC (rev 890)
+++ pkg/pomp/R/pomp-methods.R 2014-03-16 21:57:33 UTC (rev 891)
@@ -47,33 +47,19 @@
partrans.internal(object=object,params=params,dir=dir,...)
)
-## a simple method to extract the data array
-setMethod(
- "data.array",
- "pomp",
- function (object, vars, ...) {
- varnames <- rownames(object at data)
- if (missing(vars))
- vars <- varnames
- else if (!all(vars%in%varnames))
- stop("some elements of ",sQuote("vars")," correspond to no observed variable")
- object at data[vars,,drop=FALSE]
- }
- )
+obs.internal <- function (object, vars, ...) {
+ varnames <- rownames(object at data)
+ if (missing(vars))
+ vars <- varnames
+ else if (!all(vars%in%varnames))
+ stop("some elements of ",sQuote("vars")," correspond to no observed variable")
+ object at data[vars,,drop=FALSE]
+}
+
## a simple method to extract the data array
-setMethod(
- "obs",
- "pomp",
- function (object, vars, ...) {
- varnames <- rownames(object at data)
- if (missing(vars))
- vars <- varnames
- else if (!all(vars%in%varnames))
- stop("some elements of ",sQuote("vars")," correspond to no observed variable")
- object at data[vars,,drop=FALSE]
- }
- )
+setMethod("obs","pomp",obs.internal)
+setMethod("data.array","pomp",obs.internal)
## a simple method to extract the array of states
setMethod(
@@ -293,6 +279,10 @@
show(object at rmeasure)
cat("measurement model density, dmeasure = \n")
show(object at dmeasure)
+ cat("prior simulator, rprior = \n")
+ show(object at rprior)
+ cat("prior density, dprior = \n")
+ show(object at dprior)
if (!is.na(object at skeleton.type)) {
cat("skeleton (",object at skeleton.type,") = \n")
show(object at skeleton)
Modified: pkg/pomp/R/pomp.R
===================================================================
--- pkg/pomp/R/pomp.R 2014-03-11 13:23:39 UTC (rev 890)
+++ pkg/pomp/R/pomp.R 2014-03-16 21:57:33 UTC (rev 891)
@@ -1,9 +1,12 @@
## basic constructor of the pomp class
pomp.constructor <- function (data, times, t0, ..., rprocess, dprocess,
rmeasure, dmeasure, measurement.model,
- skeleton = NULL, skeleton.type = c("map","vectorfield"),
+ skeleton = NULL,
+ skeleton.type = c("map","vectorfield"),
skelmap.delta.t = 1,
- initializer, params, covar, tcovar,
+ initializer,
+ rprior, dprior,
+ params, covar, tcovar,
obsnames, statenames, paramnames, covarnames, zeronames,
PACKAGE, parameter.transform, parameter.inv.transform) {
@@ -36,14 +39,16 @@
if (!is.numeric(times) || !all(diff(times)>0))
stop("pomp error: ",sQuote("times")," must be an increasing numeric vector",call.=TRUE)
if (length(times)!=ncol(data))
- stop("pomp error: the length of ",sQuote("times")," does not equal the number of columns in ",sQuote("data"),call.=TRUE)
+ stop("pomp error: the length of ",sQuote("times")," does not equal the number of columns in ",
+ sQuote("data"),call.=TRUE)
storage.mode(times) <- 'double'
## check t0
if (!is.numeric(t0) || length(t0) > 1)
stop("pomp error: the zero-time ",sQuote("t0")," must be a single number",call.=TRUE)
if (t0 > times[1L])
- stop("pomp error: the zero-time ",sQuote("t0")," must occur no later than the first observation",call.=TRUE)
+ stop("pomp error: the zero-time ",sQuote("t0")," must occur no later than the first observation",
+ call.=TRUE)
storage.mode(t0) <- 'double'
if (missing(PACKAGE)) PACKAGE <- ""
@@ -57,7 +62,8 @@
if (!(missing(dmeasure)&&missing(rmeasure))) {
warning(
"specifying ",sQuote("measurement.model"),
- " overrides specification of ",sQuote("rmeasure")," and ",sQuote("dmeasure")
+ " overrides specification of ",
+ sQuote("rmeasure")," and ",sQuote("dmeasure")
)
}
mm <- measform2pomp(measurement.model)
@@ -71,9 +77,10 @@
stop(sQuote("skelmap.delta.t")," must be positive")
if (is.null(skeleton)) {
- skeleton <- pomp.fun(f=function(x,t,params,covars,...)stop(sQuote("skeleton")," not specified"))
+ skeleton <- pomp.fun()
} else {
- skeleton <- pomp.fun(f=skeleton,PACKAGE=PACKAGE,proto=quote(skeleton(x,t,params,...)))
+ skeleton <- pomp.fun(f=skeleton,PACKAGE=PACKAGE,
+ proto=quote(skeleton(x,t,params,...)))
}
if (missing(initializer)) {
@@ -94,13 +101,29 @@
if (missing(rmeasure))
rmeasure <- pomp.fun()
else
- rmeasure <- pomp.fun(f=rmeasure,PACKAGE=PACKAGE,proto=quote(rmeasure(x,t,params,...)))
+ rmeasure <- pomp.fun(f=rmeasure,PACKAGE=PACKAGE,
+ proto=quote(rmeasure(x,t,params,...)))
if (missing(dmeasure))
dmeasure <- pomp.fun()
else
- dmeasure <- pomp.fun(f=dmeasure,PACKAGE=PACKAGE,proto=quote(dmeasure(y,x,t,params,log,...)))
+ dmeasure <- pomp.fun(f=dmeasure,PACKAGE=PACKAGE,
+ proto=quote(dmeasure(y,x,t,params,log,...)))
+ if (missing(rprior)) { ## no default rprior
+ rprior <- pomp.fun()
+ } else {
+ rprior <- pomp.fun(f=rprior,PACKAGE=PACKAGE,
+ proto=quote(rprior(params,...)))
+ }
+
+ if (missing(dprior)) { ## by default, use flat improper prior
+ dprior <- pomp.fun(f="_pomp_default_dprior")
+ } else {
+ dprior <- pomp.fun(f=dprior,PACKAGE=PACKAGE,
+ proto=quote(dprior(params,log,...)))
+ }
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 891
More information about the pomp-commits
mailing list