[Pomp-commits] r303 - in pkg: . data inst inst/data-R inst/doc man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Aug 24 20:06:44 CEST 2010
Author: kingaa
Date: 2010-08-24 20:06:44 +0200 (Tue, 24 Aug 2010)
New Revision: 303
Added:
pkg/data/ricker.rda
pkg/inst/data-R/ricker.R
pkg/man/ricker.Rd
pkg/src/ricker.c
pkg/tests/ricker.R
pkg/tests/ricker.Rout.save
Modified:
pkg/DESCRIPTION
pkg/data/dacca.rda
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/data/ou2.rda
pkg/inst/ChangeLog
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.Rnw
pkg/inst/doc/intro_to_pomp.pdf
Log:
- include a Poisson-Ricker example
- update the vignettes
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-08-24 18:06:08 UTC (rev 302)
+++ pkg/DESCRIPTION 2010-08-24 18:06:44 UTC (rev 303)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.32-1
-Date: 2010-08-23
+Version: 0.32-2
+Date: 2010-08-24
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/data/dacca.rda
===================================================================
(Binary files differ)
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Added: pkg/data/ricker.rda
===================================================================
(Binary files differ)
Property changes on: pkg/data/ricker.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2010-08-24 18:06:08 UTC (rev 302)
+++ pkg/inst/ChangeLog 2010-08-24 18:06:44 UTC (rev 303)
@@ -1,3 +1,14 @@
+2010-08-24 kingaa
+
+ * [r299] inst/doc/pomp.bib: - update bibliography
+
+2010-08-23 kingaa
+
+ * [r297] R/basic-probes.R: - another minor bug fix
+ * [r296] R/basic-probes.R: - bug fix in 'probe.acf'
+ * [r295] DESCRIPTION:
+ * [r293] inst/ChangeLog, inst/NEWS: - update NEWS and ChangeLog
+
2010-08-19 kingaa
* [r291] DESCRIPTION, NAMESPACE, R/bsmc.R, tests/ou2-probe-match.R,
Added: pkg/inst/data-R/ricker.R
===================================================================
--- pkg/inst/data-R/ricker.R (rev 0)
+++ pkg/inst/data-R/ricker.R 2010-08-24 18:06:44 UTC (rev 303)
@@ -0,0 +1,27 @@
+require(pomp)
+
+simulate(
+ pomp(
+ data=data.frame(time=seq(0,50,by=1),y=NA),
+ times="time",
+ t0=0,
+ rprocess=discrete.time.sim(
+ step.fun="ricker_simulator"
+ ),
+ rmeasure="poisson_rmeasure",
+ dmeasure="poisson_dmeasure",
+ skeleton.map="ricker_skeleton",
+ paramnames=c("log.r","log.sigma","log.phi"),
+ statenames=c("N","e"),
+ obsnames=c("y")
+ ),
+ params=c(
+ log.r=3.8,
+ log.sigma=log(0.3),
+ log.phi=log(10),
+ N.0=7,
+ e.0=0
+ ),
+ nsim=1,
+ seed=73691676L
+ ) -> ricker
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.Rnw
===================================================================
--- pkg/inst/doc/intro_to_pomp.Rnw 2010-08-24 18:06:08 UTC (rev 302)
+++ pkg/inst/doc/intro_to_pomp.Rnw 2010-08-24 18:06:44 UTC (rev 303)
@@ -236,8 +236,7 @@
),
times="time",
rprocess=discrete.time.sim(
- step.fun=ou2.proc.sim,
- delta.t=1
+ step.fun=ou2.proc.sim
),
rmeasure=ou2.meas.sim,
t0=0
@@ -249,8 +248,8 @@
The second argument (\code{times}) specifies which of the columns of \code{data} is the time variable.
The third argument (\code{rprocess}) specifies that the process model simulator will be in discrete-time, one step at a time.
The function \code{discrete.time.sim} belongs to the package \pomp.
-It takes several arguments, one of which, \code{step.fun}, is the particular function that actually takes the step.
-Its second argument, \code{delta.t}, sets the size of the time-step.
+It takes the argument \code{step.fun}, which specifies the particular function that actually takes the step.
+The step is assumed to be a unit interval of time.
The argument \code{rmeasure} specifies the measurement model simulator function.
\code{t0} fixes $t_0$ for this model; here we have chosen this to be one time unit before the first observation.
@@ -350,7 +349,7 @@
ou2 <- pomp(
data=dat[c("time","y1","y2")],
times="time",
- rprocess=discrete.time.sim(ou2.proc.sim,delta.t=1),
+ rprocess=discrete.time.sim(ou2.proc.sim),
rmeasure=ou2.meas.sim,
dmeasure=ou2.meas.dens,
t0=0
@@ -704,7 +703,7 @@
new.ou2 <- pomp(
data=dat[c("time","y1","y2")],
times="time",
- rprocess=discrete.time.sim(ou2.proc.sim,delta.t=1),
+ rprocess=discrete.time.sim(ou2.proc.sim),
rmeasure=ou2.meas.sim,
dmeasure=ou2.meas.dens,
skeleton.map=ou2.skel,
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Added: pkg/man/ricker.Rd
===================================================================
--- pkg/man/ricker.Rd (rev 0)
+++ pkg/man/ricker.Rd 2010-08-24 18:06:44 UTC (rev 303)
@@ -0,0 +1,19 @@
+\name{ricker}
+\alias{ricker}
+\docType{data}
+\title{Ricker model with Poisson observations.}
+\description{
+ \code{ricker} is a \code{pomp} object encoding a stochastic Ricker model with Poisson measurement error.
+}
+\usage{data(ricker)}
+\details{
+ The state process is \eqn{N_{t+1} = r N_{t} exp(-N_{t}+e_{t})}, where the \eqn{e_t} are i.i.d. normal random deviates with variance \eqn{sigma^2}.
+ The observed variables \eqn{y_t} are distributed as \eqn{Poisson(\phi N_t)}.
+}
+\examples{
+data(ricker)
+plot(ricker)
+coef(ricker)
+}
+\seealso{\code{\link{pomp-class}} and the vignettes}
+\keyword{datasets}
Added: pkg/src/ricker.c
===================================================================
--- pkg/src/ricker.c (rev 0)
+++ pkg/src/ricker.c 2010-08-24 18:06:44 UTC (rev 303)
@@ -0,0 +1,55 @@
+// dear emacs, please treat this as -*- C++ -*-
+
+#include <Rmath.h>
+
+#include "pomp.h"
+
+#define LOG_R (p[parindex[0]]) // growth rate
+#define LOG_SIGMA (p[parindex[1]]) // process noise level
+#define LOG_PHI (p[parindex[2]]) // measurement scale parameter
+
+#define N (x[stateindex[0]]) // population size
+#define E (x[stateindex[1]]) // process noise
+
+#define Y (y[obsindex[0]]) // observed population size
+
+void poisson_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) {
+ *lik = dpois(Y,exp(LOG_PHI)*N,give_log);
+}
+
+void poisson_rmeasure (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t) {
+ Y = rpois(exp(LOG_PHI)*N);
+}
+
+#undef Y
+
+// Ricker model with log-normal process noise
+void ricker_simulator (double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar,
+ double t, double dt)
+{
+ double sigma = exp(LOG_SIGMA);
+ double e = (sigma > 0.0) ? rnorm(0,sigma) : 0.0;
+ N = exp(LOG_R+log(N)-N+e);
+ E = e;
+}
+
+void ricker_skeleton (double *f, double *x, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar, double t)
+{
+ f[0] = exp(LOG_R+log(N)-N);
+ f[1] = 0.0;
+}
+
+#undef N
+#undef E
+
+#undef LOG_R
+#undef LOG_SIGMA
+#undef LOG_PHI
Property changes on: pkg/src/ricker.c
___________________________________________________________________
Added: svn:eol-style
+ native
Added: pkg/tests/ricker.R
===================================================================
--- pkg/tests/ricker.R (rev 0)
+++ pkg/tests/ricker.R 2010-08-24 18:06:44 UTC (rev 303)
@@ -0,0 +1,37 @@
+library(pomp)
+
+data(ricker)
+
+set.seed(64857673L)
+
+po <- ricker
+coef(po) <- c(log(2),log(1),log(20),7,0)
+
+pf.tr <- pfilter(ricker,Np=1000,max.fail=50)
+pf.po <- pfilter(po,Np=1000,max.fail=50)
+
+mf <- mif(
+ po,
+ Nmif=100,
+ Np=1000,
+ cooling.factor=0.99,
+ var.factor=2,
+ ic.lag=3,
+ max.fail=50,
+ rw.sd=c(log.r=0.1,log.sigma=0.1,log.phi=0.1)
+ )
+plot(mf <- mif(mf,Nmif=500,max.fail=20))
+pf.mf <- pfilter(mf,Np=1000)
+
+res <- rbind(
+ cbind(truth=coef(ricker),MLE=coef(mf),guess=coef(po)),
+ loglik=c(pf.tr$loglik,pf.mf$loglik,pf.po$loglik)
+ )
+
+print(res,digits=3)
+
+tj.1 <- trajectory(ricker)
+plot(time(ricker),tj.1[1,,-1],type='l')
+tj.2 <- trajectory(ricker,times=c(0,30:50))
+lines(30:50,tj.2[1,,-1],col='red',lwd=2)
+max(abs(tj.1[,,time(ricker,t0=T)>=30]-tj.2[,,-1]))
Added: pkg/tests/ricker.Rout.save
===================================================================
--- pkg/tests/ricker.Rout.save (rev 0)
+++ pkg/tests/ricker.Rout.save 2010-08-24 18:06:44 UTC (rev 303)
@@ -0,0 +1,63 @@
+
+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.
+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)
+>
+> data(ricker)
+>
+> set.seed(64857673L)
+>
+> po <- ricker
+> coef(po) <- c(log(2),log(1),log(20),7,0)
+>
+> pf.tr <- pfilter(ricker,Np=1000,max.fail=50)
+> pf.po <- pfilter(po,Np=1000,max.fail=50)
+>
+> mf <- mif(
++ po,
++ Nmif=100,
++ Np=1000,
++ cooling.factor=0.99,
++ var.factor=2,
++ ic.lag=3,
++ max.fail=50,
++ rw.sd=c(log.r=0.1,log.sigma=0.1,log.phi=0.1)
++ )
+> plot(mf <- mif(mf,Nmif=500,max.fail=20))
+> pf.mf <- pfilter(mf,Np=1000)
+>
+> res <- rbind(
++ cbind(truth=coef(ricker),MLE=coef(mf),guess=coef(po)),
++ loglik=c(pf.tr$loglik,pf.mf$loglik,pf.po$loglik)
++ )
+>
+> print(res,digits=3)
+ truth MLE guess
+log.r 3.80 3.73 0.693
+log.sigma -1.20 -1.38 0.000
+log.phi 2.30 2.31 2.996
+N.0 7.00 7.00 7.000
+e.0 0.00 0.00 0.000
+loglik -142.74 -143.25 -266.548
+>
+> tj.1 <- trajectory(ricker)
+> plot(time(ricker),tj.1[1,,-1],type='l')
+> tj.2 <- trajectory(ricker,times=c(0,30:50))
+> lines(30:50,tj.2[1,,-1],col='red',lwd=2)
+> max(abs(tj.1[,,time(ricker,t0=T)>=30]-tj.2[,,-1]))
+[1] 0
+>
More information about the pomp-commits
mailing list