[Pomp-commits] r159 - in pkg: R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Sep 21 19:52:19 CEST 2009
Author: kingaa
Date: 2009-09-21 19:52:18 +0200 (Mon, 21 Sep 2009)
New Revision: 159
Modified:
pkg/R/aaa.R
pkg/R/pomp-methods.R
pkg/man/pomp-methods.Rd
pkg/tests/rw2.R
pkg/tests/rw2.Rout.save
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
- add a new "time<-" method for manipulating/subsetting time series
Modified: pkg/R/aaa.R
===================================================================
--- pkg/R/aaa.R 2009-07-08 16:58:19 UTC (rev 158)
+++ pkg/R/aaa.R 2009-09-21 17:52:18 UTC (rev 159)
@@ -64,6 +64,10 @@
stop("function ",sQuote("data.array")," is undefined for objects of class ",sQuote(class(object)))
setGeneric('data.array')
+"time<-" <- function (object, ..., value)
+ stop("function ",sQuote("time<-")," is undefined for objects of class ",sQuote(class(object)))
+setGeneric("time<-")
+
rprocess <- function (object, xstart, times, params, ...)
stop("function ",sQuote("rprocess")," is undefined for objects of class ",sQuote(class(object)))
setGeneric('rprocess')
@@ -122,9 +126,9 @@
stop("function ",sQuote("continue")," is undefined for objects of class ",sQuote(class(object)))
setGeneric('continue')
-'coef<-' <- function (object, pars, ..., value)
- stop("function ",sQuote("coef<")," is undefined for objects of class ",sQuote(class(object)))
-setGeneric('coef<-')
+"coef<-" <- function (object, pars, ..., value)
+ stop("function ",sQuote("coef<-")," is undefined for objects of class ",sQuote(class(object)))
+setGeneric("coef<-")
states <- function (object, ...)
stop("function ",sQuote("states")," is undefined for objects of class ",sQuote(class(object)))
Modified: pkg/R/pomp-methods.R
===================================================================
--- pkg/R/pomp-methods.R 2009-07-08 16:58:19 UTC (rev 158)
+++ pkg/R/pomp-methods.R 2009-09-21 17:52:18 UTC (rev 159)
@@ -47,6 +47,52 @@
}
)
+## modify the time
+setMethod(
+ "time<-",
+ "pomp",
+ function (object, include.t0 = FALSE, ..., value) {
+ if (!is.numeric(value))
+ stop(sQuote("value")," must be a numeric vector",call.=TRUE)
+ storage.mode(value) <- "double"
+ tt0 <- object at t0
+ tt <- object at times
+ dd <- object at data
+ ss <- object at states
+ if (include.t0) {
+ object at t0 <- value[1]
+ object at times <- value[-1]
+ } else {
+ object at times <- value
+ }
+ if (!all(diff(object at times)>0))
+ stop("the times specified must be an increasing sequence",call.=TRUE)
+ if (object at t0>object at times[1])
+ stop("the zero-time ",sQuote("t0")," must occur no later than the first observation",call.=TRUE)
+ object at data <- array(
+ data=NA,
+ dim=c(nrow(dd),length(object at times)),
+ dimnames=list(rownames(dd),NULL)
+ )
+ object at data[,object at times%in%tt] <- dd[,tt%in%object at times]
+ if (length(ss)>0) {
+ object at states <- array(
+ data=NA,
+ dim=c(nrow(ss),length(object at times)+1),
+ dimnames=list(rownames(ss),NULL)
+ )
+ if (ncol(ss)>length(tt)) tt <- c(tt0,tt)
+ nt <- c(object at t0,object at times)
+ for (kt in seq(length=length(nt))) {
+ wr <- which(nt[kt]==tt)
+ if (length(wr)>0)
+ object at states[,kt] <- ss[,wr[1]]
+ }
+ }
+ object
+ }
+ )
+
## extract the coefficients
setMethod(
'coef',
Modified: pkg/man/pomp-methods.Rd
===================================================================
--- pkg/man/pomp-methods.Rd 2009-07-08 16:58:19 UTC (rev 158)
+++ pkg/man/pomp-methods.Rd 2009-09-21 17:52:18 UTC (rev 159)
@@ -9,6 +9,9 @@
\alias{show-pomp}
\alias{time,pomp-method}
\alias{time-pomp}
+\alias{time<-}
+\alias{time<-,pomp-method}
+\alias{time<--pomp}
\alias{data.array,pomp-method}
\alias{data.array-pomp}
\alias{data.array}
@@ -30,6 +33,7 @@
\S4method{data.array}{pomp}(object, vars, \dots)
\S4method{states}{pomp}(object, vars, \dots)
\S4method{time}{pomp}(x, t0 = FALSE, \dots)
+\S4method{time}{pomp}(object, include.t0 = FALSE, \dots) <- value
\S4method{show}{pomp}(object)
\S4method{as}{pomp}(object, class)
\S4method{coerce}{pomp,data.frame}(from, to = "data.frame", strict = TRUE)
@@ -57,6 +61,10 @@
logical;
if TRUE, the zero time is prepended to the time vector.
}
+ \item{include.t0}{
+ logical;
+ if TRUE, the first element in \code{value} is taken to be the initial time.
+ }
\item{class}{
character;
name of the class to which \code{object} should be coerced.
@@ -126,6 +134,12 @@
\code{time(object)} returns the vector of observation times.
\code{time(object,t0=TRUE)} returns the vector of observation times with the zero-time \code{t0} prepended.
}
+ \item{time<-}{
+ \code{time(object) <- value} replaces the observation times slot (\code{times}) of \code{object} with \code{value}.
+ \code{time(object,include.t0 = TRUE) <- value} has the same effect, but the first element in \code{value} is taken to be the initial time.
+ 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{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 2009-07-08 16:58:19 UTC (rev 158)
+++ pkg/tests/rw2.R 2009-09-21 17:52:18 UTC (rev 159)
@@ -177,3 +177,20 @@
max(abs(b2-b3))
max(abs(d2-d3),na.rm=T)
max(abs(e2-e3),na.rm=T)
+
+time(rw2) <- seq(1,1000,by=20)
+x <- simulate(rw2)
+states(x)[,1:5]
+try(
+ time(rw2) <- seq(-20,1000,by=20)
+ )
+try(
+ time(rw2) <- c(0,5,10,15,12,20)
+ )
+time(rw2,include.t0=TRUE) <- seq(-20,1000,by=20)
+x <- simulate(rw2)
+time(rw2) <- c(0,20,25.8,50,60)
+time(rw2,include.t0=TRUE) <- c(0,20,25.8,50,60)
+time(rw2,include.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 2009-07-08 16:58:19 UTC (rev 158)
+++ pkg/tests/rw2.Rout.save 2009-09-21 17:52:18 UTC (rev 159)
@@ -1,6 +1,6 @@
-R version 2.8.1 (2008-12-22)
-Copyright (C) 2008 The R Foundation for Statistical Computing
+R version 2.9.1 (2009-06-26)
+Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -251,7 +251,7 @@
}
y
}
-<environment: 0x17d8bd8>
+<environment: 0x27c8440>
measurement model density, dmeasure =
function (y, x, t, params, log, covars, ...)
{
@@ -264,7 +264,7 @@
f
else exp(f)
}
-<environment: 0x17d8bd8>
+<environment: 0x27c8440>
initializer =
function (params, t0, ...)
{
@@ -282,7 +282,7 @@
+ )
Error in try(.Call(do_init_state, object, as.matrix(params), t0), silent = FALSE) :
a state variable and a parameter share a single name: 'x1.0'
-Error : init.state error: error in user ‘initializer’
+Error : init.state error: error in user 'initializer'
>
> rw2 <- pomp(
+ rprocess = rw.rprocess,
@@ -386,3 +386,27 @@
> max(abs(e2-e3),na.rm=T)
[1] 0
>
+> 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
+> try(
++ time(rw2) <- seq(-20,1000,by=20)
++ )
+Error in .local(object, ..., value = value) :
+ the zero-time 't0' must occur no later than the first observation
+> try(
++ time(rw2) <- c(0,5,10,15,12,20)
++ )
+Error in .local(object, ..., value = value) :
+ the times specified must be an increasing sequence
+> time(rw2,include.t0=TRUE) <- seq(-20,1000,by=20)
+> x <- simulate(rw2)
+> time(rw2) <- c(0,20,25.8,50,60)
+> time(rw2,include.t0=TRUE) <- c(0,20,25.8,50,60)
+> time(rw2,include.t0=TRUE) <- c(0,0,20,25.8,50,60)
+> time(rw2) <- c(0,20,25.8,50,60)
+>
+>
Modified: pkg/tests/sir.R
===================================================================
--- pkg/tests/sir.R 2009-07-08 16:58:19 UTC (rev 158)
+++ pkg/tests/sir.R 2009-09-21 17:52:18 UTC (rev 159)
@@ -246,4 +246,11 @@
print(max(abs(g2-g1),na.rm=T),digits=4)
print(max(abs(h2-h1),na.rm=T),digits=4)
+data(euler.sir)
+states(euler.sir)[,1:2]
+time(euler.sir) <- seq(0,1,by=1/52)
+states(euler.sir)[,1:3]
+states(simulate(euler.sir))[,1:3]
+
dev.off()
+
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2009-07-08 16:58:19 UTC (rev 158)
+++ pkg/tests/sir.Rout.save 2009-09-21 17:52:18 UTC (rev 159)
@@ -1,6 +1,6 @@
-R version 2.8.1 (2008-12-22)
-Copyright (C) 2008 The R Foundation for Statistical Computing
+R version 2.9.1 (2009-06-26)
+Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -423,7 +423,7 @@
}
y
}
-<environment: 0x167da18>
+<environment: 0x1203c28>
measurement model density, dmeasure =
function (y, x, t, params, log, covars, ...)
{
@@ -436,7 +436,7 @@
f
else exp(f)
}
-<environment: 0x167da18>
+<environment: 0x1203c28>
initializer =
function (params, t0, ...)
{
@@ -507,7 +507,7 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 4.143093 secs
+Time difference of 1.975832 secs
>
> pdf(file='sir.pdf')
>
@@ -524,7 +524,7 @@
> X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 9.751354 secs
+Time difference of 4.516915 secs
> plot(t3,X3['I',1,],type='l')
>
> f1 <- dprocess(
@@ -1097,7 +1097,7 @@
> x <- simulate(euler.sir,nsim=100)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 2.693998 secs
+Time difference of 1.300607 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t3 <- seq(0,20,by=1/52)
@@ -1105,7 +1105,7 @@
> X4 <- trajectory(euler.sir,times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 7.67392 secs
+Time difference of 3.551238 secs
> plot(t3,X4['I',1,],type='l')
>
> f2 <- dprocess(
@@ -1162,7 +1162,53 @@
> print(max(abs(h2-h1),na.rm=T),digits=4)
[1] 8.731e-11
>
+> data(euler.sir)
+> states(euler.sir)[,1:2]
+ [,1] [,2]
+S 45500 4.531800e+04
+I 2100 2.038000e+03
+R 2052400 2.052647e+06
+cases 0 1.045000e+03
+W 0 -2.334331e-01
+B 0 5.200000e+01
+SI 0 4.900000e+01
+SD 0 0.000000e+00
+IR 0 5.800000e+01
+ID 0 0.000000e+00
+RD 0 4.000000e+01
+dW 0 9.808045e-04
+> time(euler.sir) <- seq(0,1,by=1/52)
+> states(euler.sir)[,1:3]
+ [,1] [,2] [,3]
+S 45500 45500 4.531800e+04
+I 2100 2100 2.038000e+03
+R 2052400 2052400 2.052647e+06
+cases 0 0 1.045000e+03
+W 0 0 -2.334331e-01
+B 0 0 5.200000e+01
+SI 0 0 4.900000e+01
+SD 0 0 0.000000e+00
+IR 0 0 5.800000e+01
+ID 0 0 0.000000e+00
+RD 0 0 4.000000e+01
+dW 0 0 9.808045e-04
+> states(simulate(euler.sir))[,1:3]
+ [,1] [,2] [,3]
+S 45500 45500 4.527200e+04
+I 2100 2100 2.030000e+03
+R 2052400 2052400 2.052644e+06
+cases 0 0 1.034000e+03
+W 0 0 1.415062e-01
+B 0 0 3.700000e+01
+SI 0 0 4.300000e+01
+SD 0 0 0.000000e+00
+IR 0 0 5.300000e+01
+ID 0 0 0.000000e+00
+RD 0 0 3.700000e+01
+dW 0 0 9.553724e-04
+>
> dev.off()
null device
1
>
+>
More information about the pomp-commits
mailing list