[Pomp-commits] r1016 - in pkg/pomp: . inst/examples src

[noreply at r-forge.r-project.org]

Author: kingaa
Date: 2014-12-10 18:12:08 +0100 (Wed, 10 Dec 2014)
New Revision: 1016

Added:
   pkg/pomp/inst/examples/parus.R
   pkg/pomp/src/parus.c
Modified:
   pkg/pomp/DESCRIPTION
Log:
- add Parus major example codes

Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION	2014-12-05 17:44:32 UTC (rev 1015)
+++ pkg/pomp/DESCRIPTION	2014-12-10 17:12:08 UTC (rev 1016)
@@ -1,8 +1,8 @@
 Package: pomp
 Type: Package
 Title: Statistical Inference for Partially Observed Markov Processes
-Version: 0.55-1
-Date: 2014-12-03
+Version: 0.55-2
+Date: 2014-12-10
 Authors at R: c(person(given=c("Aaron","A."),family="King",
 		role=c("aut","cre"),email="kingaa at umich.edu"),
 	  person(given=c("Edward","L."),family="Ionides",role=c("aut")),

Added: pkg/pomp/inst/examples/parus.R
===================================================================
--- pkg/pomp/inst/examples/parus.R	                        (rev 0)
+++ pkg/pomp/inst/examples/parus.R	2014-12-10 17:12:08 UTC (rev 1016)
@@ -0,0 +1,74 @@
+require(pomp)
+
+dat <- 'year,pop
+1960,148
+1961,258
+1962,185
+1963,170
+1964,267
+1965,239
+1966,196
+1967,132
+1968,167
+1969,186
+1970,128
+1971,227
+1972,174
+1973,177
+1974,137
+1975,172
+1976,119
+1977,226
+1978,166
+1979,161
+1980,199
+1981,306
+1982,206
+1983,350
+1984,214
+1985,175
+1986,211
+'
+
+dat <- read.csv(text=dat)
+
+pomp(
+     data=dat,
+     times="year",
+     t0=1960,
+     params=c(K=190,r=2.7,sigma=0.2,tau=0.05,N.0=148),
+     rprocess=discrete.time.sim(
+       step.fun="_parus_gompertz_simulator"
+       ),
+     rmeasure="_parus_lognormal_rmeasure",
+     dmeasure="_parus_lognormal_dmeasure",
+     skeleton="_parus_gompertz_skeleton",
+     skeleton.type="map",
+     paramnames=c("r","K","sigma","tau"),
+     statenames=c("N"),
+     obsnames=c("pop"),
+     parameter.transform=function(params,...){
+       exp(params)
+     },
+     parameter.inv.transform=function(params,...){
+       log(params)
+     }
+     ) -> parusG
+
+pomp(
+     parusG,
+     rprocess=discrete.time.sim(
+       step.fun="_parus_ricker_simulator"
+       ),
+     rmeasure="_parus_poisson_rmeasure",
+     dmeasure="_parus_poisson_dmeasure",
+     skeleton="_parus_ricker_skeleton",
+     skeleton.type="map",
+#     paramnames=c("r","K","sigma","tau"),
+#     statenames=c("N"),
+#     obsnames=c("pop"),
+     PACKAGE="pomp"
+     ) -> parusR
+
+c("parusG","parusR")
+

Added: pkg/pomp/src/parus.c
===================================================================
--- pkg/pomp/src/parus.c	                        (rev 0)
+++ pkg/pomp/src/parus.c	2014-12-10 17:12:08 UTC (rev 1016)
@@ -0,0 +1,75 @@
+// dear emacs, please treat this as -*- C++ -*-
+
+#include <Rmath.h>
+
+#include "pomp.h"
+
+#define R       (p[parindex[0]]) // growth rate
+#define K       (p[parindex[1]]) // carrying capacity
+#define SIGMA   (p[parindex[2]]) // process noise level
+#define TAU     (p[parindex[3]]) // measurement noise level
+
+#define POP         (y[obsindex[0]])
+#define N           (x[stateindex[0]])
+#define NPRIME      (f[stateindex[0]])
+
+void _parus_lognormal_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 = dlnorm(POP,log(N),TAU,give_log);
+}
+
+void _parus_lognormal_rmeasure (double *y, double *x, double *p, 
+				int *obsindex, int *stateindex, int *parindex, int *covindex,
+				int ncovars, double *covars, double t) {
+  POP = rlnorm(log(N),TAU);
+}
+
+void _parus_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(POP,N,give_log);
+}
+
+void _parus_poisson_rmeasure (double *y, double *x, double *p, 
+			      int *obsindex, int *stateindex, int *parindex, int *covindex,
+			      int ncovars, double *covars, double t) {
+  POP = rpois(N);
+}
+
+void _parus_gompertz_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 S = exp(-R*dt);
+  double eps = (SIGMA > 0.0) ? exp(rnorm(0,SIGMA)) : 1.0;
+  N = pow(K,(1-S))*pow(N,S)*eps;
+}
+
+// the deterministic skeleton
+void _parus_gompertz_skeleton (double *f, double *x, const double *p, 
+			       const int *stateindex, const int *parindex, const int *covindex,
+			       int covdim, const double *covar, double t) 
+{
+  double dt = 1.0;
+  double S = exp(-R*dt);
+  NPRIME = pow(K,(1-S))*pow(N,S);
+}
+
+// Ricker model with log-normal process noise
+void _parus_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 e = (SIGMA > 0.0) ? rnorm(0,SIGMA) : 0.0;
+  N = exp(log(N)+R*(1-N/K)+e);
+}
+
+void _parus_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) 
+{
+  NPRIME = exp(log(N)+R*(1-N/K));
+}