[Pomp-commits] r83 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Apr 10 18:51:30 CEST 2009
Author: kingaa
Date: 2009-04-10 18:51:30 +0200 (Fri, 10 Apr 2009)
New Revision: 83
Modified:
pkg/DESCRIPTION
pkg/R/mif.R
pkg/R/pfilter.R
pkg/man/mif.Rd
pkg/man/pfilter.Rd
Log:
The handling of the argument 'pars' and 'ivps' in mif has been changed to eliminate redundancy and a possible source of errors. Now 'pars' and 'ivps' are checked to make sure that all names are in 'names(start)' and 'names(rw.sd)' and that they are mutually disjoint. Further, all the entries in 'rw.sd[pars]' and 'rw.sd[ivps]' must be positive or an error is generated. If 'rw.sd' has positive values in slots not in the union of 'ivps' and 'pars', these are ignored and a warning is generated. The 'random.walk.sd' slot of the resulting 'mif' object contains only nonzero entries corresponding to the union of 'pars' and 'ivps'.
If 'pars' is not specified in the arguments, then it is taken to be
'names(rw.sd)[!(names(rw.sd)%in%ivps)]'
Thus it is no longer necessary to specify 'pars' explicitly, provided one takes care that all and only nonzero entries in 'rw.sd' correspond to parameters that one wishes to estimate. One can be a bit more careless with 'rw.sd' if one uses 'pars'. In any case, entries in 'rw.sd' that are zero are ignored.
The computation of prediction means and variances and filtering means in 'pfilter' is handled differently. When used as a toplevel function (i.e., when the .rw.sd parameter used internally by 'mif' is unset), then the prediction means and variances and the filtering means are only computed for the state variables. This eliminates an annoying and unnecessary behavior that arose when parameters were set to nonfinite values.
When .rw.sd is set (as when 'mif' calls 'pfilter'), the prediction means and variances are computed only for all state variables, but for only those parameters which are undergoing the MIF random walk. Filtering means are computed for all state variables and parameters, since they are needed for the fixed-lag smoothing on the initial-value parameters.
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2009-04-10 16:37:44 UTC (rev 82)
+++ pkg/DESCRIPTION 2009-04-10 16:51:30 UTC (rev 83)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.22-6
-Date: 2009-03-11
+Version: 0.22-7
+Date: 2009-04-09
Author: Aaron A. King, Edward L. Ionides, Carles Martinez Breto, Steve Ellner, Bruce Kendall
Maintainer: Aaron A. King <kingaa at umich.edu>
Description: Inference methods for partially-observed Markov processes
Modified: pkg/R/mif.R
===================================================================
--- pkg/R/mif.R 2009-04-10 16:37:44 UTC (rev 82)
+++ pkg/R/mif.R 2009-04-10 16:51:30 UTC (rev 83)
@@ -24,10 +24,8 @@
particles,
rw.sd, alg.pars,
weighted = TRUE, tol = 1e-17, warn = TRUE, max.fail = 0,
- verbose = FALSE, .ndone = 0)
+ verbose = FALSE)
{
- if (missing(pars))
- stop("mif error: ",sQuote("pars")," must be specified",call.=FALSE)
if (missing(rw.sd))
stop("mif error: ",sQuote("rw.sd")," must be specified",call.=FALSE)
if (missing(alg.pars))
@@ -62,6 +60,16 @@
if (is.null(start.names))
stop("mif error: ",sQuote("start")," must be a named vector",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)
+
+ 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 (
@@ -70,33 +78,33 @@
!all(pars%in%start.names) ||
!all(ivps%in%start.names) ||
any(pars%in%ivps) ||
- any(ivps%in%pars)
+ 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 elements of ",
+ " 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(rw.names))
- stop("mif error: ",sQuote("rw.sd")," must be a named vector",call.=FALSE)
- if (any(!(rw.names%in%start.names)))
- stop("mif error: all the names of ",sQuote("rw.sd")," must be names of ",sQuote("start"),call.=FALSE)
- if (any(rw.sd[c(pars,ivps)]<=0)) {
- zero.pars <- names(which(rw.sd[c(pars,ivps)]<=0))
- stop(
- "mif error: for every parameter you wish to estimate, you must specify a positive ",
- sQuote("rw.sd"),
- ". ",
- sQuote("rw.sd"),
- " is non-positive for ",
- paste(zero.pars,collapse=", "),
- call.=FALSE
- )
- }
if (!all(c('Np','cooling.factor','ic.lag','var.factor')%in%names(alg.pars)))
stop(
@@ -171,6 +179,16 @@
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 (
@@ -179,33 +197,34 @@
!all(pars%in%start.names) ||
!all(ivps%in%start.names) ||
any(pars%in%ivps) ||
- any(ivps%in%pars)
+ 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 elements of ",
+ " must be mutually disjoint subsets of ",
sQuote("names(start)"),
+ " and must have a positive random-walk SDs specified in ",
+ sQuote("rw.sd"),
call.=FALSE
)
- sigma <- rep(0,length(start))
- names(sigma) <- start.names
+ 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 (any(!(rw.names%in%start.names)))
- stop("mif error: all the names of ",sQuote("rw.sd")," must be names of ",sQuote("start"),call.=FALSE)
- if (any(rw.sd[c(pars,ivps)]<=0)) {
- zero.pars <- names(which(rw.sd[c(pars,ivps)]<=0))
- stop(
- "mif error: for every parameter you wish to estimate, you must specify a positive ",
- sQuote("rw.sd"),
- ". ",
- sQuote("rw.sd"),
- " is non-positive for ",
- paste(zero.pars,collapse=", "),
- call.=FALSE
- )
- }
+
sigma[rw.names] <- rw.sd
if (!all(c('Np','cooling.factor','ic.lag','var.factor')%in%names(alg.pars)))
@@ -264,7 +283,7 @@
warn=warn,
max.fail=max.fail,
pred.mean=(n==Nmif),
- pred.var=TRUE,
+ pred.var=(weighted||(n==Nmif)),
filter.mean=TRUE,
.rw.sd=sigma.n[pars]
),
@@ -305,7 +324,7 @@
Nmif=as.integer(Nmif),
particles=object at particles,
alg.pars=alg.pars,
- random.walk.sd=sigma,
+ random.walk.sd=sigma[rw.names],
pred.mean=x$pred.mean,
pred.var=x$pred.var,
filter.mean=x$filter.mean,
Modified: pkg/R/pfilter.R
===================================================================
--- pkg/R/pfilter.R 2009-04-10 16:37:44 UTC (rev 82)
+++ pkg/R/pfilter.R 2009-04-10 16:51:30 UTC (rev 83)
@@ -31,7 +31,6 @@
)
)
}
- npars <- nrow(params)
paramnames <- rownames(params)
if (is.null(paramnames))
stop("pfilter error: ",sQuote("params")," must have rownames",call.=FALSE)
@@ -48,18 +47,21 @@
if (any(!(rw.names%in%paramnames)))
stop("pfilter error: the rownames of ",sQuote("params")," must include all of the names of ",sQuote(".rw.sd"),"",call.=FALSE)
sigma <- .rw.sd
+ } else {
+ rw.names <- character(0)
}
loglik <- rep(NA,ntimes)
eff.sample.size <- rep(NA,ntimes)
nfail <- 0
+ npars <- length(rw.names)
pred.m <- if (pred.mean)
matrix(
data=0,
nrow=nvars+npars,
ncol=ntimes,
- dimnames=list(c(statenames,paramnames),NULL)
+ dimnames=list(c(statenames,rw.names),NULL)
)
else NULL
@@ -68,17 +70,25 @@
data=0,
nrow=nvars+npars,
ncol=ntimes,
- dimnames=list(c(statenames,paramnames),NULL)
+ dimnames=list(c(statenames,rw.names),NULL)
)
else NULL
filt.m <- if (filter.mean)
- matrix(
- data=0,
- nrow=nvars+npars,
- ncol=ntimes,
- dimnames=list(c(statenames,paramnames),NULL)
- )
+ if (random.walk)
+ matrix(
+ data=0,
+ nrow=nvars+length(paramnames),
+ ncol=ntimes,
+ dimnames=list(c(statenames,paramnames),NULL)
+ )
+ else
+ matrix(
+ data=0,
+ nrow=nvars,
+ ncol=ntimes,
+ dimnames=list(statenames,NULL)
+ )
else NULL
times <- time(object,t0=TRUE)
@@ -106,7 +116,7 @@
xx <- try(
c(
apply(x,1,mean),
- apply(params,1,mean)
+ apply(params[rw.names,,drop=FALSE],1,mean)
),
silent=FALSE
)
@@ -127,18 +137,20 @@
call.=FALSE
)
}
- problem.indices <- unique(which(!is.finite(params),arr.ind=TRUE)[,1])
- if (length(problem.indices)>0) {
- stop(
- "pfilter error: non-finite parameter(s): ",
- paste(rownames(params)[problem.indices],collapse=', '),
- call.=FALSE
- )
+ if (random.walk) {
+ problem.indices <- unique(which(!is.finite(params[rw.names,,drop=FALSE]),arr.ind=TRUE)[,1])
+ if (length(problem.indices)>0) {
+ stop(
+ "pfilter error: non-finite parameter(s): ",
+ paste(rw.names[problem.indices],collapse=', '),
+ call.=FALSE
+ )
+ }
}
xx <- try(
c(
apply(x,1,var),
- apply(params,1,var)
+ apply(params[rw.names,,drop=FALSE],1,var)
),
silent=FALSE
)
@@ -192,7 +204,8 @@
## compute filtering means
if (filter.mean) {
filt.m[statenames,nt] <- x %*% weights
- filt.m[paramnames,nt] <- params %*% weights
+ if (random.walk)
+ filt.m[paramnames,nt] <- params %*% weights
}
## Matrix with samples (columns) from filtering distribution theta.t | Y.t
Modified: pkg/man/mif.Rd
===================================================================
--- pkg/man/mif.Rd 2009-04-10 16:37:44 UTC (rev 82)
+++ pkg/man/mif.Rd 2009-04-10 16:51:30 UTC (rev 83)
@@ -16,7 +16,7 @@
mif(object, \dots)
\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, .ndone = 0)
+ 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)
@@ -34,12 +34,13 @@
This must be a named vector.
}
\item{pars}{
- character vector;
- names of ordinary parameters to be estimated.
+ optional character vector naming the ordinary parameters to be estimated.
+ Every parameter named in \code{pars} must have a positive random-walk standard deviation specified in \code{rw.sd}.
+ Leaving \code{pars} unspecified is equivalent to setting it equal to the names of all parameters with a positive value of \code{rw.sd} that are not \code{ivps}.
}
\item{ivps}{
- character vector;
- names of initial-value parameters to be estimated.
+ optional character vector naming the initial-value parameters to be estimated.
+ Every parameter named in \code{ivps} must have a positive random-walk standard deviation specified in \code{rw.sd}.
}
\item{particles}{
Function of prototype \code{particles(Np,center,sd,...)} which sets up the initial particle matrix by drawing a sample of size \code{Np} from the initial particle distribution centered at \code{center} and of width \code{sd}.
@@ -47,11 +48,14 @@
}
\item{rw.sd}{
numeric vector with names; the intensity of the random walk to be applied to parameters.
- The random walk is only applied to parameters named in \code{pars}.
- The algorithm requires that the random walk be nontrivial.
- Thus, each element in \code{rw.sd[pars]} must be positive.
+ The random walk is only applied to parameters named in \code{pars} (i.e., not to those named in \code{ivps}).
+ The algorithm requires that the random walk be nontrivial, so each element in \code{rw.sd[pars]} must be positive.
\code{rw.sd} is also used to scale the initial-value parameters (via the \code{particles} function).
Therefore, each element of \code{rw.sd[ivps]} must be positive.
+ The following must be satisfied:
+ \code{names(rw.sd)} must be a subset of \code{names(start)},
+ \code{rw.sd} must be non-negative (zeros are simply ignored),
+ the name of every positive element of \code{rw.sd} must be in either \code{pars} or \code{ivps}.
}
\item{alg.pars}{
A named list of algorithm parameters.
@@ -96,7 +100,7 @@
}
\item{.ndone}{
for internal use by \code{continue}.
- do not meddle with this!
+ Do not meddle with this!
}
\item{\dots}{
Additional arguments that can be used to override the defaults.
Modified: pkg/man/pfilter.Rd
===================================================================
--- pkg/man/pfilter.Rd 2009-04-10 16:37:44 UTC (rev 82)
+++ pkg/man/pfilter.Rd 2009-04-10 16:51:30 UTC (rev 83)
@@ -54,7 +54,6 @@
}
\item{.rw.sd}{
For internal use with the MIF algorithm.
- If \code{TRUE}, the specified random walk SD is used.
}
\item{verbose}{
If \code{TRUE}, progress information is reported as \code{pfilter} works.
@@ -66,8 +65,8 @@
\value{
A list with the following elements:
\item{pred.mean}{
- The \code{nvars+npars} x \code{ntimes} matrix of prediction means, where \code{ntimes} is the length of the time series contained in \code{object}.
- The rows correspond to states and parameters, in that order.
+ The matrix of prediction means.
+ The rows correspond to states and parameters (if appropriate), in that order, the columns to successive observations in the time series contained in \code{object}.
}
\item{pred.variance}{
The matrix of prediction variances, in the same format as \code{pred.mean}.
More information about the pomp-commits
mailing list