[Pomp-commits] r830 - in branches/mif2: . R inst man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 28 06:48:57 CET 2013
Author: nxdao2000
Date: 2013-02-28 06:48:56 +0100 (Thu, 28 Feb 2013)
New Revision: 830
Modified:
branches/mif2/DESCRIPTION
branches/mif2/NAMESPACE
branches/mif2/R/aaa.R
branches/mif2/R/basic-probes.R
branches/mif2/R/bsmc.R
branches/mif2/R/dmeasure-pomp.R
branches/mif2/R/dprocess-pomp.R
branches/mif2/R/init-state-pomp.R
branches/mif2/R/mif-class.R
branches/mif2/R/mif-methods.R
branches/mif2/R/particles-mif.R
branches/mif2/R/plot-pomp.R
branches/mif2/R/plugins.R
branches/mif2/R/pmcmc.R
branches/mif2/R/pomp-fun.R
branches/mif2/R/pomp-methods.R
branches/mif2/R/pomp.R
branches/mif2/R/probe.R
branches/mif2/R/rmeasure-pomp.R
branches/mif2/R/rprocess-pomp.R
branches/mif2/R/simulate-pomp.R
branches/mif2/R/skeleton-pomp.R
branches/mif2/R/sobol.R
branches/mif2/R/spect-match.R
branches/mif2/R/spect.R
branches/mif2/R/traj-match.R
branches/mif2/R/trajectory-pomp.R
branches/mif2/inst/NEWS
branches/mif2/man/mif-class.Rd
branches/mif2/man/mif.Rd
branches/mif2/man/pmcmc.Rd
branches/mif2/man/probe.Rd
branches/mif2/src/SSA_wrapper.c
branches/mif2/src/dmeasure.c
branches/mif2/src/dprocess.c
branches/mif2/src/euler.c
branches/mif2/src/partrans.c
branches/mif2/src/pomp_fun.c
branches/mif2/src/pomp_internal.h
branches/mif2/src/rmeasure.c
branches/mif2/src/rprocess.c
branches/mif2/src/simulate.c
branches/mif2/src/skeleton.c
branches/mif2/src/trajectory.c
branches/mif2/tests/bbs-trajmatch.Rout.save
branches/mif2/tests/bbs.Rout.save
branches/mif2/tests/blowflies.Rout.save
branches/mif2/tests/dacca.Rout.save
branches/mif2/tests/dimchecks.Rout.save
branches/mif2/tests/fhn.Rout.save
branches/mif2/tests/filtfail.Rout.save
branches/mif2/tests/gillespie.Rout.save
branches/mif2/tests/gompertz.R
branches/mif2/tests/gompertz.Rout.save
branches/mif2/tests/logistic.Rout.save
branches/mif2/tests/ou2-bsmc.Rout.save
branches/mif2/tests/ou2-forecast.R
branches/mif2/tests/ou2-forecast.Rout.save
branches/mif2/tests/ou2-icfit.R
branches/mif2/tests/ou2-icfit.Rout.save
branches/mif2/tests/ou2-kalman.Rout.save
branches/mif2/tests/ou2-mif-fp.R
branches/mif2/tests/ou2-mif-fp.Rout.save
branches/mif2/tests/ou2-mif.R
branches/mif2/tests/ou2-mif.Rout.save
branches/mif2/tests/ou2-mif2.R
branches/mif2/tests/ou2-mif2.Rout.save
branches/mif2/tests/ou2-nlf.Rout.save
branches/mif2/tests/ou2-pmcmc.R
branches/mif2/tests/ou2-pmcmc.Rout.save
branches/mif2/tests/ou2-probe.Rout.save
branches/mif2/tests/ou2-procmeas.Rout.save
branches/mif2/tests/ou2-simulate.Rout.save
branches/mif2/tests/ou2-trajmatch.Rout.save
branches/mif2/tests/pfilter.Rout.save
branches/mif2/tests/pomppomp.Rout.save
branches/mif2/tests/ricker-bsmc.Rout.save
branches/mif2/tests/ricker-probe.Rout.save
branches/mif2/tests/ricker-spect.Rout.save
branches/mif2/tests/ricker.Rout.save
branches/mif2/tests/rw2.Rout.save
branches/mif2/tests/sir.Rout.save
branches/mif2/tests/skeleton.Rout.save
branches/mif2/tests/steps.Rout.save
branches/mif2/tests/synlik.Rout.save
branches/mif2/tests/verhulst.Rout.save
Log:
commit the whole directory for compatibility
Modified: branches/mif2/DESCRIPTION
===================================================================
--- branches/mif2/DESCRIPTION 2013-02-28 05:47:07 UTC (rev 829)
+++ branches/mif2/DESCRIPTION 2013-02-28 05:48:56 UTC (rev 830)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.44-1
-Date: 2013-01-15
+Date: 2013-02-26
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: branches/mif2/NAMESPACE
===================================================================
--- branches/mif2/NAMESPACE 2013-02-28 05:47:07 UTC (rev 829)
+++ branches/mif2/NAMESPACE 2013-02-28 05:48:56 UTC (rev 830)
@@ -1,6 +1,5 @@
useDynLib(
pomp,
- get_pomp_fun,
bspline_basis,
periodic_bspline_basis,
bspline_basis_function,
Modified: branches/mif2/R/aaa.R
===================================================================
--- branches/mif2/R/aaa.R 2013-02-28 05:47:07 UTC (rev 829)
+++ branches/mif2/R/aaa.R 2013-02-28 05:48:56 UTC (rev 830)
@@ -1,7 +1,7 @@
## .onAttach <- function (...) {
-## version <- library(help=pomp)$info[[1]]
-## version <- strsplit(version[pmatch("Version",version)]," ")[[1]]
-## version <- version[nchar(version)>0][2]
+## version <- library(help=pomp)$info[[1L]]
+## version <- strsplit(version[pmatch("Version",version)]," ")[[1L]]
+## version <- version[nchar(version)>0][2L]
## packageStartupMessage("This is pomp version ",version,"\n")
## }
Modified: branches/mif2/R/basic-probes.R
===================================================================
--- branches/mif2/R/basic-probes.R 2013-02-28 05:47:07 UTC (rev 829)
+++ branches/mif2/R/basic-probes.R 2013-02-28 05:48:56 UTC (rev 830)
@@ -53,9 +53,9 @@
method <- match.arg(method)
lag <- as.integer(lag)
transform <- match.fun(transform)
- var1 <- vars[1]
+ var1 <- vars[1L]
if (length(vars)>1)
- var2 <- vars[2]
+ var2 <- vars[2L]
else
var2 <- var1
function (y) {
@@ -85,9 +85,9 @@
method <- match.arg(method)
lag <- as.integer(lag)
transform <- match.fun(transform)
- var1 <- vars[1]
+ var1 <- vars[1L]
if (length(vars)>1)
- var2 <- vars[2]
+ var2 <- vars[2L]
else
var2 <- var1
function (y) {
@@ -134,8 +134,8 @@
lags <- as.integer(lags)
function (y) .Call(
probe_ccf,
- x=transform(y[vars[1],,drop=TRUE]),
- y=transform(y[vars[2],,drop=TRUE]),
+ x=transform(y[vars[1L],,drop=TRUE]),
+ y=transform(y[vars[2L],,drop=TRUE]),
lags=lags,
corr=corr
)
Modified: branches/mif2/R/bsmc.R
===================================================================
--- branches/mif2/R/bsmc.R 2013-02-28 05:47:07 UTC (rev 829)
+++ branches/mif2/R/bsmc.R 2013-02-28 05:48:56 UTC (rev 830)
@@ -31,327 +31,368 @@
setGeneric("bsmc",function(object,...)standardGeneric("bsmc"))
-setMethod(
- "bsmc",
- "pomp",
- function (object, params, Np, est,
- smooth = 0.1,
- ntries = 1,
- tol = 1e-17,
- lower = -Inf, upper = Inf,
- seed = NULL,
- verbose = getOption("verbose"),
- max.fail = 0,
- transform = FALSE,
- ...) {
+bsmc.internal <- function (object, params, Np, est,
+ smooth = 0.1,
+ ntries = 1,
+ tol = 1e-17,
+ lower = -Inf, upper = Inf,
+ seed = NULL,
+ verbose = getOption("verbose"),
+ max.fail = 0,
+ transform = FALSE,
+ .getnativesymbolinfo = TRUE,
+ ...) {
- transform <- as.logical(transform)
+ 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(1) ## need to initialize the RNG
- save.seed <- get(".Random.seed",pos=.GlobalEnv)
- set.seed(seed)
- }
+ 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("bsmc"),"error: ")
+ error.prefix <- paste(sQuote("bsmc"),"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(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 (transform)
- params <- partrans(object,params,dir="inverse")
+ 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 (transform)
+ params <- partrans(object,params,dir="inverse",
+ .getnativesymbolinfo=ptsi.inv)
+ 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
+ )
+ )
+ }
- 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
- 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 (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 (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)
- 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(3.6) p10
+ shrink <- sqrt(1-hsq)
- hsq <- smooth^2 # see Liu & West eq(3.6) p10
- shrink <- sqrt(1-hsq)
+ if (
+ ((length(lower)>1)&&(length(lower)!=npars.est))||
+ ((length(upper)>1)&&(length(upper)!=npars.est))
+ ) {
+ stop(
+ error.prefix,
+ sQuote("lower")," and ",sQuote("upper"),
+ " must each have length 1 or length equal to that of ",sQuote("est"),
+ call.=FALSE
+ )
+ }
- if (
- ((length(lower)>1)&&(length(lower)!=npars.est))||
- ((length(upper)>1)&&(length(upper)!=npars.est))
- ) {
- stop(
- error.prefix,
- sQuote("lower")," and ",sQuote("upper"),
- " must each have length 1 or length equal to that of ",sQuote("est"),
- call.=FALSE
- )
- }
+ for (j in seq_len(Np)) {
+ if (any((params[estind,j]<lower)|(params[estind,j]>upper))) {
+ ind <- which((params[estind,j]<lower)|(params[estind,j]>upper))
+ stop(
+ error.prefix,
+ "parameter(s) ",paste(paramnames[estind[ind]],collapse=","),
+ " in column ",j," in ",sQuote("params"),
+ " is/are outside the box defined by ",
+ sQuote("lower")," and ",sQuote("upper"),
+ call.=FALSE
+ )
+ }
+ }
- for (j in seq_len(Np)) {
- if (any((params[estind,j]<lower)|(params[estind,j]>upper))) {
- ind <- which((params[estind,j]<lower)|(params[estind,j]>upper))
- stop(
- error.prefix,
- "parameter(s) ",paste(paramnames[estind[ind]],collapse=","),
- " in column ",j," in ",sQuote("params"),
- " is/are outside the box defined by ",
- sQuote("lower")," and ",sQuote("upper"),
- call.=FALSE
- )
- }
- }
+ 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
- xstart <- init.state(
- object,
- params=if (transform) {
- partrans(object,params,dir="forward")
- } else {
- params
- }
- )
- statenames <- rownames(xstart)
- nvars <- nrow(xstart)
-
- times <- time(object,t0=TRUE)
- x <- xstart
+ evidence <- rep(NA,ntimes)
+ eff.sample.size <- rep(NA,ntimes)
+ nfail <- 0
+
+ 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]))
- evidence <- rep(NA,ntimes)
- eff.sample.size <- rep(NA,ntimes)
- nfail <- 0
-
- 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))
+ )
+ )
+ }
- if (verbose) {
- cat("at step",nt,"(time =",times[nt+1],")\n")
- print(
- rbind(
- prior.mean=params.mean[estind],
- prior.sd=sqrt(diag(params.var))
- )
+ ## update mean of states at time nt as per L&W AGM (1)
+ tries <- rprocess(
+ object,
+ xstart=parmat(x,nrep=ntries),
+ times=times[c(nt,nt+1)],
+ params=if (transform) {
+ partrans(object,params,dir="forward",
+ .getnativesymbolinfo=ptsi.for)
+ } else {
+ params
+ },
+ offset=1,
+ .getnativesymbolinfo=gnsi.rproc
)
- }
+ dim(tries) <- c(nvars,Np,ntries,1)
+ mu <- apply(tries,c(1,2,4),mean)
+ rownames(mu) <- statenames
+ ## shrink parameters towards mean as per Liu & West eq (3.3) and L&W AGM (1)
+ m <- shrink*params+(1-shrink)*params.mean
+ gnsi.rproc <- FALSE
+
+ ## evaluate probability of obervation given mean value of parameters and states (used in L&W AGM (5) below)
+ g <- dmeasure(
+ object,
+ y=object at data[,nt,drop=FALSE],
+ x=mu,
+ times=times[nt+1],
+ params=if (transform) {
+ partrans(object,m,dir="forward",
+ .getnativesymbolinfo=ptsi.for)
+ } else {
+ m
+ },
+ .getnativesymbolinfo=gnsi.dmeas
+ )
+ gnsi.dmeas <- FALSE
+ storeForEvidence1 <- log(sum(g))
+ ## sample indices -- From L&W AGM (2)
+ ## k <- .Call(systematic_resampling,g)
+ k <- sample.int(n=Np,size=Np,replace=TRUE,prob=g)
+ params <- params[,k]
+ m <- m[,k]
+ g <- g[k]
- ## update mean of states at time nt as per L&W AGM (1)
- tries <- rprocess(
- object,
- xstart=parmat(x,nrep=ntries),
- times=times[c(nt,nt+1)],
- params=if (transform) {
- partrans(object,params,dir="forward")
- } else {
- params
- },
- offset=1
- )
- dim(tries) <- c(nvars,Np,ntries,1)
- mu <- apply(tries,c(1,2,4),mean)
- rownames(mu) <- statenames
- ## shrink parameters towards mean as per Liu & West eq (3.3) and L&W AGM (1)
- m <- shrink*params+(1-shrink)*params.mean
-
- ## evaluate probability of obervation given mean value of parameters and states (used in L&W AGM (5) below)
- g <- dmeasure(
- object,
- y=object at data[,nt,drop=FALSE],
- x=mu,
- times=times[nt+1],
- params=if (transform) {
- partrans(object,m,dir="forward")
- } else {
- m
- }
- )
- storeForEvidence1 <- log(sum(g))
- ## sample indices -- From L&W AGM (2)
-## k <- .Call(systematic_resampling,g)
- k <- sample.int(n=Np,size=Np,replace=TRUE,prob=g)
- params <- params[,k]
- m <- m[,k]
- g <- g[k]
+ ## sample new parameter vector as per L&W AGM (3) and Liu & West eq(3.2)
+ pvec <- try(
+ mvtnorm::rmvnorm(
+ n=Np,
+ mean=rep(0,npars.est),
+ sigma=hsq*params.var,
+ method="svd"
+ ),
+ silent=FALSE
+ )
+ if (inherits(pvec,"try-error"))
+ stop(error.prefix,"error in ",sQuote("rmvnorm"),call.=FALSE)
+ if (any(!is.finite(pvec)))
+ stop(error.prefix,"extreme particle depletion",call.=FALSE)
+ params[estind,] <- m[estind,]+t(pvec)
- ## sample new parameter vector as per L&W AGM (3) and Liu & West eq(3.2)
- pvec <- try(
- mvtnorm::rmvnorm(
- n=Np,
- mean=rep(0,npars.est),
- sigma=hsq*params.var,
- method="svd"
- ),
- silent=FALSE
- )
- if (inherits(pvec,"try-error"))
- stop(error.prefix,"error in ",sQuote("rmvnorm"),call.=FALSE)
- if (any(!is.finite(pvec)))
- stop(error.prefix,"extreme particle depletion",call.=FALSE)
- params[estind,] <- m[estind,]+t(pvec)
+ if (transform)
+ tparams <- partrans(object,params,dir="forward",
+ .getnativesymbolinfo=ptsi.for)
+
+ ## sample current state vector x^(g)_(t+1) as per L&W AGM (4)
+ X <- rprocess(
+ object,
+ xstart=x[,k,drop=FALSE],
+ times=times[c(nt,nt+1)],
+ params=if (transform) {
+ tparams
+ } else {
+ params
+ },
+ offset=1,
+ .getnativesymbolinfo=gnsi.rproc
+ )
- if (transform)
- tparams <- partrans(object,params,dir="forward")
-
- ## sample current state vector x^(g)_(t+1) as per L&W AGM (4)
- X <- rprocess(
- object,
- xstart=x[,k,drop=FALSE],
- times=times[c(nt,nt+1)],
- params=if (transform) {
- tparams
- } else {
- params
- },
- offset=1
- )
+ ## evaluate likelihood of observation given X (from L&W AGM (4))
+ numer <- dmeasure(
+ object,
+ y=object at data[,nt,drop=FALSE],
+ x=X,
+ times=times[nt+1],
+ params=if (transform) {
+ tparams
+ } else {
+ params
+ },
+ .getnativesymbolinfo=gnsi.dmeas
+ )
+ ## evaluate weights as per L&W AGM (5)
- ## evaluate likelihood of observation given X (from L&W AGM (4))
- numer <- dmeasure(
- object,
- y=object at data[,nt,drop=FALSE],
- x=X,
- times=times[nt+1],
- params=if (transform) {
- tparams
- } else {
- params
- }
- )
- ## evaluate weights as per L&W AGM (5)
+ weights <- numer/g
+ storeForEvidence2 <- log(mean(weights))
+
+ ## apply box constraints as per the priors
+ for (j in seq_len(Np)) {
+ ## the following seems problematic: will it tend to make the boundaries repellors
+ if (any((params[estind,j]<lower)|(params[estind,j]>upper))) {
+ weights[j] <- 0
+ }
+ ## might this rejection method be preferable?
+ ## while (any((params[estind,j]<lower)|(params[estind,j]>upper))) {
+ ## ## rejection method
+ ## pvec <- try(
+ ## mvtnorm::rmvnorm(
+ ## n=1,
+ ## mean=rep(0,npars.est),
+ ## sigma=hsq*params.var,
+ ## method="eigen"
+ ## ),
+ ## silent=FALSE
+ ## )
+ ## if (inherits(pvec,"try-error"))
+ ## stop(error.prefix,"error in ",sQuote("rmvnorm"),call.=FALSE)
+ ## if (any(!is.finite(pvec)))
+ ## stop(error.prefix,"extreme particle depletion",call.=FALSE)
+ ## params[estind,j] <- m[estind,j]+pvec[1,]
+ ## }
+ }
- weights <- numer/g
- storeForEvidence2 <- log(mean(weights))
-
- ## apply box constraints as per the priors
- for (j in seq_len(Np)) {
- ## the following seems problematic: will it tend to make the boundaries repellors
- if (any((params[estind,j]<lower)|(params[estind,j]>upper))) {
- weights[j] <- 0
- }
- ## might this rejection method be preferable?
- ## while (any((params[estind,j]<lower)|(params[estind,j]>upper))) {
- ## ## rejection method
- ## pvec <- try(
- ## mvtnorm::rmvnorm(
- ## n=1,
- ## mean=rep(0,npars.est),
- ## sigma=hsq*params.var,
- ## method="eigen"
- ## ),
- ## silent=FALSE
- ## )
- ## if (inherits(pvec,"try-error"))
- ## stop(error.prefix,"error in ",sQuote("rmvnorm"),call.=FALSE)
- ## if (any(!is.finite(pvec)))
- ## stop(error.prefix,"extreme particle depletion",call.=FALSE)
- ## params[estind,j] <- m[estind,j]+pvec[1,]
- ## }
- }
+ x[,] <- X
+
+ ## test for failure to filter
+ dim(weights) <- NULL
+ 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] <- storeForEvidence1+storeForEvidence2
+ weights[failures] <- 0
+ weights <- weights/sum(weights)
+ ## compute effective sample-size
+ eff.sample.size[nt] <- 1/crossprod(weights)
+ }
- x[,] <- X
-
- ## test for failure to filter
- dim(weights) <- NULL
- 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] <- storeForEvidence1+storeForEvidence2
- 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")
+ }
- 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]
+ }
- ## 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]
- }
-
- }
-
- if (!is.null(seed)) {
- assign(".Random.seed",save.seed,pos=.GlobalEnv)
- seed <- save.seed
- }
-
- ## if (transform) {
- ## params <- partrans(object,params,dir="forward")
- ## prior <- partrans(object,prior,dir="forward")
- ## }
+ .getnativesymbolinfo <- FALSE
+
+ }
- ## replace parameters with point estimate (posterior median)
- coef(object,transform=transform) <- apply(params,1,median)
+ 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),
- weights=weights
- )
+ 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),
+ weights=weights
+ )
+}
+
+setMethod(
+ "bsmc",
+ "pomp",
+ function (object, params, Np, est,
+ smooth = 0.1,
+ ntries = 1,
+ tol = 1e-17,
+ lower = -Inf, upper = Inf,
+ seed = NULL,
+ verbose = getOption("verbose"),
+ max.fail = 0,
+ transform = FALSE,
+ ...) {
+ bsmc.internal(
+ object=object,
+ params=params,
+ Np=Np,
+ est=est,
+ smooth=smooth,
+ ntries=ntries,
+ tol=tol,
+ lower=lower,
+ upper=upper,
+ seed=seed,
+ verbose=verbose,
+ max.fail=max.fail,
+ transform=transform,
+ ...
+ )
}
)
@@ -374,17 +415,17 @@
panel=function (x, y, ...) { ## prior, posterior pairwise scatterplot
op <- par(new=TRUE)
on.exit(par(op))
- i <- which(x[1]==all[1,])
- j <- which(y[1]==all[1,])
+ i <- which(x[1L]==all[1L,])
+ j <- which(y[1L]==all[1L,])
points(prior[p1,i],prior[p1,j],pch=20,col=rgb(0.85,0.85,0.85,0.1),xlim=range(all[,i]),ylim=range(all[,j]))
points(post[p2,i],post[p2,j],pch=20,col=rgb(0,0,1,0.01))
},
diag.panel=function (x, ...) { ## marginal posterior histogram
- i <- which(x[1]==all[1,])
+ i <- which(x[1L]==all[1L,])
d1 <- density(prior[,i])
d2 <- density(post[,i])
usr <- par('usr')
- op <- par(usr=c(usr[1:2],0,1.5*max(d1$y,d2$y)))
+ op <- par(usr=c(usr[c(1L,2L)],0,1.5*max(d1$y,d2$y)))
on.exit(par(op))
polygon(d1,col=rgb(0.85,0.85,0.85,0.5))
polygon(d2,col=rgb(0,0,1,0.5))
Modified: branches/mif2/R/dmeasure-pomp.R
===================================================================
--- branches/mif2/R/dmeasure-pomp.R 2013-02-28 05:47:07 UTC (rev 829)
+++ branches/mif2/R/dmeasure-pomp.R 2013-02-28 05:48:56 UTC (rev 830)
@@ -1,9 +1,11 @@
## evaluate the measurement model density function
setGeneric("dmeasure",function(object,...)standardGeneric("dmeasure"))
-dmeasure.internal <- function (object, y, x, times, params, log = FALSE, ...) {
- fun <- get.pomp.fun(object at dmeasure)
- .Call(do_dmeasure,object,y,x,times,params,log,fun)
+dmeasure.internal <- function (object, y, x, times, params, log = FALSE, .getnativesymbolinfo = TRUE, ...) {
+ .Call(do_dmeasure,object,y,x,times,params,log,.getnativesymbolinfo)
}
-setMethod("dmeasure","pomp",dmeasure.internal)
+setMethod("dmeasure","pomp",
+ function (object, y, x, times, params, log = FALSE, ...)
+ dmeasure.internal(object=object,y=y,x=x,times=times,params=params,log=log,...)
+ )
Modified: branches/mif2/R/dprocess-pomp.R
===================================================================
--- branches/mif2/R/dprocess-pomp.R 2013-02-28 05:47:07 UTC (rev 829)
+++ branches/mif2/R/dprocess-pomp.R 2013-02-28 05:48:56 UTC (rev 830)
@@ -1,7 +1,10 @@
## evaluate the process model density function
setGeneric("dprocess",function(object,...)standardGeneric("dprocess"))
-dprocess.internal <- function (object, x, times, params, log = FALSE, ...)
- .Call(do_dprocess,object,x,times,params,log)
+dprocess.internal <- function (object, x, times, params, log = FALSE, .getnativesymbolinfo = TRUE, ...)
+ .Call(do_dprocess,object,x,times,params,log,.getnativesymbolinfo)
-setMethod("dprocess","pomp",dprocess.internal)
+setMethod("dprocess","pomp",
+ function (object, x, times, params, log = FALSE, ...)
+ dprocess.internal(object=object,x=x,times=times,params=params,log=log,...)
+ )
Modified: branches/mif2/R/init-state-pomp.R
===================================================================
--- branches/mif2/R/init-state-pomp.R 2013-02-28 05:47:07 UTC (rev 829)
+++ branches/mif2/R/init-state-pomp.R 2013-02-28 05:48:56 UTC (rev 830)
@@ -1,12 +1,16 @@
+## initialize the state variables of the process model
+init.state.internal <- function (object, params, t0, ...) {
+ if (missing(t0)) t0 <- object at t0
+ if (missing(params)) params <- coef(object)
+ .Call(do_init_state,object,params,t0)
+}
+
setGeneric("init.state",function(object,...)standardGeneric("init.state"))
-## initialize the process model
setMethod(
'init.state',
'pomp',
- function (object, params, t0, ...) { # the package algorithms will only use these arguments
- if (missing(t0)) t0 <- object at t0
- if (missing(params)) params <- coef(object)
- .Call(do_init_state,object,params,t0)
+ function (object, params, t0, ...) {
+ init.state.internal(object=object,params=params,t0=t0,...)
}
)
Modified: branches/mif2/R/mif-class.R
===================================================================
--- branches/mif2/R/mif-class.R 2013-02-28 05:47:07 UTC (rev 829)
+++ branches/mif2/R/mif-class.R 2013-02-28 05:48:56 UTC (rev 830)
@@ -10,7 +10,7 @@
particles = 'function',
var.factor='numeric',
ic.lag='integer',
- cooling.factor='numeric',
+ cooling.type='character',
cooling.fraction='numeric',
method='character',
random.walk.sd = 'numeric',
Modified: branches/mif2/R/mif-methods.R
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 830
More information about the pomp-commits
mailing list