[Pomp-commits] r439 - in pkg: . R inst inst/doc man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 13 00:29:46 CEST 2011
Author: kingaa
Date: 2011-04-13 00:29:44 +0200 (Wed, 13 Apr 2011)
New Revision: 439
Modified:
pkg/DESCRIPTION
pkg/R/plugins.R
pkg/R/probe-match.R
pkg/inst/ChangeLog
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.Rnw
pkg/inst/doc/intro_to_pomp.pdf
pkg/man/plugins.Rd
Log:
- add 'sannbox' constrained simulated annealing method for probe matching
- re-implement the 'discrete.time.sim' plugin
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2011-04-12 22:07:16 UTC (rev 438)
+++ pkg/DESCRIPTION 2011-04-12 22:29:44 UTC (rev 439)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.36-2
-Date: 2011-02-16
+Date: 2011-04-12
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>
Modified: pkg/R/plugins.R
===================================================================
--- pkg/R/plugins.R 2011-04-12 22:07:16 UTC (rev 438)
+++ pkg/R/plugins.R 2011-04-12 22:29:44 UTC (rev 439)
@@ -29,37 +29,6 @@
}
}
-discrete.time.sim <- function (step.fun, PACKAGE) {
- efun <- pomp.fun(
- f=step.fun,
- PACKAGE=PACKAGE,
- proto=quote(step.fun(x,t,params,...))
- )
- function (xstart, times, params, ...,
- statenames = character(0),
- paramnames = character(0),
- covarnames = character(0),
- zeronames = character(0),
- tcovar, covar) {
- .Call(
- euler_model_simulator,
- func=efun,
- xstart=xstart,
- times=times,
- params=params,
- dt=1,
- method=0L,
- statenames=statenames,
- paramnames=paramnames,
- covarnames=covarnames,
- zeronames=zeronames,
- tcovar=tcovar,
- covar=covar,
- args=pairlist(...)
- )
- }
-}
-
euler.sim <- function (step.fun, delta.t, PACKAGE) {
efun <- pomp.fun(
f=step.fun,
@@ -91,6 +60,9 @@
}
}
+discrete.time.sim <- function (step.fun, delta.t = 1, PACKAGE)
+ euler.sim(step.fun=step.fun,delta.t=delta.t,PACKAGE=PACKAGE)
+
onestep.dens <- function (dens.fun, PACKAGE) {
efun <- pomp.fun(
f=dens.fun,
Modified: pkg/R/probe-match.R
===================================================================
--- pkg/R/probe-match.R 2011-04-12 22:07:16 UTC (rev 438)
+++ pkg/R/probe-match.R 2011-04-12 22:29:44 UTC (rev 439)
@@ -151,6 +151,23 @@
fail.value=fail.value,
control=list(...)
)
+
+ } else if (method=="sannbox") {
+ opt <- sannbox(
+ par=guess,
+ fn=obj.fn,
+ est=par.index,
+ object=object,
+ probes=probes,
+ params=params,
+ nsim=nsim,
+ seed=seed,
+ weights=weights,
+ datval=datval,
+ fail.value=fail.value,
+ control=list(...)
+ )
+
} else {
opt <- optim(
par=guess,
@@ -202,7 +219,7 @@
function(object, start, est = character(0),
probes, weights,
nsim, seed = NULL,
- method = c("subplex","Nelder-Mead","SANN"),
+ method = c("subplex","Nelder-Mead","SANN","BFGS","sannbox"),
verbose = getOption("verbose"),
eval.only = FALSE, fail.value = NA, ...) {
@@ -241,7 +258,7 @@
function(object, start, est = character(0),
probes, weights,
nsim, seed = NULL,
- method = c("subplex","Nelder-Mead","SANN"),
+ method = c("subplex","Nelder-Mead","SANN","BFGS","sannbox"),
verbose = getOption("verbose"),
eval.only = FALSE, fail.value = NA, ...) {
@@ -280,7 +297,7 @@
function(object, start, est,
probes, weights,
nsim, seed = NULL,
- method = c("subplex","Nelder-Mead","SANN"),
+ method = c("subplex","Nelder-Mead","SANN","BFGS","sannbox"),
verbose = getOption("verbose"),
eval.only = FALSE, fail.value, ...) {
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2011-04-12 22:07:16 UTC (rev 438)
+++ pkg/inst/ChangeLog 2011-04-12 22:29:44 UTC (rev 439)
@@ -1,3 +1,17 @@
+2011-02-22 kingaa
+
+ * [r437] inst/doc/intro_to_pomp.Rnw, inst/doc/intro_to_pomp.pdf,
+ inst/doc/ou2-multi-mif.rda: - tinker with mif example in intro
+ vignette
+ * [r436] R/pomp.R: - do a check in 'pomp' to make sure tcovar
+ embraces times; warn else
+
+2011-02-19 kingaa
+
+ * [r435] inst/ChangeLog, inst/doc/advanced_topics_in_pomp.pdf,
+ inst/doc/intro_to_pomp.Rnw, inst/doc/intro_to_pomp.pdf: -
+ improvements to the intro vignette
+
2011-02-18 kingaa
* [r434] inst/ChangeLog, inst/doc/ou2-first-mif.rda,
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.Rnw
===================================================================
--- pkg/inst/doc/intro_to_pomp.Rnw 2011-04-12 22:07:16 UTC (rev 438)
+++ pkg/inst/doc/intro_to_pomp.Rnw 2011-04-12 22:29:44 UTC (rev 439)
@@ -329,7 +329,7 @@
Since some parameter estimation algorithms in the \pomp\ package only require simulations of the full process, we are already in a position to use them.
At present, the only such algorithm is nonlinear forecasting (\code{nlf}).
-In the near future, probe-matching algorithms that require only \code{rprocess} and \code{rmeasure} will be included in the package.
+Probe-matching algorithms that require only \code{rprocess} and \code{rmeasure} are also included in the package.
If we want to work with likelihood-based methods, however, we will need to be able to compute the likelihood of the data $Y_t$ given the states $X_t$.
To implement this in \pomp, we write another function:
<<ou2-meas-dens-def>>=
@@ -373,8 +373,8 @@
To compute the likelihood of the data, we can use the function \code{pfilter}.
This runs a plain vanilla particle filter (AKA sequential Monte Carlo) algorithm and results in an unbiased estimate of the likelihood.
-See \citet{Arulampalam2002} for an excellent tutorial on particle filtering.
-To do this, we must decide how many concurrrent realizations (\emph{particles}) to use: the larger the number of particles, the smaller the Monte Carlo error but the greater the computational effort.
+See \citet{Arulampalam2002} for an excellent tutorial on particle filtering and \citet{Ionides2006} for a pseudocode description of the algorithm implemented in \pomp.
+We must decide how many concurrent realizations (\emph{particles}) to use: the larger the number of particles, the smaller the Monte Carlo error but the greater the computational effort.
Let's run \code{pfilter} with 1000 particles to estimate the likelihood at the true parameters:
<<ou2-pfilter-truth,eval=F>>=
pf <- pfilter(ou2,params=theta,Np=1000)
@@ -521,8 +521,10 @@
\section{Estimating parameters using iterated filtering: \code{mif}}
Iterated filtering is a technique for maximizing the likelihood obtained by filtering.
-In \pomp, the particle filter is used as a basis for iterated filtering.
+In \pomp, it is the particle filter that is iterated.
Iterated filtering is implemented in the \code{mif} function.
+For a description of the algorithm and a description of its theoretical basis, see \citet{Ionides2006}.
+A more complete set of proofs is provided in a forthcoming paper.
The key idea of iterated filtering is to replace the model we are interested in fitting---which has time-invariant parameters---with a model that is just the same except that its parameters take a random walk in time.
As the intensity of this random walk approaches zero, the modified model approaches the original model.
@@ -940,7 +942,7 @@
\clearpage
\section{Nonlinear forecasting: \code{nlf}}
-To be added.
+\citet{Ellner1998,Kendall1999,Kendall2005}. To be added.
<<first-nlf,echo=F,eval=F>>=
estnames <- c('alpha.2','alpha.3','tau')
@@ -971,6 +973,14 @@
)
@
+\section{Bayesian sequential Monte Carlo: \code{bsmc}}
+
+\citet{Liu2001b}. To be added.
+
+\section{Particle Markov chain Monte Carlo: \code{pmcmc}}
+
+\citet{Andrieu2010}. To be added.
+
\clearpage
\section{A more complex example: a seasonal epidemic model}
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/man/plugins.Rd
===================================================================
--- pkg/man/plugins.Rd 2011-04-12 22:07:16 UTC (rev 438)
+++ pkg/man/plugins.Rd 2011-04-12 22:29:44 UTC (rev 439)
@@ -13,7 +13,7 @@
\usage{
onestep.sim(step.fun, PACKAGE)
euler.sim(step.fun, delta.t, PACKAGE)
-discrete.time.sim(step.fun, PACKAGE)
+discrete.time.sim(step.fun, delta.t = 1, PACKAGE)
gillespie.sim(rate.fun, v, d, PACKAGE)
onestep.dens(dens.fun, PACKAGE)
}
@@ -68,8 +68,7 @@
\code{euler.sim} is appropriate when one cannot do this but can compute the change in state via a sequence of smaller steps.
This is desirable, for example, if one is simulating a continuous time process but is willing to approximate it using an Euler approach.
\code{discrete.time.sim} is appropriate when the process evolves in discrete time.
- In this case, it is assumed that the intervals between observations are integers.
- Beware, however: this assumption is not checked.
+ In this case, by default, the intervals between observations are integers.
To use \code{euler.sim} or \code{discrete.time.sim}, you must write a function \code{step.fun} that will take a single Euler step, of size at most \code{delta.t}.
\code{euler.sim} and \code{discrete.time.sim} will create simulators that take as many steps as needed to get from one time to another.
More information about the pomp-commits
mailing list