[Pomp-commits] r806 - in pkg/pomp: R inst man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jan 10 00:01:27 CET 2013
Author: kingaa
Date: 2013-01-10 00:01:26 +0100 (Thu, 10 Jan 2013)
New Revision: 806
Added:
pkg/pomp/tests/ou2-mif2.R
pkg/pomp/tests/ou2-mif2.Rout.save
Modified:
pkg/pomp/R/mif-class.R
pkg/pomp/R/mif.R
pkg/pomp/R/pfilter.R
pkg/pomp/R/pmcmc.R
pkg/pomp/inst/NEWS
pkg/pomp/man/mif-class.Rd
pkg/pomp/man/mif.Rd
pkg/pomp/man/pfilter.Rd
Log:
- integrate MIF2
Modified: pkg/pomp/R/mif-class.R
===================================================================
--- pkg/pomp/R/mif-class.R 2013-01-09 22:59:28 UTC (rev 805)
+++ pkg/pomp/R/mif-class.R 2013-01-09 23:01:26 UTC (rev 806)
@@ -11,6 +11,8 @@
var.factor='numeric',
ic.lag='integer',
cooling.factor='numeric',
+ cooling.fraction='numeric',
+ method='character',
random.walk.sd = 'numeric',
conv.rec = 'matrix'
)
Modified: pkg/pomp/R/mif.R
===================================================================
--- pkg/pomp/R/mif.R 2013-01-09 22:59:28 UTC (rev 805)
+++ pkg/pomp/R/mif.R 2013-01-09 23:01:26 UTC (rev 806)
@@ -16,11 +16,18 @@
)
}
-mif.cooling <- function (factor, n) { # default cooling schedule
+mif.cooling <- function (factor, n) { # default geometric cooling schedule
alpha <- factor^(n-1)
list(alpha=alpha,gamma=alpha^2)
}
+mif2.cooling <- function (frac, nt, m, n) { # cooling schedule for mif2
+ ## frac is the fraction of cooling after 50 iterations
+ cooling.scalar <- (50*n*frac-1)/(1-frac)
+ alpha <- (1+cooling.scalar)/(cooling.scalar+nt+n*(m-1))
+ list(alpha=alpha)
+}
+
powerlaw.cooling <- function (init = 1, delta = 0.1, eps = (1-delta)/2, n) {
m <- init
if (n <= m) { # linear cooling regime
@@ -36,45 +43,45 @@
mif.internal <- function (object, Nmif,
start, pars, ivps,
particles,
- rw.sd,
+ rw.sd,
Np, cooling.factor, var.factor, ic.lag,
+ cooling.fraction,
method,
tol, max.fail,
- verbose, transform, .ndone) {
-
+ verbose, transform, .ndone = 0,
+ paramMatrix) {
+
transform <- as.logical(transform)
-
+
if (length(start)==0)
stop(
"mif error: ",sQuote("start")," must be specified if ",
sQuote("coef(object)")," is NULL",
call.=FALSE
)
-
+
if (transform)
start <- partrans(object,start,dir="inverse")
-
+
start.names <- names(start)
- if (missing(start.names))
+ if (is.null(start.names))
stop("mif error: ",sQuote("start")," must be a named vector",call.=FALSE)
-
+
if (missing(rw.sd))
stop("mif error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
-
+
rw.names <- names(rw.sd)
- if (missing(rw.names) || any(rw.sd<0))
+ if (is.null(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)
rw.names <- names(rw.sd[rw.sd>0])
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 (missing(pars))
- stop("mif error: ",sQuote("pars")," must be specified",call.=FALSE)
- if (missing(ivps))
- stop("mif error: ",sQuote("ivps")," must be specified",call.=FALSE)
-
+
+ if (missing(pars)) stop("mif error: ",sQuote("pars")," 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) ||
@@ -94,7 +101,7 @@
sQuote("rw.sd"),
call.=FALSE
)
-
+
if (!all(rw.names%in%c(pars,ivps))) {
extra.rws <- rw.names[!(rw.names%in%c(pars,ivps))]
warning(
@@ -108,13 +115,12 @@
}
rw.sd <- rw.sd[c(pars,ivps)]
rw.names <- names(rw.sd)
-
+
if (missing(particles))
stop("mif error: ",sQuote("particles")," must be specified",call.=FALSE)
-
+
ntimes <- length(time(object))
- if (missing(Np))
- stop("mif error: ",sQuote("Np")," must be specified",call.=FALSE)
+ if (missing(Np)) stop("mif error: ",sQuote("Np")," must be specified",call.=FALSE)
if (is.function(Np)) {
Np <- try(
vapply(seq.int(from=0,to=ntimes,by=1),Np,numeric(1)),
@@ -132,11 +138,10 @@
if (!is.numeric(Np))
stop(sQuote("Np")," must be a number, a vector of numbers, or a function")
Np <- as.integer(Np)
-
- if (missing(ic.lag))
- 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))
+ if ((length(ic.lag)!=1)||(ic.lag<1))
stop("mif error: ",sQuote("ic.lag")," must be a positive integer",call.=FALSE)
if (ic.lag>ntimes) {
warning(
@@ -155,17 +160,32 @@
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 (method=="mif2") {
+ if (missing(cooling.fraction) || is.na(cooling.fraction))
+ stop("mif error: ",sQuote("cooling.fraction")," must be specified for method = ",sQuote("mif2"),call.=FALSE)
+ cooling.fraction <- as.numeric(cooling.fraction)
+ if ((length(cooling.fraction)!=1)||(cooling.fraction<0)||(cooling.fraction>1))
+ stop("mif error: ",sQuote("cooling.fraction")," must be a number between 0 and 1",call.=FALSE)
+ if (!missing(cooling.factor) && !(is.na(cooling.factor)))
+ warning(sQuote("cooling.factor")," ignored for method = ",sQuote("mif2"),call.=FALSE)
+ cooling.factor <- as.numeric(NA)
+ } else {
+ if (missing(cooling.factor) || is.na(cooling.factor))
+ stop("mif error: ",sQuote("cooling.factor")," must be specified",call.=FALSE)
+ cooling.factor <- as.numeric(cooling.factor)
+ 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 (!missing(cooling.fraction) && !(is.na(cooling.fraction)))
+ warning(sQuote("cooling.fraction")," ignored for method != ",sQuote("mif2"),call.=FALSE)
+ cooling.fraction <- as.numeric(NA)
+ }
+
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)
@@ -173,15 +193,15 @@
stop("mif error: ",sQuote("Nmif")," must be a positive integer",call.=FALSE)
theta <- start
-
+
sigma <- rep(0,length(start))
names(sigma) <- start.names
-
+
rw.sd <- rw.sd[c(pars,ivps)]
rw.names <- names(rw.sd)
-
+
sigma[rw.names] <- rw.sd
-
+
conv.rec <- matrix(
data=NA,
nrow=Nmif+1,
@@ -192,7 +212,7 @@
)
)
conv.rec[1,] <- c(NA,NA,theta)
-
+
if (!all(is.finite(theta[c(pars,ivps)]))) {
stop(
sQuote("mif")," cannot estimate non-finite parameters.\n",
@@ -206,85 +226,102 @@
}
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
+
+ if (Nmif>0) {
+ tmp.mif <- new("mif",object,particles=particles,Np=Np)
+ } else {
pfp <- obj
+ }
+
+ have.parmat <- !(missing(paramMatrix) || length(paramMatrix)==0)
- for (n in seq_len(Nmif)) { # main loop
+ for (n in seq_len(Nmif)) { ## iterate the filtering
- ## compute the cooled sigma
+ ## get the intensity of artificial noise from the cooling schedule
cool.sched <- try(
- mif.cooling(cooling.factor,.ndone+n),
+ switch(
+ method,
+ mif2=mif2.cooling(frac=cooling.fraction,nt=1,m=.ndone+n,n=ntimes),
+ mif.cooling(factor=cooling.factor,n=.ndone+n)
+ ),
silent=FALSE
)
- if (inherits(cool.sched,'try-error'))
+ if (inherits(cool.sched,"try-error"))
stop("mif error: cooling schedule error",call.=FALSE)
sigma.n <- sigma*cool.sched$alpha
-
- ## initialize the particles' parameter portion...
+
+ ## initialize the parameter portions of the particles
P <- try(
- particles(tmp.mif,Np=Np[1],center=theta,sd=sigma.n*var.factor),
- silent=FALSE
+ particles(
+ tmp.mif,
+ Np=Np[1],
+ center=theta,
+ sd=sigma.n*var.factor
+ ),
+ silent = FALSE
)
- if (inherits(P,'try-error'))
+ if (inherits(P,"try-error"))
stop("mif error: error in ",sQuote("particles"),call.=FALSE)
- ## run the particle filter
+ if ((method=="mif2") && ((n>1) || have.parmat)) {
+ ## use pre-existing particle matrix
+ P[pars,] <- paramMatrix[pars,]
+ }
+
pfp <- try(
pfilter.internal(
object=obj,
- params=P,
+ params=P,
tol=tol,
max.fail=max.fail,
pred.mean=(n==Nmif),
pred.var=((method=="mif")||(n==Nmif)),
filter.mean=TRUE,
- save.states=FALSE,
- save.params=FALSE,
+ cooling.fraction=cooling.fraction,
+ cooling.m=.ndone+n,
+ .mif2=(method=="mif2"),
.rw.sd=sigma.n[pars],
- verbose=verbose,
- transform=transform
+ .transform=transform,
+ save.states=FALSE,
+ save.params=FALSE,
+ verbose=verbose
),
silent=FALSE
)
- if (inherits(pfp,'try-error'))
+ if (inherits(pfp,"try-error"))
stop("mif error: error in ",sQuote("pfilter"),call.=FALSE)
+ ## update parameters
switch(
method,
- mif={ # mif update rule
- 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],pfp$filter.mean[pars,,drop=FALSE])
- theta[pars] <- theta[pars]+colSums(apply(theta.hat,1,diff)/t(v))*v1
+ mif={ # original Ionides et al. (2006) average
+ v <- pfp at pred.var[pars,,drop=FALSE]
+ v1 <- cool.sched$gamma * (1 + var.factor^2) * sigma[pars]^2
+ theta.hat <- cbind(theta[pars],pfp at filter.mean[pars,,drop=FALSE])
+ theta[pars] <- theta[pars] + colSums(apply(theta.hat,1,diff)/t(v))*v1
},
- unweighted={ # unweighted (flat) average
- theta.hat <- pfp$filter.mean[pars,,drop=FALSE]
- theta[pars] <- rowMeans(theta.hat)
+ unweighted={ # unweighted average
+ theta[pars] <- rowMeans(pfp at filter.mean[pars,,drop=FALSE])
},
- fp={
- theta.hat <- pfp$filter.mean[pars,ntimes,drop=FALSE]
- theta[pars] <- theta.hat
+ fp={ # fixed-point iteration
+ theta[pars] <- pfp at filter.mean[pars,ntimes,drop=FALSE]
},
- stop("unrecognized method",sQuote(method))
+ mif2={ # "efficient" iterated filtering
+ paramMatrix <- pfp at paramMatrix
+ theta[pars] <- rowMeans(paramMatrix[pars,,drop=FALSE])
+ },
+ stop("unrecognized method ",sQuote(method))
)
-
- ## update the IVPs using fixed-lag smoothing
- theta[ivps] <- pfp$filter.mean[ivps,ic.lag]
-
- ## store a record of this iteration
+ theta[ivps] <- pfp at filter.mean[ivps,ic.lag]
conv.rec[n+1,-c(1,2)] <- theta
- conv.rec[n,c(1,2)] <- c(pfp$loglik,pfp$nfail)
-
+ conv.rec[n,c(1,2)] <- c(pfp at loglik,pfp at nfail)
+
if (verbose) cat("MIF iteration ",n," of ",Nmif," completed\n")
+
+ } ### end of main loop
- }
-
## back transform the parameter estimate if necessary
- if (transform)
- theta <- partrans(pfp,theta,dir="forward")
+ if (transform) theta <- partrans(pfp,theta,dir="forward")
new(
"mif",
@@ -297,10 +334,13 @@
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
+ conv.rec=conv.rec,
+ method=method,
+ cooling.factor=cooling.factor,
+ cooling.fraction=cooling.fraction,
+ paramMatrix=if (method=="mif2") paramMatrix else array(data=numeric(0),dim=c(0,0))
)
}
@@ -314,13 +354,15 @@
pars, ivps = character(0),
particles, rw.sd,
Np, ic.lag, var.factor, cooling.factor,
- weighted, method = c("mif","unweighted","fp"),
+ cooling.fraction,
+ method = c("mif","unweighted","fp","mif2"),
tol = 1e-17, max.fail = 0,
verbose = getOption("verbose"),
transform = FALSE, ...) {
-
+
transform <- as.logical(transform)
-
+ method <- match.arg(method)
+
if (missing(start)) start <- coef(object)
if (missing(rw.sd))
stop("mif error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
@@ -330,30 +372,12 @@
}
if (missing(Np))
stop("mif error: ",sQuote("Np")," must be specified",call.=FALSE)
- if (missing(ic.lag))
- stop("mif error: ",sQuote("ic.lag")," must be specified",call.=FALSE)
+ if (missing(ic.lag) && length(ivps)>0)
+ stop("mif error: ",sQuote("ic.lag")," must be specified if ",sQuote("ivps")," are",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)
- method <- match.arg(method)
- if (!missing(weighted)) {
- warning(sQuote("mif")," warning: ",sQuote("weighted")," flag is deprecated, use ",sQuote("method"))
- if (weighted) {
- if (method!="mif") {
- warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
- }
- method <- "mif"
- } else {
- if (method!="unweighted") {
- warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
- }
- method <- "unweighted"
- }
- }
-
- if (missing(particles)) { # use default: normal distribution
+ if (missing(particles)) { # use default: normal distribution
particles <- default.pomp.particles.fun
} else {
particles <- match.fun(particles)
@@ -377,16 +401,16 @@
rw.sd=rw.sd,
Np=Np,
cooling.factor=cooling.factor,
+ cooling.fraction=cooling.fraction,
var.factor=var.factor,
ic.lag=ic.lag,
method=method,
tol=tol,
max.fail=max.fail,
verbose=verbose,
- transform=transform,
- .ndone=0
+ transform=transform
)
-
+
}
)
@@ -395,82 +419,22 @@
"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, method = c("mif","unweighted","fp"),
- tol = 1e-17, max.fail = 0,
+ Np, tol, max.fail = 0,
verbose = getOption("verbose"),
- transform = FALSE, ...) {
-
- transform <- as.logical(transform)
-
- 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)
- method <- match.arg(method)
- if (!missing(weighted)) {
- warning(sQuote("mif")," warning: ",sQuote("weighted")," flag is deprecated, use ",sQuote("method"))
- if (weighted) {
- if (method!="mif") {
- warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
- }
- method <- "mif"
- } else {
- if (method!="unweighted") {
- warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
- }
- method <- "unweighted"
- }
- }
-
- 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,
- method=method,
- tol=tol,
- max.fail=max.fail,
- verbose=verbose,
- transform=transform,
- .ndone=0
- )
+ mif(
+ object=as(object,"pomp"),
+ Nmif=Nmif,
+ Np=Np,
+ tol=tol,
+ max.fail=max.fail,
+ verbose=verbose,
+ ...
+ )
}
)
@@ -482,60 +446,49 @@
pars, ivps,
particles, rw.sd,
Np, ic.lag, var.factor, cooling.factor,
- weighted, method = c("mif","unweighted","fp"),
- tol = 1e-17, max.fail = 0,
+ cooling.fraction,
+ method,
+ tol, max.fail = 0,
verbose = getOption("verbose"),
transform, ...) {
-
+
if (missing(Nmif)) Nmif <- object at 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
+ if (missing(cooling.fraction)) cooling.fraction <- object at cooling.fraction
+ if (missing(method)) method <- object at method
if (missing(transform)) transform <- object at transform
transform <- as.logical(transform)
- method <- match.arg(method)
- if (!missing(weighted)) {
- warning(sQuote("mif")," warning: ",sQuote("weighted")," flag is deprecated, use ",sQuote("method"))
- if (weighted) {
- if (method!="mif") {
- warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
- }
- method <- "mif"
- } else {
- if (method!="unweighted") {
- warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
- }
- method <- "unweighted"
- }
- }
+ if (missing(Np)) Np <- object at Np
+ 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,
- method=method,
- tol=tol,
- max.fail=max.fail,
- verbose=verbose,
- transform=transform,
- .ndone=0
- )
+ mif(
+ object=as(object,"pomp"),
+ Nmif=Nmif,
+ start=start,
+ pars=pars,
+ ivps=ivps,
+ particles=particles,
+ rw.sd=rw.sd,
+ Np=Np,
+ cooling.factor=cooling.factor,
+ cooling.fraction=cooling.fraction,
+ var.factor=var.factor,
+ ic.lag=ic.lag,
+ method=method,
+ tol=tol,
+ max.fail=max.fail,
+ verbose=verbose,
+ transform=transform,
+ ...
+ )
}
)
@@ -543,114 +496,74 @@
'continue',
signature=signature(object='mif'),
function (object, Nmif = 1,
- start,
- pars, ivps,
- particles, rw.sd,
- Np, ic.lag, var.factor, cooling.factor,
- weighted, method = c("mif","unweighted","fp"),
- tol = 1e-17, max.fail = 0,
+ max.fail = 0,
verbose = getOption("verbose"),
- transform, ...) {
-
+ ...) {
+
ndone <- object at 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
- if (missing(transform)) transform <- object at transform
- transform <- as.logical(transform)
-
- method <- match.arg(method)
- if (!missing(weighted)) {
- warning(sQuote("mif")," warning: ",sQuote("weighted")," flag is deprecated, use ",sQuote("method"))
- if (weighted) {
- if (method!="mif") {
- warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
- }
- method <- "mif"
- } else {
- if (method!="unweighted") {
- warning(sQuote("mif")," warning: use of ",sQuote("weighted")," argument overrides choice of ",sQuote("method"))
- }
- method <- "unweighted"
- }
- }
-
- 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,
- method=method,
- tol=tol,
- max.fail=max.fail,
- verbose=verbose,
- transform=transform,
- .ndone=ndone
- )
-
+
+ obj <- mif(
+ object=object,
+ Nmif=Nmif,
+ max.fail=max.fail,
+ verbose=verbose,
+ .ndone=ndone,
+ paramMatrix=object at paramMatrix,
+ ...
+ )
+
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
}
)
mif.profileDesign <- function (object, profile, lower, upper, nprof, ivps,
- rw.sd, Np, ic.lag, var.factor, cooling.factor, ...)
- {
- if (missing(profile)) profile <- list()
- if (missing(lower)) lower <- numeric(0)
- if (missing(upper)) upper <- lower
- if (length(lower)!=length(upper))
- stop(sQuote("lower")," and ",sQuote("upper")," must be of the same length")
- pars <- names(lower)
- if (missing(ivps)) ivps <- character(0)
- Np <- as.integer(Np)
-
- pd <- do.call(profileDesign,c(profile,list(lower=lower,upper=upper,nprof=nprof)))
-
- object <- as(object,"pomp")
-
- 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(
- mf=mif(
- object,
- Nmif=0,
- start=unlist(pd[k,]),
- pars=pars,
- ivps=ivps,
- rw.sd=rw.sd,
- Np=Np,
- ic.lag=ic.lag,
- var.factor=var.factor,
- cooling.factor=cooling.factor,
- ...
- )
+ rw.sd, Np, ic.lag, var.factor, cooling.factor,option, cooling.fraction, paramMatrix, ...)
+{
+ if (missing(profile)) profile <- list()
+ if (missing(lower)) lower <- numeric(0)
+ if (missing(upper)) upper <- lower
+ if (length(lower)!=length(upper))
+ stop(sQuote("lower")," and ",sQuote("upper")," must be of the same length")
+ pars <- names(lower)
+ if (missing(ivps)) ivps <- character(0)
+ Np <- as.integer(Np)
+
+ pd <- do.call(profileDesign,c(profile,list(lower=lower,upper=upper,nprof=nprof)))
+
+ object <- as(object,"pomp")
+
+ 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(
+ mf=mif(
+ object,
+ Nmif=0,
+ start=unlist(pd[k,]),
+ pars=pars,
+ ivps=ivps,
+ rw.sd=rw.sd,
+ Np=Np,
+ ic.lag=ic.lag,
+ var.factor=var.factor,
+ cooling.factor=cooling.factor,
+ option=option,
+ cooling.fraction=cooling.fraction,
+ paramMatrix=paramMatrix,
+ ...
)
- }
-
- ans
+ )
}
+
+ ans
+}
Modified: pkg/pomp/R/pfilter.R
===================================================================
--- pkg/pomp/R/pfilter.R 2013-01-09 22:59:28 UTC (rev 805)
+++ pkg/pomp/R/pfilter.R 2013-01-09 23:01:26 UTC (rev 806)
@@ -35,21 +35,16 @@
)
)
-## 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?
-
pfilter.internal <- function (object, params, Np,
tol, max.fail,
pred.mean, pred.var, filter.mean,
- cooling.fraction, cooling.m, mif2 = FALSE,
+ cooling.fraction, cooling.m, .mif2 = FALSE,
.rw.sd, seed, verbose,
save.states, save.params,
- transform,
- .ndone = 0) {
+ .transform) {
- mif2 <- as.logical(mif2)
- transform <- as.logical(transform)
+ mif2 <- as.logical(.mif2)
+ transform <- as.logical(.transform)
if (missing(seed)) seed <- NULL
if (!is.null(seed)) {
@@ -78,16 +73,16 @@
silent=FALSE
)
if (inherits(Np,"try-error"))
- stop("if ",sQuote("Np")," is a function, it must return a single positive integer")
+ stop("if ",sQuote("Np")," is a function, it must return a single positive integer",call.=FALSE)
}
if (length(Np)==1)
Np <- rep(Np,times=ntimes+1)
else if (length(Np)!=(ntimes+1))
- stop(sQuote("Np")," must have length 1 or length ",ntimes+1)
+ stop(sQuote("Np")," must have length 1 or length ",ntimes+1,call.=FALSE)
if (any(Np<=0))
- stop("number of particles, ",sQuote("Np"),", must always be positive")
+ stop("number of particles, ",sQuote("Np"),", must always be positive",call.=FALSE)
if (!is.numeric(Np))
- stop(sQuote("Np")," must be a number, a vector of numbers, or a function")
+ stop(sQuote("Np")," must be a number, a vector of numbers, or a function",call.=FALSE)
Np <- as.integer(Np)
if (is.null(dim(params))) {
@@ -199,7 +194,7 @@
if (mif2) {
cool.sched <- try(
- mif2.cooling(cooling.fraction,nt,cooling.m+.ndone,ntimes),
+ mif2.cooling(frac=cooling.fraction,nt=nt,m=cooling.m,n=ntimes),
silent=FALSE
)
if (inherits(cool.sched,"try-error"))
@@ -373,7 +368,7 @@
save.params=save.params,
seed=seed,
verbose=verbose,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 806
More information about the pomp-commits
mailing list