[Pomp-commits] r1008 - in pkg/pomp: . R inst man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 23 21:22:56 CEST 2014
Author: kingaa
Date: 2014-09-23 21:22:56 +0200 (Tue, 23 Sep 2014)
New Revision: 1008
Added:
pkg/pomp/R/bsmc2.R
pkg/pomp/man/bsmc2.Rd
pkg/pomp/tests/ou2-bsmc2.R
pkg/pomp/tests/ou2-bsmc2.Rout.save
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/NAMESPACE
pkg/pomp/R/bsmc.R
pkg/pomp/R/generics.R
pkg/pomp/inst/NEWS
pkg/pomp/inst/NEWS.Rd
pkg/pomp/man/bsmc.Rd
pkg/pomp/tests/bbs.R
Log:
- added modified version of Liu & West algorithm, 'bsmc'. Seems better all-round. Some testing still to do.
- small modifications (of no effect) to 'bsmc.R'
- mv 'methods' back into 'Depends' field
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2014-09-03 11:06:32 UTC (rev 1007)
+++ pkg/pomp/DESCRIPTION 2014-09-23 19:22:56 UTC (rev 1008)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.53-6
-Date: 2014-09-02
+Version: 0.54-1
+Date: 2014-09-23
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")),
@@ -18,8 +18,8 @@
)
URL: http://pomp.r-forge.r-project.org
Description: Inference methods for partially-observed Markov processes
-Depends: R(>= 3.0.0), subplex, nloptr
-Imports: stats, graphics, methods, mvtnorm, deSolve, coda
+Depends: R(>= 3.0.0), methods, subplex, nloptr
+Imports: stats, graphics, mvtnorm, deSolve, coda
License: GPL(>= 2)
LazyData: 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.
@@ -32,7 +32,8 @@
dmeasure-pomp.R dprocess-pomp.R skeleton-pomp.R
dprior-pomp.R rprior-pomp.R
simulate-pomp.R trajectory-pomp.R plot-pomp.R
- pfilter.R pfilter-methods.R minim.R traj-match.R bsmc.R
+ pfilter.R pfilter-methods.R minim.R traj-match.R
+ bsmc.R bsmc2.R
mif.R mif-methods.R
pmcmc.R pmcmc-methods.R
nlf-funcs.R nlf-guts.R nlf-objfun.R nlf.R
Modified: pkg/pomp/NAMESPACE
===================================================================
--- pkg/pomp/NAMESPACE 2014-09-03 11:06:32 UTC (rev 1007)
+++ pkg/pomp/NAMESPACE 2014-09-23 19:22:56 UTC (rev 1008)
@@ -64,7 +64,7 @@
eff.sample.size,cond.logLik,
particles,mif,continue,states,trajectory,
pred.mean,pred.var,filter.mean,conv.rec,
- bsmc,pmcmc,abc,nlf,
+ bsmc2,bsmc,pmcmc,abc,nlf,
traj.match.objfun,
probe.match.objfun,
spect,probe,probe.match,
Modified: pkg/pomp/R/bsmc.R
===================================================================
--- pkg/pomp/R/bsmc.R 2014-09-03 11:06:32 UTC (rev 1007)
+++ pkg/pomp/R/bsmc.R 2014-09-23 19:22:56 UTC (rev 1008)
@@ -24,8 +24,7 @@
seed="integer",
nfail="integer",
cond.log.evidence="numeric",
- log.evidence="numeric",
- weights="numeric"
+ log.evidence="numeric"
)
)
@@ -76,18 +75,6 @@
ptsi.inv <- FALSE
ntimes <- length(time(object))
- if (is.null(dim(params))) {
- params <- matrix(
- params,
- nrow=length(params),
- ncol=Np,
- dimnames=list(
- names(params),
- NULL
- )
- )
- }
-
npars <- nrow(params)
paramnames <- rownames(params)
prior <- params
@@ -360,8 +347,7 @@
seed=as.integer(seed),
nfail=as.integer(nfail),
cond.log.evidence=evidence,
- log.evidence=sum(evidence),
- weights=weights
+ log.evidence=sum(evidence)
)
}
@@ -443,7 +429,7 @@
function (x, y, ..., pars, thin) {
if (!missing(y))
warning(sQuote("y")," is ignored")
- if (missing(pars)) pars <- names(coef(x,transform=!x at transform))
+ if (missing(pars)) pars <- x at est
if (missing(thin)) thin <- Inf
bsmc.plot(
prior=if (x at transform) partrans(x,x at prior,dir="forward") else x at prior,
Added: pkg/pomp/R/bsmc2.R
===================================================================
--- pkg/pomp/R/bsmc2.R (rev 0)
+++ pkg/pomp/R/bsmc2.R 2014-09-23 19:22:56 UTC (rev 1008)
@@ -0,0 +1,257 @@
+## Bayesian particle filtering codes
+##
+## in annotation L&W AGM == Liu & West "A General Algorithm"
+##
+## params = the initial particles for the parameter values;
+## these should be drawn from the prior distribution for the parameters
+## est = names of parameters to estimate; other parameters are not updated.
+## smooth = parameter 'h' from AGM
+## ntries = number of samplesto draw from x_{t+1} | x(k)_{t} to estimate
+## mean of mu(k)_t+1 as in sect 2.2 Liu & West
+
+bsmc2.internal <- function (object, params, Np, est,
+ smooth, tol, seed = NULL,
+ verbose = getOption("verbose"),
+ max.fail, transform, .getnativesymbolinfo = TRUE,
+ ...) {
+
+ gnsi.rproc <- gnsi.dmeas <- as.logical(.getnativesymbolinfo)
+ ptsi.inv <- ptsi.for <- TRUE
+ transform <- as.logical(transform)
+
+ if (missing(seed)) seed <- NULL
+ if (!is.null(seed)) {
+ if (!exists(".Random.seed",where=.GlobalEnv))
+ runif(n=1L) ## need to initialize the RNG
+ save.seed <- get(".Random.seed",pos=.GlobalEnv)
+ set.seed(seed)
+ }
+
+ error.prefix <- paste(sQuote("bsmc2"),"error: ")
+
+ if (missing(params)) {
+ if (length(coef(object))>0) {
+ params <- coef(object)
+ } else {
+ stop(error.prefix,sQuote("params")," must be supplied",call.=FALSE)
+ }
+ }
+
+ if (missing(Np)) Np <- NCOL(params)
+ else if (is.matrix(params)&&(Np!=ncol(params)))
+ warning(sQuote("Np")," is ignored when ",sQuote("params")," is a matrix")
+
+ if ((!is.matrix(params)) && (Np > 1))
+ params <- rprior(object,params=parmat(params,Np))
+
+ if (transform)
+ params <- partrans(object,params,dir="inverse",
+ .getnativesymbolinfo=ptsi.inv)
+ ptsi.inv <- FALSE
+
+ ntimes <- length(time(object))
+ npars <- nrow(params)
+ paramnames <- rownames(params)
+ prior <- params
+
+ if (missing(est))
+ est <- paramnames[apply(params,1,function(x)diff(range(x))>0)]
+ estind <- match(est,paramnames)
+ npars.est <- length(estind)
+
+ if (npars.est<1)
+ stop(error.prefix,"no parameters to estimate",call.=FALSE)
+
+ if (is.null(paramnames))
+ stop(error.prefix,sQuote("params")," must have rownames",call.=FALSE)
+
+ if ((length(smooth)!=1)||(smooth>1)||(smooth<=0))
+ stop(error.prefix,sQuote("smooth")," must be a scalar in [0,1)",call.=FALSE)
+
+ hsq <- smooth^2 # see Liu & West eq(10.3.12)
+ shrink <- sqrt(1-hsq) # 'a' parameter of Liu & West
+
+ xstart <- init.state(
+ object,
+ params=if (transform) {
+ partrans(object,params,dir="forward",
+ .getnativesymbolinfo=ptsi.for)
+ } else {
+ params
+ }
+ )
+ statenames <- rownames(xstart)
+ nvars <- nrow(xstart)
+ ptsi.for <- FALSE
+
+ times <- time(object,t0=TRUE)
+ x <- xstart
+
+ evidence <- as.numeric(rep(NA,ntimes))
+ eff.sample.size <- as.numeric(rep(NA,ntimes))
+ nfail <- 0L
+
+ mu <- array(data=NA,dim=c(nvars,Np,1))
+ rownames(mu) <- rownames(xstart)
+ m <- array(data=NA,dim=c(npars,Np))
+ rownames(m) <- rownames(params)
+
+ for (nt in seq_len(ntimes)) {
+
+ ## calculate particle means ; as per L&W AGM (1)
+ params.mean <- apply(params,1,mean)
+ ## calculate particle covariances : as per L&W AGM (1)
+ params.var <- cov(t(params[estind,,drop=FALSE]))
+
+ if (verbose) {
+ cat("at step",nt,"(time =",times[nt+1],")\n")
+ print(
+ rbind(
+ prior.mean=params.mean[estind],
+ prior.sd=sqrt(diag(params.var))
+ )
+ )
+ }
+
+ m <- shrink*params+(1-shrink)*params.mean
+
+ ## sample new parameter vector as per L&W AGM (3) and Liu & West eq(3.2)
+ pert <- try(
+ mvtnorm::rmvnorm(
+ n=Np,
+ mean=rep(0,npars.est),
+ sigma=hsq*params.var,
+ method="svd"
+ ),
+ silent=FALSE
+ )
+ if (inherits(pert,"try-error"))
+ stop(error.prefix,"error in ",sQuote("rmvnorm"),call.=FALSE)
+ if (any(!is.finite(pert)))
+ stop(error.prefix,"extreme particle depletion",call.=FALSE)
+ params[estind,] <- m[estind,]+t(pert)
+
+ if (transform)
+ tparams <- partrans(object,params,dir="forward",
+ .getnativesymbolinfo=ptsi.for)
+
+ xpred <- rprocess(
+ object,
+ xstart=x,
+ times=times[c(nt,nt+1)],
+ params=if (transform) {
+ tparams
+ } else {
+ params
+ },
+ offset=1,
+ .getnativesymbolinfo=gnsi.rproc
+ )
+ gnsi.rproc <- FALSE
+
+ ## evaluate likelihood of observation given xpred (from L&W AGM (4))
+ weights <- dmeasure(
+ object,
+ y=object at data[,nt,drop=FALSE],
+ x=xpred,
+ times=times[nt+1],
+ params=if (transform) {
+ tparams
+ } else {
+ params
+ },
+ .getnativesymbolinfo=gnsi.dmeas
+ )
+ gnsi.dmeas <- FALSE
+ ## evaluate weights as per L&W AGM (5)
+
+ storeForEvidence <- log(mean(weights))
+
+ x[,] <- xpred
+
+ ## test for failure to filter
+ dim(weights) <- NULL ### needed?
+ failures <- ((weights<tol)|(!is.finite(weights))) # test for NA weights
+ all.fail <- all(failures)
+ if (all.fail) { # all particles are lost
+ if (verbose) {
+ message("filtering failure at time t = ",times[nt+1])
+ }
+ nfail <- nfail+1
+ if (nfail > max.fail)
+ stop(error.prefix,"too many filtering failures",call.=FALSE)
+ evidence[nt] <- log(tol) # worst log-likelihood
+ weights <- rep(1/Np,Np)
+ eff.sample.size[nt] <- 0
+ } else { # not all particles are lost
+ ## compute log-likelihood
+ evidence[nt] <- storeForEvidence
+ weights[failures] <- 0
+ weights <- weights/sum(weights)
+ ## compute effective sample-size
+ eff.sample.size[nt] <- 1/crossprod(weights)
+ }
+
+ if (verbose) {
+ cat("effective sample size =",round(eff.sample.size[nt],1),"\n")
+ }
+
+ ## Matrix with samples (columns) from filtering distribution theta.t | Y.t
+ if (!all.fail) {
+ smp <- .Call(systematic_resampling,weights)
+### smp <- sample.int(n=Np,size=Np,replace=TRUE,prob=weights)
+ x <- x[,smp,drop=FALSE]
+ params[estind,] <- params[estind,smp,drop=FALSE]
+ }
+
+ .getnativesymbolinfo <- FALSE
+
+ }
+
+ if (!is.null(seed)) {
+ assign(".Random.seed",save.seed,pos=.GlobalEnv)
+ seed <- save.seed
+ }
+
+ ## replace parameters with point estimate (posterior median)
+ coef(object,transform=transform) <- apply(params,1,median)
+
+ new(
+ "bsmcd.pomp",
+ object,
+ transform=transform,
+ post=params,
+ prior=prior,
+ est=as.character(est),
+ eff.sample.size=eff.sample.size,
+ smooth=smooth,
+ seed=as.integer(seed),
+ nfail=as.integer(nfail),
+ cond.log.evidence=evidence,
+ log.evidence=sum(evidence)
+ )
+}
+
+setMethod(
+ "bsmc2",
+ signature=signature(object="pomp"),
+ definition = function (object, params, Np, est,
+ smooth = 0.1, tol = 1e-17, seed = NULL,
+ verbose = getOption("verbose"),
+ max.fail = 0, transform = FALSE,
+ ...) {
+ bsmc2.internal(
+ object=object,
+ params=params,
+ Np=Np,
+ est=est,
+ smooth=smooth,
+ tol=tol,
+ seed=seed,
+ verbose=verbose,
+ max.fail=max.fail,
+ transform=transform,
+ ...
+ )
+ }
+ )
Modified: pkg/pomp/R/generics.R
===================================================================
--- pkg/pomp/R/generics.R 2014-09-03 11:06:32 UTC (rev 1007)
+++ pkg/pomp/R/generics.R 2014-09-23 19:22:56 UTC (rev 1008)
@@ -59,6 +59,7 @@
## Bayesian SMC (Liu & West)
setGeneric("bsmc",function(object,...)standardGeneric("bsmc"))
+setGeneric("bsmc2",function(object,...)standardGeneric("bsmc2"))
## basic SMC (particle filter)
setGeneric("pfilter",function(object,...)standardGeneric("pfilter"))
Modified: pkg/pomp/inst/NEWS
===================================================================
--- pkg/pomp/inst/NEWS 2014-09-03 11:06:32 UTC (rev 1007)
+++ pkg/pomp/inst/NEWS 2014-09-23 19:22:56 UTC (rev 1008)
@@ -1,5 +1,13 @@
_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._5_4-_1:
+
+ • A modified version of the Liu and West (2001) algorithm is
+ included as ‘bsmc2’.
+
+ • By default, when ‘B’ is of class ‘bsmcd.pomp’, ‘plot(B)’
+ ignores fixed parameters.
+
_C_h_a_n_g_e_s _i_n '_p_o_m_p' _v_e_r_s_i_o_n _0._5_3-_6:
• A bug having to do with paths to the temporary C files,
Modified: pkg/pomp/inst/NEWS.Rd
===================================================================
--- pkg/pomp/inst/NEWS.Rd 2014-09-03 11:06:32 UTC (rev 1007)
+++ pkg/pomp/inst/NEWS.Rd 2014-09-23 19:22:56 UTC (rev 1008)
@@ -1,5 +1,11 @@
\name{NEWS}
\title{News for package `pomp'}
+\section{Changes in \pkg{pomp} version 0.54-1}{
+ \itemize{
+ \item A modified version of the Liu and West (2001) algorithm is included as \code{bsmc2}.
+ \item By default, when \code{B} is of class \code{bsmcd.pomp}, \code{plot(B)} ignores fixed parameters.
+ }
+}
\section{Changes in \pkg{pomp} version 0.53-6}{
\itemize{
\item A bug having to do with paths to the temporary C files, encountered when using \code{Csnippet}s under windows, has been fixed.
Modified: pkg/pomp/man/bsmc.Rd
===================================================================
--- pkg/pomp/man/bsmc.Rd 2014-09-03 11:06:32 UTC (rev 1007)
+++ pkg/pomp/man/bsmc.Rd 2014-09-23 19:22:56 UTC (rev 1008)
@@ -2,10 +2,6 @@
\alias{bsmc}
\alias{bsmc-pomp}
\alias{bsmc,pomp-method}
-\alias{$,bsmcd.pomp-method}
-\alias{$-bsmcd.pomp}
-\alias{plot,bsmcd.pomp-method}
-\alias{plot-bsmcd.pomp}
\title{Liu and West Bayesian Particle Filter}
\description{
Generates draws from the posterior distribution for the parameters using the Liu and West algorithm.
Copied: pkg/pomp/man/bsmc2.Rd (from rev 1007, pkg/pomp/man/bsmc.Rd)
===================================================================
--- pkg/pomp/man/bsmc2.Rd (rev 0)
+++ pkg/pomp/man/bsmc2.Rd 2014-09-23 19:22:56 UTC (rev 1008)
@@ -0,0 +1,119 @@
+\name{bsmc2}
+\alias{bsmc2}
+\alias{bsmc2-pomp}
+\alias{bsmc2,pomp-method}
+\alias{$,bsmcd.pomp-method}
+\alias{$-bsmcd.pomp}
+\alias{plot,bsmcd.pomp-method}
+\alias{plot-bsmcd.pomp}
+\title{Modified Liu and West Bayesian Particle Filter}
+\description{
+ Generates draws from the posterior distribution for the parameters using a plug-and-play variant of the Liu and West algorithm.
+}
+\usage{
+\S4method{bsmc2}{pomp}(object, params, Np, est, smooth = 0.1,
+ tol = 1e-17, seed = NULL,
+ verbose = getOption("verbose"), max.fail = 0,
+ transform = FALSE, \dots)
+}
+\arguments{
+\item{object}{
+ An object of class \code{pomp} or inheriting class \code{pomp}.
+ }
+ \item{params, Np}{
+ Specifications for the prior distribution of particles.
+ See details below.
+ }
+ \item{est}{
+ Names of the rows of \code{params} that are to be estimated.
+ No updates will be made to the other parameters.
+ If \code{est} is not specified, all parameters for which there is variation in \code{params} will be estimated.
+ }
+ \item{smooth}{
+ Kernel density smoothing parameters.
+ The compensating shrinkage factor will be \code{sqrt(1-smooth^2)}.
+ Thus, \code{smooth=0} means that no noise will be added to parameters.
+ Generally, the value of \code{smooth} should be chosen close to 0 (i.e., \code{shrink~0.1}).
+ }
+ \item{tol}{
+ 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.
+ When all particles are lost, the conditional log likelihood at that time point is set to be \code{log(tol)}.
+ }
+ \item{seed}{
+ optional; an object specifying if and how the random number generator should be initialized (\sQuote{seeded}).
+ If \code{seed} is an integer, it is passed to \code{set.seed} prior to any simulation and is returned as the \dQuote{seed} element of the return list.
+ By default, the state of the random number generator is not changed and the value of \code{.Random.seed} on the call is stored in the \dQuote{seed} element of the return list.
+ }
+ \item{verbose}{
+ logical; if \code{TRUE}, print diagnostic messages.
+ }
+ \item{max.fail}{
+ The maximum number of filtering failures allowed.
+ If the number of filtering failures exceeds this number, execution will terminate with an error.
+ }
+ \item{transform}{
+ logical;
+ if \code{TRUE}, the algorithm operates on the transformed scale.
+ }
+ \item{\dots}{
+ currently ignored.
+ }
+}
+\value{
+ An object of class \dQuote{bsmcd.pomp}.
+ The \dQuote{params} slot of this object will hold the parameter posterior medians.
+ The slots of this class include:
+ \item{post}{
+ A matrix containing draws from the approximate posterior distribution.
+ }
+ \item{prior}{
+ A matrix containing draws from the prior distribution (identical to \code{params} on call).
+ }
+ \item{eff.sample.size}{
+ A vector containing the effective number of particles at each time point.
+ }
+ \item{smooth}{
+ The smoothing parameter used (see above).
+ }
+ \item{seed}{
+ The state of the random number generator at the time \code{bsmc} was called.
+ If the argument \code{seed} was specified, this is a copy;
+ if not, this is the internal state of the random number generator at the time of call.
+ }
+ \item{nfail}{
+ The number of filtering failures encountered.
+ }
+ \item{cond.log.evidence}{
+ A vector containing the conditional log evidence scores at each time point.
+ }
+ \item{log.evidence}{
+ The estimated log evidence.
+ }
+ \item{weights}{
+ The resampling weights for each particle.
+ }
+}
+\details{
+ There are two ways to specify the prior distribution of particles.
+ If \code{params} is unspecified or is a named vector, \code{Np} draws are made from the prior distribution, as specified by \code{\link{rprior}}.
+ Alternatively, \code{params} can be specified as an \code{npars} x \code{Np} matrix (with rownames).
+}
+\author{
+ Michael Lavine (lavine at math dot umass dot edu),
+ Matthew Ferrari (mferrari at psu dot edu),
+ Aaron A. King
+ Edward L. Ionides
+}
+\references{
+ Liu, J. and M. West.
+ Combining Parameter and State Estimation in Simulation-Based Filtering.
+ In A. Doucet, N. de Freitas, and N. J. Gordon, editors,
+ Sequential Monte Carlo Methods in Practice, pages 197-224.
+ Springer, New York, 2001.
+}
+\examples{
+## See the "Introducton to pomp" document for examples.
+}
+\seealso{\link{pomp-class}}
+\keyword{ts}
Modified: pkg/pomp/tests/bbs.R
===================================================================
--- pkg/pomp/tests/bbs.R 2014-09-03 11:06:32 UTC (rev 1007)
+++ pkg/pomp/tests/bbs.R 2014-09-23 19:22:56 UTC (rev 1008)
@@ -15,7 +15,8 @@
}
)
-fit1 <- bsmc(bbs,params=coef(bbs),Np=1000,transform=TRUE,est=c("beta","sigma"),smooth=0.2)
+fit1 <- bsmc2(bbs,params=coef(bbs),Np=5000,transform=TRUE,
+ est=c("beta","sigma"),smooth=0.2)
signif(coef(fit1),3)
fit2 <- traj.match(bbs,est=c("beta","sigma"),transform=TRUE)
Added: pkg/pomp/tests/ou2-bsmc2.R
===================================================================
--- pkg/pomp/tests/ou2-bsmc2.R (rev 0)
+++ pkg/pomp/tests/ou2-bsmc2.R 2014-09-23 19:22:56 UTC (rev 1008)
@@ -0,0 +1,59 @@
+library(pomp)
+
+set.seed(398585L)
+pompExample(ou2)
+
+time(ou2) <- 1:10
+
+Np <- 50000
+
+prior.bounds <- rbind(
+ alpha.2=c(-0.55,-0.45),
+ alpha.3=c(0.25,0.35)
+ )
+colnames(prior.bounds) <- c("lower","upper")
+
+estnames <- rownames(prior.bounds)
+
+prior <- matrix(data=coef(ou2),nrow=length(coef(ou2)),ncol=Np)
+rownames(prior) <- names(coef(ou2))
+for (n in estnames) {
+ prior[n,] <- runif(n=Np,min=prior.bounds[n,1],max=prior.bounds[n,2])
+}
+
+##Run Liu & West particle filter
+tic <- Sys.time()
+smc <- bsmc2(
+ ou2,
+ params=prior,
+ est=estnames,
+ smooth=0.02
+ )
+toc <- Sys.time()
+
+prior <- smc$prior
+post <- smc$post
+
+print(etime <- toc-tic)
+
+print(
+ cbind(
+ prior.mean=apply(prior,1,mean),
+ posterior.mean=apply(post,1,mean),
+ truth=coef(ou2),
+ t(apply(post,1,quantile,c(0.025,0.5,0.975)))
+ )
+ )
+
+print(min(smc$eff.sample.size))
+print(smc$log.evidence)
+
+ou2 <- pomp(ou2,
+ rprior=function(params,...){
+ params
+ }
+ )
+
+smc <- bsmc2(ou2,Np=25000,smooth=0.1,est=estnames,seed=648651945L)
+print(smc$eff.sample.size)
+print(smc$log.evidence)
Added: pkg/pomp/tests/ou2-bsmc2.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-bsmc2.Rout.save (rev 0)
+++ pkg/pomp/tests/ou2-bsmc2.Rout.save 2014-09-23 19:22:56 UTC (rev 1008)
@@ -0,0 +1,104 @@
+
+R version 3.1.1 (2014-07-10) -- "Sock it to Me"
+Copyright (C) 2014 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.
+
+> library(pomp)
+Loading required package: subplex
+Loading required package: nloptr
+>
+> set.seed(398585L)
+> pompExample(ou2)
+newly created pomp object(s):
+ ou2
+>
+> time(ou2) <- 1:10
+>
+> Np <- 10000
+>
+> prior.bounds <- rbind(
++ alpha.2=c(-0.55,-0.45),
++ alpha.3=c(0.25,0.35)
++ )
+> colnames(prior.bounds) <- c("lower","upper")
+>
+> estnames <- rownames(prior.bounds)
+>
+> prior <- matrix(data=coef(ou2),nrow=length(coef(ou2)),ncol=Np)
+> rownames(prior) <- names(coef(ou2))
+> for (n in estnames) {
++ prior[n,] <- runif(n=Np,min=prior.bounds[n,1],max=prior.bounds[n,2])
++ }
+>
+> ##Run Liu & West particle filter
+> tic <- Sys.time()
+> smc <- bsmc(
++ ou2,
++ params=prior,
++ est=estnames,
++ ntries=5,
++ smooth=0.02,
++ lower=prior.bounds[estnames,"lower"],
++ upper=prior.bounds[estnames,"upper"]
++ )
+> toc <- Sys.time()
+>
+> prior <- smc$prior
+> post <- smc$post
+>
+> print(etime <- toc-tic)
+Time difference of 2.931616 secs
+>
+> print(
++ cbind(
++ prior.mean=apply(prior,1,mean),
++ posterior.mean=apply(post,1,mean),
++ truth=coef(ou2),
++ t(apply(post,1,quantile,c(0.025,0.5,0.975)))
++ )
++ )
+ prior.mean posterior.mean truth 2.5% 50% 97.5%
+alpha.1 0.8000000 0.8000000 0.8 0.8000000 0.8000000 0.8000000
+alpha.2 -0.4999287 -0.5105023 -0.5 -0.5402483 -0.4993459 -0.4536930
+alpha.3 0.2996065 0.3148637 0.3 0.2823821 0.3260754 0.3388949
+alpha.4 0.9000000 0.9000000 0.9 0.9000000 0.9000000 0.9000000
+sigma.1 3.0000000 3.0000000 3.0 3.0000000 3.0000000 3.0000000
+sigma.2 -0.5000000 -0.5000000 -0.5 -0.5000000 -0.5000000 -0.5000000
+sigma.3 2.0000000 2.0000000 2.0 2.0000000 2.0000000 2.0000000
+tau 1.0000000 1.0000000 1.0 1.0000000 1.0000000 1.0000000
+x1.0 -3.0000000 -3.0000000 -3.0 -3.0000000 -3.0000000 -3.0000000
+x2.0 4.0000000 4.0000000 4.0 4.0000000 4.0000000 4.0000000
+>
+> print(min(smc$eff.sample.size))
+[1] 22.94863
+> print(smc$log.evidence)
+[1] 45.47584
+>
+> ou2 <- pomp(ou2,
++ rprior=function(params,...){
++ params
++ }
++ )
+>
+> smc <- bsmc(ou2,ntries=5,Np=5000,smooth=0.1,est=estnames,seed=648651945L)
+> print(smc$eff.sample.size)
+ [1] 186.40437 36.29100 57.56951 29.30424 180.23722 34.63366 156.94264
+ [8] 24.49006 178.39269 125.05970
+> print(smc$log.evidence)
+[1] 37.68127
+>
+> proc.time()
+ user system elapsed
+ 4.908 0.068 5.011
More information about the pomp-commits
mailing list