[Pomp-commits] r275 - in pkg: . R inst inst/doc man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jun 30 21:30:22 CEST 2010
Author: kingaa
Date: 2010-06-30 21:30:22 +0200 (Wed, 30 Jun 2010)
New Revision: 275
Added:
pkg/R/compare.pmcmc.R
pkg/R/pmcmc-methods.R
pkg/R/pmcmc.R
pkg/man/pmcmc-class.Rd
pkg/man/pmcmc-methods.Rd
pkg/man/pmcmc.Rd
pkg/tests/ou2-pmcmc.R
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/simulate-pomp.R
pkg/inst/NEWS
pkg/inst/doc/pomp.bib
Log:
- add PMCMC method from Andrieu et al. (2010)
- add error check on 'nsim' in 'simulate' method for 'pomp'
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-06-24 15:01:26 UTC (rev 274)
+++ pkg/DESCRIPTION 2010-06-30 19:30:22 UTC (rev 275)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.30-1
-Date: 2010-06-23
+Version: 0.31-1
+Date: 2010-06-30
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing,
Matthew J. Ferrari, Michael Lavine
Maintainer: Aaron A. King <kingaa at umich.edu>
@@ -19,4 +19,5 @@
trajectory-pomp.R plot-pomp.R init.state-pomp.R
pfilter.R trajmatch.R bsmc.R
mif-class.R particles-mif.R pfilter-mif.R mif.R mif-methods.R compare.mif.R
+ pmcmc.R pmcmc-methods.R compare.pmcmc.R
nlf-funcs.R nlf-guts.R nlf-lql.R nlf-objfun.R nlf.R
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-06-24 15:01:26 UTC (rev 274)
+++ pkg/NAMESPACE 2010-06-30 19:30:22 UTC (rev 275)
@@ -23,16 +23,17 @@
importFrom(deSolve,lsoda)
importFrom(subplex,subplex)
-exportClasses('pomp','mif')
+exportClasses('pomp','mif','pmcmc')
exportMethods(
'plot','show','print','coerce',
'dprocess','rprocess','rmeasure','dmeasure','init.state','skeleton',
- 'simulate','data.array','pfilter',
- 'coef','logLik','time','time<-',
+ 'data.array','coef','logLik','time','time<-',
+ 'simulate','pfilter',
+ 'particles','mif','continue','coef<-','states','trajectory',
'pred.mean','pred.var','filter.mean','conv.rec',
- 'particles','mif','continue','coef<-','states','trajectory',
- 'bsmc'
+ 'bsmc',
+ 'pmcmc','dprior'
)
export(
Added: pkg/R/compare.pmcmc.R
===================================================================
--- pkg/R/compare.pmcmc.R (rev 0)
+++ pkg/R/compare.pmcmc.R 2010-06-30 19:30:22 UTC (rev 275)
@@ -0,0 +1,103 @@
+compare.pmcmc <- function (z) {
+ ## assumes that x is a list of pmcmcs with identical structure
+ if (!is.list(z)) z <- list(z)
+ if (!all(sapply(z,function(x)is(x,'pmcmc'))))
+ stop("compare.pmcmc error: ",sQuote("z")," must be a pmcmc object or a list of pmcmc objects",call.=FALSE)
+ mar.multi <- c(0,5.1,0,2.1)
+ oma.multi <- c(6,0,5,0)
+ xx <- z[[1]]
+ estnames <- xx at pars
+ parnames <- names(coef(xx))
+ unestnames <- parnames[-match(estnames,parnames)]
+
+ ## plot filter means
+ filt.diag <- rbind("eff. sample size"=xx at eff.sample.size,filter.mean(xx))
+ filtnames <- rownames(filt.diag)
+# plotnames <- if(length(unestnames)>0) filtnames[-match(unestnames,filtnames)] else filtnames
+ plotnames <- filtnames
+ lognames <- filtnames[1] # eff. sample size
+ nplots <- length(plotnames)
+ n.per.page <- min(nplots,10)
+ if(n.per.page<=4) nc <- 1 else nc <- 2
+ nr <- ceiling(n.per.page/nc)
+ oldpar <- par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc),ask=dev.interactive(orNone=TRUE))
+ on.exit(par(oldpar))
+ low <- 1
+ hi <- 0
+ time <- time(xx)
+ while (hi<nplots) {
+ hi <- min(low+n.per.page-1,nplots)
+ for (i in seq(from=low,to=hi,by=1)) {
+ n <- i-low+1
+ logplot <- if (plotnames[i]%in%lognames) "y" else ""
+ dat <- sapply(
+ z,
+ function(po, label) {
+ if (label=="eff. sample size")
+ po at eff.sample.size
+ else
+ filter.mean(po,label)
+ },
+ label=plotnames[i]
+ )
+ matplot(
+ y=dat,
+ x=time,
+ axes = FALSE,
+ xlab = "",
+ log=logplot,
+ ylab = "",
+ type = "l"
+ )
+ box()
+ y.side <- 2
+ axis(y.side, xpd = NA)
+ mtext(plotnames[i], y.side, line = 3)
+ do.xax <- (n%%nr==0||n==n.per.page)
+ if (do.xax) axis(1,xpd=NA)
+ if (do.xax) mtext("time",side=1,line=3)
+ }
+ low <- hi+1
+ mtext("Filter diagnostics (last iteration)",3,line=2,outer=TRUE)
+ }
+
+ ## plot pmcmc convergence diagnostics
+ other.diagnostics <- c("loglik", "log.prior","nfail")
+ plotnames <- c(other.diagnostics,estnames)
+ nplots <- length(plotnames)
+ n.per.page <- min(nplots,10)
+ nc <- if (n.per.page<=4) 1 else 2
+ nr <- ceiling(n.per.page/nc)
+ par(mar=mar.multi,oma=oma.multi,mfcol=c(nr,nc))
+ ## on.exit(par(oldpar))
+ low <- 1
+ hi <- 0
+ iteration <- seq(0,xx at Nmcmc)
+ while (hi<nplots) {
+ hi <- min(low+n.per.page-1,nplots)
+ for (i in seq(from=low,to=hi,by=1)) {
+ n <- i-low+1
+ dat <- sapply(z,function(po,label) conv.rec(po,label),label=plotnames[i])
+ matplot(
+ y=dat,
+ x=iteration,
+ axes = FALSE,
+ xlab = "",
+ ylab = "",
+ type = "l"
+ )
+ box()
+ y.side <- 2
+ axis(y.side,xpd=NA)
+ mtext(plotnames[i],y.side,line=3)
+ do.xax <- (n%%nr==0||n==n.per.page)
+ if (do.xax) axis(1,xpd=NA)
+ if (do.xax) mtext("PMCMC iteration",side=1,line=3)
+ }
+ low <- hi+1
+ mtext("PMCMC convergence diagnostics",3,line=2,outer=TRUE)
+ }
+ invisible(NULL)
+}
+
+
Added: pkg/R/pmcmc-methods.R
===================================================================
--- pkg/R/pmcmc-methods.R (rev 0)
+++ pkg/R/pmcmc-methods.R 2010-06-30 19:30:22 UTC (rev 275)
@@ -0,0 +1,34 @@
+## this file contains short definitions of methods for the 'pmcmc' class
+
+## extract the estimated log likelihood
+setMethod('logLik','pmcmc',function(object,...)object at loglik)
+
+## extract the filtering means
+setMethod(
+ 'filter.mean',
+ 'pmcmc',
+ function (object, pars, ...) {
+ if (missing(pars)) pars <- rownames(object at filter.mean)
+ object at filter.mean[pars,]
+ }
+ )
+
+## extract the convergence record
+setMethod(
+ 'conv.rec',
+ 'pmcmc',
+ function (object, pars, ...) {
+ if (missing(pars)) pars <- colnames(object at conv.rec)
+ object at conv.rec[,pars]
+ }
+ )
+
+## plot pmcmc object
+setMethod(
+ "plot",
+ "pmcmc",
+ function (x, y = NULL, ...) {
+ compare.pmcmc(x)
+ }
+ )
+
Added: pkg/R/pmcmc.R
===================================================================
--- pkg/R/pmcmc.R (rev 0)
+++ pkg/R/pmcmc.R 2010-06-30 19:30:22 UTC (rev 275)
@@ -0,0 +1,354 @@
+## define the pmcmc class
+setClass(
+ 'pmcmc',
+ representation(
+ pars = 'character',
+ Nmcmc = 'integer',
+ dprior = 'function',
+ hyperparams = 'list',
+ Np = 'integer',
+ random.walk.sd = 'numeric',
+ filter.mean = 'matrix',
+ conv.rec = 'matrix',
+ eff.sample.size = 'numeric',
+ cond.loglik = 'numeric',
+ loglik = 'numeric',
+ log.prior = 'numeric'
+ ),
+ contains='pomp'
+ )
+
+## PMCMC algorithm functions
+pmcmc <- function (object, ... )
+ stop("function ",sQuote("pmcmc")," is undefined for objects of class ",sQuote(class(object)))
+setGeneric('pmcmc')
+
+dprior <- function (object, params, log)
+ stop("function ",sQuote("dprior")," is undefined for objects of class ",sQuote(class(object)))
+setGeneric('dprior')
+
+setMethod(
+ "dprior",
+ "pmcmc",
+ function (object, params, log = FALSE) {
+ do.call(object at dprior,list(params=params,hyperparams=object at hyperparams,log=log))
+ }
+ )
+
+pmcmc.internal <- function (object, Nmcmc = 1,
+ start = NULL,
+ pars = NULL,
+ dprior.fun = NULL,
+ rw.sd = NULL,
+ Np = NULL,
+ hyperparams = NULL,
+ tol = 1e-17, max.fail = 0,
+ verbose = FALSE,
+ .ndone = 0, .prevcomp = NULL
+ ) {
+ is.pmcmc <- is(object,"pmcmc")
+ if (is.null(start)) {
+ start <- coef(object)
+ if (length(start)==0)
+ stop("pmcmc error: ",sQuote("start")," must be specified if ",
+ sQuote("coef(object)")," is NULL",
+ call.=FALSE)
+ } else if (length(start)==0)
+ stop("pmcmc error: ",sQuote("start")," must be a named vector",call.=FALSE)
+ start.names <- names(start)
+ if (is.null(start.names))
+ stop("pmcmc error: ",sQuote("start")," must be a named vector",call.=FALSE)
+
+ if (is.null(rw.sd)) {
+ if (is.pmcmc)
+ rw.sd <- object at random.walk.sd
+ else
+ 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.null(pars)) {
+ if (is.pmcmc)
+ pars <- object at pars
+ else
+ 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
+ )
+
+ 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)
+
+ if (is.null(dprior.fun)) {
+ if (is.pmcmc)
+ dprior.fun <- object at dprior
+ else
+ stop("pmcmc error: ",sQuote("dprior")," must be specified",call.=FALSE)
+ }
+
+ if (is.null(Np)) {
+ if (is.pmcmc)
+ Np <- object at Np
+ else
+ stop("pmcmc error: ",sQuote("Np")," must be specified",call.=FALSE)
+ }
+ Np <- as.integer(Np)
+ if ((length(Np)!=1)||(Np < 1))
+ stop("pmcmc error: ",sQuote("Np")," must be a positive integer",call.=FALSE)
+
+ if (is.null(hyperparams)) {
+ if (is.pmcmc)
+ hyperparams <- object at hyperparams
+ else
+ stop("pmcmc error: ",sQuote("hyperparams")," must be specified",call.=FALSE)
+ }
+
+ Nmcmc <- as.integer(Nmcmc)
+ if (Nmcmc<0)
+ 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")
+ }
+
+ theta <- start
+
+ conv.rec <- matrix(
+ data=NA,
+ nrow=Nmcmc+1,
+ ncol=length(theta)+3,
+ dimnames=list(
+ rownames=seq(from=0,to=Nmcmc,by=1),
+ colnames=c('loglik','log.prior','nfail',names(theta))
+ )
+ )
+
+ if (!all(is.finite(theta[pars]))) {
+ stop(
+ sQuote("pmcmc"),
+ " error: cannot estimate non-finite parameters: ",
+ paste(
+ pars[!is.finite(theta[pars])],
+ collapse=","
+ ),
+ call.=FALSE
+ )
+ }
+
+ obj <- new(
+ "pmcmc",
+ as(object,"pomp"),
+ pars=pars,
+ Nmcmc=0L,
+ dprior=dprior.fun,
+ Np=Np,
+ hyperparams=hyperparams,
+ random.walk.sd=rw.sd
+ )
+
+ if (.ndone==0) { ## compute prior and likelihood on initial parameter vector
+ x <- try(
+ pfilter.internal(
+ object=obj,
+ params=theta,
+ Np=Np,
+ tol=tol,
+ max.fail=max.fail,
+ pred.mean=FALSE,
+ pred.var=FALSE,
+ filter.mean=TRUE,
+ save.states=FALSE,
+ verbose=verbose
+ ),
+ silent=FALSE
+ )
+ if (inherits(x,'try-error'))
+ stop("pmcmc error: error in ",sQuote("pfilter"),call.=FALSE)
+ x$log.prior <- dprior(object=obj,params=theta,log=TRUE)
+ } else { ## has been computed previously
+ x <- .prevcomp
+ }
+ conv.rec[1,names(theta)] <- theta
+ conv.rec[1,c(1,2,3)] <- c(x$loglik,x$log.prior,x$nfail)
+
+ for (n in seq_len(Nmcmc)) { # main loop
+
+ theta.prop <- theta
+ theta.prop[pars] <- rnorm(n=length(pars),mean=theta.prop[pars],sd=rw.sd)
+
+ ## run the particle filter on the proposed new parameter values
+ x.prop <- try(
+ pfilter.internal(
+ object=obj,
+ params=theta.prop,
+ Np=Np,
+ tol=tol,
+ max.fail=max.fail,
+ pred.mean=FALSE,
+ pred.var=FALSE,
+ filter.mean=TRUE,
+ save.states=FALSE,
+ verbose=verbose
+ ),
+ silent=FALSE
+ )
+ if (inherits(x.prop,'try-error'))
+ stop("pmcmc error: error in ",sQuote("pfilter"),call.=FALSE)
+ x.prop$log.prior <- dprior(object=obj,params=theta.prop,log=TRUE)
+
+ ## PMCMC update rule (OK because proposal is symmetric)
+ if (runif(1) < exp(x.prop$loglik-x$loglik+x.prop$log.prior-x$log.prior)) {
+ theta <- theta.prop
+ x <- x.prop
+ }
+
+ ## store a record of this iteration
+ conv.rec[n+1,names(theta)] <- theta
+ conv.rec[n+1,c(1,2,3)] <- c(x$loglik,x$log.prior,x$nfail)
+
+ if (verbose) cat("PMCMC iteration ",n," of ",Nmcmc," completed\n")
+
+ }
+
+ coef(obj) <- theta
+
+ if (Nmcmc>0) {
+ obj at Nmcmc <- as.integer(Nmcmc)
+ obj at filter.mean <- x$filter.mean
+ obj at conv.rec <- conv.rec
+ obj at eff.sample.size <- x$eff.sample.size
+ obj at cond.loglik <- x$cond.loglik
+ obj at loglik <- x$loglik
+ obj at log.prior <- x$log.prior
+ }
+
+ obj
+}
+
+setMethod(
+ "pmcmc",
+ "pomp",
+ function (
+ object, Nmcmc = 1,
+ start, pars, rw.sd,
+ dprior, Np, hyperparams,
+ tol = 1e-17, max.fail = 0, verbose = getOption("verbose")
+ )
+ {
+ if (missing(start)) start <- NULL
+ 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(hyperparams))
+ stop("pmcmc 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(
+ "pmcmc error: ",
+ sQuote("dprior"),
+ " must be a function of prototype ",
+ sQuote("dprior(params,hyperparams,log)"),
+ call.=FALSE
+ )
+ }
+
+ pmcmc.internal(
+ object=object,
+ Nmcmc=Nmcmc,
+ start=start,
+ pars=pars,
+ dprior.fun=dprior,
+ rw.sd=rw.sd,
+ Np=Np,
+ hyperparams=hyperparams,
+ tol=tol,
+ max.fail=max.fail,
+ verbose=verbose,
+ .ndone=0
+ )
+ }
+ )
+
+
+setMethod(
+ "pmcmc",
+ "pmcmc",
+ function (object, Nmcmc, ...)
+ {
+ if (missing(Nmcmc)) Nmcmc <- object at Nmcmc
+ pmcmc.internal(object,Nmcmc=Nmcmc,...)
+ }
+ )
+
+setMethod(
+ 'continue',
+ 'pmcmc',
+ function (object, Nmcmc = 1, ...) {
+ ndone <- object at Nmcmc
+ obj <- pmcmc.internal(
+ object=object,
+ Nmcmc=Nmcmc,
+ .ndone=ndone,
+ .prevcomp=list(
+ log.prior=object at log.prior,
+ loglik=object at loglik,
+ nfail=object at conv.rec[ndone+1,"nfail"],
+ filter.mean=object at filter.mean,
+ eff.sample.size=object at eff.sample.size,
+ cond.loglik=object at cond.loglik
+ ),
+ ...
+ )
+ obj at conv.rec <- rbind(
+ object at conv.rec[,colnames(obj at conv.rec)],
+ obj at conv.rec[-1,]
+ )
+ obj at Nmcmc <- as.integer(ndone+Nmcmc)
+ obj
+ }
+ )
Modified: pkg/R/simulate-pomp.R
===================================================================
--- pkg/R/simulate-pomp.R 2010-06-24 15:01:26 UTC (rev 274)
+++ pkg/R/simulate-pomp.R 2010-06-30 19:30:22 UTC (rev 275)
@@ -24,6 +24,9 @@
save.seed <- get('.Random.seed',envir=.GlobalEnv)
set.seed(seed)
}
+ if (!is.numeric(nsim)||(length(nsim)>1)||nsim<1)
+ stop(sQuote("nsim")," must be a positive integer")
+ nsim <- as.integer(nsim)
xstart <- init.state(object,params=params,t0=times[1])
nreps <- npars*nsim # total number of replicates
## we will do the process model simulations with single calls to the user functions
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2010-06-24 15:01:26 UTC (rev 274)
+++ pkg/inst/NEWS 2010-06-30 19:30:22 UTC (rev 275)
@@ -1,5 +1,10 @@
+NEW IN VERSION 0.31-1:
+ o pomp now includes an implementation of the particle Markov chain Monte Carlo algorithm of Andrieu et al. (Andrieu, C., Doucet, A., & Holenstein, R. (2010) Particle Markov chain Monte Carlo methods. J. R. Stat. Soc. B 72:1--33).
+ This algorithm is implemented in the method 'pmcmc': see 'help(pmcmc)' for details.
+ Implementation primarily due to Ed Ionides.
+
NEW IN VERSION 0.30-1:
- o pomp now includes an implementation of the approximate Bayesian Sequential Monte Carlo algorithm of Liu & West (Liu, J. & West, M. (2001) Combining Parameter and State Estimation in Simulation-Based Filtering in Doucet, A.; de Freitas, N. & Gordon, N. J. (eds.) Sequential Monte Carlo Methods in Practice, Springer, New York, pp. 197--224).
+ o pomp now includes an implementation of the approximate Bayesian Sequential Monte Carlo algorithm of Liu & West (Liu, J. & West, M. (2001) Combining Parameter and State Estimation in Simulation-Based Filtering. In Doucet, A.; de Freitas, N. & Gordon, N. J. (eds.) Sequential Monte Carlo Methods in Practice, Springer, New York, pp. 197--224).
This algorithm is implemented in the method 'bsmc': see 'help(bsmc)' for details.
Thanks to Matt Ferrari and Michael Lavine for the implementation.
Modified: pkg/inst/doc/pomp.bib
===================================================================
--- pkg/inst/doc/pomp.bib 2010-06-24 15:01:26 UTC (rev 274)
+++ pkg/inst/doc/pomp.bib 2010-06-30 19:30:22 UTC (rev 275)
@@ -1,6 +1,34 @@
% This file was created with JabRef 2.6.
% Encoding: ISO8859_1
+ at ARTICLE{Andrieu2010,
+ author = {Andrieu, Christophe and Doucet, Arnaud and Holenstein, Roman},
+ title = {Particle Markov chain Monte Carlo methods},
+ journal = {Journal of the Royal Statistical Society, Series B},
+ year = {2010},
+ volume = {72},
+ pages = {269--342},
+ number = {3},
+ abstract = {Markov chain Monte Carlo and sequential Monte Carlo methods have emerged
+ as the two main tools to sample from high dimensional probability
+ distributions. Although asymptotic convergence of Markov chain Monte
+ Carlo algorithms is ensured under weak assumptions, the performance
+ of these algorithms is unreliable when the proposal distributions
+ that are used to explore the space are poorly chosen and/or if highly
+ correlated variables are updated independently. We show here how
+ it is possible to build efficient high dimensional proposal distributions
+ by using sequential Monte Carlo methods. This allows us not only
+ to improve over standard Markov chain Monte Carlo schemes but also
+ to make Bayesian inference feasible for a large class of statistical
+ models where this was not previously so. We demonstrate these algorithms
+ on a non-linear state space model and a Levy-driven stochastic volatility
+ model.},
+ doi = {10.1111/j.1467-9868.2009.00736.x},
+ file = {Andrieu2010.pdf:Andrieu2010.pdf:PDF},
+ owner = {kingaa},
+ timestamp = {2010.06.30}
+}
+
@ARTICLE{Arulampalam2002,
author = {Arulampalam, M. S. and Maskell, S. and Gordon, N. and Clapp, T.},
title = {A Tutorial on Particle Filters for Online Nonlinear, Non-{G}aussian
Added: pkg/man/pmcmc-class.Rd
===================================================================
--- pkg/man/pmcmc-class.Rd (rev 0)
+++ pkg/man/pmcmc-class.Rd 2010-06-30 19:30:22 UTC (rev 275)
@@ -0,0 +1,85 @@
+\name{pmcmc-class}
+\docType{class}
+\alias{pmcmc-class}
+\title{The "pmcmc" class}
+\description{
+ The \code{pmcmc} class holds a fitted model and is created by a call to \code{\link{pmcmc}}.
+ See \code{\link{pmcmc}} for usage.
+}
+\section{Objects from the Class}{
+ Objects can be created by calls to the \code{\link{pmcmc}} method on an \code{\link{pomp}} object.
+ Such a call uses the Particle MCMC algorithm to fit the model parameters.
+}
+\section{Slots}{
+ A \code{pmcmc} object is derived from a \code{pomp} object and therefore has all the slots of such an object.
+ See \code{\link{pomp-class}} for details.
+ A full description of slots in a \code{pmcmc} object follows.
+ \describe{
+ \item{pars}{
+ A character vector containing the names of parameters to be estimated using PMCMC.
+ }
+ \item{Np}{
+ Number of particles to use in each filtering operation.
+ }
+ \item{Nmcmc}{
+ Number of PMCMC iterations that have been completed.
+ }
+ \item{dprior}{
+ A function of prototype \code{dprior(params,hyperparams,log)} that evaluates the prior density. This defaults to an improper uniform prior.
+ }
+ \item{hyperparams}{
+ A list passed as an argument to \code{dprior}.
+ }
+ \item{random.walk.sd}{
+ A named vector containing the random-walk variances used to parameterize a Gaussian random walk MCMC proposal.
+ }
+ \item{filter.mean}{
+ Matrix of filtering means.
+ See \code{\link{pfilter}}.
+ }
+ \item{eff.sample.size}{
+ A vector containing the effective number of particles at each time point.
+ See \code{\link{pfilter}}.
+ }
+ \item{cond.loglik}{
+ A vector containing the conditional log likelihoods at each time point.
+ See \code{\link{pfilter}}.
+ }
+ \item{conv.rec}{
+ The \dQuote{convergence record}: a matrix containing a record of the parameter values, log likelihoods, log prior probabilities, and other pertinent information, with one row for each PMCMC iteration.
+ }
+ \item{loglik}{
+ A numeric value containing the value of the log likelihood, as evaluated for the random-parameter model.
+ Note that this will not be equal to the log likelihood for the fixed-parameter model.
+ }
+ \item{log.prior}{
+ A numeric value containing the log of the prior density evaluated at the parameter vector in the params slot.
+ }
+ \item{data, times, t0, rprocess, dprocess, dmeasure, rmeasure,
+ skeleton.type, skeleton, initializer, states, params,
+ statenames, paramnames, covarnames, tcovar, covar,
+ PACKAGE, userdata}{
+ Inherited from the \code{pomp} class.
+ }
+ }
+}
+\section{Extends}{
+ Class \code{pomp}, directly.
+ See \code{\link{pomp-class}}.
+}
+\section{Methods}{
+ See \code{\link{pmcmc}}, \link{pmcmc-methods}.
+}
+\references{
+ C. Andrieu, A. Doucet and R. Holenstein,
+ Particle Markov chain Monte Carlo methods,
+ J. Roy. Stat. Soc B, to appear, 2010.
+
+ C. Andrieu and G.O. Roberts,
+ The pseudo-marginal approach for efficient computation,
+ Ann Stat 37:697-725, 2009.
+}
+\author{Edward L. Ionides \email{ionides at umich dot edu}, Aaron A. King \email{kingaa at umich dot edu}}
+\seealso{\code{\link{pmcmc}}, \link{pmcmc-methods}, \code{\link{pomp}}, \link{pomp-class}}
+\keyword{models}
+\keyword{ts}
Added: pkg/man/pmcmc-methods.Rd
===================================================================
--- pkg/man/pmcmc-methods.Rd (rev 0)
+++ pkg/man/pmcmc-methods.Rd 2010-06-30 19:30:22 UTC (rev 275)
@@ -0,0 +1,89 @@
+\name{pmcmc-methods}
+\docType{methods}
+\alias{pmcmc-methods}
+\alias{dprior,pmcmc-method}
+\alias{dprior-pmcmc}
+\alias{logLik,pmcmc-method}
+\alias{logLik-pmcmc}
+\alias{conv.rec,pmcmc-method}
+\alias{conv.rec-pmcmc}
+\alias{filter.mean,pmcmc-method}
+\alias{filter.mean-pmcmc}
+\alias{plot-pmcmc}
+\alias{plot,pmcmc-method}
+\alias{compare.pmcmc}
+\alias{dprior}
+\title{Methods of the "pmcmc" class}
+\description{Methods of the "pmcmc" class.}
+\usage{
+\S4method{logLik}{pmcmc}(object, \dots)
+\S4method{conv.rec}{pmcmc}(object, pars, \dots)
+\S4method{filter.mean}{pmcmc}(object, pars, \dots)
+\S4method{plot}{pmcmc}(x, y = NULL, \dots)
+\S4method{dprior}{pmcmc}(object, params, log)
+compare.pmcmc(z)
+}
+\arguments{
+ \item{object, x}{The \code{pmcmc} object.}
+ \item{pars}{Names of parameters.}
+ \item{y}{Ignored.}
+ \item{z}{A \code{pmcmc} object or list of \code{pmcmc} objects.}
+ \item{params}{
+ 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{logLik}{
+ Returns the value in the \code{loglik} slot.
+ }
+ \item{dprior}{
+ \code{dprior(object,params,log)} evaluates the prior density at \code{params} with values of the hyperparameters given by \code{object at hyperparams}.
+ }
+ \item{pmcmc}{
+ Re-runs the PMCMC iterations.
+ See the documentation for \code{\link{pmcmc}}.
+ }
+ \item{compare.pmcmc}{
+ Given a \code{pmcmc} object or a list of \code{pmcmc} objects, \code{compare.pmcmc} produces a set of diagnostic plots.
+ }
+ \item{plot}{
+ Plots a series of diagnostic plots.
+ When \code{x} is a \code{pmcmc} object, \code{plot(x)} is equivalent to \code{compare.pmcmc(list(x))}.
+ }
+ \item{filter.mean}{
+ \code{filter.mean(object, pars = NULL)} returns the rows of the filtering-mean matrix corresponding to the names in \code{pars}.
+ By default, all rows are returned.
+ }
+ \item{print}{
+ Prints a summary of the \code{pmcmc} object.
+ }
+ \item{show}{
+ Displays the \code{pmcmc} object.
+ }
+ \item{pfilter}{
+ See \code{\link{pfilter}}.
+ }
+ }
+}
+\references{
+ C. Andrieu, A. Doucet and R. Holenstein,
+ Particle Markov chain Monte Carlo methods,
+ J. Roy. Stat. Soc B, to appear, 2010.
+
+ C. Andrieu and G.O. Roberts,
+ The pseudo-marginal approach for efficient computation,
+ Ann Stat 37:697-725, 2009.
+}
+\author{Edward L. Ionides \email{ionides at umich dot edu}, Aaron A. King \email{kingaa at umich dot edu}}
+\seealso{\code{\link{pmcmc}}, \code{\link{pomp}}, \code{\link{pomp-class}}, \code{\link{pfilter}}}
+\keyword{models}
+\keyword{ts}
Added: pkg/man/pmcmc.Rd
===================================================================
--- pkg/man/pmcmc.Rd (rev 0)
+++ pkg/man/pmcmc.Rd 2010-06-30 19:30:22 UTC (rev 275)
@@ -0,0 +1,102 @@
+\name{pmcmc}
+\docType{methods}
+\alias{pmcmc}
+\alias{pmcmc,pmcmc-method}
+\alias{pmcmc-pmcmc}
+\alias{pmcmc,pomp-method}
+\alias{pmcmc-pomp}
+\alias{continue,pmcmc-method}
+\alias{continue-pmcmc}
+\title{The PMCMC algorithm}
+\description{
+ The Particle MCMC algorithm for estimating the parameters of a partially-observed Markov process.
+}
+\usage{
+pmcmc(object, \dots)
+\S4method{pmcmc}{pomp}(object, Nmcmc = 1, start, pars,
+ rw.sd, dprior, Np, hyperparams, tol = 1e-17, max.fail = 0,
+ verbose = getOption("verbose"))
+\S4method{pmcmc}{pmcmc}(object, Nmcmc, \dots)
+\S4method{continue}{pmcmc}(object, Nmcmc = 1, \dots)
+}
+\arguments{
+ \item{object}{
+ An object of class \code{pomp}.
+ }
+ \item{Nmcmc}{
+ The number of PMCMC 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{dprior}{
+ Function of prototype \code{dprior(params,hyperparams,...,log)} that evaluates the prior density.
+ This defaults to an improper uniform prior.
+ }
+ \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{Np}{
+ a positive integer;
+ the number of particles to use in each filtering operation.
+ }
+ \item{hyperparams}{
+ optional list; parameters to be passed to \code{dprior}.
+ }
+ \item{tol}{
+ numeric scalar; particles with log likelihood below \code{tol} are considered to be \dQuote{lost}.
+ A filtering failure occurs when, at some time point, all particles are lost.
+ }
+ \item{max.fail}{
+ integer; maximum number of filtering failures permitted.
+ If the number of failures exceeds this number, execution will terminate with an error.
+ }
+ \item{verbose}{
+ logical; if TRUE, print progress reports.
+ }
+ \item{\dots}{
+ Additional arguments that can be used to override the defaults.
+ }
+}
+\section{Re-running PMCMC Iterations}{
+ To re-run a sequence of PMCMC iterations, one can use the \code{pmcmc} method on a \code{pmcmc} object.
+ By default, the same parameters used for the original PMCMC 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 PMCMC Iterations}{
+ One can continue a series of PMCMC iterations from where one left off using the \code{continue} method.
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 275
More information about the pomp-commits
mailing list