[Pomp-commits] r414 - in pkg: . R data inst/data-R man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Nov 8 21:57:15 CET 2010
Author: kingaa
Date: 2010-11-08 21:57:14 +0100 (Mon, 08 Nov 2010)
New Revision: 414
Modified:
pkg/DESCRIPTION
pkg/R/traj-match.R
pkg/data/euler.sir.rda
pkg/inst/data-R/euler.sir.R
pkg/man/traj-match.Rd
pkg/src/sir.c
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
- 'traj.match' is now implemented as a generic with methods for 'pomp' and 'traj.matched.pomp' objects
- the 'euler.sir' example has a different skeleton, one for which 'cases' gives the approximate weekly number of new cases
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-11-08 11:54:34 UTC (rev 413)
+++ pkg/DESCRIPTION 2010-11-08 20:57:14 UTC (rev 414)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.36-1
-Date: 2010-11-07
+Date: 2010-11-08
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing,
Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
Maintainer: Aaron A. King <kingaa at umich.edu>
Modified: pkg/R/traj-match.R
===================================================================
--- pkg/R/traj-match.R 2010-11-08 11:54:34 UTC (rev 413)
+++ pkg/R/traj-match.R 2010-11-08 20:57:14 UTC (rev 414)
@@ -2,6 +2,7 @@
"traj.matched.pomp",
contains="pomp",
representation=representation(
+ est="character",
evals="integer",
convergence="integer",
msg="character",
@@ -29,62 +30,15 @@
}
)
-trajectory.nll <- function (par, est, object, params, t0, fail.value = NA, ...) {
- if (missing(par)) par <- numeric(0)
- if (missing(est)) est <- integer(0)
- if (missing(params)) params <- coef(object)
- if (missing(t0)) t0 <- timezero(object)
- params[est] <- par
- x <- trajectory(object,params=params,t0=t0)
- d <- dmeasure(
- object,
- y=obs(object),
- x=x,
- times=time(object),
- params=as.matrix(params),
- log=TRUE
- )
- -sum(d)
-}
-
-trajectory.ls <- function (par, est, object, params, t0, fail.value = NA, transform, ...) {
- if (missing(transform)) transform <- identity
- if (missing(par)) par <- numeric(0)
- if (missing(est)) est <- integer(0)
- if (missing(params)) params <- coef(object)
- if (missing(t0)) t0 <- timezero(object)
- params[est] <- par
- x <- trajectory(object,params=params,t0=t0)
- d <- dmeasure(
- object,
- y=obs(object),
- x=x,
- times=time(object),
- params=as.matrix(params),
- log=TRUE
- )
- -sum(d)
-}
-
-traj.match <- function (object, start, est,
- method = c("Nelder-Mead","SANN","subplex"),
- gr = NULL, eval.only = FALSE, ...) {
+traj.match.internal <- function (object, start, est, method, gr, eval.only, ...) {
- if (!is(object,'pomp'))
- stop("traj.match error: ",sQuote("object")," must be a ",sQuote("pomp")," object",call.=FALSE)
-
- if (missing(start)) start <- coef(object)
-
- method <- match.arg(method)
-
if (eval.only) {
par.est <- integer(0)
} else {
- if (missing(est))
- stop("traj.match error: ",sQuote("est")," must be specified")
if (!is.character(est)) stop(sQuote("est")," must be a vector of parameter names")
if (!all(est%in%names(start)))
- stop("traj.match error: parameters named in ",sQuote("est")," must exist in ",sQuote("start"),call.=FALSE)
+ stop(sQuote("traj.match")," error: parameters named in ",sQuote("est"),
+ " must exist in ",sQuote("start"),call.=FALSE)
par.est <- which(names(start)%in%est)
guess <- start[par.est]
}
@@ -109,7 +63,6 @@
if (eval.only) {
coef(obj,names(start)) <- unname(start)
- obj at states[,] <- trajectory(obj,t0=t0)[,1,]
val <- obj.fn(start)
conv <- NA
evals <- c(1,0)
@@ -137,7 +90,6 @@
}
coef(obj,names(opt$par)) <- unname(opt$par)
- obj at states[,] <- trajectory(obj,t0=t0)[,1,]
msg <- if (is.null(opt$message)) character(0) else opt$message
conv <- opt$convergence
evals <- opt$counts
@@ -145,12 +97,62 @@
}
- ans <- new(
- "traj.matched.pomp",
- obj,
- evals=as.integer(evals),
- convergence=as.integer(conv),
- msg=msg,
- value=as.numeric(-val)
- )
+ obj at states <- trajectory(obj,t0=t0)[,1,]
+
+ new(
+ "traj.matched.pomp",
+ obj,
+ est=as.character(est),
+ evals=as.integer(evals),
+ convergence=as.integer(conv),
+ msg=msg,
+ value=as.numeric(-val)
+ )
}
+
+setGeneric("traj.match",function(object,...)standardGeneric("traj.match"))
+
+setMethod(
+ "traj.match",
+ signature=signature(object="pomp"),
+ function (object, start, est,
+ method = c("Nelder-Mead","SANN","subplex"),
+ gr = NULL, eval.only = FALSE, ...) {
+ if (missing(start)) start <- coef(object)
+ if (missing(est)) {
+ est <- character(0)
+ eval.only <- TRUE
+ }
+ method <- match.arg(method)
+ traj.match.internal(
+ object=object,
+ start=start,
+ est=est,
+ method=method,
+ gr=gr,
+ eval.only=eval.only,
+ ...
+ )
+ }
+ )
+
+setMethod(
+ "traj.match",
+ signature=signature(object="traj.matched.pomp"),
+ function (object, start, est,
+ method = c("Nelder-Mead","SANN","subplex"),
+ gr = NULL, eval.only = FALSE, ...) {
+ if (missing(start)) start <- coef(object)
+ if (missing(est)) est <- object at est
+ method <- match.arg(method)
+ traj.match.internal(
+ object=object,
+ start=start,
+ est=est,
+ method=method,
+ gr=gr,
+ eval.only=eval.only,
+ ...
+ )
+ }
+ )
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/data-R/euler.sir.R
===================================================================
--- pkg/inst/data-R/euler.sir.R 2010-11-08 11:54:34 UTC (rev 413)
+++ pkg/inst/data-R/euler.sir.R 2010-11-08 20:57:14 UTC (rev 414)
@@ -33,6 +33,10 @@
x0 <- numeric(length(snames))
names(x0) <- snames
x0[comp.names] <- round(p['pop']*fracs/sum(fracs))
+ ## since 'cases' is in 'zeronames' above, the IC for this variable
+ ## will only matter in trajectory computations
+ ## In trajectory computations, however, 'cases' will be roughly the weekly number of new cases
+ x0["cases"] <- x0["I"]*exp(params["gamma"])/52
x0
}
),
Modified: pkg/man/traj-match.Rd
===================================================================
--- pkg/man/traj-match.Rd 2010-11-08 11:54:34 UTC (rev 413)
+++ pkg/man/traj-match.Rd 2010-11-08 20:57:14 UTC (rev 414)
@@ -1,5 +1,9 @@
\name{traj.match}
\alias{traj.match}
+\alias{traj.match-pomp}
+\alias{traj.match,pomp-method}
+\alias{traj.match-traj.matched.pomp}
+\alias{traj.match,traj.matched.pomp-method}
\alias{logLik,traj.matched.pomp-method}
\alias{logLik-traj.matched.pomp}
\alias{$,traj.matched.pomp-method}
@@ -12,13 +16,17 @@
Match trajectories to data.
}
\usage{
-traj.match(object, start, est, method = c("Nelder-Mead", "SANN", "subplex"),
- gr = NULL, eval.only = FALSE, \dots)
+ \S4method{traj.match}{pomp}(object, start, est,
+ method = c("Nelder-Mead", "SANN", "subplex"),
+ gr = NULL, eval.only = FALSE, \dots)
+ \S4method{traj.match}{traj.matched.pomp}(object, start, est,
+ method = c("Nelder-Mead", "SANN", "subplex"),
+ gr = NULL, eval.only = FALSE, \dots)
}
\arguments{
\item{object}{A \code{pomp} object.}
- \item{start}{Initial guess for parameters.}
- \item{est}{Character vector containing the names of parameters to be estimated.}
+ \item{start}{initial guess for parameters.}
+ \item{est}{character vector containing the names of parameters to be estimated.}
\item{method}{
Optimization method.
Choices are \code{\link[subplex]{subplex}} and any of the methods used by \code{\link{optim}}.
@@ -27,7 +35,8 @@
Passed to \code{\link{optim}}.
}
\item{eval.only}{
- logical; if \code{TRUE}, no optimization is attempted and the quasi-loglikelihood value is evaluated at the \code{start} parameters.
+ logical;
+ if \code{TRUE}, no optimization is attempted and the log likelihood value is evaluated at the \code{start} parameters.
}
\item{\dots}{
Arguments that will be passed to \code{\link{optim}} in its \code{control} list.
@@ -46,9 +55,12 @@
number of function and gradient evaluations by the optimizer.
See \code{\link{optim}}.
}
- \item{value}{Value of the objective function.}
+ \item{value}{
+ value of the objective function.
+ Larger values indicate better fit (i.e., \code{traj.match} attempts to maximize this quantity.
+ }
\item{convergence, msg}{
- Convergence code and message from the optimizer.
+ convergence code and message from the optimizer.
See \code{\link{optim}}.
}
}
@@ -84,10 +96,10 @@
lines(x2~time,data=as(res,"data.frame"),col='red')
}
\seealso{
- \code{\link{optim}},
\code{\link{trajectory}},
\code{\link{pomp}},
- \code{\link{pomp-class}}
+ \code{\link{optim}},
+ \code{\link[subplex]{subplex}}
}
\keyword{models}
\keyword{ts}
Modified: pkg/src/sir.c
===================================================================
--- pkg/src/sir.c 2010-11-08 11:54:34 UTC (rev 413)
+++ pkg/src/sir.c 2010-11-08 20:57:14 UTC (rev 414)
@@ -182,7 +182,7 @@
DSDT = term[0]-term[1]-term[2];
DIDT = term[1]-term[3]-term[4];
DRDT = term[3]-term[5];
- DCDT = term[3]; // cases are cumulative recoveries
+ DCDT = DIDT*gamma/52; // roughly the weekly number of new cases
}
Modified: pkg/tests/sir.R
===================================================================
--- pkg/tests/sir.R 2010-11-08 11:54:34 UTC (rev 413)
+++ pkg/tests/sir.R 2010-11-08 20:57:14 UTC (rev 414)
@@ -101,7 +101,7 @@
terms[1]-terms[2]-terms[3],
terms[2]-terms[4]-terms[5],
terms[4]-terms[6],
- terms[4]
+ gamma/52*(terms[2]-terms[4]-terms[5])
)
xdot
}
@@ -119,6 +119,7 @@
rep(0,9) # zeros for 'cases', 'W', and the transition numbers
)
names(x0) <- c("S","I","R","cases","W","B","SI","SD","IR","ID","RD","dW")
+ x0["cases"] <- gamma/52*x0["I"]
x0
}
)
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2010-11-08 11:54:34 UTC (rev 413)
+++ pkg/tests/sir.Rout.save 2010-11-08 20:57:14 UTC (rev 414)
@@ -118,7 +118,7 @@
+ terms[1]-terms[2]-terms[3],
+ terms[2]-terms[4]-terms[5],
+ terms[4]-terms[6],
-+ terms[4]
++ gamma/52*(terms[2]-terms[4]-terms[5])
+ )
+ xdot
+ }
@@ -136,6 +136,7 @@
+ rep(0,9) # zeros for 'cases', 'W', and the transition numbers
+ )
+ names(x0) <- c("S","I","R","cases","W","B","SI","SD","IR","ID","RD","dW")
++ x0["cases"] <- gamma/52*x0["I"]
+ x0
+ }
+ )
@@ -406,7 +407,8 @@
terms <- c(mu * pop, foi * S, mu * S, gamma * I,
mu * I, mu * R)
xdot[1:4] <- c(terms[1] - terms[2] - terms[3], terms[2] -
- terms[4] - terms[5], terms[4] - terms[6], terms[4])
+ terms[4] - terms[5], terms[4] - terms[6], gamma/52 *
+ (terms[2] - terms[4] - terms[5]))
xdot
})
}, initializer = function(params, t0, ...) {
@@ -416,6 +418,7 @@
x0 <- c(round(pop * fracs/sum(fracs)), rep(0, 9))
names(x0) <- c("S", "I", "R", "cases", "W", "B",
"SI", "SD", "IR", "ID", "RD", "dW")
+ x0["cases"] <- gamma/52 * x0["I"]
x0
})
}, covar = basis, tcovar = tbasis)
@@ -432,7 +435,7 @@
zeronames = zeronames, tcovar = tcovar, covar = covar,
args = pairlist(...))
}
-<environment: 0x2633e80>
+<environment: 0x34f04e8>
process model density, dprocess =
function (x, times, params, ..., statenames = character(0), paramnames = character(0),
covarnames = character(0), tcovar, covar, log = FALSE)
@@ -442,7 +445,7 @@
covarnames = covarnames, tcovar = tcovar, covar = covar,
log = log, args = pairlist(...))
}
-<environment: 0x2ca6dd0>
+<environment: 0x318a708>
measurement model simulator, rmeasure =
function (x, t, params, covars, ...)
{
@@ -454,7 +457,7 @@
}
y
}
-<environment: 0x2c3f9d0>
+<environment: 0x359df90>
measurement model density, dmeasure =
function (y, x, t, params, log, covars, ...)
{
@@ -467,7 +470,7 @@
f
else exp(f)
}
-<environment: 0x2c3f9d0>
+<environment: 0x359df90>
skeleton ( vectorfield ) =
function (x, t, params, covars, ...)
{
@@ -479,7 +482,8 @@
terms <- c(mu * pop, foi * S, mu * S, gamma * I, mu *
I, mu * R)
xdot[1:4] <- c(terms[1] - terms[2] - terms[3], terms[2] -
- terms[4] - terms[5], terms[4] - terms[6], terms[4])
+ terms[4] - terms[5], terms[4] - terms[6], gamma/52 *
+ (terms[2] - terms[4] - terms[5]))
xdot
})
}
@@ -492,6 +496,7 @@
x0 <- c(round(pop * fracs/sum(fracs)), rep(0, 9))
names(x0) <- c("S", "I", "R", "cases", "W", "B", "SI",
"SD", "IR", "ID", "RD", "dW")
+ x0["cases"] <- gamma/52 * x0["I"]
x0
})
}
@@ -506,7 +511,7 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 2.066529 secs
+Time difference of 2.096126 secs
>
> pdf(file='sir.pdf')
>
@@ -538,7 +543,7 @@
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 4.560865 secs
+Time difference of 4.700443 secs
> plot(t3,X3['I',1,],type='l')
>
> f1 <- dprocess(
@@ -598,7 +603,7 @@
> x <- simulate(po,nsim=100)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 2.543359 secs
+Time difference of 2.518889 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t3 <- seq(0,20,by=1/52)
@@ -611,7 +616,7 @@
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.726177 secs
+Time difference of 0.71978 secs
> plot(t3,X4['I',1,],type='l')
>
> g2 <- dmeasure(
More information about the pomp-commits
mailing list