[Pomp-commits] r676 - in pkg/pomp: . R data inst inst/data-R src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 23 13:43:09 CEST 2012
Author: kingaa
Date: 2012-04-23 13:43:09 +0200 (Mon, 23 Apr 2012)
New Revision: 676
Added:
pkg/pomp/tests/blowflies.R
pkg/pomp/tests/blowflies.Rout.save
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/R/builder.R
pkg/pomp/data/bbs.rda
pkg/pomp/data/blowflies.rda
pkg/pomp/data/dacca.rda
pkg/pomp/data/euler.sir.rda
pkg/pomp/data/gillespie.sir.rda
pkg/pomp/data/gompertz.rda
pkg/pomp/data/ou2.rda
pkg/pomp/data/ricker.rda
pkg/pomp/data/rw2.rda
pkg/pomp/data/verhulst.rda
pkg/pomp/inst/NEWS
pkg/pomp/inst/data-R/blowflies.R
pkg/pomp/src/blowfly.c
pkg/pomp/tests/sir.Rout.save
Log:
- remove the delay 'tau' from the blowflies example
- make the 'pompBuilder' work under 2.14.2
- update the data()-loadable examples
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2012-04-21 14:29:05 UTC (rev 675)
+++ pkg/pomp/DESCRIPTION 2012-04-23 11:43:09 UTC (rev 676)
@@ -1,14 +1,14 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.42-1
+Version: 0.41-7
Date: 2012-04-21
Revision: $Rev$
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>
URL: http://pomp.r-forge.r-project.org
Description: Inference methods for partially-observed Markov processes
-Depends: R(>= 2.15.0), stats, methods, graphics, mvtnorm, subplex, deSolve
+Depends: R(>= 2.14.2), stats, methods, graphics, mvtnorm, subplex, deSolve
License: GPL(>= 2)
LazyLoad: true
LazyData: false
Modified: pkg/pomp/R/builder.R
===================================================================
--- pkg/pomp/R/builder.R 2012-04-21 14:29:05 UTC (rev 675)
+++ pkg/pomp/R/builder.R 2012-04-23 11:43:09 UTC (rev 676)
@@ -203,5 +203,5 @@
}
retval[[i]] <- tpl
}
- do.call(paste0,retval)
+ do.call(function(...)paste(...,sep='.'),retval)
}
Modified: pkg/pomp/data/bbs.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/blowflies.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/dacca.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/gompertz.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/ricker.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/rw2.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/data/verhulst.rda
===================================================================
(Binary files differ)
Modified: pkg/pomp/inst/NEWS
===================================================================
--- pkg/pomp/inst/NEWS 2012-04-21 14:29:05 UTC (rev 675)
+++ pkg/pomp/inst/NEWS 2012-04-23 11:43:09 UTC (rev 676)
@@ -1,12 +1,9 @@
NEWS
-0.42-1
- o An experimental facility for constructing pomp objects with native C routines is now included.
- This requires R version 2.15.0 or greater.
-
-0.41-6
+0.41-7
o A bug in the 'blowflies' example has been fixed.
Thanks to Greg Minshall for discovering it.
+0.41-6
o The 'gompertz' example parameter transformations have been changed.
No longer are the names of the parameter vector changed in the transformation.
This change is not backward-compatible, but only codes using this example are affected.
Modified: pkg/pomp/inst/data-R/blowflies.R
===================================================================
--- pkg/pomp/inst/data-R/blowflies.R 2012-04-21 14:29:05 UTC (rev 675)
+++ pkg/pomp/inst/data-R/blowflies.R 2012-04-23 11:43:09 UTC (rev 676)
@@ -18,28 +18,26 @@
times="day",
t0=14,
rprocess=discrete.time.sim(
- step.fun="_blowfly_model_simulator",
+ step.fun="_blowfly_simulator_one",
delta.t=1,
PACKAGE="pomp"
),
- paramnames=c("P","N0","delta","sigma.P","sigma.d","tau","sigma.y"),
+ rmeasure="_blowfly_rmeasure",
+ dmeasure="_blowfly_dmeasure",
+ PACKAGE="pomp",
+ paramnames=c("P","N0","delta","sigma.P","sigma.d","sigma.y"),
statenames=c("N1","R","S","e","eps"),
obsnames=c("y"),
- measurement.model=y~nbinom(mu=N1,size=1/(sigma.y^2)),
y.init=with( ## initial data
raw.data,
approx(x=day,y=y,xout=seq(from=0,to=14,by=1),rule=2)$y
),
# y.init=c(948, 948, 942, 930, 911, 885, 858, 833.7, 801, 748.3, 676, 589.8, 504, 434.9, 397),
parameter.inv.transform=function(params,...) {
- p <- c(log(params[c("delta","sigma.P","sigma.d","sigma.y","P","N0")]),params["tau"])
- names(p) <- c(paste("log",c("delta","sigma.P","sigma.d","sigma.y","P","N0"),sep="."),"tau")
- p
+ log(params)
},
parameter.transform=function(params,...) {
- p <- c(exp(params[paste("log",c("delta","sigma.P","sigma.d","sigma.y","P","N0"),sep=".")]),params["tau"])
- names(p) <- c(c("delta","sigma.P","sigma.d","sigma.y","P","N0"),"tau")
- p
+ exp(params)
},
initializer=function (params, t0, y.init, ...) {
ntau <- length(y.init)
@@ -52,7 +50,7 @@
pomp(
blowflies1,
rprocess=discrete.time.sim(
- step.fun="_blowfly_model_simulator",
+ step.fun="_blowfly_simulator_two",
delta.t=2,
PACKAGE="pomp"
),
@@ -71,24 +69,22 @@
## mle from search to date
coef(blowflies1,transform=TRUE) <- c(
- log.P = 1.189,
- log.delta = -1.828,
- log.N0 = 6.522,
- log.sigma.P = 0.301,
- log.sigma.d = -0.292,
- log.sigma.y = -3.625,
- tau = 14
+ P = 1.189,
+ delta = -1.828,
+ N0 = 6.522,
+ sigma.P = 0.301,
+ sigma.d = -0.292,
+ sigma.y = -3.625
)
## mle from search to date
coef(blowflies2,transform=TRUE) <- c(
- log.P = 1.005,
- log.delta = -1.75,
- log.N0 = 6.685,
- log.sigma.P = 0.366,
- log.sigma.d = -0.274,
- log.sigma.y = -4.524,
- tau = 7
+ P = 1.005,
+ delta = -1.75,
+ N0 = 6.685,
+ sigma.P = 0.366,
+ sigma.d = -0.274,
+ sigma.y = -4.524
)
test <- FALSE
Modified: pkg/pomp/src/blowfly.c
===================================================================
--- pkg/pomp/src/blowfly.c 2012-04-21 14:29:05 UTC (rev 675)
+++ pkg/pomp/src/blowfly.c 2012-04-23 11:43:09 UTC (rev 676)
@@ -9,7 +9,7 @@
#define DELTA (p[parindex[2]]) // survival parameter
#define SIGMAP (p[parindex[3]]) // recruitment noise SD
#define SIGMAD (p[parindex[4]]) // survivorship noise SD
-#define TAU (p[parindex[5]]) // delay
+#define SIGMAY (p[parindex[5]]) // survivorship noise SD
#define N (&x[stateindex[0]]) // total population
#define R (x[stateindex[1]]) // recruits
@@ -17,21 +17,18 @@
#define E (x[stateindex[3]]) // recruitment noise
#define EPS (x[stateindex[4]]) // survival noise
+#define Y (y[obsindex[0]]) // observation
+
// Ricker model with log-normal process noise
-void _blowfly_model_simulator (double *x, const double *p,
- const int *stateindex, const int *parindex, const int *covindex,
- int covdim, const double *covar,
- double t, double dt)
+static void _blowfly_simulator (double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar,
+ double t, double dt, int tau)
{
double e = rgammawn(SIGMAP,dt)/dt;
double eps = rgammawn(SIGMAD,dt)/dt;
- int tau;
int k;
- if (!(R_FINITE(TAU)))
- error("non-finite value of 'tau'");
-
- tau = nearbyint(TAU);
R = rpois(P*N[tau]*exp(-N[tau]/NZERO)*dt*e);
S = rbinom(N[0],exp(-DELTA*dt*eps));
E = e;
@@ -40,16 +37,46 @@
N[0] = R+S;
}
+void _blowfly_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ double size = 1.0/SIGMAY/SIGMAY;
+ double prob = size/(size+N[0]);
+ *lik = dnbinom(Y,size,prob,give_log);
+}
+
+void _blowfly_rmeasure (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ double size = 1.0/SIGMAY/SIGMAY;
+ double prob = size/(size+N[0]);
+ Y = rnbinom(size,prob);
+}
+
#undef N
#undef R
#undef S
#undef E
#undef EPS
-#undef LOG_P
-#undef LOG_NZERO
-#undef LOG_DELTA
-#undef LOG_SIGMAP
-#undef LOG_SIGMAD
-#undef TAU
+#undef P
+#undef NZERO
+#undef DELTA
+#undef SIGMAP
+#undef SIGMAD
+#undef SIGMAY
+void _blowfly_simulator_one (double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar,
+ double t, double dt) {
+ _blowfly_simulator(x, p, stateindex, parindex, covindex, covdim, covar, t, dt, 14);
+}
+
+void _blowfly_simulator_two (double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar,
+ double t, double dt) {
+ _blowfly_simulator(x, p, stateindex, parindex, covindex, covdim, covar, t, dt, 7);
+}
+
Added: pkg/pomp/tests/blowflies.R
===================================================================
--- pkg/pomp/tests/blowflies.R (rev 0)
+++ pkg/pomp/tests/blowflies.R 2012-04-23 11:43:09 UTC (rev 676)
@@ -0,0 +1,13 @@
+library(pomp)
+
+data(blowflies)
+
+init.state(blowflies1)
+x1 <- simulate(blowflies1)
+f1 <- pfilter(blowflies1,Np=1000,seed=599688L)
+logLik(f1)
+
+init.state(blowflies2)
+x2 <- simulate(blowflies2)
+f2 <- pfilter(blowflies2,Np=1000,seed=599688L)
+logLik(f2)
Added: pkg/pomp/tests/blowflies.Rout.save
===================================================================
--- pkg/pomp/tests/blowflies.Rout.save (rev 0)
+++ pkg/pomp/tests/blowflies.Rout.save 2012-04-23 11:43:09 UTC (rev 676)
@@ -0,0 +1,70 @@
+
+R version 2.14.2 (2012-02-29)
+Copyright (C) 2012 The R Foundation for Statistical Computing
+ISBN 3-900051-07-0
+Platform: x86_64-unknown-linux-gnu (64-bit)
+
+R is free software and comes with ABSOLUTELY NO WARRANTY.
+You are welcome to redistribute it under certain conditions.
+Type 'license()' or 'licence()' for distribution details.
+
+R is a collaborative project with many contributors.
+Type 'contributors()' for more information and
+'citation()' on how to cite R or R packages in publications.
+
+Type 'demo()' for some demos, 'help()' for on-line help, or
+'help.start()' for an HTML browser interface to help.
+Type 'q()' to quit R.
+
+> library(pomp)
+Loading required package: mvtnorm
+Loading required package: subplex
+Loading required package: deSolve
+>
+> data(blowflies)
+>
+> init.state(blowflies1)
+ [,1]
+N1 397.0
+N2 450.5
+N3 504.0
+N4 590.0
+N5 676.0
+N6 738.5
+N7 801.0
+N8 829.5
+N9 858.0
+N10 884.5
+N11 911.0
+N12 926.5
+N13 942.0
+N14 945.0
+N15 948.0
+R 0.0
+S 0.0
+e 0.0
+eps 0.0
+> x1 <- simulate(blowflies1)
+> f1 <- pfilter(blowflies1,Np=1000,seed=599688L)
+> logLik(f1)
+[1] -1466.694
+>
+> init.state(blowflies2)
+ [,1]
+N1 397
+N2 504
+N3 676
+N4 801
+N5 858
+N6 911
+N7 942
+N8 948
+R 0
+S 0
+e 0
+eps 0
+> x2 <- simulate(blowflies2)
+> f2 <- pfilter(blowflies2,Np=1000,seed=599688L)
+> logLik(f2)
+[1] -1471.781
+>
Modified: pkg/pomp/tests/sir.Rout.save
===================================================================
--- pkg/pomp/tests/sir.Rout.save 2012-04-21 14:29:05 UTC (rev 675)
+++ pkg/pomp/tests/sir.Rout.save 2012-04-23 11:43:09 UTC (rev 676)
@@ -389,7 +389,7 @@
zeronames = zeronames, tcovar = tcovar, covar = covar,
args = pairlist(...))
}
-<environment: 0x25543e0>
+<environment: 0x26bb410>
process model density, dprocess =
function (x, times, params, ..., statenames = character(0), paramnames = character(0),
covarnames = character(0), tcovar, covar, log = FALSE)
@@ -399,7 +399,7 @@
covarnames = covarnames, tcovar = tcovar, covar = covar,
log = log, args = pairlist(...))
}
-<environment: 0x24e7238>
+<environment: 0x264e268>
measurement model simulator, rmeasure =
function (x, t, params, covars, ...)
{
@@ -463,7 +463,7 @@
> x <- simulate(po,params=params,nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 1.024599 secs
+Time difference of 1.033125 secs
>
> pdf(file='sir.pdf')
>
@@ -480,7 +480,7 @@
> X3 <- trajectory(po,params=params,times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.7358234 secs
+Time difference of 0.7432306 secs
> plot(t3,X3['I',1,],type='l')
>
> f1 <- dprocess(
@@ -530,7 +530,7 @@
> x <- simulate(po,nsim=100)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 1.200051 secs
+Time difference of 1.196313 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t3 <- seq(0,20,by=1/52)
@@ -538,7 +538,7 @@
> X4 <- trajectory(po,times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.01725101 secs
+Time difference of 0.01725578 secs
> plot(t3,X4['I',1,],type='l')
>
> g2 <- dmeasure(
@@ -575,18 +575,18 @@
>
> states(po)[,1:2]
[,1] [,2]
-S 1.361010e+05 1.357700e+05
-I 2.064000e+03 2.065000e+03
-R 1.961799e+06 1.962133e+06
-cases 1.062000e+03 1.060000e+03
+S 1.361160e+05 1.357770e+05
+I 2.066000e+03 2.074000e+03
+R 1.961786e+06 1.962120e+06
+cases 1.049000e+03 1.060000e+03
W 6.596952e-02 9.166551e-02
> time(po) <- seq(0,1,by=1/52)
> states(po)[,1:3]
[,1] [,2] [,3]
-S NA 1.361010e+05 1.357700e+05
-I NA 2.064000e+03 2.065000e+03
-R NA 1.961799e+06 1.962133e+06
-cases NA 1.062000e+03 1.060000e+03
+S NA 1.361160e+05 1.357770e+05
+I NA 2.066000e+03 2.074000e+03
+R NA 1.961786e+06 1.962120e+06
+cases NA 1.049000e+03 1.060000e+03
W NA 6.596952e-02 9.166551e-02
> states(simulate(po))[,1:3]
[,1] [,2] [,3]
@@ -610,4 +610,4 @@
>
> proc.time()
user system elapsed
- 3.920 0.036 4.195
+ 3.960 0.020 4.143
More information about the pomp-commits
mailing list