[Pomp-commits] r826 - in pkg/pomp: . R man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 26 14:20:27 CET 2013
Author: kingaa
Date: 2013-02-26 14:20:27 +0100 (Tue, 26 Feb 2013)
New Revision: 826
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/NAMESPACE
pkg/pomp/R/bsmc.R
pkg/pomp/R/dmeasure-pomp.R
pkg/pomp/R/dprocess-pomp.R
pkg/pomp/R/init-state-pomp.R
pkg/pomp/R/mif-methods.R
pkg/pomp/R/mif.R
pkg/pomp/R/particles-mif.R
pkg/pomp/R/pfilter.R
pkg/pomp/R/plot-pomp.R
pkg/pomp/R/plugins.R
pkg/pomp/R/pmcmc.R
pkg/pomp/R/pomp-fun.R
pkg/pomp/R/pomp-methods.R
pkg/pomp/R/probe.R
pkg/pomp/R/rmeasure-pomp.R
pkg/pomp/R/rprocess-pomp.R
pkg/pomp/R/simulate-pomp.R
pkg/pomp/R/skeleton-pomp.R
pkg/pomp/R/spect-match.R
pkg/pomp/R/trajectory-pomp.R
pkg/pomp/man/pmcmc.Rd
pkg/pomp/man/probe.Rd
pkg/pomp/src/SSA_wrapper.c
pkg/pomp/src/dmeasure.c
pkg/pomp/src/dprocess.c
pkg/pomp/src/euler.c
pkg/pomp/src/partrans.c
pkg/pomp/src/pomp_fun.c
pkg/pomp/src/pomp_internal.h
pkg/pomp/src/rmeasure.c
pkg/pomp/src/rprocess.c
pkg/pomp/src/simulate.c
pkg/pomp/src/skeleton.c
pkg/pomp/src/trajectory.c
pkg/pomp/tests/dimchecks.Rout.save
pkg/pomp/tests/ou2-mif.Rout.save
pkg/pomp/tests/ou2-mif2.Rout.save
pkg/pomp/tests/ou2-pmcmc.R
pkg/pomp/tests/ou2-pmcmc.Rout.save
pkg/pomp/tests/ricker-probe.Rout.save
pkg/pomp/tests/sir.Rout.save
Log:
- new mechanism to prevent unnecessary symbol-table lookups when native routines are used
- named functions are now used for most methods, to make error messages more sensible
- some changes to the implementation of pmcmc
- fix bug in pmcmc method on pmcmc objects
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/DESCRIPTION 2013-02-26 13:20:27 UTC (rev 826)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.44-1
-Date: 2013-02-04
+Date: 2013-02-14
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: pkg/pomp/NAMESPACE
===================================================================
--- pkg/pomp/NAMESPACE 2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/NAMESPACE 2013-02-26 13:20:27 UTC (rev 826)
@@ -1,6 +1,5 @@
useDynLib(
pomp,
- get_pomp_fun,
bspline_basis,
periodic_bspline_basis,
bspline_basis_function,
Modified: pkg/pomp/R/bsmc.R
===================================================================
--- pkg/pomp/R/bsmc.R 2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/R/bsmc.R 2013-02-26 13:20:27 UTC (rev 826)
@@ -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(n=1L) ## 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,
+ ...
+ )
}
)
Modified: pkg/pomp/R/dmeasure-pomp.R
===================================================================
--- pkg/pomp/R/dmeasure-pomp.R 2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/R/dmeasure-pomp.R 2013-02-26 13:20:27 UTC (rev 826)
@@ -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: pkg/pomp/R/dprocess-pomp.R
===================================================================
--- pkg/pomp/R/dprocess-pomp.R 2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/R/dprocess-pomp.R 2013-02-26 13:20:27 UTC (rev 826)
@@ -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: pkg/pomp/R/init-state-pomp.R
===================================================================
--- pkg/pomp/R/init-state-pomp.R 2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/R/init-state-pomp.R 2013-02-26 13:20:27 UTC (rev 826)
@@ -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: pkg/pomp/R/mif-methods.R
===================================================================
--- pkg/pomp/R/mif-methods.R 2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/R/mif-methods.R 2013-02-26 13:20:27 UTC (rev 826)
@@ -1,48 +1,48 @@
## this file contains short definitions of methods for the 'mif' class
-conv.rec <- function (object, ...)
- stop("function ",sQuote("conv.rec")," is undefined for objects of class ",sQuote(class(object)))
-setGeneric('conv.rec')
-
## extract the estimated log likelihood
setMethod('logLik','mif',function(object,...)object at loglik)
## extract the convergence record
-setMethod(
- 'conv.rec',
- 'mif',
+conv.rec.internal <- function (object, pars, transform = FALSE, ...) {
+ if (transform) {
+ pars.improper <- c("loglik","nfail")
+ pars.proper <- setdiff(colnames(object at conv.rec),pars.improper)
+ retval <- cbind(
+ t(
+ partrans(
+ object,
+ params=t(object at conv.rec[,pars.proper]),
+ dir="forward"
+ )
+ ),
+ object at conv.rec[,pars.improper]
+ )
+ } else {
+ retval <- object at conv.rec
+ }
+ if (missing(pars))
+ retval
+ else {
+ bad.pars <- setdiff(pars,colnames(retval))
+ if (length(bad.pars)>0)
+ stop(
+ "in ",sQuote("conv.rec"),": name(s) ",
+ paste(sQuote(bad.pars),collapse=","),
+ " correspond to no parameter(s) in ",
+ if (transform) sQuote("conv.rec(object,transform=TRUE)")
+ else sQuote("conv.rec(object,transform=FALSE)"),
+ call.=FALSE
+ )
+ retval[,pars]
+ }
+}
+
+setGeneric("conv.rec",function(object,...)standardGeneric("conv.rec"))
+
+setMethod('conv.rec','mif',
function (object, pars, transform = FALSE, ...) {
- if (transform) {
- pars.improper <- c("loglik","nfail")
- pars.proper <- setdiff(colnames(object at conv.rec),pars.improper)
- retval <- cbind(
- t(
- partrans(
- object,
- params=t(object at conv.rec[,pars.proper]),
- dir="forward"
- )
- ),
- object at conv.rec[,pars.improper]
- )
- } else {
- retval <- object at conv.rec
- }
- if (missing(pars))
- retval
- else {
- bad.pars <- setdiff(pars,colnames(retval))
- if (length(bad.pars)>0)
- stop(
- "in ",sQuote("conv.rec"),": name(s) ",
- paste(sQuote(bad.pars),collapse=","),
- " correspond to no parameter(s) in ",
- if (transform) sQuote("conv.rec(object,transform=TRUE)")
- else sQuote("conv.rec(object,transform=FALSE)"),
- call.=FALSE
- )
- retval[,pars]
- }
+ conv.rec.internal(object=object,pars=pars,transform=transform,...)
}
)
Modified: pkg/pomp/R/mif.R
===================================================================
--- pkg/pomp/R/mif.R 2013-02-05 22:49:05 UTC (rev 825)
+++ pkg/pomp/R/mif.R 2013-02-26 13:20:27 UTC (rev 826)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 826
More information about the pomp-commits
mailing list