[Pomp-commits] r90 - in pkg: R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 13 14:20:33 CEST 2009
Author: kingaa
Date: 2009-04-13 14:20:33 +0200 (Mon, 13 Apr 2009)
New Revision: 90
Modified:
pkg/R/mif.R
pkg/man/mif.Rd
pkg/tests/ou2-mif.R
Log:
rewrite mif methods to use 'mif.internal'.
this obviates the need for a user-visible '.ndone' argument.
add some reporting when verbose=TRUE.
the 'particles' function can now be altered after the mif object is created.
some new tests have been added to tests/ou2-mif.R
Modified: pkg/R/mif.R
===================================================================
--- pkg/R/mif.R 2009-04-13 12:15:53 UTC (rev 89)
+++ pkg/R/mif.R 2009-04-13 12:20:33 UTC (rev 90)
@@ -15,6 +15,219 @@
list(alpha=alpha,gamma=gamma)
}
+mif.internal <- function (object, Nmif = 1,
+ start = NULL,
+ pars = NULL, ivps = NULL,
+ particles = NULL,
+ rw.sd = NULL, alg.pars = NULL,
+ weighted = TRUE, tol = 1e-17, warn = TRUE, max.fail = 0,
+ verbose = FALSE, .ndone = 0) {
+ 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)
+ start.names <- names(start)
+ if (is.null(start.names))
+ stop("mif error: ",sQuote("start")," must be a named vector",call.=FALSE)
+
+ if (is.null(rw.sd))
+ rw.sd <- object at random.walk.sd
+ rw.names <- names(rw.sd)
+ 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 (is.null(pars)) {
+ pars <- object at pars
+ }
+ if (length(pars)==0)
+ stop("mif error: at least one ordinary (non-IVP) parameter must be estimated",call.=FALSE)
+ if (is.null(ivps))
+ ivps <- object at ivps
+ if (
+ !is.character(pars) ||
+ !is.character(ivps) ||
+ !all(pars%in%start.names) ||
+ !all(ivps%in%start.names) ||
+ any(pars%in%ivps) ||
+ any(ivps%in%pars) ||
+ !all(pars%in%rw.names) ||
+ !all(ivps%in%rw.names)
+ )
+ stop(
+ "mif error: ",
+ sQuote("pars")," and ",sQuote("ivps"),
+ " must be mutually disjoint subsets of ",
+ sQuote("names(start)"),
+ " and must have a positive random-walk SDs specified in ",
+ sQuote("rw.sd"),
+ call.=FALSE
+ )
+
+ if (!all(rw.names%in%c(pars,ivps))) {
+ extra.rws <- rw.names[!(rw.names%in%c(pars,ivps))]
+ warning(
+ "mif warning: the variable(s) ",
+ paste(extra.rws,collapse=", "),
+ " have positive random-walk SDs specified, but are included in neither ",
+ sQuote("pars")," nor ",sQuote("ivps"),
+ ". These random walk SDs are ignored.",
+ call.=FALSE
+ )
+ }
+ rw.sd <- rw.sd[c(pars,ivps)]
+ rw.names <- names(rw.sd)
+
+ if (is.null(particles))
+ particles <- object at particles
+
+ if (is.null(alg.pars))
+ alg.pars <- object at alg.pars
+ if (!all(c('Np','cooling.factor','ic.lag','var.factor')%in%names(alg.pars)))
+ stop(
+ "mif error: ",sQuote("alg.pars"),
+ " must be a named list with elements ",
+ sQuote("Np"),",",sQuote("cooling.factor"),",",sQuote("ic.lag"),
+ ",and ",sQuote("var.factor"),
+ call.=FALSE
+ )
+
+ if (verbose) {
+ cat("performing",Nmif,"MIF iteration(s) to estimate parameter(s)",
+ paste(pars,collapse=", "))
+ if (length(ivps)>0)
+ cat(" and IVP(s)",paste(ivps,collapse=", "))
+ cat(" using random-walk with SD\n")
+ print(rw.sd)
+ cat(
+ "using",alg.pars$Np,"particles, variance factor",alg.pars$var.factor,
+ "\ninitial condition smoothing lag",alg.pars$ic.lag,
+ "and cooling factor",alg.pars$cooling.factor,"\n"
+ )
+ }
+
+ 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,
+ ncol=length(theta)+2,
+ dimnames=list(
+ seq(.ndone,.ndone+Nmif),
+ c('loglik','nfail',names(theta))
+ )
+ )
+ conv.rec[1,] <- c(NA,NA,theta)
+
+ for (n in seq(length=Nmif)) { # main loop
+
+ ## compute the cooled sigma
+ cool.sched <- try(
+ mif.cooling(alg.pars$cooling.factor,.ndone+n),
+ silent=FALSE
+ )
+ 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...
+ P <- try(
+ particles(object,Np=alg.pars$Np,center=theta,sd=sigma.n*alg.pars$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=object,
+ params=P,
+ tol=tol,
+ warn=warn,
+ max.fail=max.fail,
+ pred.mean=(n==Nmif),
+ pred.var=(weighted||(n==Nmif)),
+ filter.mean=TRUE,
+ .rw.sd=sigma.n[pars],
+ verbose=verbose
+ ),
+ silent=FALSE
+ )
+ if (inherits(x,'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
+ v1 <- cool.sched$gamma*(1+alg.pars$var.factor^2)*sigma[pars]^2
+ theta.hat <- cbind(theta[pars],x$filter.mean[pars,,drop=FALSE])
+ theta[pars] <- theta[pars]+apply(apply(theta.hat,1,diff)/t(v),2,sum)*v1
+ } else { # unweighted (flat) average
+ theta.hat <- x$filter.mean[pars,,drop=FALSE]
+ theta[pars] <- apply(theta.hat,1,mean)
+ }
+
+ ## update the IVPs using fixed-lag smoothing
+ theta[ivps] <- x$filter.mean[ivps,alg.pars$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)
+
+ if (verbose) cat("MIF iteration ",n," of ",Nmif," completed\n")
+
+ }
+
+ coef(object) <- theta
+
+ if (Nmif>0)
+ new(
+ "mif",
+ as(object,"pomp"),
+ ivps=ivps,
+ pars=pars,
+ Nmif=as.integer(Nmif),
+ particles=particles,
+ alg.pars=alg.pars,
+ random.walk.sd=sigma[rw.names],
+ pred.mean=x$pred.mean,
+ pred.var=x$pred.var,
+ filter.mean=x$filter.mean,
+ conv.rec=conv.rec,
+ eff.sample.size=x$eff.sample.size,
+ cond.loglik=x$cond.loglik,
+ loglik=x$loglik
+ )
+ else
+ new(
+ "mif",
+ as(object,"pomp"),
+ ivps=ivps,
+ pars=pars,
+ Nmif=0L,
+ particles=particles,
+ alg.pars=alg.pars,
+ random.walk.sd=sigma[rw.names],
+ conv.rec=conv.rec
+ )
+}
+
setMethod(
"mif",
"pomp",
@@ -26,10 +239,13 @@
weighted = TRUE, tol = 1e-17, warn = TRUE, max.fail = 0,
verbose = FALSE)
{
+ if (missing(start)) start <- NULL
if (missing(rw.sd))
stop("mif error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
- if (missing(alg.pars))
- stop("mif error: ",sQuote("alg.pars")," 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(particles)) { # use default: normal distribution
particles <- function (Np, center, sd, ...) {
matrix(
@@ -48,121 +264,22 @@
}
} 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
+ )
}
- 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
- )
- if (missing(start)) {
- start <- coef(object)
- if (length(start)==0)
- stop("mif error: ",sQuote("start")," must be specified",call.=FALSE)
- }
- start.names <- names(start)
- if (is.null(start.names))
- stop("mif error: ",sQuote("start")," must be a named vector",call.=FALSE)
+ if (missing(alg.pars))
+ stop("mif error: ",sQuote("alg.pars")," must be specified",call.=FALSE)
- rw.names <- names(rw.sd)
- 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)
+ mif.internal(object,Nmif=Nmif,start=start,pars=pars,ivps=ivps,particles=particles,
+ rw.sd=rw.sd,alg.pars=alg.pars,weighted=weighted,tol=tol,warn=warn,
+ max.fail=max.fail,verbose=verbose,.ndone=0)
- rw.names <- names(rw.sd[rw.sd>0])
- if (missing(pars)) {
- pars <- rw.names[!(rw.names%in%ivps)]
- }
- if (length(pars) == 0)
- stop("mif error: ",sQuote("pars")," must be a nonempty character vector",call.=FALSE)
- if (
- !is.character(pars) ||
- !is.character(ivps) ||
- !all(pars%in%start.names) ||
- !all(ivps%in%start.names) ||
- any(pars%in%ivps) ||
- any(ivps%in%pars) ||
- !all(pars%in%rw.names) ||
- !all(ivps%in%rw.names)
- )
- stop(
- "mif error: ",
- sQuote("pars")," and ",sQuote("ivps"),
- " must be mutually disjoint subsets of ",
- sQuote("names(start)"),
- " and must have a positive random-walk SDs specified in ",
- sQuote("rw.sd"),
- call.=FALSE
- )
-
- if (!all(rw.names%in%c(pars,ivps))) {
- extra.rws <- rw.names[!(rw.names%in%c(pars,ivps))]
- warning(
- "mif warning: the variable(s) ",
- paste(extra.rws,collapse=", "),
- " have positive random-walk SDs specified, but are included in neither ",
- sQuote("pars")," nor ",sQuote("ivps"),
- ". These random walk SDs are ignored.",
- call.=FALSE
- )
- }
- rw.sd <- rw.sd[c(pars,ivps)]
- rw.names <- names(rw.sd)
-
- if (!all(c('Np','cooling.factor','ic.lag','var.factor')%in%names(alg.pars)))
- stop(
- "mif error: ",sQuote("alg.pars"),
- " must be a named list with elements ",
- sQuote("Np"),",",sQuote("cooling.factor"),",",sQuote("ic.lag"),
- ",and ",sQuote("var.factor"),
- call.=FALSE
- )
-
- coef(object) <- start
-
- newmif <- new(
- "mif",
- object,
- ivps=ivps,
- pars=pars,
- Nmif=as.integer(0),
- particles=particles,
- alg.pars=alg.pars,
- random.walk.sd=rw.sd,
- pred.mean=matrix(NA,0,0),
- pred.var=matrix(NA,0,0),
- filter.mean=matrix(NA,0,0),
- conv.rec=matrix(
- c(NA,NA,start),
- nrow=1,
- ncol=2+length(start),
- dimnames=list(
- 0,
- c('loglik','nfail',start.names)
- )
- ),
- cond.loglik=numeric(0),
- eff.sample.size=numeric(0),
- loglik=as.numeric(NA)
- )
-
- if (Nmif > 0) {
- mif(
- newmif,
- Nmif=Nmif,
- weighted=weighted,
- tol=tol,
- warn=warn,
- max.fail=max.fail,
- verbose=verbose,
- .ndone=0
- )
- } else {
- newmif
- }
}
)
@@ -170,175 +287,19 @@
setMethod(
"mif",
"mif",
- function (object, Nmif, start, pars, ivps, rw.sd, alg.pars,
- weighted = TRUE, tol = 1e-17, warn = TRUE, max.fail = 0,
- verbose = FALSE, .ndone = 0)
+ function (object, Nmif, ...)
{
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(rw.sd)) rw.sd <- object at random.walk.sd
- if (missing(alg.pars)) alg.pars <- object at alg.pars
- theta <- start
- start.names <- names(start)
- if (is.null(start.names))
- stop("mif error: ",sQuote("start")," must be a named vector",call.=FALSE)
-
- sigma <- rep(0,length(start))
- names(sigma) <- start.names
-
- rw.names <- names(rw.sd)
- 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(pars) == 0)
- stop("mif error: ",sQuote("pars")," must be a nonempty character vector",call.=FALSE)
- if (
- !is.character(pars) ||
- !is.character(ivps) ||
- !all(pars%in%start.names) ||
- !all(ivps%in%start.names) ||
- any(pars%in%ivps) ||
- any(ivps%in%pars) ||
- !all(pars%in%rw.names) ||
- !all(ivps%in%rw.names)
- )
- stop(
- "mif error: ",
- sQuote("pars")," and ",sQuote("ivps"),
- " must be mutually disjoint subsets of ",
- sQuote("names(start)"),
- " and must have a positive random-walk SDs specified in ",
- sQuote("rw.sd"),
- call.=FALSE
- )
-
- if (!all(rw.names%in%c(pars,ivps))) {
- extra.rws <- rw.names[!(rw.names%in%c(pars,ivps))]
- warning(
- "mif warning: the variable(s) ",
- paste(extra.rws,collapse=", "),
- " have positive random-walk SDs specified, but are included in neither ",
- sQuote("pars")," nor ",sQuote("ivps"),
- ". These random walk SDs are ignored.",
- call.=FALSE
- )
- }
- rw.sd <- rw.sd[c(pars,ivps)]
- rw.names <- names(rw.sd)
-
- sigma[rw.names] <- rw.sd
-
- if (!all(c('Np','cooling.factor','ic.lag','var.factor')%in%names(alg.pars)))
- stop(
- "mif error: ",sQuote("alg.pars")," must be a named list with elements ",
- sQuote("Np"),",",sQuote("cooling.factor"),",",sQuote("ic.lag"),
- ",and ",sQuote("var.factor"),
- call.=FALSE
- )
-
- conv.rec <- matrix(
- data=NA,
- nrow=Nmif+1,
- ncol=length(theta)+2,
- dimnames=list(
- seq(.ndone,.ndone+Nmif),
- c('loglik','nfail',names(theta))
- )
- )
- conv.rec[1,] <- c(NA,NA,theta)
-
- for (n in seq(length=Nmif)) { # main loop
-
- ## compute the cooled sigma
- cool.sched <- try(
- mif.cooling(alg.pars$cooling.factor,.ndone+n),
- silent=FALSE
- )
- 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...
- P <- try(
- particles(object,Np=alg.pars$Np,center=theta,sd=sigma.n*alg.pars$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=object,
- params=P,
- tol=tol,
- warn=warn,
- max.fail=max.fail,
- pred.mean=(n==Nmif),
- pred.var=(weighted||(n==Nmif)),
- filter.mean=TRUE,
- .rw.sd=sigma.n[pars],
- verbose=verbose
- ),
- silent=FALSE
- )
- if (inherits(x,'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
- v1 <- cool.sched$gamma*(1+alg.pars$var.factor^2)*sigma[pars]^2
- theta.hat <- cbind(theta[pars],x$filter.mean[pars,,drop=FALSE])
- theta[pars] <- theta[pars]+apply(apply(theta.hat,1,diff)/t(v),2,sum)*v1
- } else { # unweighted (flat) average
- theta.hat <- x$filter.mean[pars,,drop=FALSE]
- theta[pars] <- apply(theta.hat,1,mean)
- }
-
- ## update the IVPs using fixed-lag smoothing
- theta[ivps] <- x$filter.mean[ivps,alg.pars$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)
-
- if (verbose) cat("MIF iteration ",n," of ",Nmif," completed\n")
-
- }
-
- coef(object) <- theta
-
- new(
- "mif",
- object,
- ivps=ivps,
- pars=pars,
- Nmif=as.integer(Nmif),
- particles=object at particles,
- alg.pars=alg.pars,
- random.walk.sd=sigma[rw.names],
- pred.mean=x$pred.mean,
- pred.var=x$pred.var,
- filter.mean=x$filter.mean,
- conv.rec=conv.rec,
- cond.loglik=x$cond.loglik,
- eff.sample.size=x$eff.sample.size,
- loglik=x$loglik
- )
+ mif.internal(object,Nmif=Nmif,...)
}
)
setMethod(
'continue',
'mif',
- function (object, Nmif, ...) {
+ function (object, Nmif = 1, ...) {
ndone <- object at Nmif
- obj <- mif(object,Nmif=Nmif,.ndone=ndone,...)
+ obj <- mif.internal(object,Nmif=Nmif,.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,
@@ -348,4 +309,3 @@
obj
}
)
-
Modified: pkg/man/mif.Rd
===================================================================
--- pkg/man/mif.Rd 2009-04-13 12:15:53 UTC (rev 89)
+++ pkg/man/mif.Rd 2009-04-13 12:20:33 UTC (rev 90)
@@ -17,10 +17,8 @@
\S4method{mif}{pomp}(object, Nmif = 1, start, pars, ivps = character(0),
particles, rw.sd, alg.pars, weighted = TRUE,
tol = 1e-17, warn = TRUE, max.fail = 0, verbose = FALSE)
-\S4method{mif}{mif}(object, Nmif, start, pars, ivps, rw.sd, alg.pars,
- weighted = TRUE, tol = 1e-17, warn = TRUE, max.fail = 0,
- verbose = FALSE, .ndone = 0)
-\S4method{continue}{mif}(object, Nmif, \dots)
+\S4method{mif}{mif}(object, Nmif, \dots)
+\S4method{continue}{mif}(object, Nmif = 1, \dots)
}
\arguments{
\item{object}{
@@ -98,10 +96,6 @@
\item{verbose}{
logical; if TRUE, print progress reports.
}
- \item{.ndone}{
- for internal use by \code{continue}.
- Do not meddle with this!
- }
\item{\dots}{
Additional arguments that can be used to override the defaults.
}
@@ -111,15 +105,15 @@
The call sequence is \code{mif(object)}.
By default, the same parameters used for the original MIF run are re-used.
If one does specify additional arguments, these will override the defaults.
- An exception is that one cannot override the \code{particles} function.
}
\section{Continuing MIF Iterations}{
One can continue a series of MIF iterations from where one left off.
The call sequence is \code{continue(object, Nmif)}.
This will perform \code{Nmif} additional MIF iterations on the \code{mif} object \code{object}.
A call to \code{mif} to perform \code{Nmif=m} iterations followed by a call to \code{continue} to perform \code{Nmif=n} iterations will produce precisely the same effect as a single call to \code{mif} to perform \code{Nmif=m+n} iterations.
- Additional arguments are passed to \code{mif}.
- This feature can be used to change any of the parameters (except the \code{particles} function).
+ By default, all the algorithmic parameters are the same as used in the original call to \code{mif}.
+ Additional arguments will override the defaults.
+ This feature can be used to change any of the parameters.
}
\section{Details}{
\strong{It is the user's responsibility to ensure that, if the optional \code{particles} argument is given, that the \code{particles} function satisfies the following conditions:}
Modified: pkg/tests/ou2-mif.R
===================================================================
--- pkg/tests/ou2-mif.R 2009-04-13 12:15:53 UTC (rev 89)
+++ pkg/tests/ou2-mif.R 2009-04-13 12:20:33 UTC (rev 90)
@@ -81,3 +81,9 @@
print(ff$loglik)
fit <- mif(fit,rw.sd=c(x1.0=5,x2.0=5,alpha.1=0.1,alpha.4=0.1))
fit <- continue(fit,Nmif=2,ivps=c("x1.0"),pars=c("alpha.1"))
+s <- coef(fit)
+s[2] <- 0.01
+fit <- mif(fit,Nmif=10,start=s)
+fit <- mif(ou2,Nmif=3,rw.sd=c(alpha.1=0.1,alpha.4=0.1),alg.pars=list(Np=1000,cooling.factor=0.98,var.factor=1,ic.lag=2))
+fit <- continue(fit,Nmif=5,alg.pars=list(Np=2000,cooling.factor=0.98,var.factor=1,ic.lag=2))
+fit <- continue(fit,ivps=c("x1.0"),rw.sd=c(alpha.1=0.1,alpha.4=0.1,x1.0=5,x2.0=5),Nmif=3)
More information about the pomp-commits
mailing list