[Pomp-commits] r407 - in pkg: . R data inst inst/doc man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Nov 5 22:12:09 CET 2010
Author: kingaa
Date: 2010-11-05 22:12:09 +0100 (Fri, 05 Nov 2010)
New Revision: 407
Removed:
pkg/R/pfilter-mif.R
Modified:
pkg/DESCRIPTION
pkg/R/mif-class.R
pkg/R/mif.R
pkg/R/particles-mif.R
pkg/R/pfilter.R
pkg/R/pmcmc.R
pkg/data/dacca.rda
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/data/ou2.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/ChangeLog
pkg/inst/NEWS
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/doc/ou2-first-mif.rda
pkg/inst/doc/ou2-trajmatch.rda
pkg/man/bsmc.Rd
pkg/man/mif-class.Rd
pkg/man/mif-methods.Rd
pkg/man/mif.Rd
pkg/man/pfilter.Rd
pkg/man/pmcmc-class.Rd
pkg/man/pmcmc-methods.Rd
pkg/man/pmcmc.Rd
pkg/man/probe.Rd
pkg/man/probed-pomp-methods.Rd
pkg/man/spect.Rd
pkg/man/traj-match.Rd
pkg/man/trajectory-pomp.Rd
pkg/tests/ou2-forecast.R
pkg/tests/ou2-forecast.Rout.save
pkg/tests/ou2-mif.R
pkg/tests/ou2-mif.Rout.save
pkg/tests/ou2-pmcmc.R
pkg/tests/ou2-pmcmc.Rout.save
Log:
- create new class 'pfilterd.pomp' to hold results of 'pfilter'
- the classes 'mif' and 'pmcmc' now inherit from 'pfilterd.pomp'
- changed the internal representation of 'mif' object (expanded 'alg.pars' list to separate slots)
- the arguments to the generic 'dprior' have changed
- rearrange internal structure of 'mif' algorithm
- rearrange internal structure of 'pmcmc' algorithm
- add 'mif' and 'pmcmc' methods for 'pfilterd.pomp' objects
- remove 'pfilter' method for 'mif' objects: the 'pfilter' method for 'pfilterd.pomp' objects does the job now
- change output of 'pfilter' when 'save.states=TRUE': now states at final time are stored in slot 'last.states'
- it is now allowed to call 'pfilter' with 'save.states=TRUE' on a 'mif' object.
- the internal man pages for various of the classes have yet to be updated
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-11-05 20:51:47 UTC (rev 406)
+++ pkg/DESCRIPTION 2010-11-05 21:12:09 UTC (rev 407)
@@ -8,7 +8,7 @@
Maintainer: Aaron A. King <kingaa at umich.edu>
URL: http://pomp.r-forge.r-project.org
Description: Inference methods for partially-observed Markov processes
-Depends: R(>= 2.9.1), stats, methods, graphics
+Depends: R(>= 2.11.1), stats, methods, graphics
Imports: mvtnorm, subplex, deSolve
License: GPL(>= 2)
LazyLoad: true
@@ -18,7 +18,7 @@
pomp-fun.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 traj-match.R bsmc.R
- mif-class.R particles-mif.R pfilter-mif.R mif.R mif-methods.R compare-mif.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
Modified: pkg/R/mif-class.R
===================================================================
--- pkg/R/mif-class.R 2010-11-05 20:51:47 UTC (rev 406)
+++ pkg/R/mif-class.R 2010-11-05 21:12:09 UTC (rev 407)
@@ -1,20 +1,16 @@
## define the mif class
setClass(
'mif',
- representation(
- ivps = 'character',
- pars = 'character',
- Nmif = 'integer',
- particles = 'function',
- alg.pars = 'list',
- random.walk.sd = 'numeric',
- pred.mean = 'matrix',
- pred.var = 'matrix',
- filter.mean = 'matrix',
- conv.rec = 'matrix',
- eff.sample.size = 'numeric',
- cond.loglik = 'numeric',
- loglik = 'numeric'
- ),
- contains='pomp'
+ contains='pfilterd.pomp',
+ representation=representation(
+ ivps = 'character',
+ pars = 'character',
+ Nmif = 'integer',
+ particles = 'function',
+ var.factor='numeric',
+ ic.lag='integer',
+ cooling.factor='numeric',
+ random.walk.sd = 'numeric',
+ conv.rec = 'matrix'
+ )
)
Modified: pkg/R/mif.R
===================================================================
--- pkg/R/mif.R 2010-11-05 20:51:47 UTC (rev 406)
+++ pkg/R/mif.R 2010-11-05 21:12:09 UTC (rev 407)
@@ -1,6 +1,22 @@
## MIF algorithm functions
-mif.cooling <- function (factor, n) {
+default.pomp.particles.fun <- function (Np, center, sd, ...) {
+ matrix(
+ data=rnorm(
+ n=Np*length(center),
+ mean=center,
+ sd=sd
+ ),
+ nrow=length(center),
+ ncol=Np,
+ dimnames=list(
+ names(center),
+ NULL
+ )
+ )
+}
+
+mif.cooling <- function (factor, n) { # default cooling schedule
alpha <- factor^(n-1)
list(alpha=alpha,gamma=alpha^2)
}
@@ -17,35 +33,29 @@
list(alpha=alpha,gamma=gamma)
}
-mif.internal <- function (object, Nmif = 1,
- start = NULL,
- pars = NULL, ivps = NULL,
- particles = NULL,
- rw.sd = NULL,
- Np = NULL, cooling.factor = NULL, var.factor = NULL, ic.lag = NULL,
- weighted = TRUE, tol = 1e-17, max.fail = 0,
- verbose = FALSE, .ndone = 0) {
- is.mif <- is(object,"mif")
- if (is.null(start)) {
- start <- coef(object)
- if (length(start)==0)
- stop("mif error: ",sQuote("start")," must be specified if ",
- sQuote("coef(object)")," is NULL",
- call.=FALSE)
- } else if (length(start)==0)
- stop("mif error: ",sQuote("start")," must be a named vector",call.=FALSE)
+mif.internal <- function (object, Nmif,
+ start, pars, ivps,
+ particles,
+ rw.sd,
+ Np, cooling.factor, var.factor, ic.lag,
+ weighted, tol, max.fail,
+ verbose, .ndone) {
+
+ if (length(start)==0)
+ stop(
+ "mif error: ",sQuote("start")," must be specified if ",
+ sQuote("coef(object)")," is NULL",
+ call.=FALSE
+ )
start.names <- names(start)
- if (is.null(start.names))
+ if (missing(start.names))
stop("mif error: ",sQuote("start")," must be a named vector",call.=FALSE)
- if (is.null(rw.sd)) {
- if (is.mif)
- rw.sd <- object at random.walk.sd
- else
- stop("mif error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
- }
+ if (missing(rw.sd))
+ stop("mif error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
+
rw.names <- names(rw.sd)
- if (is.null(rw.names) || any(rw.sd<0))
+ if (missing(rw.names) || any(rw.sd<0))
stop("mif error: ",sQuote("rw.sd")," must be a named non-negative numerical vector",call.=FALSE)
if (!all(rw.names%in%start.names))
stop("mif error: all the names of ",sQuote("rw.sd")," must be names of ",sQuote("start"),call.=FALSE)
@@ -53,20 +63,14 @@
if (length(rw.names) == 0)
stop("mif error: ",sQuote("rw.sd")," must have one positive entry for each parameter to be estimated",call.=FALSE)
- if (is.null(pars)) {
- if (is.mif)
- pars <- object at pars
- else
- stop("mif error: ",sQuote("pars")," must be specified",call.=FALSE)
- }
+ if (missing(pars))
+ stop("mif error: ",sQuote("pars")," must be specified",call.=FALSE)
if (length(pars)==0)
stop("mif error: at least one ordinary (non-IVP) parameter must be estimated",call.=FALSE)
- if (is.null(ivps)) {
- if (is.mif)
- ivps <- object at ivps
- else
- stop("mif error: ",sQuote("ivps")," must be specified",call.=FALSE)
- }
+
+ if (missing(ivps))
+ stop("mif error: ",sQuote("ivps")," must be specified",call.=FALSE)
+
if (
!is.character(pars) ||
!is.character(ivps) ||
@@ -101,48 +105,33 @@
rw.sd <- rw.sd[c(pars,ivps)]
rw.names <- names(rw.sd)
- if (is.null(particles)) {
- if (is.mif)
- particles <- object at particles
- else
- stop("mif error: ",sQuote("particles")," must be specified",call.=FALSE)
- }
-
- if (is.null(Np)) {
- if (is.mif)
- Np <- object at alg.pars$Np
- else
- stop("mif error: ",sQuote("Np")," must be specified",call.=FALSE)
- }
+ if (missing(particles))
+ stop("mif error: ",sQuote("particles")," must be specified",call.=FALSE)
+
+ if (missing(Np))
+ stop("mif error: ",sQuote("Np")," must be specified",call.=FALSE)
Np <- as.integer(Np)
if ((length(Np)!=1)||(Np < 1))
stop("mif error: ",sQuote("Np")," must be a positive integer",call.=FALSE)
- if (is.null(ic.lag)) {
- if (is.mif)
- ic.lag <- object at alg.pars$ic.lag
- else
- stop("mif error: ",sQuote("ic.lag")," must be specified",call.=FALSE)
- }
+
+ if (missing(ic.lag))
+ stop("mif error: ",sQuote("ic.lag")," must be specified",call.=FALSE)
ic.lag <- as.integer(ic.lag)
if ((length(ic.lag)!=1)||(ic.lag < 1))
stop("mif error: ",sQuote("ic.lag")," must be a positive integer",call.=FALSE)
- if (is.null(cooling.factor)) {
- if (is.mif)
- cooling.factor <- object at alg.pars$cooling.factor
- else
- stop("mif error: ",sQuote("cooling.factor")," must be specified",call.=FALSE)
- }
+
+ if (missing(cooling.factor))
+ stop("mif error: ",sQuote("cooling.factor")," must be specified",call.=FALSE)
if ((length(cooling.factor)!=1)||(cooling.factor < 0)||(cooling.factor>1))
stop("mif error: ",sQuote("cooling.factor")," must be a number between 0 and 1",call.=FALSE)
- if (is.null(var.factor)) {
- if (is.mif)
- var.factor <- object at alg.pars$var.factor
- else
- stop("mif error: ",sQuote("var.factor")," must be specified",call.=FALSE)
- }
+
+ if (missing(var.factor))
+ stop("mif error: ",sQuote("var.factor")," must be specified",call.=FALSE)
if ((length(var.factor)!=1)||(var.factor < 0))
stop("mif error: ",sQuote("var.factor")," must be a positive number",call.=FALSE)
+ if (missing(Nmif))
+ stop("mif error: ",sQuote("Nmif")," must be specified",call.=FALSE)
Nmif <- as.integer(Nmif)
if (Nmif<0)
stop("mif error: ",sQuote("Nmif")," must be a positive integer",call.=FALSE)
@@ -194,23 +183,13 @@
)
}
- obj <- new(
- "mif",
- as(object,"pomp"),
- ivps=ivps,
- pars=pars,
- Nmif=0L,
- particles=particles,
- alg.pars=list(
- Np=Np,
- cooling.factor=cooling.factor,
- var.factor=var.factor,
- ic.lag=ic.lag
- ),
- random.walk.sd=sigma[rw.names],
- conv.rec=conv.rec
- )
+ obj <- as(object,"pomp")
+ if (Nmif>0)
+ tmp.mif <- new("mif",object,particles=particles,Np=Np) # only needed so that we can use the 'particles' method below
+ else
+ pfp <- obj
+
for (n in seq_len(Nmif)) { # main loop
## compute the cooled sigma
@@ -224,88 +203,84 @@
## initialize the particles' parameter portion...
P <- try(
- particles(obj,Np=Np,center=theta,sd=sigma.n*var.factor),
+ particles(tmp.mif,Np=Np,center=theta,sd=sigma.n*var.factor),
silent=FALSE
)
if (inherits(P,'try-error'))
stop("mif error: error in ",sQuote("particles"),call.=FALSE)
## run the particle filter
- x <- try(
- pfilter.internal(
- object=obj,
- params=P,
- tol=tol,
- max.fail=max.fail,
- pred.mean=(n==Nmif),
- pred.var=(weighted||(n==Nmif)),
- filter.mean=TRUE,
- save.states=FALSE,
- .rw.sd=sigma.n[pars],
- verbose=verbose
- ),
- silent=FALSE
- )
- if (inherits(x,'try-error'))
+ pfp <- try(
+ pfilter.internal(
+ object=obj,
+ params=P,
+ tol=tol,
+ max.fail=max.fail,
+ pred.mean=(n==Nmif),
+ pred.var=(weighted||(n==Nmif)),
+ filter.mean=TRUE,
+ save.states=FALSE,
+ .rw.sd=sigma.n[pars],
+ verbose=verbose
+ ),
+ silent=FALSE
+ )
+ if (inherits(pfp,'try-error'))
stop("mif error: error in ",sQuote("pfilter"),call.=FALSE)
if (weighted) { # MIF update rule
- v <- x$pred.var[pars,,drop=FALSE] # the prediction variance
+ v <- pfp$pred.var[pars,,drop=FALSE] # the prediction variance
v1 <- cool.sched$gamma*(1+var.factor^2)*sigma[pars]^2
- theta.hat <- cbind(theta[pars],x$filter.mean[pars,,drop=FALSE])
+ theta.hat <- cbind(theta[pars],pfp$filter.mean[pars,,drop=FALSE])
theta[pars] <- theta[pars]+colSums(apply(theta.hat,1,diff)/t(v))*v1
} else { # unweighted (flat) average
- theta.hat <- x$filter.mean[pars,,drop=FALSE]
+ theta.hat <- pfp$filter.mean[pars,,drop=FALSE]
theta[pars] <- rowMeans(theta.hat)
}
## update the IVPs using fixed-lag smoothing
- theta[ivps] <- x$filter.mean[ivps,ic.lag]
+ theta[ivps] <- pfp$filter.mean[ivps,ic.lag]
## store a record of this iteration
conv.rec[n+1,-c(1,2)] <- theta
- conv.rec[n,c(1,2)] <- c(x$loglik,x$nfail)
+ conv.rec[n,c(1,2)] <- c(pfp$loglik,pfp$nfail)
if (verbose) cat("MIF iteration ",n," of ",Nmif," completed\n")
}
- coef(obj) <- theta
-
- if (Nmif>0) {
- obj at Nmif <- as.integer(Nmif)
- obj at conv.rec <- conv.rec
- obj at pred.mean <- x$pred.mean
- obj at pred.var <- x$pred.var
- obj at filter.mean <- x$filter.mean
- obj at eff.sample.size <- x$eff.sample.size
- obj at cond.loglik <- x$cond.loglik
- obj at loglik <- x$loglik
- }
-
- obj
+ new(
+ "mif",
+ pfp,
+ params=theta,
+ ivps=ivps,
+ pars=pars,
+ Nmif=Nmif,
+ particles=particles,
+ var.factor=var.factor,
+ ic.lag=ic.lag,
+ cooling.factor=cooling.factor,
+ random.walk.sd=sigma[rw.names],
+ tol=tol,
+ conv.rec=conv.rec
+ )
}
-mif <- function (object, ... )
- stop("function ",sQuote("mif")," is undefined for objects of class ",sQuote(class(object)))
-setGeneric('mif')
+setGeneric('mif',function(object,...)standardGeneric("mif"))
+setGeneric('continue',function(object,...)standardGeneric("continue"))
-continue <- function (object, ... )
- stop("function ",sQuote("continue")," is undefined for objects of class ",sQuote(class(object)))
-setGeneric('continue')
-
setMethod(
"mif",
- "pomp",
+ signature=signature(object="pomp"),
function (object, Nmif = 1,
start,
pars, ivps = character(0),
particles, rw.sd,
Np, ic.lag, var.factor, cooling.factor,
weighted = TRUE, tol = 1e-17, max.fail = 0,
- verbose = getOption("verbose"))
- {
- if (missing(start)) start <- NULL
+ verbose = getOption("verbose"), ...) {
+
+ if (missing(start)) start <- coef(object)
if (missing(rw.sd))
stop("mif error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
if (missing(pars)) {
@@ -320,7 +295,7 @@
stop("mif error: ",sQuote("var.factor")," must be specified",call.=FALSE)
if (missing(cooling.factor))
stop("mif error: ",sQuote("cooling.factor")," must be specified",call.=FALSE)
-
+
if (missing(particles)) { # use default: normal distribution
particles <- function (Np, center, sd, ...) {
matrix(
@@ -348,39 +323,185 @@
call.=FALSE
)
}
-
- mif.internal(object,Nmif=Nmif,start=start,pars=pars,ivps=ivps,particles=particles,
- rw.sd=rw.sd,Np=Np,cooling.factor=cooling.factor,
- var.factor=var.factor,ic.lag=ic.lag,
- weighted=weighted,tol=tol,max.fail=max.fail,
- verbose=verbose,.ndone=0)
+
+ mif.internal(
+ object=object,
+ Nmif=Nmif,
+ start=start,
+ pars=pars,
+ ivps=ivps,
+ particles=particles,
+ rw.sd=rw.sd,
+ Np=Np,
+ cooling.factor=cooling.factor,
+ var.factor=var.factor,
+ ic.lag=ic.lag,
+ weighted=weighted,
+ tol=tol,
+ max.fail=max.fail,
+ verbose=verbose,
+ .ndone=0
+ )
}
)
-
+
setMethod(
"mif",
+ signature=signature(object="pfilterd.pomp"),
+ function (object, Nmif = 1,
+ start,
+ pars, ivps = character(0),
+ particles, rw.sd,
+ Np, ic.lag, var.factor, cooling.factor,
+ weighted = TRUE, tol, max.fail = 0,
+ verbose = getOption("verbose"), ...) {
+
+ if (missing(start)) start <- coef(object)
+ if (missing(rw.sd))
+ stop("mif error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
+ if (missing(pars)) {
+ rw.names <- names(rw.sd)[rw.sd>0]
+ pars <- rw.names[!(rw.names%in%ivps)]
+ }
+ if (missing(Np)) Np <- object at Np
+ if (missing(tol)) tol <- object at tol
+ if (missing(ic.lag))
+ stop("mif error: ",sQuote("ic.lag")," must be specified",call.=FALSE)
+ if (missing(var.factor))
+ stop("mif error: ",sQuote("var.factor")," must be specified",call.=FALSE)
+ if (missing(cooling.factor))
+ stop("mif error: ",sQuote("cooling.factor")," must be specified",call.=FALSE)
+
+ if (missing(particles)) { # use default: normal distribution
+ particles <- default.pomp.particles.fun
+ } else {
+ particles <- match.fun(particles)
+ if (!all(c('Np','center','sd','...')%in%names(formals(particles))))
+ stop(
+ "mif error: ",
+ sQuote("particles"),
+ " must be a function of prototype ",
+ sQuote("particles(Np,center,sd,...)"),
+ call.=FALSE
+ )
+ }
+
+ mif.internal(
+ object=as(object,"pomp"),
+ Nmif=Nmif,
+ start=start,
+ pars=pars,
+ ivps=ivps,
+ particles=particles,
+ rw.sd=rw.sd,
+ Np=Np,
+ cooling.factor=cooling.factor,
+ var.factor=var.factor,
+ ic.lag=ic.lag,
+ weighted=weighted,
+ tol=tol,
+ max.fail=max.fail,
+ verbose=verbose,
+ .ndone=0
+ )
+ }
+ )
+
+setMethod(
"mif",
- function (object, Nmif, ...)
- {
+ signature=signature(object="mif"),
+ function (object, Nmif,
+ start,
+ pars, ivps,
+ particles, rw.sd,
+ Np, ic.lag, var.factor, cooling.factor,
+ weighted = TRUE, tol, max.fail = 0,
+ verbose = getOption("verbose"), ...) {
+
if (missing(Nmif)) Nmif <- object at Nmif
- mif.internal(object,Nmif=Nmif,...)
+ if (missing(start)) start <- coef(object)
+ if (missing(pars)) pars <- object at pars
+ if (missing(ivps)) ivps <- object at ivps
+ if (missing(particles)) particles <- object at particles
+ if (missing(rw.sd)) rw.sd <- object at random.walk.sd
+ if (missing(Np)) Np <- object at Np
+ if (missing(ic.lag)) ic.lag <- object at ic.lag
+ if (missing(var.factor)) var.factor <- object at var.factor
+ if (missing(cooling.factor)) cooling.factor <- object at cooling.factor
+ if (missing(tol)) tol <- object at tol
+
+ mif.internal(
+ object=as(object,"pomp"),
+ Nmif=Nmif,
+ start=start,
+ pars=pars,
+ ivps=ivps,
+ particles=particles,
+ rw.sd=rw.sd,
+ Np=Np,
+ cooling.factor=cooling.factor,
+ var.factor=var.factor,
+ ic.lag=ic.lag,
+ weighted=weighted,
+ tol=tol,
+ max.fail=max.fail,
+ verbose=verbose,
+ .ndone=0
+ )
}
)
setMethod(
'continue',
- 'mif',
- function (object, Nmif = 1, ...) {
+ signature=signature(object='mif'),
+ function (object, Nmif = 1,
+ start,
+ pars, ivps,
+ particles, rw.sd,
+ Np, ic.lag, var.factor, cooling.factor,
+ weighted = TRUE, tol, max.fail = 0,
+ verbose = getOption("verbose"), ...) {
+
ndone <- object at Nmif
- obj <- mif.internal(object,Nmif=Nmif,.ndone=ndone,...)
+ if (missing(start)) start <- coef(object)
+ if (missing(pars)) pars <- object at pars
+ if (missing(ivps)) ivps <- object at ivps
+ if (missing(particles)) particles <- object at particles
+ if (missing(rw.sd)) rw.sd <- object at random.walk.sd
+ if (missing(Np)) Np <- object at Np
+ if (missing(ic.lag)) ic.lag <- object at ic.lag
+ if (missing(var.factor)) var.factor <- object at var.factor
+ if (missing(cooling.factor)) cooling.factor <- object at cooling.factor
+ if (missing(tol)) tol <- object at tol
+
+ obj <- mif.internal(
+ object=as(object,"pomp"),
+ Nmif=Nmif,
+ start=start,
+ pars=pars,
+ ivps=ivps,
+ particles=particles,
+ rw.sd=rw.sd,
+ Np=Np,
+ cooling.factor=cooling.factor,
+ var.factor=var.factor,
+ ic.lag=ic.lag,
+ weighted=weighted,
+ tol=tol,
+ max.fail=max.fail,
+ verbose=verbose,
+ .ndone=ndone
+ )
+
object at conv.rec[ndone+1,c('loglik','nfail')] <- obj at conv.rec[1,c('loglik','nfail')]
obj at conv.rec <- rbind(
object at conv.rec,
obj at conv.rec[-1,colnames(object at conv.rec)]
)
obj at Nmif <- as.integer(ndone+Nmif)
+
obj
}
)
@@ -404,7 +525,7 @@
pp <- coef(object)
idx <- !(names(pp)%in%names(pd))
if (any(idx)) pd <- cbind(pd,as.list(pp[idx]))
-
+
ans <- vector(mode="list",length=nrow(pd))
for (k in seq_len(nrow(pd))) {
ans[[k]] <- list(
Modified: pkg/R/particles-mif.R
===================================================================
--- pkg/R/particles-mif.R 2010-11-05 20:51:47 UTC (rev 406)
+++ pkg/R/particles-mif.R 2010-11-05 21:12:09 UTC (rev 407)
@@ -1,8 +1,6 @@
## draw a set of Np particles from the user-specified distribution
-particles <- function (object, ...)
- stop("function ",sQuote("particles")," is undefined for objects of class ",sQuote(class(object)))
-setGeneric('particles')
+setGeneric('particles',function(object,...)standardGeneric("particles"))
setMethod(
"particles",
Deleted: pkg/R/pfilter-mif.R
===================================================================
--- pkg/R/pfilter-mif.R 2010-11-05 20:51:47 UTC (rev 406)
+++ pkg/R/pfilter-mif.R 2010-11-05 21:12:09 UTC (rev 407)
@@ -1,33 +0,0 @@
-setMethod(
- "pfilter",
- "mif",
- function (object, params, Np,
- tol = 1e-17, max.fail = 0,
- pred.mean = FALSE,
- pred.var = FALSE,
- filter.mean = FALSE,
- ...) {
- if (missing(params))
- params <- coef(object)
- if (missing(Np))
- Np <- object at alg.pars$Np
- if ("save.states"%in%names(list(...)))
- stop(
- "pfilter error: you cannot set ",sQuote("save.states"),
- " when ",sQuote("object")," is of class mif",
- call.=FALSE
- )
- pfilter(
- object=as(object,"pomp"),
- params=params,
- Np=Np,
- tol=tol,
- max.fail=max.fail,
- pred.mean=pred.mean,
- pred.var=pred.var,
- filter.mean=filter.mean,
- save.states=FALSE,
- ...
- )
- }
- )
Modified: pkg/R/pfilter.R
===================================================================
--- pkg/R/pfilter.R 2010-11-05 20:51:47 UTC (rev 406)
+++ pkg/R/pfilter.R 2010-11-05 21:12:09 UTC (rev 407)
@@ -1,18 +1,31 @@
## particle filtering codes
-## generic particle filter
-pfilter <- function (object, ...)
- stop("function ",sQuote("pfilter")," is undefined for objects of class ",sQuote(class(object)))
-setGeneric('pfilter')
+setClass(
+ "pfilterd.pomp",
+ contains="pomp",
+ representation=representation(
+ pred.mean="array",
+ pred.var="array",
+ filter.mean="array",
+ eff.sample.size="numeric",
+ cond.loglik="numeric",
+ last.states="array",
+ seed="integer",
+ Np="integer",
+ tol="numeric",
+ nfail="integer",
+ loglik="numeric"
+ )
+ )
## question: when pfilter.internal is called by mif, do we need to compute the prediction means and variances of the state variables each time, or only at the end?
-## question: how much efficiency would be realized by eliminating the calls to 'apply' with something else?
pfilter.internal <- function (object, params, Np,
tol, max.fail,
pred.mean, pred.var, filter.mean,
.rw.sd, seed, verbose,
save.states) {
+
if (missing(seed)) seed <- NULL
if (!is.null(seed)) {
if (!exists(".Random.seed",where=.GlobalEnv)) { # need to initialize the RNG
@@ -22,14 +35,15 @@
set.seed(seed)
}
- if (missing(params)) {
- params <- coef(object)
- if (length(params)==0) {
- stop(sQuote("pfilter")," error: ",sQuote("params")," must be supplied",call.=FALSE)
- }
- }
+ if (length(params)==0)
+ stop(sQuote("pfilter")," error: ",sQuote("params")," must be specified",call.=FALSE)
+
if (missing(Np))
Np <- NCOL(params)
+
+ if (missing(tol))
+ stop(sQuote("pfilter")," error: ",sQuote("tol")," must be specified",call.=FALSE)
+
times <- time(object,t0=TRUE)
ntimes <- length(times)-1
if (is.null(dim(params))) {
@@ -50,15 +64,14 @@
x <- init.state(object,params=params)
statenames <- rownames(x)
nvars <- nrow(x)
- if (save.states) {
+ if (save.states)
xparticles <- array(
data=NA,
dim=c(nvars,Np,ntimes),
dimnames=list(statenames,NULL,NULL)
)
- } else {
- xparticles <- NULL
- }
+ else
+ xparticles <- array(dim=c(0,0,0))
random.walk <- !missing(.rw.sd)
if (random.walk) {
@@ -81,10 +94,6 @@
nfail <- 0
npars <- length(rw.names)
- pred.m <- NULL
- pred.v <- NULL
- filt.m <- NULL
-
## set up storage for prediction means, variances, etc.
if (pred.mean)
pred.m <- matrix(
@@ -93,7 +102,9 @@
ncol=ntimes,
dimnames=list(c(statenames,rw.names),NULL)
)
-
+ else
+ pred.m <- array(dim=c(0,0))
+
if (pred.var)
pred.v <- matrix(
data=0,
@@ -101,24 +112,26 @@
ncol=ntimes,
dimnames=list(c(statenames,rw.names),NULL)
)
+ else
+ pred.v <- array(dim=c(0,0))
- if (filter.mean) {
- if (random.walk) {
+ if (filter.mean)
+ if (random.walk)
filt.m <- matrix(
data=0,
nrow=nvars+length(paramnames),
ncol=ntimes,
dimnames=list(c(statenames,paramnames),NULL)
)
- } else {
+ else
filt.m <- matrix(
data=0,
nrow=nvars,
ncol=ntimes,
dimnames=list(statenames,NULL)
)
- }
- }
+ else
+ filt.m <- array(dim=c(0,0))
for (nt in seq_len(ntimes)) {
@@ -232,24 +245,29 @@
seed <- save.seed
}
- list(
- pred.mean=pred.m,
- pred.var=pred.v,
- filter.mean=filt.m,
- eff.sample.size=eff.sample.size,
- cond.loglik=loglik,
- states=xparticles,
- seed=seed,
- Np=Np,
- tol=tol,
- nfail=nfail,
- loglik=sum(loglik)
- )
+ new(
+ "pfilterd.pomp",
+ object,
+ pred.mean=pred.m,
+ pred.var=pred.v,
+ filter.mean=filt.m,
+ eff.sample.size=eff.sample.size,
+ cond.loglik=loglik,
+ last.states=xparticles,
+ seed=as.integer(seed),
+ Np=as.integer(Np),
+ tol=tol,
+ nfail=as.integer(nfail),
+ loglik=sum(loglik)
+ )
}
+## generic particle filter
+setGeneric("pfilter",function(object,...)standardGeneric("pfilter"))
+
setMethod(
"pfilter",
- "pomp",
+ signature=signature(object="pomp"),
function (object, params, Np,
tol = 1e-17,
max.fail = 0,
@@ -260,6 +278,7 @@
seed = NULL,
verbose = getOption("verbose"),
...) {
+ if (missing(params)) params <- coef(object)
pfilter.internal(
object=object,
params=params,
@@ -275,3 +294,39 @@
)
}
)
+
+setMethod(
+ "pfilter",
+ signature=signature(object="pfilterd.pomp"),
+ function (object, params, Np,
+ tol,
+ max.fail = 0,
+ pred.mean = FALSE,
+ pred.var = FALSE,
+ filter.mean = FALSE,
+ save.states = FALSE,
+ seed = NULL,
+ verbose = getOption("verbose"),
+ ...) {
+ if (missing(params)) params <- coef(object)
+ if (missing(Np)) Np <- object at Np
+ if (missing(tol)) tol <- object at tol
+ pfilter.internal(
+ object=as(object,"pomp"),
+ params=params,
+ Np=Np,
+ tol=tol,
+ max.fail=max.fail,
+ pred.mean=pred.mean,
+ pred.var=pred.var,
+ filter.mean=filter.mean,
+ save.states=save.states,
+ seed=seed,
+ verbose=verbose
+ )
+ }
+ )
+
+setMethod("$",signature(x="pfilterd.pomp"),function (x,name) slot(x,name))
+setMethod("logLik",signature(object="pfilterd.pomp"),function(object,...)object at loglik)
+
Modified: pkg/R/pmcmc.R
===================================================================
--- pkg/R/pmcmc.R 2010-11-05 20:51:47 UTC (rev 406)
+++ pkg/R/pmcmc.R 2010-11-05 21:12:09 UTC (rev 407)
@@ -1,70 +1,53 @@
## define the pmcmc class
setClass(
'pmcmc',
+ contains='pfilterd.pomp',
representation(
pars = 'character',
Nmcmc = 'integer',
dprior = 'function',
hyperparams = 'list',
- Np = 'integer',
random.walk.sd = 'numeric',
- filter.mean = 'matrix',
conv.rec = 'matrix',
- eff.sample.size = 'numeric',
- cond.loglik = 'numeric',
- loglik = 'numeric',
log.prior = 'numeric'
- ),
- contains='pomp'
+ )
)
## PMCMC algorithm functions
-pmcmc <- function (object, ... )
- stop("function ",sQuote("pmcmc")," is undefined for objects of class ",sQuote(class(object)))
-setGeneric('pmcmc')
+setGeneric('pmcmc',function(object,...)standardGeneric("pmcmc"))
-dprior <- function (object, params, log)
- stop("function ",sQuote("dprior")," is undefined for objects of class ",sQuote(class(object)))
-setGeneric('dprior')
+setGeneric('dprior', function(object,...)standardGeneric("dprior"))
setMethod(
"dprior",
- "pmcmc",
- function (object, params, log = FALSE) {
+ signature=signature(object="pmcmc"),
+ function (object, params, log = FALSE, ...) {
do.call(object at dprior,list(params=params,hyperparams=object at hyperparams,log=log))
}
)
-pmcmc.internal <- function (object, Nmcmc = 1,
- start = NULL,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 407
More information about the pomp-commits
mailing list