[Pomp-commits] r652 - in pkg: . demo inst inst/examples man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Apr 10 17:47:46 CEST 2012
Author: kingaa
Date: 2012-04-10 17:47:46 +0200 (Tue, 10 Apr 2012)
New Revision: 652
Added:
pkg/demo/
pkg/demo/00Index
pkg/demo/gompertz.R
pkg/demo/logistic.R
pkg/demo/rw2.R
pkg/demo/sir.R
Removed:
pkg/inst/examples/gompertz.R
pkg/inst/examples/logistic.R
pkg/inst/examples/rw2.R
pkg/inst/examples/sir.R
Modified:
pkg/inst/NEWS
pkg/inst/examples/gompertz.c
pkg/inst/examples/ou2.c
pkg/inst/examples/sir.c
pkg/man/trajectory-pomp.Rd
Log:
- add demos
Added: pkg/demo/00Index
===================================================================
--- pkg/demo/00Index (rev 0)
+++ pkg/demo/00Index 2012-04-10 15:47:46 UTC (rev 652)
@@ -0,0 +1,4 @@
+logistic The stochastic logistic (Verhulst-Pearl) model.
+gompertz The Gompertz model using both R code and native C code.
+sir A Seasonally-forced SIR model using both R and native C.
+rw2 Trivial 2-d random walk model.
Copied: pkg/demo/gompertz.R (from rev 651, pkg/inst/examples/gompertz.R)
===================================================================
--- pkg/demo/gompertz.R (rev 0)
+++ pkg/demo/gompertz.R 2012-04-10 15:47:46 UTC (rev 652)
@@ -0,0 +1,141 @@
+require(pomp)
+
+## First, code up the Gompertz example in R:
+
+pomp(
+ data=data.frame(time=1:100,Y=NA),
+ times="time",
+ t0=0,
+ rprocess=discrete.time.sim( # a discrete-time process (see ?plugins)
+ step.fun=function (x, t, params, delta.t, ...) { # this function takes one step t -> t+delta.t
+ ## unpack the parameters:
+ r <- params["r"]
+ K <- params["K"]
+ sigma <- params["sigma"]
+ ## the state at time t:
+ X <- x["X"]
+ ## generate a log-normal random variable:
+ eps <- exp(rnorm(n=1,mean=0,sd=sigma))
+ ## compute the state at time t+delta.t:
+ S <- exp(-r*delta.t)
+ xnew <- c(X=unname(K^(1-S)*X^S*eps))
+ return(xnew)
+ },
+ delta.t=1 # the size of the discrete time-step
+ ),
+ rmeasure=function (x, t, params, ...) {# the measurement model simulator
+ ## unpack the parameters:
+ tau <- params["tau"]
+ ## state at time t:
+ X <- x["X"]
+ ## generate a simulated observation:
+ y <- c(Y=unname(rlnorm(n=1,meanlog=log(X),sdlog=tau)))
+ return(y)
+ },
+ dmeasure=function (y, x, t, params, log, ...) { # measurement model density
+ ## unpack the parameters:
+ tau <- params["tau"]
+ ## state at time t:
+ X <- x["X"]
+ ## observation at time t:
+ Y <- y["Y"]
+ ## compute the likelihood of Y|X,tau
+ f <- dlnorm(x=Y,meanlog=log(X),sdlog=tau,log=log)
+ return(f)
+ },
+ parameter.inv.transform=function(params,...){
+ params <- log(params[c("X.0","r","K","tau","sigma")])
+ names(params) <- c("log.X.0","log.r","log.K","log.tau","log.sigma")
+ params
+ },
+ parameter.transform=function(params,...){
+ params <- exp(params[c("log.X.0","log.r","log.K","log.tau","log.sigma")])
+ names(params) <- c("X.0","r","K","tau","sigma")
+ params
+ }
+ ) -> gompertz
+
+## Now code up the Gompertz example using native routines results in much faster computations.
+## The C codes are included in the "examples" directory (file "gompertz.c")
+
+if (Sys.info()['sysname']=='Linux') { # only run this under linux
+
+ model <- "gompertz"
+ ## name of the file holding the native codes for this model:
+ modelfile <- system.file("examples",paste(model,".c",sep=""),package="pomp")
+ ## name of the shared-object library
+ solib <- paste(model,.Platform$dynlib.ext,sep="")
+ ## environment variables needed to locate the pomp header file:
+ cflags <- paste("PKG_CFLAGS=\"-I",system.file("include/",package="pomp"),"\"")
+
+ ## compile the shared-object library containing the model codes:
+ rv <- system2(
+ command=R.home("bin/R"),
+ args=c("CMD","SHLIB","-o",solib,modelfile),
+ env=cflags
+ )
+ ## compile the shared-object library containing the model codes:
+ if (rv!=0)
+ stop("cannot compile shared-object library ",sQuote(solib))
+
+ dyn.load(solib) ## load the shared-object library
+
+ po <- pomp(
+ data=data.frame(time=seq(0,100,by=1),Y=NA),
+ times="time",
+ t0=0,
+ rprocess=discrete.time.sim(
+ step.fun="gompertz_simulator",
+ PACKAGE="gompertz"
+ ),
+ rmeasure="gompertz_normal_rmeasure",
+ dmeasure="gompertz_normal_dmeasure",
+ skeleton.type="map",
+ skeleton="gompertz_skeleton",
+ PACKAGE="gompertz",
+ paramnames=c("r","K","sigma","tau"),
+ statenames=c("X"),
+ obsnames=c("Y"),
+ parameter.transform=function(params,...){
+ params <- log(params[c("X.0","r","K","tau","sigma")])
+ names(params) <- c("log.X.0","log.r","log.K","log.tau","log.sigma")
+ params
+ },
+ parameter.inv.transform=function(params,...){
+ params <- exp(params[c("log.X.0","log.r","log.K","log.tau","log.sigma")])
+ names(params) <- c("X.0","r","K","tau","sigma")
+ params
+ }
+ )
+
+ ## set the parameters
+ coef(po) <- c(K=1,r=0.1,sigma=0.1,tau=0.1,X.0=1)
+
+ p <- parmat(coef(po),nrep=4)
+ p["X.0",] <- c(0.5,0.9,1.1,1.5)
+
+ ## compute a trajectory of the deterministic skeleton
+ tic <- Sys.time()
+ X <- trajectory(po,params=p,as.data.frame=TRUE)
+ toc <- Sys.time()
+ print(toc-tic)
+ X <- reshape(X,dir="wide",v.names="X",timevar="traj",idvar="time")
+ matplot(X$time,X[-1],type='l',lty=1,bty='l',xlab="time",ylab="X",
+ main="Gompertz model\ndeterministic trajectories")
+
+ ## simulate from the model
+ tic <- Sys.time()
+ x <- simulate(po,params=p,as.data.frame=TRUE)
+ toc <- Sys.time()
+ print(toc-tic)
+
+ x <- reshape(x,dir="wide",v.names=c("Y","X"),timevar="sim",idvar="time")
+ op <- par(mfrow=c(2,1),mgp=c(2,1,0),mar=c(3,3,0,0),bty='l')
+ matplot(x$time,x[c("X.1","X.2","X.3")],lty=1,type='l',xlab="time",ylab="X",
+ main="Gompertz model\nstochastic simulations")
+ matplot(x$time,x[c("Y.1","Y.2","Y.3")],lty=1,type='l',xlab="time",ylab="Y")
+ par(op)
+
+ dyn.unload(solib)
+
+}
Copied: pkg/demo/logistic.R (from rev 651, pkg/inst/examples/logistic.R)
===================================================================
--- pkg/demo/logistic.R (rev 0)
+++ pkg/demo/logistic.R 2012-04-10 15:47:46 UTC (rev 652)
@@ -0,0 +1,60 @@
+require(pomp)
+
+## a stochastic version of the Verhulst-Pearl logistic model
+## this evolves in continuous time, but we approximate the
+## stochastic dynamics using an Euler approximation
+## (plugin 'euler.sim') with fixed step-size
+
+po <- pomp(
+ data=rbind(obs=rep(0,1000)),
+ times=seq(0.1,by=0.1,length=1000),
+ t0=0,
+ rprocess=euler.sim(
+ step.fun=function(x,t,params,delta.t,...){
+ with(
+ as.list(c(x,params)),
+ rnorm(
+ n=1,
+ mean=n+r*n*(1-n/K)*delta.t,
+ sd=sigma*n*sqrt(delta.t)
+ )
+ )
+ },
+ delta.t=0.01
+ ),
+ dprocess=onestep.dens(
+ dens.fun=function(x1,x2,t1,t2,params,log,...){
+ delta.t <- t2-t1
+ with(
+ as.list(c(x1,params)),
+ dnorm(
+ x=x2['n'],
+ mean=n+r*n*(1-n/K)*delta.t,
+ sd=sigma*n*sqrt(delta.t),
+ log=log
+ )
+ )
+ }
+ ),
+ measurement.model=obs~lnorm(meanlog=log(n),sdlog=log(1+tau)),
+ skeleton.type="vectorfield",
+ skeleton=function(x,t,params,...){
+ with(
+ as.list(c(x,params)),
+ r*n*(1-n/K)
+ )
+ }
+ )
+
+params <- c(n.0=10000,K=10000,r=0.9,sigma=0.4,tau=0.1)
+set.seed(73658676)
+po <- simulate(po,params=params)
+plot(po)
+
+params <- cbind(
+ c(n.0=100,K=10000,r=0.2,sigma=0.4,tau=0.1),
+ c(n.0=1000,K=11000,r=0.1,sigma=0.4,tau=0.1)
+ )
+x <- trajectory(po,params=params,as.data.frame=TRUE)
+x <- reshape(x,dir="wide",idvar="time",timevar="traj")
+matplot(x$time,x[-1],type='l',bty='l',lty=1,xlab="time",ylab="n")
Copied: pkg/demo/rw2.R (from rev 651, pkg/inst/examples/rw2.R)
===================================================================
--- pkg/demo/rw2.R (rev 0)
+++ pkg/demo/rw2.R 2012-04-10 15:47:46 UTC (rev 652)
@@ -0,0 +1,43 @@
+require(pomp)
+
+## a simple two-dimensional random walk
+## this makes use of the 'onestep.sim' plugin
+## which we can use since we can simulate the
+## increment of a random walk over any time
+
+rw2 <- pomp(
+ 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)
+ ),
+ times=1:100,
+ data=rbind(
+ y1=rep(0,100),
+ y2=rep(0,100)
+ ),
+ t0=0
+ )
+
+rw2 <- simulate(rw2,params=c(s1=3,s2=1,x1.0=0,x2.0=1,tau=10))
+plot(rw2)
Copied: pkg/demo/sir.R (from rev 651, pkg/inst/examples/sir.R)
===================================================================
--- pkg/demo/sir.R (rev 0)
+++ pkg/demo/sir.R 2012-04-10 15:47:46 UTC (rev 652)
@@ -0,0 +1,112 @@
+require(pomp)
+
+## Coding up the SIR example using native routines results in much faster computations.
+## The C codes "sir_euler_simulator", "sir_ODE", "binom_rmeasure", and "binom_dmeasure"
+## are included in the "examples" directory (file "sir.c")
+
+if (Sys.info()['sysname']=='Linux') { # only run this under linux
+
+ model <- "sir"
+ ## name of the file holding the model codes:
+ modelfile <- system.file("examples",paste(model,".c",sep=""),package="pomp")
+ ## name of the shared-object library:
+ solib <- paste(model,.Platform$dynlib.ext,sep="")
+ ## environment variables needed to locate the pomp header file:
+ cflags <- paste("PKG_CFLAGS=\"-I",system.file("include/",package="pomp"),"\"")
+
+ ## compile the shared-object library containing the model codes:
+ rv <- system2(
+ command=R.home("bin/R"),
+ args=c("CMD","SHLIB","-o",solib,modelfile),
+ env=cflags
+ )
+ if (rv!=0)
+ stop("cannot compile shared-object library ",sQuote(solib))
+
+ po <- pomp(
+ data=data.frame(
+ time=seq(from=1/52,to=4,by=1/52),
+ reports=NA
+ ),
+ times="time",
+ t0=0,
+ rprocess=euler.sim(
+ step.fun="sir_euler_simulator", # native routine for the simulation step
+ delta.t=1/52/20 # Euler stepsize: 1/20 wk
+ ),
+ skeleton.type="vectorfield",
+ skeleton="sir_ODE", # native routine for the skeleton
+ rmeasure="binomial_rmeasure", # binomial measurement model
+ dmeasure="binomial_dmeasure", # binomial measurement model
+ PACKAGE="sir", ## name of the shared-object library
+ ## the order of the observables assumed in the native routines:
+ obsnames="reports",
+ ## the order of the state variables assumed in the native routines:
+ statenames=c("S","I","R","cases","W"),
+ ## the order of the parameters assumed in the native routines:
+ paramnames=c(
+ "gamma","mu","iota","beta1","beta.sd","pop","rho",
+ "nbasis","degree","period"
+ ),
+ ## reset cases to zero after each new observation:
+ zeronames=c("cases"),
+ to.log.transform=c(
+ "gamma","mu","iota",
+ "beta1","beta2","beta3","beta.sd",
+ "S.0","I.0","R.0"
+ ),
+ to.logit.transform="rho",
+ comp.names=c("S","I","R"),
+ ic.names=c("S.0","I.0","R.0"),
+ parameter.transform=function (params, to.log.transform, to.logit.transform, ic.names, ...) {
+ params[log.transformed] <- exp(params[log.transformed])
+ params[logit.transformed] <- plogis(params[logit.transformed])
+ params[ic.names] <- params[ic.names]/sum(params[ic.names])
+ params
+ },
+ parameter.inv.transform=function (params, log.transformed, logit.transformed, ic.names, ...) {
+ params[ic.names] <- params[ic.names]/sum(params[ic.names])
+ params[log.transformed] <- log(params[log.transformed])
+ params[logit.transformed] <- qlogis(params[logit.transformed])
+ params
+ },
+ initializer=function(params, t0, comp.names, ic.names, ...) {
+ x0 <- numeric(5)
+ names(x0) <- c("S","I","R","cases","W")
+ fracs <- params[ic.names]
+ x0[comp.names] <- round(params['pop']*fracs/sum(fracs))
+ x0
+ }
+ )
+
+ coef(po) <- c(
+ gamma=26,mu=0.02,iota=0.01,
+ beta1=1200,beta2=1800,beta3=600,
+ beta.sd=1e-3,
+ pop=2.1e6,
+ rho=0.6,
+ S.0=26/1200,I.0=0.001,R.0=1-26/1200,
+ nbasis=3,degree=3,period=1
+ )
+
+ dyn.load(solib) ## load the shared-object library
+
+ ## compute a trajectory of the deterministic skeleton
+ tic <- Sys.time()
+ X <- trajectory(po,hmax=1/52,as.data.frame=TRUE)
+ toc <- Sys.time()
+ print(toc-tic)
+
+ plot(cases~time,data=X,type='l')
+
+ ## simulate from the model
+ tic <- Sys.time()
+ x <- simulate(po,nsim=3,as.data.frame=TRUE)
+ toc <- Sys.time()
+ print(toc-tic)
+
+ plot(cases~time,data=x,col=as.factor(x$sim),pch=16)
+
+ dyn.unload(solib)
+
+}
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2012-04-10 14:16:20 UTC (rev 651)
+++ pkg/inst/NEWS 2012-04-10 15:47:46 UTC (rev 652)
@@ -3,6 +3,8 @@
o In 'trajectory', all vectorfield and map evaluation is now done in C for speed.
For continuous-time dynamical systems, this gives an approximately 4-fold speedup.
+ o Demos have been created to show some examples.
+
0.41-2
o Fix bug in 'bbs' example.
Deleted: pkg/inst/examples/gompertz.R
===================================================================
--- pkg/inst/examples/gompertz.R 2012-04-10 14:16:20 UTC (rev 651)
+++ pkg/inst/examples/gompertz.R 2012-04-10 15:47:46 UTC (rev 652)
@@ -1,125 +0,0 @@
-require(pomp)
-
-## First code up the Gompertz example in R:
-
-pomp(
- data=data.frame(time=1:100,Y=NA),
- times="time",
- t0=0,
- rprocess=discrete.time.sim( # a discrete-time process (see ?plugins)
- step.fun=function (x, t, params, delta.t, ...) { # this function takes one step t -> t+delta.t
- ## unpack the parameters:
- r <- params["r"]
- K <- params["K"]
- sigma <- params["sigma"]
- ## the state at time t:
- X <- x["X"]
- ## generate a log-normal random variable:
- eps <- exp(rnorm(n=1,mean=0,sd=sigma))
- ## compute the state at time t+delta.t:
- S <- exp(-r*delta.t)
- xnew <- c(X=unname(K^(1-S)*X^S*eps))
- return(xnew)
- },
- delta.t=1 # the size of the discrete time-step
- ),
- rmeasure=function (x, t, params, ...) {# the measurement model simulator
- ## unpack the parameters:
- tau <- params["tau"]
- ## state at time t:
- X <- x["X"]
- ## generate a simulated observation:
- y <- c(Y=unname(rlnorm(n=1,meanlog=log(X),sdlog=tau)))
- return(y)
- },
- dmeasure=function (y, x, t, params, log, ...) { # measurement model density
- ## unpack the parameters:
- tau <- params["log.tau"]
- ## state at time t:
- X <- x["X"]
- ## observation at time t:
- Y <- y["Y"]
- ## compute the likelihood of Y|X,tau
- f <- dlnorm(x=Y,meanlog=log(X),sdlog=tau,log=log)
- return(f)
- },
- parameter.inv.transform=function(params,...){
- params <- log(params[c("X.0","r","K","tau","sigma")])
- names(params) <- c("log.X.0","log.r","log.K","log.tau","log.sigma")
- params
- },
- parameter.transform=function(params,...){
- params <- exp(params[c("log.X.0","log.r","log.K","log.tau","log.sigma")])
- names(params) <- c("X.0","r","K","tau","sigma")
- params
- }
- ) -> gompertz
-
-## Now code up the Gompertz example using native routines results in much faster computations.
-## The C codes are included in the "examples" directory (file "gompertz.c")
-
-po <- pomp(
- data=data.frame(time=seq(0,100,by=1),Y=NA),
- times="time",
- t0=0,
- rprocess=discrete.time.sim(
- step.fun="gompertz_simulator",
- PACKAGE="gompertz"
- ),
- rmeasure="gompertz_normal_rmeasure",
- dmeasure="gompertz_normal_dmeasure",
- skeleton.type="map",
- skeleton="gompertz_skeleton",
- PACKAGE="gompertz",
- paramnames=c("r","K","sigma","tau"),
- statenames=c("X"),
- obsnames=c("Y"),
- parameter.transform=function(params,...){
- params <- log(params[c("X.0","r","K","tau","sigma")])
- names(params) <- c("log.X.0","log.r","log.K","log.tau","log.sigma")
- params
- },
- parameter.inv.transform=function(params,...){
- params <- exp(params[c("log.X.0","log.r","log.K","log.tau","log.sigma")])
- names(params) <- c("X.0","r","K","tau","sigma")
- params
- }
- )
-
-## set the parameters
-coef(po) <- c(K=1,r=0.1,sigma=0.1,tau=0.1,X.0=1)
-
-if (Sys.info()['sysname']=='Linux') { # only run this under linux
-
- model <- "gompertz"
- pkg <- "pomp"
- modelfile <- paste(model,".c",sep="")
- solib <- paste(model,.Platform$dynlib.ext,sep="")
-
- ## compile the model into a shared-object library
- if (!file.copy(from=system.file(paste("examples/",modelfile,sep=""),package=pkg),to=getwd()))
- stop("cannot copy source code ",modelfile," to ",getwd())
- if (!file.copy(from=system.file("include/pomp.h",package=pkg),to=getwd()))
- stop("cannot copy header file ",modelfile," to ",getwd())
- cmd <- paste(R.home("bin/R"),"CMD SHLIB -o",solib,modelfile)
- rv <- system(cmd)
- if (rv!=0)
- stop("cannot compile shared-object library ",solib)
-
- dyn.load(solib) ## load the shared-object library
-
- ## compute a trajectory of the deterministic skeleton
- tic <- Sys.time()
- X <- trajectory(po)
- toc <- Sys.time()
- print(toc-tic)
-
- ## simulate from the model
- tic <- Sys.time()
- x <- simulate(po,nsim=3)
- toc <- Sys.time()
- print(toc-tic)
-
- dyn.unload(solib)
-
-}
Modified: pkg/inst/examples/gompertz.c
===================================================================
--- pkg/inst/examples/gompertz.c 2012-04-10 14:16:20 UTC (rev 651)
+++ pkg/inst/examples/gompertz.c 2012-04-10 15:47:46 UTC (rev 652)
@@ -1,27 +1,30 @@
// dear emacs, please treat this as -*- C++ -*-
+// Gompertz example as described in the "Introduction to pomp" vignette.
+// for a demonstration of how to compile, link, and run this example,
+// do 'demo("gompertz",package="pomp")' at the R prompt.
+
#include <Rmath.h>
-#include "pomp.h" // in R, do 'system.file("include/pomp.h",package="pomp")' to find this header file
+#include "pomp.h" // in R, do 'system.file("include","pomp.h",package="pomp")' to find this header file
// define some macros to make the code easier to read
-#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 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 Y (y[obsindex[0]]) // observed population size
+#define X (x[stateindex[0]]) // actual population size
+#define XPRIME (f[stateindex[0]]) // new population size (for skeleton function only)
-#define Y (y[obsindex[0]]) // observed population size
-#define X (x[stateindex[0]]) // actual population size
-#define XPRIME (f[stateindex[0]]) // new population size (for skeleton function only)
-
-// normal measurement model density
+// log-normal measurement model density
void gompertz_normal_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(Y,log(X),TAU,give_log);
}
-// normal measurement model simulator
+// log-normal measurement model simulator
void gompertz_normal_rmeasure (double *y, double *x, double *p,
int *obsindex, int *stateindex, int *parindex, int *covindex,
int ncovars, double *covars, double t) {
Deleted: pkg/inst/examples/logistic.R
===================================================================
--- pkg/inst/examples/logistic.R 2012-04-10 14:16:20 UTC (rev 651)
+++ pkg/inst/examples/logistic.R 2012-04-10 15:47:46 UTC (rev 652)
@@ -1,54 +0,0 @@
-require(pomp)
-
-## a stochastic version of the Verhulst-Pearl logistic model
-## this evolves in continuous time, but we approximate the
-## stochastic dynamics using an Euler approximation
-## (plugin 'euler.sim') with fixed step-size
-
-po <- pomp(
- data=rbind(obs=rep(0,1000)),
- times=seq(0.1,by=0.1,length=1000),
- t0=0,
- rprocess=euler.sim(
- step.fun=function(x,t,params,delta.t,...){
- with(
- as.list(c(x,params)),
- rnorm(
- n=1,
- mean=n+r*n*(1-n/K)*delta.t,
- sd=sigma*n*sqrt(delta.t)
- )
- )
- },
- delta.t=0.01
- ),
- dprocess=onestep.dens(
- dens.fun=function(x1,x2,t1,t2,params,log,...){
- delta.t <- t2-t1
- with(
- as.list(c(x1,params)),
- dnorm(
- x=x2['n'],
- mean=n+r*n*(1-n/K)*delta.t,
- sd=sigma*n*sqrt(delta.t),
- log=log
- )
- )
- }
- ),
- measurement.model=obs~lnorm(meanlog=log(n),sdlog=log(1+tau)),
- skeleton.type="vectorfield",
- skeleton=function(x,t,params,...){
- with(
- as.list(c(x,params)),
- r*n*(1-n/K)
- )
- }
- )
-
-params <- c(n.0=10000,K=10000,r=0.9,sigma=0.4,tau=0.1)
-set.seed(73658676)
-po <- simulate(po,params=params)
-
-params <- cbind(c(n.0=100,K=10000,r=0.2,sigma=0.4,tau=0.1),c(n.0=1000,K=11000,r=0.1,sigma=0.4,tau=0.1))
-x <- trajectory(po,params=params)
Modified: pkg/inst/examples/ou2.c
===================================================================
--- pkg/inst/examples/ou2.c 2012-04-10 14:16:20 UTC (rev 651)
+++ pkg/inst/examples/ou2.c 2012-04-10 15:47:46 UTC (rev 652)
@@ -1,5 +1,7 @@
// -*- C++ -*-
+// 2-D Ornstein-Uhlenbeck examples as discussed in the "Advanced Topics in pomp" vignette.
+
#include <R.h>
#include <Rmath.h>
#include <Rdefines.h>
Deleted: pkg/inst/examples/rw2.R
===================================================================
--- pkg/inst/examples/rw2.R 2012-04-10 14:16:20 UTC (rev 651)
+++ pkg/inst/examples/rw2.R 2012-04-10 15:47:46 UTC (rev 652)
@@ -1,40 +0,0 @@
-require(pomp)
-
-## a simple two-dimensional random walk
-## this makes use of the 'onestep.sim' plugin
-## which we can use since we can simulate the
-## increment of a random walk over any time
-
-rw2 <- pomp(
- 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)
- ),
- times=1:100,
- data=rbind(
- y1=rep(0,100),
- y2=rep(0,100)
- ),
- t0=0
- )
Deleted: pkg/inst/examples/sir.R
===================================================================
--- pkg/inst/examples/sir.R 2012-04-10 14:16:20 UTC (rev 651)
+++ pkg/inst/examples/sir.R 2012-04-10 15:47:46 UTC (rev 652)
@@ -1,108 +0,0 @@
-require(pomp)
-
-## Coding up the SIR example using native routines results in much faster computations.
-## The C codes "sir_euler_simulator", "sir_ODE", "binom_rmeasure", and "binom_dmeasure"
-## are included in the "examples" directory (file "sir.c")
-
-if (Sys.info()['sysname']=='Linux') { # only run this under linux
-
- model <- "sir"
- ## name of the file holding the model codes:
- modelfile <- system.file(file.path("examples",paste(model,".c",sep="")),package="pomp")
- ## name of the shared-object library:
- solib <- paste(model,.Platform$dynlib.ext,sep="")
- ## environment variables needed to locate the pomp header file:
- cflags <- paste("PKG_CFLAGS=\"-I",system.file("include/",package="pomp"),"\"")
-
- ## compile the shared-object library containing the model codes:
- rv <- system2(
- command=R.home("bin/R"),
- args=c("CMD","SHLIB","-o",solib,modelfile),
- env=cflags
- )
- if (rv!=0)
- stop("cannot compile shared-object library ",sQuote(solib))
-
- po <- pomp(
- data=data.frame(
- time=seq(from=1/52,to=4,by=1/52),
- reports=NA
- ),
- times="time",
- t0=0,
- rprocess=euler.sim(
- step.fun="sir_euler_simulator", # native routine for the simulation step
- delta.t=1/52/20 # Euler stepsize: 1/20 wk
- ),
- skeleton.type="vectorfield",
- skeleton="sir_ODE", # native routine for the skeleton
- rmeasure="binomial_rmeasure", # binomial measurement model
- dmeasure="binomial_dmeasure", # binomial measurement model
- PACKAGE="sir", ## name of the shared-object library
- ## the order of the observables assumed in the native routines:
- obsnames="reports",
- ## the order of the state variables assumed in the native routines:
- statenames=c("S","I","R","cases","W"),
- ## the order of the parameters assumed in the native routines:
- paramnames=c(
- "gamma","mu","iota","beta1","beta.sd","pop","rho",
- "nbasis","degree","period"
- ),
- ## reset cases to zero after each new observation:
- zeronames=c("cases"),
- to.log.transform=c(
- "gamma","mu","iota",
- "beta1","beta2","beta3","beta.sd",
- "S.0","I.0","R.0"
- ),
- to.logit.transform="rho",
- comp.names=c("S","I","R"),
- ic.names=c("S.0","I.0","R.0"),
- parameter.transform=function (params, to.log.transform, to.logit.transform, ic.names, ...) {
- params[log.transformed] <- exp(params[log.transformed])
- params[logit.transformed] <- plogis(params[logit.transformed])
- params[ic.names] <- params[ic.names]/sum(params[ic.names])
- params
- },
- parameter.inv.transform=function (params, log.transformed, logit.transformed, ic.names, ...) {
- params[ic.names] <- params[ic.names]/sum(params[ic.names])
- params[log.transformed] <- log(params[log.transformed])
- params[logit.transformed] <- qlogis(params[logit.transformed])
- params
- },
- initializer=function(params, t0, comp.names, ic.names, ...) {
- x0 <- numeric(5)
- names(x0) <- c("S","I","R","cases","W")
- fracs <- params[ic.names]
- x0[comp.names] <- round(params['pop']*fracs/sum(fracs))
- x0
- }
- )
-
- coef(po) <- c(
- gamma=26,mu=0.02,iota=0.01,
- beta1=1200,beta2=1800,beta3=600,
- beta.sd=1e-3,
- pop=2.1e6,
- rho=0.6,
- S.0=26/1200,I.0=0.001,R.0=1-26/1200,
- nbasis=3,degree=3,period=1
- )
-
- dyn.load(solib) ## load the shared-object library
-
- ## compute a trajectory of the deterministic skeleton
- tic <- Sys.time()
- X <- trajectory(po,hmax=1/52)
- toc <- Sys.time()
- print(toc-tic)
-
- ## simulate from the model
- tic <- Sys.time()
- x <- simulate(po,nsim=3)
- toc <- Sys.time()
- print(toc-tic)
-
- dyn.unload(solib)
-
-}
Modified: pkg/inst/examples/sir.c
===================================================================
--- pkg/inst/examples/sir.c 2012-04-10 14:16:20 UTC (rev 651)
+++ pkg/inst/examples/sir.c 2012-04-10 15:47:46 UTC (rev 652)
@@ -5,6 +5,10 @@
#include <R_ext/Rdynload.h>
#include <pomp.h>
+// SIR example as described in the "Advanced Topics in pomp" vignette.
+// for a demonstration of how to compile, link, and run this example,
+// do 'demo("sir",package="pomp")' at the R prompt.
+
// some macros to make the code easier to read.
// note how variables and parameters use the indices.
// the integer vectors parindex, stateindex, obsindex
Modified: pkg/man/trajectory-pomp.Rd
===================================================================
--- pkg/man/trajectory-pomp.Rd 2012-04-10 14:16:20 UTC (rev 651)
+++ pkg/man/trajectory-pomp.Rd 2012-04-10 15:47:46 UTC (rev 652)
@@ -4,7 +4,7 @@
\alias{trajectory,pomp-method}
\alias{trajectory-pomp}
-\title{Compute trajectories of the determinstic skeleton.}
+\title{Compute trajectories of the deterministic skeleton.}
\description{
The method \code{trajectory} computes a trajectory of the deterministic skeleton of a Markov process.
In the case of a discrete-time system, the deterministic skeleton is a map and a trajectory is obtained by iterating the map.
More information about the pomp-commits
mailing list