[Pomp-commits] r923 - in branches/mif2: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Apr 17 18:47:07 CEST 2014
Author: kingaa
Date: 2014-04-17 18:47:05 +0200 (Thu, 17 Apr 2014)
New Revision: 923
Added:
branches/mif2/R/generics.R
branches/mif2/R/mif2-methods.R
branches/mif2/R/mif2.R
branches/mif2/inst/NEWS
branches/mif2/inst/NEWS.Rd
branches/mif2/man/mif2.Rd
Removed:
branches/mif2/R/authors.R
branches/mif2/R/basic-probes.R
branches/mif2/R/bsmc.R
branches/mif2/R/bsplines.R
branches/mif2/R/builder.R
branches/mif2/R/compare-mif.R
branches/mif2/R/compare-pmcmc.R
branches/mif2/R/dmeasure-pomp.R
branches/mif2/R/dprocess-pomp.R
branches/mif2/R/eulermultinom.R
branches/mif2/R/example.R
branches/mif2/R/init-state-pomp.R
branches/mif2/R/mif-class.R
branches/mif2/R/mif-methods.R
branches/mif2/R/mif.R
branches/mif2/R/nlf-funcs.R
branches/mif2/R/nlf-guts.R
branches/mif2/R/nlf-objfun.R
branches/mif2/R/nlf.R
branches/mif2/R/parmat.R
branches/mif2/R/particles-mif.R
branches/mif2/R/pfilter-methods.R
branches/mif2/R/pfilter.R
branches/mif2/R/plot-pomp.R
branches/mif2/R/plugins.R
branches/mif2/R/pmcmc-methods.R
branches/mif2/R/pmcmc.R
branches/mif2/R/pomp-class.R
branches/mif2/R/pomp-fun.R
branches/mif2/R/pomp-methods.R
branches/mif2/R/pomp.R
branches/mif2/R/probe-match.R
branches/mif2/R/probe.R
branches/mif2/R/profile-design.R
branches/mif2/R/rmeasure-pomp.R
branches/mif2/R/rprocess-pomp.R
branches/mif2/R/sannbox.R
branches/mif2/R/simulate-pomp.R
branches/mif2/R/skeleton-pomp.R
branches/mif2/R/slice-design.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/GPL
branches/mif2/man/LondonYorke.Rd
branches/mif2/man/basic-probes.Rd
branches/mif2/man/blowflies.Rd
branches/mif2/man/bsmc.Rd
branches/mif2/man/bsplines.Rd
branches/mif2/man/builder.Rd
branches/mif2/man/dacca.Rd
branches/mif2/man/dmeasure-pomp.Rd
branches/mif2/man/dprocess-pomp.Rd
branches/mif2/man/eulermultinom.Rd
branches/mif2/man/example.Rd
branches/mif2/man/gompertz.Rd
branches/mif2/man/init.state-pomp.Rd
branches/mif2/man/mif-class.Rd
branches/mif2/man/mif-methods.Rd
branches/mif2/man/mif.Rd
branches/mif2/man/nlf.Rd
branches/mif2/man/ou2.Rd
branches/mif2/man/parmat.Rd
branches/mif2/man/particles-mif.Rd
branches/mif2/man/pfilter-methods.Rd
branches/mif2/man/pfilter.Rd
branches/mif2/man/plugins.Rd
branches/mif2/man/pmcmc-methods.Rd
branches/mif2/man/pmcmc.Rd
branches/mif2/man/pomp-class.Rd
branches/mif2/man/pomp-fun.Rd
branches/mif2/man/pomp-methods.Rd
branches/mif2/man/pomp-package.Rd
branches/mif2/man/pomp.Rd
branches/mif2/man/probe.Rd
branches/mif2/man/probed-pomp-methods.Rd
branches/mif2/man/profile-design.Rd
branches/mif2/man/ricker.Rd
branches/mif2/man/rmeasure-pomp.Rd
branches/mif2/man/rprocess-pomp.Rd
branches/mif2/man/rw2.Rd
branches/mif2/man/sannbox.Rd
branches/mif2/man/simulate-pomp.Rd
branches/mif2/man/sir.Rd
branches/mif2/man/skeleton-pomp.Rd
branches/mif2/man/slice-design.Rd
branches/mif2/man/sobol.Rd
branches/mif2/man/spect.Rd
branches/mif2/man/traj-match.Rd
branches/mif2/man/trajectory-pomp.Rd
branches/mif2/man/verhulst.Rd
branches/mif2/src/
Modified:
branches/mif2/DESCRIPTION
branches/mif2/NAMESPACE
branches/mif2/R/aaa.R
branches/mif2/R/version.R
Log:
- remove overlap with 'pomp'
Modified: branches/mif2/DESCRIPTION
===================================================================
--- branches/mif2/DESCRIPTION 2014-04-17 13:21:28 UTC (rev 922)
+++ branches/mif2/DESCRIPTION 2014-04-17 16:47:05 UTC (rev 923)
@@ -1,37 +1,17 @@
-Package: pomp
+Package: mif2
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.46-1
-Date: 2013-10-30
-Authors at R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="kingaa at umich.edu"),
+Version: 0.50-1
+Date: 2014-04-17
+Authors at R: c(person(given=c("Aaron","A."),family="King",
+ role=c("aut","cre"),email="kingaa at umich.edu"),
person(given=c("Edward","L."),family="Ionides",role=c("aut")),
- person(given=c("Carles"),family="Breto",role=c("aut")),
- person(given=c("Stephen","P."),family="Ellner",role=c("ctb")),
- person(given=c("Matthew","J."),family="Ferrari",role=c("ctb")),
- person(given=c("Bruce","E."),family="Kendall",role=c("ctb")),
- person(given=c("Michael"),family="Lavine",role=c("ctb")),
- person(given=c("Daniel","C."),family="Reuman",role=c("ctb")),
- person(given=c("Helen"),family="Wearing",role=c("ctb")),
- person(given=c("Simon","N."),family="Wood",role=c("ctb")),
- person(given="Dao",family="Nguyen",role=c("ctb")) )
-Maintainer: Aaron A. King <kingaa at umich.edu>
-Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing, Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman, Simon N. Wood, Dao Nguyen
+ person(given="Dao",family="Nguyen",role=c("ctb"))
+ )
URL: http://pomp.r-forge.r-project.org
Description: Inference methods for partially-observed Markov processes
-Depends: R(>= 2.14.1), stats, graphics, mvtnorm, subplex, deSolve
-Imports: methods
+Depends: R(>= 2.15.1), pomp(>= 0.49-1), methods
License: GPL(>= 2)
-LazyLoad: yes
-LazyData: no
-BuildVignettes: no
-Collate: aaa.R authors.R version.R eulermultinom.R plugins.R
- parmat.R slice-design.R profile-design.R sobol.R bsplines.R sannbox.R
- pomp-fun.R pomp-class.R pomp.R pomp-methods.R
- rmeasure-pomp.R rprocess-pomp.R init-state-pomp.R dmeasure-pomp.R dprocess-pomp.R skeleton-pomp.R
- simulate-pomp.R trajectory-pomp.R plot-pomp.R
- pfilter.R pfilter-methods.R traj-match.R bsmc.R
- mif-class.R particles-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-objfun.R nlf.R
- probe.R probe-match.R basic-probes.R spect.R spect-match.R
- builder.R example.R
+LazyData: true
+BuildVignettes: false
+Collate: aaa.R generics.R version.R mif2.R mif2-methods.R
Modified: branches/mif2/NAMESPACE
===================================================================
--- branches/mif2/NAMESPACE 2014-04-17 13:21:28 UTC (rev 922)
+++ branches/mif2/NAMESPACE 2014-04-17 16:47:05 UTC (rev 923)
@@ -1,98 +1,6 @@
-useDynLib(
- pomp,
- bspline_basis,
- periodic_bspline_basis,
- bspline_basis_function,
- systematic_resampling,
- euler_model_simulator,
- euler_model_density,
- lookup_in_table,
- SSA_simulator,
- R_Euler_Multinom,D_Euler_Multinom,R_GammaWN,
- mif_update,
- pfilter_computations,
- simulation_computations,
- iterate_map,traj_transp_and_copy,
- apply_probe_data,apply_probe_sim,
- probe_marginal_setup,probe_marginal_solve,
- probe_acf,probe_ccf,
- probe_nlar,
- synth_loglik,
- pomp_desolve_setup,pomp_desolve_takedown,
- pomp_vf_eval,
- do_partrans,
- do_rprocess,
- do_dprocess,
- do_rmeasure,
- do_dmeasure,
- do_skeleton,
- do_init_state
- )
+import(pomp)
+import(methods)
-importFrom(graphics,plot)
-importFrom(stats,simulate,time,coef,logLik,window)
-importFrom(mvtnorm,dmvnorm,rmvnorm)
-importFrom(subplex,subplex)
-importFrom(deSolve,ode)
-
-exportClasses(
- pomp,
- pfilterd.pomp,
- mif,
- pmcmc,
- traj.matched.pomp,
- probed.pomp,probe.matched.pomp,
- spect.pomp,spect.matched.pomp
- )
-
-exportMethods(
- pomp,
- plot,show,print,coerce,summary,logLik,window,"$",
- dprocess,rprocess,rmeasure,dmeasure,init.state,skeleton,
- data.array,obs,partrans,coef,"coef<-",time,"time<-",timezero,"timezero<-",
- simulate,pfilter,
- eff.sample.size,cond.logLik,
- particles,mif,continue,states,trajectory,
- pred.mean,pred.var,filter.mean,conv.rec,
- bsmc,
- pmcmc,dprior,
- spect,probe,
- probe.match,traj.match
- )
-
-export(
- as.data.frame.pomp,
- as.data.frame.pfilterd.pomp,
- reulermultinom,
- deulermultinom,
- rgammawn,
- euler.sim,
- discrete.time.sim,
- onestep.sim,
- onestep.dens,
- gillespie.sim,
- sobol,
- sobolDesign,
- sliceDesign,
- profileDesign,
- bspline.basis,
- periodic.bspline.basis,
- compare.mif,
- nlf,
- parmat,
- probe.mean,
- probe.median,
- probe.var,
- probe.sd,
- probe.period,
- probe.quantile,
- probe.acf,
- probe.ccf,
- probe.nlar,
- probe.marginal,
- sannbox,
- spect.match,
- traj.match.objfun,
- pompBuilder,
- pompExample
- )
+exportClasses(mif2d.pomp)
+exportMethods(mif2,continue,plot,conv.rec,logLik)
+export()
Modified: branches/mif2/R/aaa.R
===================================================================
--- branches/mif2/R/aaa.R 2014-04-17 13:21:28 UTC (rev 922)
+++ branches/mif2/R/aaa.R 2014-04-17 16:47:05 UTC (rev 923)
@@ -1,25 +1,3 @@
-## .onAttach <- function (...) {
-## 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")
-## }
-
-setGeneric("print",function(x,...)standardGeneric("print"))
-setGeneric("plot",function(x,y,...)standardGeneric("plot"))
-setGeneric("summary",function(object,...)standardGeneric("summary"))
-setGeneric("simulate",function(object,nsim=1,seed=NULL,...)standardGeneric("simulate"))
-setGeneric("time",function(x,...)standardGeneric("time"))
-setGeneric("coef",function(object,...)standardGeneric("coef"))
-setGeneric("logLik",function(object,...)standardGeneric("logLik"))
-setGeneric("window",function(x,...)standardGeneric("window"))
-setGeneric("continue",function(object,...)standardGeneric("continue"))
-setGeneric("pred.mean",function(object,...)standardGeneric("pred.mean"))
-setGeneric("pred.var",function(object,...)standardGeneric("pred.var"))
-setGeneric("filter.mean",function(object,...)standardGeneric("filter.mean"))
-setGeneric("cond.logLik",function(object,...)standardGeneric("cond.logLik"))
-setGeneric("eff.sample.size",function(object,...)standardGeneric("eff.sample.size"))
-
if (!exists("paste0",where="package:base")) {
paste0 <- function(...) paste(...,sep="")
}
Deleted: branches/mif2/R/authors.R
===================================================================
--- branches/mif2/R/authors.R 2014-04-17 13:21:28 UTC (rev 922)
+++ branches/mif2/R/authors.R 2014-04-17 16:47:05 UTC (rev 923)
@@ -1,12 +0,0 @@
-list(
- aak=person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="kingaa at umich.edu"),
- eli=person(given=c("Edward","L."),family="Ionides",role=c("ctb")),
- cb=person(given=c("Carles"),family="Breto",role=c("ctb")),
- spe=person(given=c("Stephen","P."),family="Ellner",role=c("ctb")),
- bek=person(given=c("Bruce","E."),family="Kendall",role=c("ctb")),
- mf=person(given=c("Matthew","J."),family="Ferrari",role=c("ctb")),
- ml=person(given=c("Michael"),family="Lavine",role=c("ctb")),
- dcr=person(given=c("Daniel","C."),family="Reuman",role=c("ctb")),
- hw=person(given=c("Helen"),family="Wearing",role=c("ctb")),
- snw=person(given=c("Simon","N."),family="Wood",role=c("ctb"))
- ) -> author.list
Deleted: branches/mif2/R/basic-probes.R
===================================================================
--- branches/mif2/R/basic-probes.R 2014-04-17 13:21:28 UTC (rev 922)
+++ branches/mif2/R/basic-probes.R 2014-04-17 16:47:05 UTC (rev 923)
@@ -1,176 +0,0 @@
-probe.mean <- function (var, trim = 0, transform = identity, na.rm = TRUE) {
- if (length(var)>1) stop(sQuote("probe.mean")," is a univariate probe")
- transform <- match.fun(transform)
- function(y) mean(x=transform(y[var,]),trim=trim,na.rm=na.rm)
-}
-
-probe.median <- function (var, na.rm = TRUE) {
- if (length(var)>1) stop(sQuote("probe.median")," is a univariate probe")
- function(y) median(x=as.numeric(y[var,]),na.rm=na.rm)
-}
-
-probe.var <- function (var, transform = identity, na.rm = TRUE) {
- if (length(var)>1) stop(sQuote("probe.var")," is a univariate probe")
- transform <- match.fun(transform)
- function(y) var(x=transform(y[var,]),na.rm=na.rm)
-}
-
-probe.sd <- function (var, transform = identity, na.rm = TRUE) {
- if (length(var)>1) stop(sQuote("probe.sd")," is a univariate probe")
- transform <- match.fun(transform)
- function(y) sd(x=transform(y[var,]),na.rm=na.rm)
-}
-
-probe.period <- function (var, kernel.width, transform = identity) {
- if (length(var)>1) stop(sQuote("probe.period")," is a univariate probe")
- transform <- match.fun(transform)
- function (y) {
- zz <- spec.pgram(
- x=transform(y[var,]),
- kernel=kernel("modified.daniell",m=kernel.width),
- taper=0,
- fast=FALSE,
- pad=0,
- detrend=FALSE,
- plot=FALSE
- )
- 1/zz$freq[which.max(zz$spec)]
- }
-}
-
-probe.quantile <- function (var, prob, transform = identity) {
- if (length(var)>1) stop(sQuote("probe.quantile")," is a univariate probe")
- transform <- match.fun(transform)
- function (y) quantile(transform(y[var,]),probs=prob)
-}
-
-probe.cov <- function (
- vars,
- lag,
- method = c("pearson", "kendall", "spearman"),
- transform = identity
- ) {
- method <- match.arg(method)
- lag <- as.integer(lag)
- transform <- match.fun(transform)
- var1 <- vars[1L]
- if (length(vars)>1)
- var2 <- vars[2L]
- else
- var2 <- var1
- function (y) {
- if (lag>=0) {
- val <- cov(
- x=transform(y[var1,seq(from=1+lag,to=ncol(y),by=1)]),
- y=transform(y[var2,seq(from=1,to=ncol(y)-lag,by=1)]),
- method=method
- )
- } else {
- val <- cov(
- x=transform(y[var1,seq(from=1,to=ncol(y)+lag,by=1)]),
- y=transform(y[var2,seq(from=-lag,to=ncol(y),by=1)]),
- method=method
- )
- }
- val
- }
-}
-
-probe.cor <- function (
- vars,
- lag,
- method = c("pearson", "kendall", "spearman"),
- transform = identity
- ) {
- method <- match.arg(method)
- lag <- as.integer(lag)
- transform <- match.fun(transform)
- var1 <- vars[1L]
- if (length(vars)>1)
- var2 <- vars[2L]
- else
- var2 <- var1
- function (y) {
- if (lag>=0) {
- val <- cor(
- x=transform(y[var1,seq(from=1+lag,to=ncol(y),by=1)]),
- y=transform(y[var2,seq(from=1,to=ncol(y)-lag,by=1)]),
- method=method
- )
- } else {
- val <- cor(
- x=transform(y[var1,seq(from=1,to=ncol(y)+lag,by=1)]),
- y=transform(y[var2,seq(from=-lag,to=ncol(y),by=1)]),
- method=method
- )
- }
- val
- }
-}
-
-probe.acf <- function (var, lags, type = c("covariance", "correlation"), transform = identity) {
- type <- match.arg(type)
- corr <- type=="correlation"
- transform <- match.fun(transform)
- if (corr && any(lags==0)) {
- warning("useless zero lag discarded in ",sQuote("probe.acf"))
- lags <- lags[lags!=0]
- }
- lags <- as.integer(lags)
- function (y) .Call(
- probe_acf,
- x=transform(y[var,,drop=FALSE]),
- lags=lags,
- corr=corr
- )
-}
-
-probe.ccf <- function (vars, lags, type = c("covariance", "correlation"), transform = identity) {
- type <- match.arg(type)
- corr <- type=="correlation"
- transform <- match.fun(transform)
- if (length(vars)!=2)
- stop(sQuote("vars")," must name two variables")
- lags <- as.integer(lags)
- function (y) .Call(
- probe_ccf,
- x=transform(y[vars[1L],,drop=TRUE]),
- y=transform(y[vars[2L],,drop=TRUE]),
- lags=lags,
- corr=corr
- )
-}
-
-probe.marginal <- function (var, ref, order = 3, diff = 1, transform = identity) {
- if (length(var)>1) stop(sQuote("probe.marginal")," is a univariate probe")
- transform <- match.fun(transform)
- setup <- .Call(probe_marginal_setup,transform(ref),order,diff)
- function (y) .Call(
- probe_marginal_solve,
- x=transform(y[var,,drop=TRUE]),
- setup=setup,
- diff=diff
- )
-}
-
-probe.nlar <- function (var, lags, powers, transform = identity) {
- if (length(var)>1) stop(sQuote("probe.nlar")," is a univariate probe")
- transform <- match.fun(transform)
- if (any(lags<1)||any(powers<1))
- stop(sQuote("lags")," and ",sQuote("powers")," must be positive integers")
- if (length(lags)<length(powers)) {
- if (length(lags)>1) stop(sQuote("lags")," must match ",sQuote("powers")," in length, or have length 1")
- lags <- rep(lags,length(powers))
- } else if (length(lags)>length(powers)) {
- if (length(powers)>1) stop(sQuote("powers")," must match ",sQuote("lags")," in length, or have length 1")
- powers <- rep(powers,length(lags))
- }
- lags <- as.integer(lags)
- powers <- as.integer(powers)
- function (y) .Call(
- probe_nlar,
- x=transform(y[var,,drop=TRUE]),
- lags=lags,
- powers=powers
- )
-}
Deleted: branches/mif2/R/bsmc.R
===================================================================
--- branches/mif2/R/bsmc.R 2014-04-17 13:21:28 UTC (rev 922)
+++ branches/mif2/R/bsmc.R 2014-04-17 16:47:05 UTC (rev 923)
@@ -1,450 +0,0 @@
-## 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
-## lower = lower bounds on prior
-## upper = upper bounds on prior
-
-setClass(
- "bsmcd.pomp",
- contains="pomp",
- representation=representation(
- transform="logical",
- post="array",
- prior="array",
- est="character",
- eff.sample.size="numeric",
- smooth="numeric",
- seed="integer",
- nfail="integer",
- cond.log.evidence="numeric",
- log.evidence="numeric",
- weights="numeric"
- )
- )
-
-setGeneric("bsmc",function(object,...)standardGeneric("bsmc"))
-
-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,
- ...) {
-
- 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("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(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
- )
- )
- }
-
- 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(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
- )
- }
-
- 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
-
- 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))
- )
- )
- }
-
- ## 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]
-
- ## 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
- )
-
- ## 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)
-
- 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)
- }
-
- 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),
- 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,
- ...
- )
- }
- )
-
-setMethod("$",signature(x="bsmcd.pomp"),function (x,name) slot(x,name))
-
-bsmc.plot <- function (prior, post, pars, thin, ...) {
- p1 <- sample.int(n=ncol(prior),size=min(thin,ncol(prior)))
- p2 <- sample.int(n=ncol(post),size=min(thin,ncol(post)))
- if (!all(pars%in%rownames(prior))) {
- missing <- which(!(pars%in%rownames(prior)))
- stop("unrecognized parameters: ",paste(sQuote(pars[missing]),collapse=","))
-
- }
- prior <- t(prior[pars,])
- post <- t(post[pars,])
- all <- rbind(prior,post)
- pairs(
- all,
- labels=pars,
- panel=function (x, y, ...) { ## prior, posterior pairwise scatterplot
- op <- par(new=TRUE)
- on.exit(par(op))
- 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[1L]==all[1L,])
- d1 <- density(prior[,i])
- d2 <- density(post[,i])
- usr <- par('usr')
- 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))
- }
- )
-}
-
-setMethod(
- "plot",
- signature(x="bsmcd.pomp"),
- function (x, ..., pars, thin) {
- if (missing(pars)) pars <- names(coef(x,transform=!x at transform))
- if (missing(thin)) thin <- Inf
- bsmc.plot(
- prior=if (x at transform) partrans(x,x at prior,dir="forward") else x at prior,
- post=if (x at transform) partrans(x,x at post,dir="forward") else x at post,
- pars=pars,
- thin=thin,
- ...
- )
- }
- )
Deleted: branches/mif2/R/bsplines.R
===================================================================
--- branches/mif2/R/bsplines.R 2014-04-17 13:21:28 UTC (rev 922)
+++ branches/mif2/R/bsplines.R 2014-04-17 16:47:05 UTC (rev 923)
@@ -1,33 +0,0 @@
-bspline.basis <- function (x, nbasis, degree = 3, names = NULL) {
- y <- .Call(bspline_basis,x,nbasis,degree)
- if (!is.null(names)) {
- if (length(names)==1) {
- nm <- sprintf(names,seq_len(nbasis))
- if (length(unique(nm))!=nbasis)
- nm <- paste(names,seq_len(nbasis),sep=".")
- colnames(y) <- nm
- } else if (length(names)==nbasis) {
- colnames(y) <- names
- } else {
- stop(sQuote("length(names)")," must be either 1 or ",nbasis)
- }
- }
- y
-}
-
-periodic.bspline.basis <- function (x, nbasis, degree = 3, period = 1, names = NULL) {
- y <- .Call(periodic_bspline_basis,x,nbasis,degree,period)
- if (!is.null(names)) {
- if (length(names)==1) {
- nm <- sprintf(names,seq_len(nbasis))
- if (length(unique(nm))!=nbasis)
- nm <- paste(names,seq_len(nbasis),sep=".")
- colnames(y) <- nm
- } else if (length(names)==nbasis) {
- colnames(y) <- names
- } else {
- stop(sQuote("length(names)")," must be either 1 or ",nbasis)
- }
- }
- y
-}
Deleted: branches/mif2/R/builder.R
===================================================================
--- branches/mif2/R/builder.R 2014-04-17 13:21:28 UTC (rev 922)
+++ branches/mif2/R/builder.R 2014-04-17 16:47:05 UTC (rev 923)
@@ -1,283 +0,0 @@
-setClass(
- "pompCode",
- representation=representation(
- type="character",
- slot="character",
- text="character",
- fun="function"
- ),
- prototype=prototype(
- type="ccode",
- slot=character(0),
- text=character(0),
- fun=function(...)stop("function not specified")
- )
- )
-
-
-CCode <- function (text, slot) {
- new("pompCode",type="ccode",slot=as.character(slot))
-}
-
-pompBuilder <- function (data, times, t0, name,
- statenames, paramnames, tcovar, covar,
- rmeasure, dmeasure, step.fn, step.fn.delta.t,
- skeleton, skeleton.type, skelmap.delta.t = 1,
- parameter.transform, parameter.inv.transform,
- ..., link = TRUE) {
- if (!is.data.frame(data)) stop(sQuote("data")," must be a data-frame")
- obsnames <- names(data)
- obsnames <- setdiff(obsnames,times)
- if (!missing(covar)) {
- if (!is.data.frame(covar)) stop(sQuote("covar")," must be a data-frame")
- covarnames <- colnames(covar)
- covarnames <- setdiff(covarnames,tcovar)
- } else {
- covar <- matrix(data=0,nrow=0,ncol=0)
- tcovar <- numeric(0)
- covarnames <- character(0)
- }
- solib <- pompCBuilder(
- name=name,
- statenames=statenames,
- paramnames=paramnames,
- covarnames=covarnames,
- obsnames=obsnames,
- rmeasure=rmeasure,
- dmeasure=dmeasure,
- step.fn=step.fn,
- skeleton=skeleton,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 923
More information about the pomp-commits
mailing list