[Pomp-commits] r884 - in pkg/pomp: . R inst man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Feb 2 18:11:31 CET 2014
Author: kingaa
Date: 2014-02-02 18:11:31 +0100 (Sun, 02 Feb 2014)
New Revision: 884
Added:
pkg/pomp/R/abc-methods.R
pkg/pomp/R/abc.R
pkg/pomp/man/abc-methods.Rd
pkg/pomp/man/abc.Rd
pkg/pomp/tests/abc.R
pkg/pomp/tests/abc.Rout.save
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/NAMESPACE
pkg/pomp/R/pmcmc-methods.R
pkg/pomp/inst/NEWS
pkg/pomp/man/mif.Rd
pkg/pomp/man/pmcmc.Rd
Log:
- add ABC method
- minor tweaks to documentation for mif, pmcmc
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2014-01-14 17:48:08 UTC (rev 883)
+++ pkg/pomp/DESCRIPTION 2014-02-02 17:11:31 UTC (rev 884)
@@ -1,10 +1,9 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.46-1
-Date: 2014-01-10
-Authors at R: c(
- person(given=c("Aaron","A."),family="King",
+Version: 0.47-1
+Date: 2014-02-02
+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")),
person(given=c("Carles"),family="Breto",role=c("aut")),
@@ -23,6 +22,7 @@
License: GPL(>= 2)
LazyData: true
BuildVignettes: true
+MailingList: Subscribe to pomp-announce at r-forge.r-project.org for announcements by going to http://lists.r-forge.r-project.org/mailman/listinfo/pomp-announce.
Collate: aaa.R authors.R version.R eulermultinom.R plugins.R
parmat.R logmeanexp.R slice-design.R
profile-design.R sobol.R bsplines.R sannbox.R
@@ -35,4 +35,5 @@
pmcmc.R pmcmc-methods.R compare-pmcmc.R
nlf-funcs.R nlf-guts.R nlf-objfun.R nlf.R
probe.R probe-match.R basic-probes.R spect.R spect-match.R
+ abc.R abc-methods.R
builder.R example.R
Modified: pkg/pomp/NAMESPACE
===================================================================
--- pkg/pomp/NAMESPACE 2014-01-14 17:48:08 UTC (rev 883)
+++ pkg/pomp/NAMESPACE 2014-02-02 17:11:31 UTC (rev 884)
@@ -43,7 +43,8 @@
pmcmc,
traj.matched.pomp,
probed.pomp,probe.matched.pomp,
- spect.pomp,spect.matched.pomp
+ spect.pomp,spect.matched.pomp,
+ abc
)
exportMethods(
@@ -58,7 +59,8 @@
bsmc,
pmcmc,dprior,
spect,probe,
- probe.match,traj.match
+ probe.match,abc,
+ traj.match
)
export(
Added: pkg/pomp/R/abc-methods.R
===================================================================
--- pkg/pomp/R/abc-methods.R (rev 0)
+++ pkg/pomp/R/abc-methods.R 2014-02-02 17:11:31 UTC (rev 884)
@@ -0,0 +1,40 @@
+## this file contains short definitions of methods for the 'abc' class
+
+## extract the convergence record
+setMethod(
+ 'conv.rec',
+ 'abc',
+ function (object, pars, ...) {
+ if (missing(pars)) pars <- colnames(object at conv.rec)
+ object at conv.rec[,pars]
+ }
+ )
+
+## plot pmcmc object
+setMethod(
+ "plot",
+ "abc",
+ function (x, y, pars, scatter = FALSE, ...) {
+ if (missing(pars)) pars <- x at pars
+ if (scatter) {
+ pairs(conv.rec(x, pars))
+ } else {
+ plot.ts(conv.rec(x,pars),main="Convergence record")
+ }
+ }
+ )
+
+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
+ )
+ )
+ }
+ )
Added: pkg/pomp/R/abc.R
===================================================================
--- pkg/pomp/R/abc.R (rev 0)
+++ pkg/pomp/R/abc.R 2014-02-02 17:11:31 UTC (rev 884)
@@ -0,0 +1,387 @@
+## define the abc class
+setClass(
+ 'abc',
+ contains='pomp',
+ slots=c(
+ pars = 'character',
+ transform = 'logical',
+ Nabc = 'integer',
+ dprior = 'function',
+ probes='list',
+ scale = 'numeric',
+ epsilon = 'numeric',
+ hyperparams = 'list',
+ random.walk.sd = 'numeric',
+ conv.rec = 'matrix'
+ )
+ )
+
+## ABC algorithm functions
+setGeneric('abc',function(object,...)standardGeneric("abc"))
+
+abc.internal <- function (object, Nabc,
+ start, pars, dprior.fun,
+ rw.sd, probes,
+ hyperparams,
+ epsilon, scale,
+ verbose, transform,
+ .ndone, .getnativesymbolinfo = TRUE,
+ ...) {
+
+ gnsi <- as.logical(.getnativesymbolinfo)
+
+ transform <- as.logical(transform)
+
+ Nabc <- as.integer(Nabc)
+
+ if (length(start)==0)
+ stop(
+ "abc error: ",sQuote("start")," must be specified if ",
+ sQuote("coef(object)")," is NULL",
+ call.=FALSE
+ )
+
+ start.names <- names(start)
+ 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.character(pars) ||
+ !all(pars%in%start.names) ||
+ !all(pars%in%rw.names)
+ )
+ stop(
+ "abc 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
+ )
+
+ 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")
+ if (!all(sapply(probes,function(f)length(formals(f))==1)))
+ stop("each probe must be a function of a single argument")
+
+ 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)
+ }
+
+ theta <- start
+ log.prior <- dprior.fun(params=theta,hyperparams=hyperparams,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()
+
+ conv.rec <- matrix(
+ data=NA,
+ nrow=Nabc+1,
+ ncol=length(theta),
+ dimnames=list(
+ rownames=seq(from=0,to=Nabc,by=1),
+ colnames=names(theta)
+ )
+ )
+
+ if (!all(is.finite(theta[pars]))) {
+ stop(
+ sQuote("abc"),
+ " error: cannot estimate non-finite parameters: ",
+ paste(
+ pars[!is.finite(theta[pars])],
+ collapse=","
+ ),
+ call.=FALSE
+ )
+ }
+
+ po <- as(object,"pomp")
+
+ ## apply probes to data
+ datval <- try(.Call(apply_probe_data,po,probes),silent=FALSE)
+ if (inherits(datval,'try-error'))
+ stop("abc error: error in ",sQuote("apply_probe_data"),call.=FALSE)
+
+ conv.rec[1,names(theta)] <- theta
+
+ for (n in seq_len(Nabc)) { # main loop
+
+ theta.prop <- theta
+ theta.prop[pars] <- rnorm(n=length(pars),mean=theta.prop[pars],sd=rw.sd)
+
+ ## compute the probes for the proposed new parameter values
+
+ simval <- try(
+ .Call(
+ apply_probe_sim,
+ object=po,
+ nsim=1,
+ params=theta.prop,
+ seed=NULL,
+ probes=probes,
+ datval=datval
+ ),
+ silent=FALSE
+ )
+
+ if (inherits(simval,'try-error'))
+ stop("abc error: error in ",sQuote("apply_probe_sim"),call.=FALSE)
+
+ ## 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 (runif(1) < exp(log.prior.prop-log.prior)) {
+ theta <- theta.prop
+ log.prior <- log.prior.prop
+ }
+ }
+
+ ## store a record of this iteration
+ conv.rec[n+1,names(theta)] <- theta
+ if (verbose && (n%%5==0))
+ cat("ABC iteration ",n," of ",Nabc," completed\n")
+
+ }
+
+ new(
+ 'abc',
+ po,
+ params=theta,
+ pars = pars,
+ transform = transform,
+ Nabc = Nabc,
+ dprior = dprior.fun,
+ probes=probes,
+ scale = scale,
+ epsilon = epsilon,
+ hyperparams = hyperparams,
+ random.walk.sd = rw.sd,
+ conv.rec=conv.rec
+ )
+
+}
+
+setMethod(
+ "abc",
+ signature=signature(object="pomp"),
+ function (object, Nabc = 1,
+ start, pars, rw.sd,
+ dprior, probes, scale, epsilon, hyperparams,
+ 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))
+ stop("abc error: ",sQuote("probes")," must be specified",call.=FALSE)
+
+ if (missing(scale))
+ 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)
+
+ 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
+ )
+ }
+ )
+
+setMethod(
+ "abc",
+ signature=signature(object="probed.pomp"),
+ function (object, Nabc = 1,
+ start, pars, rw.sd,
+ dprior, probes, scale, epsilon, hyperparams,
+ 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
+ )
+ }
+ )
+
+setMethod(
+ "abc",
+ signature=signature(object="abc"),
+ function (object, Nabc,
+ start, pars, rw.sd,
+ dprior, probes, scale, epsilon, hyperparams,
+ verbose = getOption("verbose"),
+ transform,
+ ...) {
+
+ 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(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(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
+ )
+ }
+ )
+
+setMethod(
+ 'continue',
+ signature=signature(object='abc'),
+ function (object, Nabc = 1, ...) {
+
+ ndone <- object at Nabc
+
+ obj <- abc(
+ object=object,
+ Nabc=Nabc,
+ .ndone=ndone,
+ ...
+ )
+
+ obj at conv.rec <- rbind(
+ object at conv.rec[,colnames(obj at conv.rec)],
+ obj at conv.rec[-1,]
+ )
+ obj at Nabc <- as.integer(ndone+Nabc)
+ obj
+ }
+ )
Modified: pkg/pomp/R/pmcmc-methods.R
===================================================================
--- pkg/pomp/R/pmcmc-methods.R 2014-01-14 17:48:08 UTC (rev 883)
+++ pkg/pomp/R/pmcmc-methods.R 2014-02-02 17:11:31 UTC (rev 884)
@@ -31,4 +31,3 @@
compare.pmcmc(x)
}
)
-
Modified: pkg/pomp/inst/NEWS
===================================================================
--- pkg/pomp/inst/NEWS 2014-01-14 17:48:08 UTC (rev 883)
+++ pkg/pomp/inst/NEWS 2014-02-02 17:11:31 UTC (rev 884)
@@ -1,4 +1,7 @@
NEWS
+0.47-1
+ o 'abc' implements Approximate Bayesian Computation for pomp models.
+
0.46-1
o 'pompExample' now has an optional argument, 'envir', determining which environment the pomp object will be loaded into.
Added: pkg/pomp/man/abc-methods.Rd
===================================================================
--- pkg/pomp/man/abc-methods.Rd (rev 0)
+++ pkg/pomp/man/abc-methods.Rd 2014-02-02 17:11:31 UTC (rev 884)
@@ -0,0 +1,43 @@
+\name{abc-methods}
+\docType{methods}
+\alias{abc-methods}
+\alias{dprior,abc-method}
+\alias{dprior-abc}
+\alias{conv.rec,abc-method}
+\alias{conv.rec-abc}
+\alias{plot-abc}
+\alias{plot,abc-method}
+\title{Methods of the "abc" class}
+\description{Methods of the "abc" class.}
+\usage{
+\S4method{conv.rec}{abc}(object, pars, \dots)
+\S4method{plot}{abc}(x, y, pars, scatter = FALSE, \dots)
+\S4method{dprior}{abc}(object, params, log = FALSE, \dots)
+}
+\arguments{
+ \item{object, x}{The \code{abc} object.}
+ \item{pars}{Names of parameters.}
+ \item{y}{Ignored.}
+ \item{scatter}{WHAT DOES THIS DO?}
+ \item{params}{
+ Named vector of parameters.
+ }
+ \item{log}{if TRUE, log probabilities are returned.}
+ \item{\dots}{
+ Further arguments (either ignored or passed to underlying functions).
+ }
+}
+\section{Methods}{
+ \describe{
+ \item{conv.rec}{
+ \code{conv.rec(object, pars = NULL)} returns the columns of the convergence-record matrix corresponding to the names in \code{pars}.
+ By default, all rows are returned.
+ }
+ \item{plot}{
+ Plots a series of diagnostic plots.
+ }
+ }
+}
+\author{Edward L. Ionides \email{ionides at umich dot edu}, Aaron A. King \email{kingaa at umich dot edu}}
+\seealso{\code{\link{abc}}, \code{\link{pomp}}}
+\keyword{ts}
Added: pkg/pomp/man/abc.Rd
===================================================================
--- pkg/pomp/man/abc.Rd (rev 0)
+++ pkg/pomp/man/abc.Rd 2014-02-02 17:11:31 UTC (rev 884)
@@ -0,0 +1,124 @@
+\name{abc}
+\docType{methods}
+\alias{abc}
+\alias{abc-abc}
+\alias{abc,probed.pomp-method}
+\alias{abc-probed.pomp}
+\alias{abc,pomp-method}
+\alias{abc-pomp}
+\alias{abc,abc-method}
+\alias{abc-abc}
+\alias{continue,abc-method}
+\alias{continue-abc}
+\alias{abc-class}
+\title{The ABC algorithm}
+\description{
+ 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, dprior, probes, scale, epsilon, hyperparams,
+ verbose = getOption("verbose"), transform = FALSE, \dots)
+\S4method{abc}{probed.pomp}(object, Nabc = 1, start, pars,
+ rw.sd, dprior, probes, scale, epsilon, hyperparams,
+ verbose = getOption("verbose"), transform = FALSE, \dots)
+\S4method{abc}{abc}(object, Nabc, start, pars,
+ rw.sd, dprior, probes, scale, epsilon, hyperparams,
+ verbose = getOption("verbose"), transform, \dots)
+\S4method{continue}{abc}(object, Nabc = 1, \dots)
+}
+\arguments{
+ \item{object}{
+ An object of class \code{pomp}.
+ }
+ \item{Nabc}{
+ The number of ABC iterations to perform.
+ }
+ \item{start}{
+ 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{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{dprior}{
+ Function of prototype \code{dprior(params,hyperparams,...,log)} that evaluates the prior density.
+ This defaults to an improper uniform prior.
+ }
+ \item{probes}{
+ }
+ \item{scale}{
+
+ }
+ \item{epsilon}{
+
+ }
+ \item{hyperparams}{
+ optional list; parameters to be passed to \code{dprior}.
+ }
+ \item{verbose}{
+ logical; if TRUE, print progress reports.
+ }
+ \item{transform}{
+ logical;
+ if \code{TRUE}, optimization is performed on the transformed scale.
+ }
+ \item{\dots}{
+ Additional arguments.
+ These are currently ignored.
+ }
+}
+\value{
+ An object of class \code{abc}.
+ This class inherits from class \code{\link[=probed.pomp-class]{probed.pomp}} and contains the following additional slots:
+ \describe{
+ \item{pars, Nabc, dprior, hyperparams, transform, 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}{
+ To re-run a sequence of ABC iterations, one can use the \code{abc} method on a \code{abc} object.
+ By default, the same parameters used for the original ABC run are re-used (except for \code{tol}, \code{max.fail}, and \code{verbose}, the defaults of which are shown above).
+ If one does specify additional arguments, these will override the defaults.
+}
+\section{Continuing ABC Iterations}{
+ One can continue a series of ABC iterations from where one left off using the \code{continue} method.
+ A call to \code{abc} to perform \code{Nabc=m} iterations followed by a call to \code{continue} to perform \code{Nabc=n} iterations will produce precisely the same effect as a single call to \code{abc} to perform \code{Nabc=m+n} iterations.
+ By default, all the algorithmic parameters are the same as used in the original call to \code{abc}.
+ Additional arguments will override the defaults.
+}
+\section{Details}{
+ TO APPEAR.
+}
+\references{
+ T. Toni and M. P. H. Stumpf,
+ Simulation-based model selection for dynamical systems in systems and population biology,
+ Bioinformatics 26:104--110, 2010.
+
+ T. Toni, D. Welch, N. Strelkowa, A. Ipsen, and M. P. H. Stumpf,
+ Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems
+ Journal of the Royal Society, Interface 6:187--202, 2009.
+}
+\author{Edward L. Ionides \email{ionides at umich dot edu}, Aaron A. King \email{kingaa at umich dot edu}}
+\seealso{
+ \link{abc-methods}, \code{\link{pomp}}, \code{\link{probe}}.
+ See the \dQuote{intro_to_pomp} vignette for an example
+}
+\keyword{ts}
Modified: pkg/pomp/man/mif.Rd
===================================================================
--- pkg/pomp/man/mif.Rd 2014-01-14 17:48:08 UTC (rev 883)
+++ pkg/pomp/man/mif.Rd 2014-02-02 17:11:31 UTC (rev 884)
@@ -175,5 +175,4 @@
\code{\link{mif-methods}}, \code{\link{pomp}}, \code{\link{pomp-class}}, \code{\link{pfilter}}.
See the \dQuote{intro_to_pomp} vignette for examples.
}
-\keyword{models}
\keyword{ts}
Modified: pkg/pomp/man/pmcmc.Rd
===================================================================
--- pkg/pomp/man/pmcmc.Rd 2014-01-14 17:48:08 UTC (rev 883)
+++ pkg/pomp/man/pmcmc.Rd 2014-02-02 17:11:31 UTC (rev 884)
@@ -124,5 +124,4 @@
See the \dQuote{intro_to_pomp} vignette for an example
[CURRENTLY, ONLY DEMONSTRATING THE MIF ALGORITHM, WHICH IS ALGORITHMICALLY VERY SIMILAR TO PMCMC SINCE THEY BOTH DEPEND CRITICALLY ON A PARTICLE FILTERING STEP].
}
-\keyword{models}
\keyword{ts}
Added: pkg/pomp/tests/abc.R
===================================================================
--- pkg/pomp/tests/abc.R (rev 0)
+++ pkg/pomp/tests/abc.R 2014-02-02 17:11:31 UTC (rev 884)
@@ -0,0 +1,48 @@
+### OU2 test of abc for pomp; nov 2013
+
+library(pomp)
+pompExample(ou2)
+plot(ou2)
+
+pdf(file='abc.pdf')
+
+set.seed(2079015564L)
+
+probes.good <- list(
+ y1.mean=probe.mean(var="y1"),
+ y2.mean=probe.mean(var="y2"),
+ probe.acf(var="y1",lags=c(0,5)),
+ probe.acf(var="y2",lags=c(0,5)),
+ probe.ccf(vars=c("y1","y2"),lags=0)
+ )
+psim <- probe(ou2,probes=probes.good,nsim=200)
+plot(psim)
+## why do simulations sometimes seem extreme with respect to these probes?
+
+scale.dat <- apply(psim at simvals,2,sd)
+
+po <- ou2
+
+theta <- coef(ou2)
+
+abc1 <- abc(po,
+ Nabc=10000,
+ start=coef(ou2),
+ pars=c("alpha.1","alpha.2"),
+ probes=probes.good,
+ scale=scale.dat,
+ epsilon=3,
+ rw.sd= c(alpha.1=0.01,alpha.2=0.01),
+ hyperparams=list(junk=0),
+ verbose=TRUE
+ )
+
+plot(abc1,scatter=TRUE)
+plot(abc1)
+
+## check how sticky the chain is:
+runs <- rle(as.vector(conv.rec(abc1)[, "alpha.1"]))
+hist(runs$lengths)
+mean(runs$length)
+
+dev.off()
Added: pkg/pomp/tests/abc.Rout.save
===================================================================
--- pkg/pomp/tests/abc.Rout.save (rev 0)
+++ pkg/pomp/tests/abc.Rout.save 2014-02-02 17:11:31 UTC (rev 884)
@@ -0,0 +1,2080 @@
+
+R version 3.0.2 (2013-09-25) -- "Frisbee Sailing"
+Copyright (C) 2013 The R Foundation for Statistical Computing
+Platform: x86_64-unknown-linux-gnu (64-bit)
+
+R is free software and comes with ABSOLUTELY NO WARRANTY.
+You are welcome to redistribute it under certain conditions.
+Type 'license()' or 'licence()' for distribution details.
+
+R is a collaborative project with many contributors.
+Type 'contributors()' for more information and
+'citation()' on how to cite R or R packages in publications.
+
+Type 'demo()' for some demos, 'help()' for on-line help, or
+'help.start()' for an HTML browser interface to help.
+Type 'q()' to quit R.
+
+> ### OU2 test of abc for pomp; nov 2013
+>
+> library(pomp)
+Loading required package: mvtnorm
+Loading required package: subplex
+Loading required package: deSolve
+> pompExample(ou2)
+newly created pomp object(s):
+ ou2
+> plot(ou2)
+>
+> pdf(file='abc.pdf')
+>
+> set.seed(2079015564L)
+>
+> probes.good <- list(
++ y1.mean=probe.mean(var="y1"),
++ y2.mean=probe.mean(var="y2"),
++ probe.acf(var="y1",lags=c(0,5)),
++ probe.acf(var="y2",lags=c(0,5)),
++ probe.ccf(vars=c("y1","y2"),lags=0)
++ )
+> psim <- probe(ou2,probes=probes.good,nsim=200)
+> plot(psim)
+> ## why do simulations sometimes seem extreme with respect to these probes?
+>
+> scale.dat <- apply(psim at simvals,2,sd)
+>
+> po <- ou2
+>
+> theta <- coef(ou2)
+>
+> abc1 <- abc(po,
++ Nabc=10000,
++ start=coef(ou2),
++ pars=c("alpha.1","alpha.2"),
++ probes=probes.good,
++ scale=scale.dat,
++ epsilon=3,
++ rw.sd= c(alpha.1=0.01,alpha.2=0.01),
++ hyperparams=list(junk=0),
++ verbose=TRUE
++ )
+performing 10000 ABC iteration(s) to estimate parameter(s) alpha.1, alpha.2 using random-walk with SD
+alpha.1 alpha.2
+ 0.01 0.01
+ABC iteration 5 of 10000 completed
+ABC iteration 10 of 10000 completed
+ABC iteration 15 of 10000 completed
+ABC iteration 20 of 10000 completed
+ABC iteration 25 of 10000 completed
+ABC iteration 30 of 10000 completed
+ABC iteration 35 of 10000 completed
+ABC iteration 40 of 10000 completed
+ABC iteration 45 of 10000 completed
+ABC iteration 50 of 10000 completed
+ABC iteration 55 of 10000 completed
+ABC iteration 60 of 10000 completed
+ABC iteration 65 of 10000 completed
+ABC iteration 70 of 10000 completed
+ABC iteration 75 of 10000 completed
+ABC iteration 80 of 10000 completed
+ABC iteration 85 of 10000 completed
+ABC iteration 90 of 10000 completed
+ABC iteration 95 of 10000 completed
+ABC iteration 100 of 10000 completed
+ABC iteration 105 of 10000 completed
+ABC iteration 110 of 10000 completed
+ABC iteration 115 of 10000 completed
+ABC iteration 120 of 10000 completed
+ABC iteration 125 of 10000 completed
+ABC iteration 130 of 10000 completed
+ABC iteration 135 of 10000 completed
+ABC iteration 140 of 10000 completed
+ABC iteration 145 of 10000 completed
+ABC iteration 150 of 10000 completed
+ABC iteration 155 of 10000 completed
+ABC iteration 160 of 10000 completed
+ABC iteration 165 of 10000 completed
+ABC iteration 170 of 10000 completed
+ABC iteration 175 of 10000 completed
+ABC iteration 180 of 10000 completed
+ABC iteration 185 of 10000 completed
+ABC iteration 190 of 10000 completed
+ABC iteration 195 of 10000 completed
+ABC iteration 200 of 10000 completed
+ABC iteration 205 of 10000 completed
+ABC iteration 210 of 10000 completed
+ABC iteration 215 of 10000 completed
+ABC iteration 220 of 10000 completed
+ABC iteration 225 of 10000 completed
+ABC iteration 230 of 10000 completed
+ABC iteration 235 of 10000 completed
+ABC iteration 240 of 10000 completed
+ABC iteration 245 of 10000 completed
+ABC iteration 250 of 10000 completed
+ABC iteration 255 of 10000 completed
+ABC iteration 260 of 10000 completed
+ABC iteration 265 of 10000 completed
+ABC iteration 270 of 10000 completed
+ABC iteration 275 of 10000 completed
+ABC iteration 280 of 10000 completed
+ABC iteration 285 of 10000 completed
+ABC iteration 290 of 10000 completed
+ABC iteration 295 of 10000 completed
+ABC iteration 300 of 10000 completed
+ABC iteration 305 of 10000 completed
+ABC iteration 310 of 10000 completed
+ABC iteration 315 of 10000 completed
+ABC iteration 320 of 10000 completed
+ABC iteration 325 of 10000 completed
+ABC iteration 330 of 10000 completed
+ABC iteration 335 of 10000 completed
+ABC iteration 340 of 10000 completed
+ABC iteration 345 of 10000 completed
+ABC iteration 350 of 10000 completed
+ABC iteration 355 of 10000 completed
+ABC iteration 360 of 10000 completed
+ABC iteration 365 of 10000 completed
+ABC iteration 370 of 10000 completed
+ABC iteration 375 of 10000 completed
+ABC iteration 380 of 10000 completed
+ABC iteration 385 of 10000 completed
+ABC iteration 390 of 10000 completed
+ABC iteration 395 of 10000 completed
+ABC iteration 400 of 10000 completed
+ABC iteration 405 of 10000 completed
+ABC iteration 410 of 10000 completed
+ABC iteration 415 of 10000 completed
+ABC iteration 420 of 10000 completed
+ABC iteration 425 of 10000 completed
+ABC iteration 430 of 10000 completed
+ABC iteration 435 of 10000 completed
+ABC iteration 440 of 10000 completed
+ABC iteration 445 of 10000 completed
+ABC iteration 450 of 10000 completed
+ABC iteration 455 of 10000 completed
+ABC iteration 460 of 10000 completed
+ABC iteration 465 of 10000 completed
+ABC iteration 470 of 10000 completed
+ABC iteration 475 of 10000 completed
+ABC iteration 480 of 10000 completed
+ABC iteration 485 of 10000 completed
+ABC iteration 490 of 10000 completed
+ABC iteration 495 of 10000 completed
+ABC iteration 500 of 10000 completed
+ABC iteration 505 of 10000 completed
+ABC iteration 510 of 10000 completed
+ABC iteration 515 of 10000 completed
+ABC iteration 520 of 10000 completed
+ABC iteration 525 of 10000 completed
+ABC iteration 530 of 10000 completed
+ABC iteration 535 of 10000 completed
+ABC iteration 540 of 10000 completed
+ABC iteration 545 of 10000 completed
+ABC iteration 550 of 10000 completed
+ABC iteration 555 of 10000 completed
+ABC iteration 560 of 10000 completed
+ABC iteration 565 of 10000 completed
+ABC iteration 570 of 10000 completed
+ABC iteration 575 of 10000 completed
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 884
More information about the pomp-commits
mailing list