[Pomp-commits] r1088 - in pkg/pomp: R inst man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Feb 22 20:43:32 CET 2015
Author: kingaa
Date: 2015-02-22 20:43:32 +0100 (Sun, 22 Feb 2015)
New Revision: 1088
Modified:
pkg/pomp/R/abc.R
pkg/pomp/inst/NEWS
pkg/pomp/inst/NEWS.Rd
pkg/pomp/man/abc.Rd
pkg/pomp/tests/ou2-abc.R
pkg/pomp/tests/ou2-abc.Rout.save
Log:
- new 'proposal' argument to 'abc'
- 'rw.sd' argument is obviated and deprecated
- 'pars' argument is deprecated and ignored
Modified: pkg/pomp/R/abc.R
===================================================================
--- pkg/pomp/R/abc.R 2015-02-22 19:43:23 UTC (rev 1087)
+++ pkg/pomp/R/abc.R 2015-02-22 19:43:32 UTC (rev 1088)
@@ -8,14 +8,22 @@
probes='list',
scale = 'numeric',
epsilon = 'numeric',
- random.walk.sd = 'numeric',
+ proposal = 'function',
conv.rec = 'matrix'
+ ),
+ prototype=prototype(
+ pars = character(0),
+ Nabc = 0L,
+ probes = list(),
+ scale = numeric(0),
+ epsilon = 1.0,
+ proposal = function (...) stop("proposal not specified"),
+ conv.rec=array(dim=c(0,0))
)
)
abc.internal <- function (object, Nabc,
- start, pars,
- rw.sd, probes,
+ start, proposal, probes,
epsilon, scale,
verbose,
.ndone = 0L,
@@ -39,46 +47,16 @@
if (is.null(start.names))
stop("abc error: ",sQuote("start")," must be a named vector",call.=FALSE)
- if (missing(rw.sd))
- stop("abc error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
- rw.names <- names(rw.sd)
- if (is.null(rw.names) || any(rw.sd<0))
- stop("abc error: ",sQuote("rw.sd")," must be a named non-negative numerical vector",call.=FALSE)
- if (!all(rw.names%in%start.names))
- stop("abc 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("abc 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 (
- !is.character(pars) ||
- !all(pars%in%start.names) ||
- !all(pars%in%rw.names)
- )
- stop(
- "abc error: ",
- sQuote("pars"),
- " must be a 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("abc error: error in proposal function",call.=FALSE)
+ if (is.null(names(theta)) || !is.numeric(theta))
+ stop("abc error: ",sQuote("proposal")," must return a named numeric vector",call.=FALSE)
- if (!all(rw.names%in%pars)) {
- extra.rws <- rw.names[!(rw.names%in%pars)]
- warning(
- "abc 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)
-
if (!is.list(probes)) probes <- list(probes)
if (!all(sapply(probes,is.function)))
stop(sQuote("probes")," must be a function or a list of functions")
@@ -88,10 +66,7 @@
ntimes <- length(time(object))
if (verbose) {
- cat("performing",Nabc,"ABC iteration(s) to estimate parameter(s)",
- paste(pars,collapse=", "))
- cat(" using random-walk with SD\n")
- print(rw.sd)
+ cat("performing",Nabc,"ABC iteration(s)\n")
}
theta <- start
@@ -113,18 +88,6 @@
)
)
- if (!all(is.finite(theta[pars]))) {
- stop(
- sQuote("abc"),
- " error: cannot estimate non-finite parameters: ",
- paste(
- pars[!is.finite(theta[pars])],
- collapse=","
- ),
- call.=FALSE
- )
- }
-
## apply probes to data
datval <- try(.Call(apply_probe_data,object,probes),silent=FALSE)
if (inherits(datval,'try-error'))
@@ -134,8 +97,7 @@
for (n in seq_len(Nabc)) { # main loop
- theta.prop <- theta
- theta.prop[pars] <- rnorm(n=length(pars),mean=theta[pars],sd=rw.sd)
+ theta.prop <- proposal(theta)
log.prior.prop <- dprior(object,params=theta.prop,log=TRUE)
if (is.finite(log.prior.prop) &&
@@ -175,6 +137,9 @@
}
+ pars <- apply(conv.rec,2,function(x)diff(range(x))>0)
+ pars <- names(pars[pars])
+
new(
'abc',
object,
@@ -184,7 +149,7 @@
probes=probes,
scale=scale,
epsilon=epsilon,
- random.walk.sd=rw.sd,
+ proposal=proposal,
conv.rec=conv.rec
)
@@ -194,7 +159,7 @@
"abc",
signature=signature(object="pomp"),
function (object, Nabc = 1,
- start, pars, rw.sd,
+ start, proposal, pars, rw.sd,
probes, scale, epsilon,
verbose = getOption("verbose"),
...) {
@@ -202,13 +167,25 @@
if (missing(start))
start <- coef(object)
- if (missing(rw.sd))
- stop("abc error: ",sQuote("rw.sd")," must be specified",
- call.=FALSE)
+ if (missing(proposal)) proposal <- NULL
- if (missing(pars))
- pars <- names(rw.sd)[rw.sd>0]
+ if (!missing(rw.sd)) {
+ warning("abc warning: ",sQuote("rw.sd")," is a deprecated argument: ",
+ "Use ",sQuote("proposal")," instead.",call.=FALSE)
+ if (is.null(proposal)) {
+ proposal <- mvn.diag.rw(rw.sd=rw.sd)
+ } else {
+ warning("abc warning: since ",sQuote("proposal"),
+ " has been specified, ",sQuote("rw.sd")," is ignored.")
+ }
+ }
+ if (is.null(proposal))
+ stop("abc error: ",sQuote("proposal")," must be specified",call.=FALSE)
+
+ if (!missing(pars))
+ warning("abc warning: ",sQuote("pars")," is a deprecated argument and will be ignored.",call.=FALSE)
+
if (missing(probes))
stop("abc error: ",sQuote("probes")," must be specified",
call.=FALSE)
@@ -225,11 +202,10 @@
object=object,
Nabc=Nabc,
start=start,
- pars=pars,
+ proposal=proposal,
probes=probes,
scale=scale,
epsilon=epsilon,
- rw.sd=rw.sd,
verbose=verbose
)
}
@@ -256,15 +232,14 @@
"abc",
signature=signature(object="abc"),
function (object, Nabc,
- start, pars, rw.sd,
+ start, proposal,
probes, scale, epsilon,
verbose = getOption("verbose"),
...) {
if (missing(Nabc)) Nabc <- object at Nabc
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(probes)) probes <- object at probes
if (missing(scale)) scale <- object at scale
if (missing(epsilon)) epsilon <- object at epsilon
@@ -275,8 +250,7 @@
object=object,
Nabc=Nabc,
start=start,
- pars=pars,
- rw.sd=rw.sd,
+ proposal=proposal,
probes=probes,
scale=scale,
epsilon=epsilon,
Modified: pkg/pomp/inst/NEWS
===================================================================
--- pkg/pomp/inst/NEWS 2015-02-22 19:43:23 UTC (rev 1087)
+++ pkg/pomp/inst/NEWS 2015-02-22 19:43:32 UTC (rev 1088)
@@ -1,5 +1,28 @@
_N_e_w_s _f_o_r _p_a_c_k_a_g_e '_p_o_m_p'
+_C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _0._6_0-_1:
+
+ • ‘pmcmc’ and ‘abc’ can now use arbitrary symmetric proposal
+ distributions via the ‘proposal’ argument. For the moment,
+ these are constrained to be symmetric. Two new functions,
+ ‘mvn.diag.rw’ and ‘mvn.rw’ generate suitable proposal
+ functions. The first generates a multivariate normal
+ random-walk proposal with diagonal variance-covariance
+ matrix; this duplicates the old behavior of both ‘abc’ and
+ ‘pmcmc’. The second, ‘mvn.rw’, corresponds to a multivariate
+ normal random-walk proposal with arbitrary
+ variance-covariance matrix.
+
+ • In ‘pmcmc’ and ‘abc’, the arguments ‘pars’ and ‘rw.sd’ are
+ now unneeded and have been deprecated. Use of ‘rw.sd’ will
+ generate a warning and result in behavior equivalent to
+ choosing ‘proposal=mvn.diag.rw(rw.sd)’. Use of ‘pars’ will
+ be ignored, with a warning.
+
+ • In ‘nlf’, the ‘transform.params’ argument is now deprecated;
+ use instead the ‘transform’ argument, as in the other
+ inference methods.
+
_C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _0._5_9-_1:
• A bug in ‘spect.match’ has been fixed. Thanks to Karsten
Modified: pkg/pomp/inst/NEWS.Rd
===================================================================
--- pkg/pomp/inst/NEWS.Rd 2015-02-22 19:43:23 UTC (rev 1087)
+++ pkg/pomp/inst/NEWS.Rd 2015-02-22 19:43:32 UTC (rev 1088)
@@ -1,5 +1,19 @@
\name{NEWS}
\title{News for package `pomp'}
+\section{Changes in \pkg{pomp} version 0.60-1}{
+ \itemize{
+ \item \code{pmcmc} and \code{abc} can now use arbitrary symmetric proposal distributions via the \code{proposal} argument.
+ For the moment, these are constrained to be symmetric.
+ Two new functions, \code{mvn.diag.rw} and \code{mvn.rw} generate suitable proposal functions.
+ The first generates a multivariate normal random-walk proposal with diagonal variance-covariance matrix; this duplicates the old behavior of both \code{abc} and \code{pmcmc}.
+ The second, \code{mvn.rw}, corresponds to a multivariate normal random-walk proposal with arbitrary variance-covariance matrix.
+ \item In \code{pmcmc} and \code{abc}, the arguments \code{pars} and \code{rw.sd} are now unneeded and have been deprecated.
+ Use of \code{rw.sd} will generate a warning and result in behavior equivalent to choosing \code{proposal=mvn.diag.rw(rw.sd)}.
+ Use of \code{pars} will be ignored, with a warning.
+ \item In \code{nlf}, the \code{transform.params} argument is now deprecated;
+ use instead the \code{transform} argument, as in the other inference methods.
+ }
+}
\section{Changes in \pkg{pomp} version 0.59-1}{
\itemize{
\item A bug in \code{spect.match} has been fixed.
Modified: pkg/pomp/man/abc.Rd
===================================================================
--- pkg/pomp/man/abc.Rd 2015-02-22 19:43:23 UTC (rev 1087)
+++ pkg/pomp/man/abc.Rd 2015-02-22 19:43:32 UTC (rev 1088)
@@ -16,13 +16,13 @@
The approximate Bayesian computation (ABC) algorithm for estimating the parameters of a partially-observed Markov process.
}
\usage{
-\S4method{abc}{pomp}(object, Nabc = 1, start, pars,
- rw.sd, probes, scale, epsilon,
+\S4method{abc}{pomp}(object, Nabc = 1, start, proposal,
+ pars, rw.sd, probes, scale, epsilon,
verbose = getOption("verbose"), \dots)
\S4method{abc}{probed.pomp}(object, probes,
verbose = getOption("verbose"), \dots)
-\S4method{abc}{abc}(object, Nabc, start, pars,
- rw.sd, probes, scale, epsilon,
+\S4method{abc}{abc}(object, Nabc, start, proposal,
+ probes, scale, epsilon,
verbose = getOption("verbose"), \dots)
\S4method{continue}{abc}(object, Nabc = 1, \dots)
}
@@ -37,27 +37,27 @@
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{pars}{Deprecated and ignored. Will be removed in a future release.}
\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}.
+ Deprecated. Will be removed in a future release.
+ Specifying \code{rw.sd} is equivalent to setting \code{proposal=mvn.diag.rw(rw.sd)}.
}
\item{probes}{
+ List of probes (AKA summary statistics).
+ See \code{\link{probe}} for details.
}
\item{scale}{
-
+ named numeric vector of scales.
}
\item{epsilon}{
-
+ ABC tolerance.
}
\item{verbose}{
logical; if TRUE, print progress reports.
@@ -74,11 +74,6 @@
\item{pars, Nabc, dprior, hyperparams, scale, epsilon}{
These slots hold the values of the corresponding arguments of the call to \code{abc}.
}
- \item{random.walk.sd}{
- a named numeric vector containing the random-walk variances used to parameterize a Gaussian random walk MCMC proposal.
- }
- \item{probes, conv.rec}{
- }
}
}
\section{Re-running ABC Iterations}{
Modified: pkg/pomp/tests/ou2-abc.R
===================================================================
--- pkg/pomp/tests/ou2-abc.R 2015-02-22 19:43:23 UTC (rev 1087)
+++ pkg/pomp/tests/ou2-abc.R 2015-02-22 19:43:32 UTC (rev 1088)
@@ -25,11 +25,10 @@
abc1 <- abc(po,
Nabc=10000,
start=coef(ou2),
- pars=c("alpha.1","alpha.2"),
probes=probes.good,
scale=scale.dat,
epsilon=1.7,
- rw.sd= c(alpha.1=0.01,alpha.2=0.01)
+ proposal=mvn.diag.rw(rw.sd=c(alpha.1=0.01,alpha.2=0.01))
)
plot(abc1,scatter=TRUE)
@@ -42,11 +41,10 @@
abc2 <- abc(po,
Nabc=2000,
- pars=c("alpha.1","alpha.2"),
probes=probes.good,
scale=scale.dat,
epsilon=1,
- rw.sd= c(alpha.1=0.01,alpha.2=0.01)
+ proposal=mvn.diag.rw(c(alpha.1=0.01,alpha.2=0.01))
)
plot(abc2)
@@ -55,16 +53,21 @@
probes=probes.good,
scale=scale.dat,
epsilon=2,
- rw.sd= c(alpha.1=0.01,alpha.2=0.01)
+ rw.sd=c(alpha.1=0.01,alpha.2=0.01)
)
abc3 <- continue(abc3,Nabc=3000)
plot(abc3)
+sig <- matrix(c(0.01,0.005,0.005,0.01),
+ 2,2,
+ dimnames=list(c("alpha.1","alpha.2"),
+ c("alpha.1","alpha.2")))
+
abc4 <- abc(probe(po,probes=probes.good,nsim=200),
Nabc=2000,
scale=scale.dat,
epsilon=2,
- rw.sd= c(alpha.1=0.01,alpha.2=0.01)
+ proposal=mvn.rw(sig)
)
plot(abc4)
@@ -89,7 +92,7 @@
probes=probes.good,
scale=scale.dat,
epsilon=1,
- rw.sd= c(alpha.1=0.01,alpha.2=0.01)
+ proposal=mvn.diag.rw(c(alpha.1=0.01,alpha.2=0.01))
)
plot(abc6)
Modified: pkg/pomp/tests/ou2-abc.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-abc.Rout.save 2015-02-22 19:43:23 UTC (rev 1087)
+++ pkg/pomp/tests/ou2-abc.Rout.save 2015-02-22 19:43:32 UTC (rev 1088)
@@ -46,11 +46,10 @@
> abc1 <- abc(po,
+ Nabc=10000,
+ start=coef(ou2),
-+ pars=c("alpha.1","alpha.2"),
+ probes=probes.good,
+ scale=scale.dat,
+ epsilon=1.7,
-+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
++ proposal=mvn.diag.rw(rw.sd=c(alpha.1=0.01,alpha.2=0.01))
+ )
>
> plot(abc1,scatter=TRUE)
@@ -60,15 +59,14 @@
> runs <- rle(as.vector(conv.rec(abc1)[, "alpha.1"]))
> hist(runs$lengths)
> mean(runs$length)
-[1] 6.949965
+[1] 13.81354
>
> abc2 <- abc(po,
+ Nabc=2000,
-+ pars=c("alpha.1","alpha.2"),
+ probes=probes.good,
+ scale=scale.dat,
+ epsilon=1,
-+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
++ proposal=mvn.diag.rw(c(alpha.1=0.01,alpha.2=0.01))
+ )
> plot(abc2)
>
@@ -77,16 +75,23 @@
+ probes=probes.good,
+ scale=scale.dat,
+ epsilon=2,
-+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
++ rw.sd=c(alpha.1=0.01,alpha.2=0.01)
+ )
+Warning message:
+abc warning: 'rw.sd' is a deprecated argument: Use 'proposal' instead.
> abc3 <- continue(abc3,Nabc=3000)
> plot(abc3)
>
+> sig <- matrix(c(0.01,0.005,0.005,0.01),
++ 2,2,
++ dimnames=list(c("alpha.1","alpha.2"),
++ c("alpha.1","alpha.2")))
+>
> abc4 <- abc(probe(po,probes=probes.good,nsim=200),
+ Nabc=2000,
+ scale=scale.dat,
+ epsilon=2,
-+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
++ proposal=mvn.rw(sig)
+ )
> plot(abc4)
>
@@ -111,8 +116,10 @@
+ probes=probes.good,
+ scale=scale.dat,
+ epsilon=1,
-+ rw.sd= c(alpha.1=0.01,alpha.2=0.01)
++ proposal=mvn.diag.rw(c(alpha.1=0.01,alpha.2=0.01))
+ )
+Warning message:
+abc warning: 'pars' is a deprecated argument and will be ignored.
> plot(abc6)
>
> try(abc7 <- c(abc2,abc3))
@@ -131,4 +138,4 @@
>
> proc.time()
user system elapsed
- 9.612 0.100 9.906
+ 12.322 0.033 12.359
More information about the pomp-commits
mailing list