[Pomp-commits] r1086 - in pkg/pomp: R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Feb 22 20:43:18 CET 2015
Author: kingaa
Date: 2015-02-22 20:43:18 +0100 (Sun, 22 Feb 2015)
New Revision: 1086
Modified:
pkg/pomp/R/pmcmc-methods.R
pkg/pomp/R/pmcmc.R
pkg/pomp/man/pmcmc.Rd
pkg/pomp/tests/ou2-pmcmc.R
pkg/pomp/tests/ou2-pmcmc.Rout.save
Log:
- new 'proposal' argument to 'pmcmc'
- 'rw.sd' argument is obviated and deprecated
- 'pars' argument is deprecated and ignored
Modified: pkg/pomp/R/pmcmc-methods.R
===================================================================
--- pkg/pomp/R/pmcmc-methods.R 2015-02-22 19:43:12 UTC (rev 1085)
+++ pkg/pomp/R/pmcmc-methods.R 2015-02-22 19:43:18 UTC (rev 1086)
@@ -127,7 +127,8 @@
mar.multi <- c(0,5.1,0,2.1)
oma.multi <- c(6,0,5,0)
xx <- z[[1]]
- estnames <- xx at pars
+ estnames <- apply(xx at conv.rec,2,function(x)diff(range(x))>0)
+ estnames <- names(estnames[estnames])
## plot pmcmc convergence diagnostics
other.diagnostics <- c("loglik", "log.prior","nfail")
Modified: pkg/pomp/R/pmcmc.R
===================================================================
--- pkg/pomp/R/pmcmc.R 2015-02-22 19:43:12 UTC (rev 1085)
+++ pkg/pomp/R/pmcmc.R 2015-02-22 19:43:18 UTC (rev 1086)
@@ -3,18 +3,22 @@
'pmcmc',
contains='pfilterd.pomp',
slots=c(
- pars = 'character',
Nmcmc = 'integer',
- random.walk.sd = 'numeric',
- conv.rec = 'matrix',
+ proposal = 'function',
+ conv.rec = 'array',
log.prior = 'numeric'
+ ),
+ prototype=prototype(
+ Nmcmc = 0L,
+ proposal = function (...) stop("proposal not specified"),
+ conv.rec=array(dim=c(0,0)),
+ log.prior=numeric(0)
)
)
pmcmc.internal <- function (object, Nmcmc,
- start, pars,
- rw.sd, Np,
- tol, max.fail,
+ start, proposal,
+ Np, tol, max.fail,
verbose,
.ndone = 0L,
.prev.pfp = NULL, .prev.log.prior = NULL,
@@ -29,57 +33,21 @@
if (length(start)==0)
stop(
sQuote("start")," must be specified if ",
- sQuote("coef(object)")," is NULL",
- call.=FALSE
+ sQuote("coef(object)")," is NULL"
)
- start.names <- names(start)
- if (is.null(start.names))
+ if (is.null(names(start)))
stop("pmcmc error: ",sQuote("start")," must be a named vector",call.=FALSE)
- if (missing(rw.sd))
- stop("pmcmc error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
- rw.names <- names(rw.sd)
- if (is.null(rw.names) || any(rw.sd<0))
- stop("pmcmc error: ",sQuote("rw.sd")," must be a named non-negative numerical vector",call.=FALSE)
- if (!all(rw.names%in%start.names))
- stop("pmcmc error: all the names of ",sQuote("rw.sd")," must be names of ",sQuote("start"),call.=FALSE)
- rw.names <- names(rw.sd[rw.sd>0])
- if (length(rw.names) == 0)
- stop("pmcmc error: ",sQuote("rw.sd")," must have one positive entry for each parameter to be estimated",call.=FALSE)
+ if (!is.function(proposal))
+ stop(sQuote("proposal")," must be a function")
- if (missing(pars))
- stop("pmcmc error: ",sQuote("pars")," must be specified",call.=FALSE)
- if (length(pars)==0)
- stop("pmcmc error: at least one parameter must be estimated",call.=FALSE)
- if (
- !is.character(pars) ||
- !all(pars%in%start.names) ||
- !all(pars%in%rw.names)
- )
- stop(
- "pmcmc error: ",
- sQuote("pars"),
- " must be a mutually disjoint subset of ",
- sQuote("names(start)"),
- " and must correspond to positive random-walk SDs specified in ",
- sQuote("rw.sd"),
- call.=FALSE
- )
+ ## test proposal distribution
+ theta <- try(proposal(start))
+ if (inherits(theta,"try-error"))
+ stop("pmcmc error: error in proposal function")
+ if (is.null(names(theta)) || !is.numeric(theta))
+ stop("pmcmc error: ",sQuote("proposal")," must return a named numeric vector")
- if (!all(rw.names%in%pars)) {
- extra.rws <- rw.names[!(rw.names%in%pars)]
- warning(
- "pmcmc warning: the variable(s) ",
- paste(extra.rws,collapse=", "),
- " have positive random-walk SDs specified, but are not included in ",
- sQuote("pars"),
- ". These random walk SDs are ignored.",
- call.=FALSE
- )
- }
- rw.sd <- rw.sd[pars]
- rw.names <- names(rw.sd)
-
ntimes <- length(time(object))
if (missing(Np))
stop("pmcmc error: ",sQuote("Np")," must be specified",call.=FALSE)
@@ -108,11 +76,7 @@
stop("pmcmc error: ",sQuote("Nmcmc")," must be a positive integer",call.=FALSE)
if (verbose) {
- cat("performing",Nmcmc,"PMCMC iteration(s) to estimate parameter(s)",
- paste(pars,collapse=", "))
- cat(" using random-walk with SD\n")
- print(rw.sd)
- cat("using",Np,"particles\n")
+ cat("performing",Nmcmc,"PMCMC iteration(s) using",Np,"particles\n")
}
theta <- start
@@ -127,18 +91,6 @@
)
)
- if (!all(is.finite(theta[pars]))) {
- stop(
- sQuote("pmcmc"),
- " error: cannot estimate non-finite parameters: ",
- paste(
- pars[!is.finite(theta[pars])],
- collapse=","
- ),
- call.=FALSE
- )
- }
-
if (.ndone==0L) { ## compute prior and likelihood on initial parameter vector
pfp <- try(
pfilter.internal(
@@ -171,8 +123,7 @@
for (n in seq_len(Nmcmc)) { # main loop
- theta.prop <- theta
- theta.prop[pars] <- rnorm(n=length(pars),mean=theta[pars],sd=rw.sd)
+ theta.prop <- proposal(theta)
## run the particle filter on the proposed new parameter values
pfp.prop <- try(
@@ -218,8 +169,7 @@
pfp,
params=theta,
Nmcmc=Nmcmc,
- pars=pars,
- random.walk.sd=rw.sd,
+ proposal=proposal,
Np=Np,
tol=tol,
conv.rec=conv.rec,
@@ -231,24 +181,34 @@
"pmcmc",
signature=signature(object="pomp"),
function (object, Nmcmc = 1,
- start, pars, rw.sd, Np,
+ start, proposal, pars, rw.sd, Np,
tol = 1e-17, max.fail = 0,
verbose = getOption("verbose"),
...) {
if (missing(start)) start <- coef(object)
- 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(rw.sd)) {
+ warning("pmcmc warning: ",sQuote("rw.sd")," is a deprecated argument.",
+ "Use ",sQuote("proposal")," instead.",call.=FALSE)
+ if (missing(proposal)) {
+ proposal <- mvn.diag.rw(rw.sd=rw.sd)
+ } else {
+ warning("pmcmc warning: since ",sQuote("proposal"),
+ " has been specified, ",sQuote("rw.sd")," is ignored.")
+ }
+ }
+
+ if (!missing(pars))
+ warning("pmcmc warning: ",sQuote("pars")," is a deprecated argument and will be ignored.",call.=FALSE)
+
pmcmc.internal(
object=object,
Nmcmc=Nmcmc,
start=start,
- pars=pars,
- rw.sd=rw.sd,
+ proposal=proposal,
Np=Np,
tol=tol,
max.fail=max.fail,
@@ -280,15 +240,14 @@
"pmcmc",
signature=signature(object="pmcmc"),
function (object, Nmcmc,
- start, pars, rw.sd,
+ start, proposal,
Np, tol, max.fail = 0,
verbose = getOption("verbose"),
...) {
if (missing(Nmcmc)) Nmcmc <- object at Nmcmc
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(proposal)) proposal <- object at proposal
if (missing(Np)) Np <- object at Np
if (missing(tol)) tol <- object at tol
@@ -296,8 +255,7 @@
object=as(object,"pomp"),
Nmcmc=Nmcmc,
start=start,
- pars=pars,
- rw.sd=rw.sd,
+ proposal=proposal,
Np=Np,
tol=tol,
max.fail=max.fail,
Modified: pkg/pomp/man/pmcmc.Rd
===================================================================
--- pkg/pomp/man/pmcmc.Rd 2015-02-22 19:43:12 UTC (rev 1085)
+++ pkg/pomp/man/pmcmc.Rd 2015-02-22 19:43:18 UTC (rev 1086)
@@ -15,11 +15,11 @@
The Particle MCMC algorithm for estimating the parameters of a partially-observed Markov process.
}
\usage{
-\S4method{pmcmc}{pomp}(object, Nmcmc = 1, start, pars, rw.sd, Np,
+\S4method{pmcmc}{pomp}(object, Nmcmc = 1, start, proposal, pars, rw.sd, Np,
tol = 1e-17, max.fail = 0, verbose = getOption("verbose"), \dots)
\S4method{pmcmc}{pfilterd.pomp}(object, Nmcmc = 1, Np, tol, \dots)
-\S4method{pmcmc}{pmcmc}(object, Nmcmc, start, pars, rw.sd, Np,
- tol, max.fail = 0, verbose = getOption("verbose"), \dots)
+\S4method{pmcmc}{pmcmc}(object, Nmcmc, start, proposal, Np, tol,
+ max.fail = 0, verbose = getOption("verbose"), \dots)
\S4method{continue}{pmcmc}(object, Nmcmc = 1, \dots)
}
\arguments{
@@ -33,18 +33,14 @@
named numeric vector;
the starting guess of the parameters.
}
- \item{pars}{
- optional character vector naming the ordinary parameters to be estimated.
- Every parameter named in \code{pars} must have a positive random-walk standard deviation specified in \code{rw.sd}.
- Leaving \code{pars} unspecified is equivalent to setting it equal to the names of all parameters with a positive value of \code{rw.sd}.
+ \item{proposal}{
+ optional function that draws from the proposal distribution.
+ Currently, the proposal distribution must be symmetric for proper inference:
+ it is the user's responsibility to ensure that it is.
+ Several functions that construct appropriate proposal function are provided:
+ see \link{MCMC proposal functions} for more information.
}
- \item{rw.sd}{
- numeric vector with names; used to parameterize a Gaussian random walk MCMC proposal. The random walk is only applied to parameters named in \code{pars}. The algorithm requires that the random walk be nontrivial, so each element in \code{rw.sd[pars]} must be positive.
- The following must be satisfied:
- \code{names(rw.sd)} must be a subset of \code{names(start)},
- \code{rw.sd} must be non-negative (zeros are simply ignored),
- the name of every positive element of \code{rw.sd} must be in \code{pars}.
- }
+ \item{pars, rw.sd}{Deprecated. Will be removed in a future release.}
\item{Np}{
a positive integer;
the number of particles to use in each filtering operation.
Modified: pkg/pomp/tests/ou2-pmcmc.R
===================================================================
--- pkg/pomp/tests/ou2-pmcmc.R 2015-02-22 19:43:12 UTC (rev 1085)
+++ pkg/pomp/tests/ou2-pmcmc.R 2015-02-22 19:43:18 UTC (rev 1086)
@@ -12,7 +12,7 @@
f1 <- pmcmc(
pomp(ou2,dprior=dprior.ou2),
Nmcmc=20,
- rw.sd=c(alpha.2=0.001,alpha.3=0.001),
+ proposal=mvn.diag.rw(c(alpha.2=0.001,alpha.3=0.001)),
Np=100,
max.fail=100,
verbose=FALSE
@@ -55,7 +55,7 @@
}
),
Nmcmc=20,
- rw.sd=c(alpha.2=0.001,alpha.3=0.001),
+ proposal=mvn.diag.rw(c(alpha.2=0.001,alpha.3=0.001)),
Np=100,
max.fail=100,
verbose=FALSE
Modified: pkg/pomp/tests/ou2-pmcmc.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-pmcmc.Rout.save 2015-02-22 19:43:12 UTC (rev 1085)
+++ pkg/pomp/tests/ou2-pmcmc.Rout.save 2015-02-22 19:43:18 UTC (rev 1086)
@@ -33,7 +33,7 @@
> f1 <- pmcmc(
+ pomp(ou2,dprior=dprior.ou2),
+ Nmcmc=20,
-+ rw.sd=c(alpha.2=0.001,alpha.3=0.001),
++ proposal=mvn.diag.rw(c(alpha.2=0.001,alpha.3=0.001)),
+ Np=100,
+ max.fail=100,
+ verbose=FALSE
@@ -49,6 +49,8 @@
+ max.fail=100,
+ verbose=FALSE
+ )
+Warning message:
+pmcmc warning: 'rw.sd' is a deprecated argument.Use 'proposal' instead.
>
> f3 <- pmcmc(f2)
> f4 <- continue(f3,Nmcmc=20)
@@ -78,7 +80,7 @@
+ }
+ ),
+ Nmcmc=20,
-+ rw.sd=c(alpha.2=0.001,alpha.3=0.001),
++ proposal=mvn.diag.rw(c(alpha.2=0.001,alpha.3=0.001)),
+ Np=100,
+ max.fail=100,
+ verbose=FALSE
@@ -114,4 +116,4 @@
>
> proc.time()
user system elapsed
- 24.125 0.068 24.493
+ 28.001 0.048 28.080
More information about the pomp-commits
mailing list