[Pomp-commits] r286 - in pkg: . R inst man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jul 19 18:12:36 CEST 2010
Author: kingaa
Date: 2010-07-19 18:12:36 +0200 (Mon, 19 Jul 2010)
New Revision: 286
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/pomp-methods.R
pkg/inst/ChangeLog
pkg/inst/NEWS
pkg/man/pomp-methods.Rd
pkg/tests/rw2.R
pkg/tests/rw2.Rout.save
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
- added new 'window', 'timezero', and 'timezero<-' methods for 'pomp' objects
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-07-06 17:47:12 UTC (rev 285)
+++ pkg/DESCRIPTION 2010-07-19 16:12:36 UTC (rev 286)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.31-1
-Date: 2010-06-30
+Version: 0.31-2
+Date: 2010-07-19
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing,
Matthew J. Ferrari, Michael Lavine
Maintainer: Aaron A. King <kingaa at umich.edu>
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-07-06 17:47:12 UTC (rev 285)
+++ pkg/NAMESPACE 2010-07-19 16:12:36 UTC (rev 286)
@@ -18,7 +18,7 @@
)
importFrom(graphics,plot)
-importFrom(stats,simulate,time,coef,logLik)
+importFrom(stats,simulate,time,coef,logLik,window)
importFrom(mvtnorm,dmvnorm,rmvnorm)
importFrom(deSolve,lsoda)
importFrom(subplex,subplex)
@@ -28,7 +28,7 @@
exportMethods(
'plot','show','print','coerce',
'dprocess','rprocess','rmeasure','dmeasure','init.state','skeleton',
- 'data.array','coef','logLik','time','time<-',
+ 'data.array','coef','logLik','time','time<-','timezero','timezero<-','window',
'simulate','pfilter',
'particles','mif','continue','coef<-','states','trajectory',
'pred.mean','pred.var','filter.mean','conv.rec',
Modified: pkg/R/pomp-methods.R
===================================================================
--- pkg/R/pomp-methods.R 2010-07-06 17:47:12 UTC (rev 285)
+++ pkg/R/pomp-methods.R 2010-07-19 16:12:36 UTC (rev 286)
@@ -9,6 +9,10 @@
setGeneric("states",function(object,...)standardGeneric("states"))
+setGeneric("timezero",function(object,...)standardGeneric("timezero"))
+
+setGeneric("timezero<-",function(object,...,value)standardGeneric("timezero<-"))
+
## 'coerce' method: allows for coercion of a "pomp" object to a data-frame
setAs(
from="pomp",
@@ -105,6 +109,38 @@
}
)
+setMethod(
+ "window",
+ "pomp",
+ function (x, start, end, ...) {
+ tm <- time(x,t0=FALSE)
+ if (missing(start))
+ start <- tm[1]
+ if (missing(end))
+ end <- tm[length(tm)]
+ tm <- tm[(tm>=start)&(tm<=end)]
+ time(x,t0=FALSE) <- tm
+ x
+ }
+ )
+
+setMethod("timezero","pomp",function(object,...)object at t0)
+
+setMethod(
+ "timezero<-",
+ "pomp",
+ function(object,...,value) {
+ if (value>object at times[1])
+ if (!is.numeric(value) || length(value) > 1)
+ stop("the zero-time ",sQuote("t0")," must be a single number",call.=TRUE)
+ if (value > object at times[1])
+ stop("the zero-time ",sQuote("t0")," must occur no later than the first observation",call.=TRUE)
+ storage.mode(value) <- "double"
+ object at t0 <- value
+ object
+ }
+ )
+
## extract the coefficients
setMethod(
"coef",
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2010-07-06 17:47:12 UTC (rev 285)
+++ pkg/inst/ChangeLog 2010-07-19 16:12:36 UTC (rev 286)
@@ -1,3 +1,23 @@
+2010-07-06 kingaa
+
+ * [r285] R/pomp-methods.R: - error checking in 'data.array'
+ * [r284] R/pfilter.R: - minor
+ * [r283] R/pmcmc.R: - unname before 'coef<-' to avoid warning
+
+2010-07-02 kingaa
+
+ * [r282] R/mif.R, R/pomp-methods.R, inst/ChangeLog,
+ inst/doc/advanced_topics_in_pomp.pdf, inst/doc/intro_to_pomp.Rnw,
+ inst/doc/intro_to_pomp.pdf, inst/doc/pomp.bib,
+ man/pomp-methods.Rd, man/simulate-pomp.Rd, tests/rw2.R,
+ tests/rw2.Rout.save: - add warning in 'coef<-' where parameter
+ values might unintentionally be mismatched
+ - redefine some of the generics using standardGeneric()
+ - change 'include.t0' argument in 'time<-()' to 't0' to match
+ 'time()'
+ - improve documentation of 'simulate' method slightly
+ - add Bhadra article to pomp.bib
+
2010-06-30 kingaa
* [r281] DESCRIPTION, R/traj-match.R, R/trajmatch.R,
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2010-07-06 17:47:12 UTC (rev 285)
+++ pkg/inst/NEWS 2010-07-19 16:12:36 UTC (rev 286)
@@ -1,3 +1,6 @@
+NEW IN VERSION 0.31-2:
+ o new "window", "timezero", and "timezero<-" methods are available for manipulating 'pomp' objects.
+
NEW IN VERSION 0.31-1:
o pomp now includes an implementation of the particle Markov chain Monte Carlo algorithm of Andrieu et al. (Andrieu, C., Doucet, A., & Holenstein, R. (2010) Particle Markov chain Monte Carlo methods. J. R. Stat. Soc. B 72:1--33).
This algorithm is implemented in the method 'pmcmc': see 'help(pmcmc)' for details.
Modified: pkg/man/pomp-methods.Rd
===================================================================
--- pkg/man/pomp-methods.Rd 2010-07-06 17:47:12 UTC (rev 285)
+++ pkg/man/pomp-methods.Rd 2010-07-19 16:12:36 UTC (rev 286)
@@ -12,6 +12,14 @@
\alias{time<-}
\alias{time<-,pomp-method}
\alias{time<--pomp}
+\alias{timezero}
+\alias{timezero,pomp-method}
+\alias{timezero-pomp}
+\alias{timezero<-}
+\alias{timezero<-,pomp-method}
+\alias{timezero<--pomp}
+\alias{window,pomp-method}
+\alias{window--pomp}
\alias{data.array,pomp-method}
\alias{data.array-pomp}
\alias{data.array}
@@ -34,6 +42,9 @@
\S4method{states}{pomp}(object, vars, \dots)
\S4method{time}{pomp}(x, t0 = FALSE, \dots)
\S4method{time}{pomp}(object, t0 = FALSE, \dots) <- value
+\S4method{timezero}{pomp}(object, \dots)
+\S4method{timezero}{pomp}(object, \dots) <- value
+\S4method{window}{pomp}(x, start, end, \dots)
\S4method{show}{pomp}(object)
\S4method{as}{pomp}(object, class)
\S4method{coerce}{pomp,data.frame}(from, to = "data.frame", strict = TRUE)
@@ -55,13 +66,16 @@
}
\item{value}{
numeric;
- values to be assigned to the parameters.
+ values to be assigned.
}
\item{t0}{
logical;
if TRUE on a call to \code{time}, the zero time is prepended to the time vector;
if TRUE on a call to \code{time<-}, the first element in \code{value} is taken to be the initial time.
}
+ \item{start, end}{
+ start and end times of the window.
+ }
\item{class}{
character;
name of the class to which \code{object} should be coerced.
@@ -137,6 +151,14 @@
The second and subsequent elements of \code{value} are taken to be the observation times.
Those data and states (if they exist) corresponding to the new times are retained.
}
+ \item{timezero, timezero<-}{
+ \code{timezero(object)} returns the zero-time \code{t0}.
+ \code{timezero(object) <- value} sets the zero-time to \code{value}.
+ }
+ \item{window}{
+ \code{window(x,start=t1,end=t2} returns a new \code{pomp} object, identical to \code{x} but with only the data in the window between times \code{t1} and \code{t2} (inclusive).
+ By default, \code{start} is the time of the first observation and \code{end} is the time of the last.
+ }
\item{show}{Displays the \code{pomp} object.}
\item{plot}{
Plots the data and state trajectories (if the latter exist).
Modified: pkg/tests/rw2.R
===================================================================
--- pkg/tests/rw2.R 2010-07-06 17:47:12 UTC (rev 285)
+++ pkg/tests/rw2.R 2010-07-19 16:12:36 UTC (rev 286)
@@ -178,6 +178,13 @@
max(abs(d2-d3),na.rm=T)
max(abs(e2-e3),na.rm=T)
+new <- window(rw2,start=20,end=30)
+new <- simulate(new)
+
+timezero(new)
+timezero(new) <- 19
+print(simulate(new))
+
time(rw2) <- seq(1,1000,by=20)
x <- simulate(rw2)
states(x)[,1:5]
@@ -194,3 +201,4 @@
time(rw2,t0=TRUE) <- c(0,0,20,25.8,50,60)
time(rw2) <- c(0,20,25.8,50,60)
+
Modified: pkg/tests/rw2.Rout.save
===================================================================
--- pkg/tests/rw2.Rout.save 2010-07-06 17:47:12 UTC (rev 285)
+++ pkg/tests/rw2.Rout.save 2010-07-19 16:12:36 UTC (rev 286)
@@ -255,7 +255,7 @@
}
y
}
-<environment: 0x2629720>
+<environment: 0x1831540>
measurement model density, dmeasure =
function (y, x, t, params, log, covars, ...)
{
@@ -268,7 +268,7 @@
f
else exp(f)
}
-<environment: 0x2629720>
+<environment: 0x1831540>
initializer =
function (params, t0, ...)
{
@@ -390,12 +390,48 @@
> max(abs(e2-e3),na.rm=T)
[1] 0
>
+> new <- window(rw2,start=20,end=30)
+> new <- simulate(new)
+>
+> timezero(new)
+[1] 0
+> timezero(new) <- 19
+> print(simulate(new))
+data and states:
+ time y1 y2 x1 x2
+1 20 1.0889842 -2.9306255 0.3928615 -3.6019244
+2 21 0.6190426 -3.5224007 0.3502208 -2.7385549
+3 22 -0.9996197 -3.9952599 -1.3872796 -2.2270139
+4 23 -1.3789056 -2.5477141 -0.4674199 -2.3980597
+5 24 0.2470824 0.0985785 0.4256256 -0.0595737
+6 25 3.5291989 3.9686175 1.2073146 3.1753928
+7 26 2.4044278 5.0563153 0.6048257 4.5525857
+8 27 1.0704905 5.3986777 0.9983714 3.9165344
+9 28 -0.3706092 1.3254624 1.0613907 2.6270702
+10 29 1.2884699 4.0593363 1.4887009 4.8853658
+11 30 0.8533498 -0.7281767 2.3528744 -1.4303754
+
+call:
+pomp(data = rbind(y1 = rep(0, 100), y2 = rep(0, 100)), times = 1:100,
+ t0 = 0, rprocess = onestep.sim(step.fun = function(x, t,
+ params, delta.t, ...) {
+ c(x1 = rnorm(n = 1, mean = x["x1"], sd = params["s1"] *
+ delta.t), x2 = rnorm(n = 1, mean = x["x2"], sd = params["s2"] *
+ delta.t))
+ }), dprocess = onestep.dens(dens.fun = function(x1, t1, x2,
+ t2, params, ...) {
+ sum(dnorm(x = x2[c("x1", "x2")], mean = x1[c("x1", "x2")],
+ sd = params[c("s1", "s2")] * (t2 - t1), log = TRUE),
+ na.rm = TRUE)
+ }), measurement.model = list(y1 ~ norm(mean = x1, sd = tau),
+ y2 ~ norm(mean = x2, sd = tau)))
+>
> time(rw2) <- seq(1,1000,by=20)
> x <- simulate(rw2)
> states(x)[,1:5]
- [,1] [,2] [,3] [,4] [,5]
-x1 0 -0.2414417 -7.85263 -33.53872 -52.31123
-x2 0 -0.9545774 -21.42633 -23.14645 -98.62724
+ [,1] [,2] [,3] [,4] [,5]
+x1 0 0.1308992 -0.8665404 -2.91768 -15.17674
+x2 0 -1.6507047 -30.2956790 -71.13064 -199.25743
> try(
+ time(rw2) <- seq(-20,1000,by=20)
+ )
@@ -414,3 +450,4 @@
> time(rw2) <- c(0,20,25.8,50,60)
>
>
+>
Modified: pkg/tests/sir.R
===================================================================
--- pkg/tests/sir.R 2010-07-06 17:47:12 UTC (rev 285)
+++ pkg/tests/sir.R 2010-07-19 16:12:36 UTC (rev 286)
@@ -238,5 +238,11 @@
states(po)[,1:3]
states(simulate(po))[,1:3]
+po <- window(euler.sir,start=1,end=2)
+plot(simulate(po))
+timezero(po)
+timezero(po)<-2*time(po)[1]-time(po)[2]
+plot(simulate(po))
+
dev.off()
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2010-07-06 17:47:12 UTC (rev 285)
+++ pkg/tests/sir.Rout.save 2010-07-19 16:12:36 UTC (rev 286)
@@ -1,6 +1,6 @@
-R version 2.10.1 (2009-12-14)
-Copyright (C) 2009 The R Foundation for Statistical Computing
+R version 2.11.1 (2010-05-31)
+Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -16,9 +16,6 @@
Type 'q()' to quit R.
> library(pomp)
-Loading required package: deSolve
-Loading required package: subplex
-Loading required package: mvtnorm
>
> tbasis <- seq(0,25,by=1/52)
> basis <- periodic.bspline.basis(tbasis,nbasis=3)
@@ -435,7 +432,7 @@
zeronames = zeronames, tcovar = tcovar, covar = covar,
args = pairlist(...))
}
-<environment: 0xcfaab8>
+<environment: 0x2dfe318>
process model density, dprocess =
function (x, times, params, ..., statenames = character(0), paramnames = character(0),
covarnames = character(0), tcovar, covar, log = FALSE)
@@ -445,7 +442,7 @@
covarnames = covarnames, tcovar = tcovar, covar = covar,
log = log, args = pairlist(...))
}
-<environment: 0xc60008>
+<environment: 0x369b420>
measurement model simulator, rmeasure =
function (x, t, params, covars, ...)
{
@@ -457,7 +454,7 @@
}
y
}
-<environment: 0xfe3fa0>
+<environment: 0x2abb470>
measurement model density, dmeasure =
function (y, x, t, params, log, covars, ...)
{
@@ -470,7 +467,7 @@
f
else exp(f)
}
-<environment: 0xfe3fa0>
+<environment: 0x2abb470>
skeleton ( vectorfield ) =
function (x, t, params, covars, ...)
{
@@ -509,7 +506,7 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 3.034172 secs
+Time difference of 2.071329 secs
>
> pdf(file='sir.pdf')
>
@@ -526,7 +523,7 @@
> X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 7.110227 secs
+Time difference of 4.534647 secs
> plot(t3,X3['I',1,],type='l')
>
> f1 <- dprocess(
@@ -586,7 +583,7 @@
> x <- simulate(po,nsim=100)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 6.164427 secs
+Time difference of 2.7706 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t3 <- seq(0,20,by=1/52)
@@ -594,7 +591,7 @@
> X4 <- trajectory(po,times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 1.135105 secs
+Time difference of 0.7220042 secs
> plot(t3,X4['I',1,],type='l')
>
> g2 <- dmeasure(
@@ -657,6 +654,13 @@
cases 0 0 1.007000e+03
W 0 0 -1.602887e-01
>
+> po <- window(euler.sir,start=1,end=2)
+> plot(simulate(po))
+> timezero(po)
+[1] 0
+> timezero(po)<-2*time(po)[1]-time(po)[2]
+> plot(simulate(po))
+>
> dev.off()
null device
1
More information about the pomp-commits
mailing list