[Pomp-commits] r608 - in pkg: . R inst man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 9 20:42:03 CET 2012
Author: kingaa
Date: 2012-02-09 20:42:03 +0100 (Thu, 09 Feb 2012)
New Revision: 608
Modified:
pkg/.Rinstignore
pkg/DESCRIPTION
pkg/R/mif.R
pkg/inst/NEWS
pkg/man/mif.Rd
pkg/src/dprocess.c
pkg/src/euler.c
pkg/src/rmeasure.c
pkg/src/skeleton.c
Log:
- more informative error messages in cases where dimensionality of state space or data space differ internally
- more flexible choice of update methods in 'mif' via the new argument 'method' which replaces the old 'weighted' flag
- the 'weighted' flag of 'mif' is now deprecated
Modified: pkg/.Rinstignore
===================================================================
--- pkg/.Rinstignore 2012-02-06 20:14:40 UTC (rev 607)
+++ pkg/.Rinstignore 2012-02-09 19:42:03 UTC (rev 608)
@@ -1 +1,2 @@
inst/doc/Makefile
+inst/doc/fullnat.bst
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-02-06 20:14:40 UTC (rev 607)
+++ pkg/DESCRIPTION 2012-02-09 19:42:03 UTC (rev 608)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.40-4
-Date: 2012-01-20
+Version: 0.40-5
+Date: 2012-02-10
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
Modified: pkg/R/mif.R
===================================================================
--- pkg/R/mif.R 2012-02-06 20:14:40 UTC (rev 607)
+++ pkg/R/mif.R 2012-02-09 19:42:03 UTC (rev 608)
@@ -38,7 +38,7 @@
particles,
rw.sd,
Np, cooling.factor, var.factor, ic.lag,
- weighted, tol, max.fail,
+ method, tol, max.fail,
verbose, .ndone) {
if (length(start)==0)
@@ -232,7 +232,7 @@
tol=tol,
max.fail=max.fail,
pred.mean=(n==Nmif),
- pred.var=(weighted||(n==Nmif)),
+ pred.var=((method=="mif")||(n==Nmif)),
filter.mean=TRUE,
save.states=FALSE,
save.params=FALSE,
@@ -244,15 +244,24 @@
if (inherits(pfp,'try-error'))
stop("mif error: error in ",sQuote("pfilter"),call.=FALSE)
- if (weighted) { # MIF update rule
- v <- pfp$pred.var[pars,,drop=FALSE] # the prediction variance
- v1 <- cool.sched$gamma*(1+var.factor^2)*sigma[pars]^2
- theta.hat <- cbind(theta[pars],pfp$filter.mean[pars,,drop=FALSE])
- theta[pars] <- theta[pars]+colSums(apply(theta.hat,1,diff)/t(v))*v1
- } else { # unweighted (flat) average
- theta.hat <- pfp$filter.mean[pars,,drop=FALSE]
- theta[pars] <- rowMeans(theta.hat)
- }
+ switch(
+ method,
+ mif={ # mif update rule
+ v <- pfp$pred.var[pars,,drop=FALSE] # the prediction variance
+ v1 <- cool.sched$gamma*(1+var.factor^2)*sigma[pars]^2
+ theta.hat <- cbind(theta[pars],pfp$filter.mean[pars,,drop=FALSE])
+ theta[pars] <- theta[pars]+colSums(apply(theta.hat,1,diff)/t(v))*v1
+ },
+ unweighted={ # unweighted (flat) average
+ theta.hat <- pfp$filter.mean[pars,,drop=FALSE]
+ theta[pars] <- rowMeans(theta.hat)
+ },
+ fp={
+ theta.hat <- pfp$filter.mean[pars,ntimes,drop=FALSE]
+ theta[pars] <- theta.hat
+ },
+ stop("unrecognized method",sQuote(method))
+ )
## update the IVPs using fixed-lag smoothing
theta[ivps] <- pfp$filter.mean[ivps,ic.lag]
@@ -292,7 +301,8 @@
pars, ivps = character(0),
particles, rw.sd,
Np, ic.lag, var.factor, cooling.factor,
- weighted = TRUE, tol = 1e-17, max.fail = 0,
+ weighted, method = c("mif","unweighted","fp"),
+ tol = 1e-17, max.fail = 0,
verbose = getOption("verbose"), ...) {
if (missing(start)) start <- coef(object)
@@ -311,6 +321,22 @@
if (missing(cooling.factor))
stop("mif error: ",sQuote("cooling.factor")," must be specified",call.=FALSE)
+ method <- match.arg(method)
+ if (!missing(weighted)) {
+ warning(sQuote("mif")," warning: ",sQuote("weighted")," flag is deprecated, use ",sQuote("method"))
+ if (weighted) {
+ if (method!="mif") {
+ warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
+ }
+ method <- "mif"
+ } else {
+ if (method!="unweighted") {
+ warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
+ }
+ method <- "unweighted"
+ }
+ }
+
if (missing(particles)) { # use default: normal distribution
particles <- function (Np, center, sd, ...) {
matrix(
@@ -351,7 +377,7 @@
cooling.factor=cooling.factor,
var.factor=var.factor,
ic.lag=ic.lag,
- weighted=weighted,
+ method=method,
tol=tol,
max.fail=max.fail,
verbose=verbose,
@@ -370,7 +396,8 @@
pars, ivps = character(0),
particles, rw.sd,
Np, ic.lag, var.factor, cooling.factor,
- weighted = TRUE, tol, max.fail = 0,
+ weighted, method = c("mif","unweighted","fp"),
+ tol = 1e-17, max.fail = 0,
verbose = getOption("verbose"), ...) {
if (missing(start)) start <- coef(object)
@@ -389,6 +416,22 @@
if (missing(cooling.factor))
stop("mif error: ",sQuote("cooling.factor")," must be specified",call.=FALSE)
+ method <- match.arg(method)
+ if (!missing(weighted)) {
+ warning(sQuote("mif")," warning: ",sQuote("weighted")," flag is deprecated, use ",sQuote("method"))
+ if (weighted) {
+ if (method!="mif") {
+ warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
+ }
+ method <- "mif"
+ } else {
+ if (method!="unweighted") {
+ warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
+ }
+ method <- "unweighted"
+ }
+ }
+
if (missing(particles)) { # use default: normal distribution
particles <- default.pomp.particles.fun
} else {
@@ -415,7 +458,7 @@
cooling.factor=cooling.factor,
var.factor=var.factor,
ic.lag=ic.lag,
- weighted=weighted,
+ method=method,
tol=tol,
max.fail=max.fail,
verbose=verbose,
@@ -432,7 +475,8 @@
pars, ivps,
particles, rw.sd,
Np, ic.lag, var.factor, cooling.factor,
- weighted = TRUE, tol, max.fail = 0,
+ weighted, method = c("mif","unweighted","fp"),
+ tol = 1e-17, max.fail = 0,
verbose = getOption("verbose"), ...) {
if (missing(Nmif)) Nmif <- object at Nmif
@@ -447,6 +491,22 @@
if (missing(cooling.factor)) cooling.factor <- object at cooling.factor
if (missing(tol)) tol <- object at tol
+ method <- match.arg(method)
+ if (!missing(weighted)) {
+ warning(sQuote("mif")," warning: ",sQuote("weighted")," flag is deprecated, use ",sQuote("method"))
+ if (weighted) {
+ if (method!="mif") {
+ warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
+ }
+ method <- "mif"
+ } else {
+ if (method!="unweighted") {
+ warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
+ }
+ method <- "unweighted"
+ }
+ }
+
mif.internal(
object=as(object,"pomp"),
Nmif=Nmif,
@@ -459,7 +519,7 @@
cooling.factor=cooling.factor,
var.factor=var.factor,
ic.lag=ic.lag,
- weighted=weighted,
+ method=method,
tol=tol,
max.fail=max.fail,
verbose=verbose,
@@ -476,7 +536,8 @@
pars, ivps,
particles, rw.sd,
Np, ic.lag, var.factor, cooling.factor,
- weighted = TRUE, tol, max.fail = 0,
+ weighted, method = c("mif","unweighted","fp"),
+ tol = 1e-17, max.fail = 0,
verbose = getOption("verbose"), ...) {
ndone <- object at Nmif
@@ -491,6 +552,22 @@
if (missing(cooling.factor)) cooling.factor <- object at cooling.factor
if (missing(tol)) tol <- object at tol
+ method <- match.arg(method)
+ if (!missing(weighted)) {
+ warning(sQuote("mif")," warning: ",sQuote("weighted")," flag is deprecated, use ",sQuote("method"))
+ if (weighted) {
+ if (method!="mif") {
+ warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
+ }
+ method <- "mif"
+ } else {
+ if (method!="unweighted") {
+ warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
+ }
+ method <- "unweighted"
+ }
+ }
+
obj <- mif.internal(
object=as(object,"pomp"),
Nmif=Nmif,
@@ -503,7 +580,7 @@
cooling.factor=cooling.factor,
var.factor=var.factor,
ic.lag=ic.lag,
- weighted=weighted,
+ method=method,
tol=tol,
max.fail=max.fail,
verbose=verbose,
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2012-02-06 20:14:40 UTC (rev 607)
+++ pkg/inst/NEWS 2012-02-09 19:42:03 UTC (rev 608)
@@ -1,4 +1,9 @@
NEWS
+0.40-5
+ o More informative error messages when dimension of state space or data space disagree internally.
+ o The 'weighted' argument to 'mif' and 'continue' is now deprecated in favor of a new argument 'method'.
+ The old 'weighted=T' corresponds to 'method="mif"' while the old 'weighted=F' corresponds to 'method="unweighted"'.
+
0.40-4
o New functions 'traj.match.objfun' and 'probe.match.objfun' have been added.
These functions construct functions of one argument suitable for use as objective functions in 'optim'-like optimizers (iincluding 'subplex' and 'sannbox').
Modified: pkg/man/mif.Rd
===================================================================
--- pkg/man/mif.Rd 2012-02-06 20:14:40 UTC (rev 607)
+++ pkg/man/mif.Rd 2012-02-09 19:42:03 UTC (rev 608)
@@ -18,19 +18,23 @@
mif(object, \dots)
\S4method{mif}{pomp}(object, Nmif = 1, start, pars, ivps = character(0),
particles, rw.sd, Np, ic.lag, var.factor, cooling.factor,
- weighted = TRUE, tol = 1e-17, max.fail = 0,
+ weighted, method = c("mif","unweighted","fp"),
+ tol = 1e-17, max.fail = 0,
verbose = getOption("verbose"), \dots)
\S4method{mif}{pfilterd.pomp}(object, Nmif = 1, start, pars, ivps = character(0),
particles, rw.sd, Np, ic.lag, var.factor, cooling.factor,
- weighted = TRUE, tol, max.fail = 0,
+ weighted, method = c("mif","unweighted","fp"),
+ tol, max.fail = 0,
verbose = getOption("verbose"), \dots)
\S4method{mif}{mif}(object, Nmif, start, pars, ivps,
particles, rw.sd, Np, ic.lag, var.factor, cooling.factor,
- weighted = TRUE, tol, max.fail = 0,
+ weighted, method = c("mif","unweighted","fp"),
+ tol, max.fail = 0,
verbose = getOption("verbose"), \dots)
\S4method{continue}{mif}(object, Nmif = 1, start, pars, ivps,
particles, rw.sd, Np, ic.lag, var.factor, cooling.factor,
- weighted = TRUE, tol, max.fail = 0,
+ weighted, method = c("mif","unweighted","fp"),
+ tol, max.fail = 0,
verbose = getOption("verbose"), \dots)
}
\arguments{
@@ -95,10 +99,15 @@
a positive number not greater than 1;
the exponential cooling factor, \code{alpha}.
}
- \item{weighted}{
- logical; if TRUE, the MIF update (a weighted average) is used.
- If FALSE, the MIF update is not used;
- instead, an unweighed average of the filtering means is used for the update.
+ \item{method, weighted}{
+ \code{method} sets the update rule used in the algorithm.
+ \code{method="mif"} uses the iterated filtering update rule (Ionides 2006, 2011);
+ \code{method="unweighted"} updates the parameter to the unweighted average of the filtering means of the parameters at each time;
+ \code{method="fp"} updates the parameter to the filtering mean at the end of the time series.
+ \code{weighted} is logical.
+ This argument is deprecated and will be removed in a future release.
+ \code{weighted=TRUE} is equivalent to \code{method="mif"};
+ \code{weighted=FALSE} is equivalent to \code{method="unweighted"}.
}
\item{tol}{
numeric scalar; particles with log likelihood below \code{tol} are considered to be \dQuote{lost}.
@@ -157,6 +166,10 @@
Inference for nonlinear dynamical systems,
Proc. Natl. Acad. Sci. U.S.A., 103:18438--18443, 2006.
+ E. L. Ionides, A. Bhadra, Y. Atchad{\\'e}, & A. A. King,
+ Iterated filtering,
+ Annals of Statistics, 39:1776--1802, 2011.
+
A. A. King, E. L. Ionides, M. Pascual, and M. J. Bouma,
Inapparent infections and cholera dynamics,
Nature, 454:877--880, 2008.
Modified: pkg/src/dprocess.c
===================================================================
--- pkg/src/dprocess.c 2012-02-06 20:14:40 UTC (rev 607)
+++ pkg/src/dprocess.c 2012-02-09 19:42:03 UTC (rev 608)
@@ -76,4 +76,3 @@
UNPROTECT(nprotect);
return X;
}
-
Modified: pkg/src/euler.c
===================================================================
--- pkg/src/euler.c 2012-02-06 20:14:40 UTC (rev 607)
+++ pkg/src/euler.c 2012-02-09 19:42:03 UTC (rev 608)
@@ -292,7 +292,8 @@
if (FIRST) {
if (LENGTH(ans) != NVAR) {
UNPROTECT(nprotect);
- error("user 'step.fun' must return a vector of length %d",NVAR);
+ error("user 'step.fun' returns a vector of %d states but %d are expected: compare initial conditions?",
+ LENGTH(ans),NVAR);
}
PROTECT(nm = GET_NAMES(ans)); nprotect++;
if (!isNull(nm)) { // match names against names from data slot
Modified: pkg/src/rmeasure.c
===================================================================
--- pkg/src/rmeasure.c 2012-02-06 20:14:40 UTC (rev 607)
+++ pkg/src/rmeasure.c 2012-02-09 19:42:03 UTC (rev 608)
@@ -101,8 +101,10 @@
PROTECT(ans = eval(FCALL,RHO)); nprotect++; // evaluate the call
if (FIRST) {
- if (LENGTH(ans) != NOBS)
- error("rmeasure error: user 'rmeasure' must return a vector of length %d",NOBS);
+ if (LENGTH(ans) != NOBS) {
+ error("user 'rmeasure' returns a vector of %d observables but %d are expected: compare 'data' slot?",
+ LENGTH(ans),NOBS);
+ }
PROTECT(nm = GET_NAMES(ans)); nprotect++;
if (!isNull(nm)) { // match names against names from data slot
PROTECT(oidx = matchnames(OBNM,nm)); nprotect++;
Modified: pkg/src/skeleton.c
===================================================================
--- pkg/src/skeleton.c 2012-02-06 20:14:40 UTC (rev 607)
+++ pkg/src/skeleton.c 2012-02-09 19:42:03 UTC (rev 608)
@@ -95,7 +95,8 @@
if (LENGTH(ans)!=NVAR) {
UNPROTECT(nprotect);
- error("skeleton error: user 'skeleton' must return a vector of length %d",NVAR);
+ error("user 'skeleton' returns a vector of %d state variables but %d are expected: compare initial conditions?",
+ LENGTH(ans),NVAR);
}
PROTECT(nm = GET_NAMES(ans)); nprotect++;
More information about the pomp-commits
mailing list