[Pomp-commits] r635 - in pkg: . data inst inst/data-R inst/examples inst/include man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 29 17:30:55 CEST 2012
Author: kingaa
Date: 2012-03-29 17:30:55 +0200 (Thu, 29 Mar 2012)
New Revision: 635
Added:
pkg/data/bbs.rda
Modified:
pkg/DESCRIPTION
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/inst/NEWS
pkg/inst/data-R/sir.R
pkg/inst/examples/gompertz.R
pkg/inst/examples/gompertz.c
pkg/inst/examples/sir.R
pkg/inst/examples/sir.c
pkg/inst/include/pomp.h
pkg/man/sir.Rd
pkg/src/eulermultinom.c
pkg/src/lookup_table.c
pkg/src/pomp.h
pkg/src/sir.c
Log:
- rejigger the 'inst/examples/*' to use the 0.41-1 parameter-transformation functionality
- put some of the utility functions into 'pomp.h' as inlines
- add the 'bbs' example
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-03-27 21:17:56 UTC (rev 634)
+++ pkg/DESCRIPTION 2012-03-29 15:30:55 UTC (rev 635)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.41-1
-Date: 2012-03-28
+Date: 2012-03-30
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>
Added: pkg/data/bbs.rda
===================================================================
(Binary files differ)
Property changes on: pkg/data/bbs.rda
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2012-03-27 21:17:56 UTC (rev 634)
+++ pkg/inst/NEWS 2012-03-29 15:30:55 UTC (rev 635)
@@ -7,10 +7,13 @@
Note that this new functionality makes it unnecessary to "un-transform" model parameters within the user-specified 'rprocess', 'dprocess', 'rmeasure', 'dmeasure', 'skeleton', and 'initializer' codes.
o The Bayesian sequential Monte Carlo command 'bsmc' now returns not a list but an object of class 'bsmcd.pomp'.
- A 'plot' method for objects of this class now exists.
+ An experimental 'plot' method for objects of this class now exists.
Also, the parameter posterior means are now stored in the 'params' slot of the 'bsmcd.pomp' object:
access them with the 'coef' command as usual.
+ o A new example, using data from an influenza outbreak in a British boarding school and an SIR model, has been included.
+ Do 'data(bbs)' to load it.
+
0.40-9
o Setting the new argument 'as.data.frame' to 'TRUE' in 'simulate' and 'trajectory' causes the results to be returned as a data-frame.
@@ -54,7 +57,7 @@
0.40-3
o The 'bsmc' method has been improved, due to the contributions of Pierre Jacob.
- Specifically, 'bsmc' no longer reports a log-likelihood (which it never really computed anyway) but a log-evidence.
+ Specifically, 'bsmc' no longer reports a log-likelihood (which it never actually computed anyway) but a log-evidence.
The latter computation was supplied by Pierre Jacob.
o A new helper function, 'exp2geom_rate_correction' has been added to the 'pomp.h' file.
@@ -62,11 +65,13 @@
This is useful in approximating a continuous-time death process by a discrete time (Euler) process.
In particular, in such a case, T is the waiting time to death in the former and N*dt is the waiting time to death in the latter.
An Euler binomial or multinomial process with rate r=exp2geom_rate_correction(R,dt) will have the same mean waiting time as the exponential process with rate R.
+ Thanks to Sebastien Ballesteros for suggesting this.
o A new helper function, 'rgammawn' has been added, with both R and C interfaces.
This function draws an increment of a Gamma white-noise process with a given intensity.
Specifically, dw=rgammwn(sigma,dt) is Gamma distributed with mean dt and variance sigma^2*dt.
In this case, mu*dw/dt is suitable for use as a random rate in an Euler-multinomial process.
+ Thanks to Sebastien Ballesteros for suggesting this.
0.40-2
o A bug to do with computation of the number of steps needed in discrete-time simulation and trajectory computations has been fixed.
Modified: pkg/inst/data-R/sir.R
===================================================================
--- pkg/inst/data-R/sir.R 2012-03-27 21:17:56 UTC (rev 634)
+++ pkg/inst/data-R/sir.R 2012-03-29 15:30:55 UTC (rev 635)
@@ -27,19 +27,19 @@
),
zeronames=c("cases"),
comp.names=c("S","I","R"),
+ ic.names=c("S.0","I.0","R.0"),
parameter.transform="_sir_par_trans",
parameter.inv.transform="_sir_par_untrans",
- initializer=function(params, t0, comp.names, ...) {
- ic.names <- paste(comp.names,".0",sep="")
+ initializer=function(params, t0, comp.names, ic.names, ...) {
snames <- c("S","I","R","cases","W")
fracs <- params[ic.names]
x0 <- numeric(length(snames))
names(x0) <- snames
x0[comp.names] <- round(params['pop']*fracs/sum(fracs))
- x0["cases"] <- 0
x0
}
)
+
coef(po) <- c(
gamma=26,mu=0.02,iota=0.01,
nbasis=3,degree=3,period=1,
@@ -54,13 +54,11 @@
save(euler.sir,file="euler.sir.rda",compress="xz")
+time(po) <- seq(from=0,to=10,by=1/52)
+
po <- pomp(
- data=data.frame(
- time=seq(from=0,to=10,by=1/52),
- reports=NA
- ),
- times="time",
- t0=0,
+ po,
+ statenames=c("S","I","R","N","cases"),
rprocess=gillespie.sim(
rate.fun="_sir_rates",
PACKAGE="pomp",
@@ -81,25 +79,13 @@
rdeath=c(0,0,1,0,0)
)
),
- obsnames="reports",
- statenames=c("S","I","R","N","cases"),
- paramnames=c(
- "gamma","mu","iota",
- "log.beta1","beta.sd","pop","rho",
- "nbasis","degree","period",
- "S.0","I.0","R.0"
- ),
- zeronames=c("cases"),
measurement.model=reports~binom(size=cases,prob=rho),
- parameter.transform="_sir_par_trans",
- parameter.inv.transform="_sir_par_untrans",
- initializer=function(params, t0, ...){
- comp.names <- c("S","I","R")
- icnames <- paste(comp.names,"0",sep=".")
- snames <- c("S","I","R","N","cases")
- fracs <- params[icnames]
- x0 <- numeric(length(snames))
- names(x0) <- snames
+ comp.names=c("S","I","R"),
+ ic.names=c("S.0","I.0","R.0"),
+ initializer=function(params, t0, comp.names, ic.names, ...) {
+ x0 <- numeric(5)
+ names(x0) <- c("S","I","R","N","cases")
+ fracs <- params[ic.names]
x0["N"] <- params["pop"]
x0[comp.names] <- round(params["pop"]*fracs/sum(fracs))
x0
@@ -107,24 +93,15 @@
)
coef(po) <- c(
- gamma=24,
- mu=1/70,
- iota=0.1,
- log.beta1=log(330),
- log.beta2=log(410),
- log.beta3=log(490),
- rho=0.1,
- S.0=0.05,
- I.0=1e-4,
- R.0=0.95,
- pop=1000000,
- nbasis=3,
- degree=3,
- period=1,
- beta.sd=0
- )
+ gamma=24,mu=1/70,iota=0.1,
+ log.beta1=log(330),log.beta2=log(410),log.beta3=log(490),
+ rho=0.1,
+ S.0=0.05,I.0=1e-4,R.0=0.95,
+ pop=1000000,
+ nbasis=3,degree=3,period=1,
+ beta.sd=0
+ )
-
simulate(
po,
nsim=1,
@@ -132,3 +109,78 @@
) -> gillespie.sir
save(gillespie.sir,file="gillespie.sir.rda",compress="xz")
+
+tc <- textConnection("
+day;flu
+1;3
+2;8
+3;28
+4;76
+5;222
+6;293
+7;257
+8;237
+9;192
+10;126
+11;70
+12;28
+13;12
+14;5
+")
+
+flu <- read.csv2(file=tc)
+close(tc)
+
+po <- pomp(
+ data=flu,
+ times="day",
+ t0=0,
+ rprocess=euler.sim(
+ step.fun="_sir_euler_simulator",
+ delta.t=1/12,
+ PACKAGE="pomp"
+ ),
+ skeleton.type="vectorfield",
+ skeleton="_sir_ODE",
+ measurement.model=flu~norm(mean=rho*cases,sd=1+sigma*cases),
+ PACKAGE="pomp",
+ obsnames = c("flu"),
+ statenames=c("S","I","R","cases","W"),
+ paramnames=c(
+ "gamma","mu","iota",
+ "beta1","beta.sd","pop","rho",
+ "nbasis","degree","period"
+ ),
+ zeronames=c("cases"),
+ comp.names=c("S","I","R"),
+ ic.names=c("S.0","I.0","R.0"),
+ log.transformed=c(
+ "gamma","mu","iota","sigma",
+ "beta1","beta.sd",
+ "S.0","I.0","R.0"
+ ),
+ logit.transformed="rho",
+ parameter.transform="_sir_par_trans",
+ parameter.inv.transform="_sir_par_untrans",
+ initializer=function(params, t0, comp.names, ic.names, ...) {
+ snames <- c("S","I","R","cases","W")
+ fracs <- params[ic.names]
+ x0 <- numeric(length(snames))
+ names(x0) <- snames
+ x0[comp.names] <- round(params['pop']*fracs/sum(fracs))
+ x0
+ }
+ )
+
+coef(po) <- c(
+ gamma=1/3,mu=0.0,iota=0.0,
+ nbasis=1,degree=0,period=1,
+ beta1=1.4,
+ beta.sd=0,
+ pop=1400,
+ rho=0.9,sigma=3.6,
+ S.0=0.999,I.0=0.001,R.0=0
+ )
+
+bbs <- po
+save(bbs,file="bbs.rda",compress="xz")
Modified: pkg/inst/examples/gompertz.R
===================================================================
--- pkg/inst/examples/gompertz.R 2012-03-27 21:17:56 UTC (rev 634)
+++ pkg/inst/examples/gompertz.R 2012-03-29 15:30:55 UTC (rev 635)
@@ -8,10 +8,10 @@
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 and untransform the parameters:
- r <- exp(params["log.r"])
- K <- exp(params["log.K"])
- sigma <- exp(params["log.sigma"])
+ ## 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:
@@ -24,8 +24,8 @@
delta.t=1 # the size of the discrete time-step
),
rmeasure=function (x, t, params, ...) {# the measurement model simulator
- ## unpack and untransform the parameters:
- tau <- exp(params["log.tau"])
+ ## unpack the parameters:
+ tau <- params["tau"]
## state at time t:
X <- x["X"]
## generate a simulated observation:
@@ -33,8 +33,8 @@
return(y)
},
dmeasure=function (y, x, t, params, log, ...) { # measurement model density
- ## unpack and untransform the parameters:
- tau <- exp(params["log.tau"])
+ ## unpack the parameters:
+ tau <- params["log.tau"]
## state at time t:
X <- x["X"]
## observation at time t:
@@ -43,13 +43,13 @@
f <- dlnorm(x=Y,meanlog=log(X),sdlog=tau,log=log)
return(f)
},
- parameter.transform=function(params,...){
- params <- c(params["X.0"],log(params[c("r","K","tau","sigma")]))
- names(params) <- c("X.0","log.r","log.K","log.tau","log.sigma")
+ 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.inv.transform=function(params,...){
- params <- c(params["X.0"],exp(params[c("log.r","log.K","log.tau","log.sigma")]))
+ 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
}
@@ -71,23 +71,23 @@
skeleton.type="map",
skeleton="gompertz_skeleton",
PACKAGE="gompertz",
- paramnames=c("log.r","log.K","log.sigma","log.tau"),
+ paramnames=c("r","K","sigma","tau"),
statenames=c("X"),
obsnames=c("Y"),
parameter.transform=function(params,...){
- params <- c(params["X.0"],log(params[c("r","K","tau","sigma")]))
- names(params) <- c("X.0","log.r","log.K","log.tau","log.sigma")
+ 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 <- c(params["X.0"],exp(params[c("log.r","log.K","log.tau","log.sigma")]))
+ 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,transform=TRUE) <- c(K=1,r=0.1,sigma=0.1,tau=0.1,X.0=1)
+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
Modified: pkg/inst/examples/gompertz.c
===================================================================
--- pkg/inst/examples/gompertz.c 2012-03-27 21:17:56 UTC (rev 634)
+++ pkg/inst/examples/gompertz.c 2012-03-29 15:30:55 UTC (rev 635)
@@ -5,10 +5,10 @@
#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 LOG_R (p[parindex[0]]) // growth rate
-#define LOG_K (p[parindex[1]]) // carrying capacity
-#define LOG_SIGMA (p[parindex[2]]) // process noise level
-#define LOG_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
@@ -18,14 +18,14 @@
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),exp(LOG_TAU),give_log);
+ *lik = dlnorm(Y,log(X),TAU,give_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) {
- Y = rlnorm(log(X),exp(LOG_TAU));
+ Y = rlnorm(log(X),TAU);
}
// stochastic Gompertz model with log-normal process noise
@@ -34,10 +34,10 @@
int covdim, const double *covar,
double t, double dt)
{
- double S = exp(-exp(LOG_R)*dt);
- double sigma = exp(LOG_SIGMA);
+ double S = exp(-R*dt);
+ double sigma = SIGMA;
double logeps = (sigma > 0.0) ? rnorm(0,sigma) : 0.0;
- X = exp((1-S)*LOG_K)*pow(X,S)*exp(logeps); // note X is over-written by this line
+ X = pow(K,(1-S))*pow(X,S)*exp(logeps); // note X is over-written by this line
}
// the deterministic skeleton
@@ -46,6 +46,6 @@
int covdim, const double *covar, double t)
{
double dt = 1.0;
- double S = exp(-exp(LOG_R)*dt);
- XPRIME = exp((1-S)*LOG_K)*pow(X,S); // X is not over-written in the skeleton function
+ double S = exp(-R*dt);
+ XPRIME = pow(K,1-S)*pow(X,S); // X is not over-written in the skeleton function
}
Modified: pkg/inst/examples/sir.R
===================================================================
--- pkg/inst/examples/sir.R 2012-03-27 21:17:56 UTC (rev 634)
+++ pkg/inst/examples/sir.R 2012-03-29 15:30:55 UTC (rev 635)
@@ -7,17 +7,21 @@
if (Sys.info()['sysname']=='Linux') { # only run this under linux
model <- "sir"
- pkg <- "pomp"
- modelfile <- paste(model,".c",sep="")
+ ## 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 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())
- cmd <- paste(R.home("bin/R"),"CMD SHLIB -o",solib,modelfile)
- rv <- system(cmd)
+ ## 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 ",solib)
+ stop("cannot compile shared-object library ",sQuote(solib))
po <- pomp(
data=data.frame(
@@ -28,7 +32,7 @@
t0=0,
rprocess=euler.sim(
step.fun="sir_euler_simulator", # native routine for the simulation step
- delta.t=1/52/20
+ delta.t=1/52/20 # Euler stepsize: 1/20 wk
),
skeleton.type="vectorfield",
skeleton="sir_ODE", # native routine for the skeleton
@@ -49,39 +53,41 @@
to.log.transform=c(
"gamma","mu","iota",
"beta1","beta2","beta3","beta.sd",
- "rho",
"S.0","I.0","R.0"
),
- parameter.transform=function (params, to.log.transform, ...) {
- params[to.log.transform] <- log(params[to.log.transform])
+ 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, to.log.transform, ...) {
- params[to.log.transform] <- exp(params[to.log.transform])
+ 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 <- c("S","I","R")
- ic.names <- c("S.0","I.0","R.0")
- snames <- c("S","I","R","cases","W")
- fracs <- exp(params[ic.names])
- x0 <- numeric(length(snames))
- names(x0) <- snames
+ 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["cases"] <- 0
x0
}
)
- coef(po,transform=TRUE) <- 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-0.001-26/1200,
- nbasis=3,degree=3,period=1
- )
+ 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
Modified: pkg/inst/examples/sir.c
===================================================================
--- pkg/inst/examples/sir.c 2012-03-27 21:17:56 UTC (rev 634)
+++ pkg/inst/examples/sir.c 2012-03-29 15:30:55 UTC (rev 635)
@@ -3,57 +3,54 @@
#include <R.h>
#include <Rmath.h>
#include <R_ext/Rdynload.h>
+#include <pomp.h>
-static double expit (double x) {
- return 1.0/(1.0 + exp(-x));
-}
-
-static double logit (double x) {
- return log(x/(1-x));
-}
-
// some macros to make the code easier to read.
// note how variables and parameters use the indices.
// the integer vectors parindex, stateindex, obsindex
// are computed by pomp by matching the names of input vectors
// against parnames, statenames, and obsnames, respectively.
-#define LOGGAMMA (p[parindex[0]]) // recovery rate
-#define LOGMU (p[parindex[1]]) // baseline birth and death rate
-#define LOGIOTA (p[parindex[2]]) // import rate
-#define LOGBETA (p[parindex[3]]) // transmission rate
-#define LOGBETA_SD (p[parindex[4]]) // environmental stochasticity SD in transmission rate
-#define POPSIZE (p[parindex[5]]) // population size
-#define LOGRHO (p[parindex[6]]) // reporting probability
-#define NBASIS (p[parindex[7]]) // number of periodic B-spline basis functions
-#define DEGREE (p[parindex[8]]) // degree of periodic B-spline basis functions
-#define PERIOD (p[parindex[9]]) // period of B-spline basis functions
+#define GAMMA (p[parindex[0]]) // recovery rate
+#define MU (p[parindex[1]]) // baseline birth and death rate
+#define IOTA (p[parindex[2]]) // import rate
+#define BETA (&p[parindex[3]]) // transmission rate
+#define BETA_SD (p[parindex[4]]) // environmental stochasticity SD in transmission rate
+#define POPSIZE (p[parindex[5]]) // population size
+#define RHO (p[parindex[6]]) // reporting probability
+#define NBASIS (p[parindex[7]]) // number of periodic B-spline basis functions
+#define DEGREE (p[parindex[8]]) // degree of periodic B-spline basis functions
+#define PERIOD (p[parindex[9]]) // period of B-spline basis functions
-#define SUSC (x[stateindex[0]]) // number of susceptibles
-#define INFD (x[stateindex[1]]) // number of infectives
-#define RCVD (x[stateindex[2]]) // number of recovereds
-#define CASE (x[stateindex[3]]) // number of cases (accumulated per reporting period)
-#define W (x[stateindex[4]]) // integrated white noise
+#define SUSC (x[stateindex[0]]) // number of susceptibles
+#define INFD (x[stateindex[1]]) // number of infectives
+#define RCVD (x[stateindex[2]]) // number of recovereds
+#define CASE (x[stateindex[3]]) // number of cases (accumulated per reporting period)
+#define W (x[stateindex[4]]) // integrated white noise
-#define REPORTS (y[obsindex[0]]) // the data
+#define REPORTS (y[obsindex[0]]) // the data
+#define DSDT (f[stateindex[0]])
+#define DIDT (f[stateindex[1]])
+#define DRDT (f[stateindex[2]])
+#define DCDT (f[stateindex[3]])
+#define DWDT (f[stateindex[4]])
+
// the measurement model.
// note that dmeasure and rmeasure are mirrors of each other.
void binomial_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 = dbinom(REPORTS,nearbyint(CASE),exp(LOGRHO),give_log);
+ *lik = dbinom(REPORTS,nearbyint(CASE),RHO,give_log);
}
void binomial_rmeasure (double *y, double *x, double *p,
int *obsindex, int *stateindex, int *parindex, int *covindex,
int ncovars, double *covars, double t) {
- REPORTS = rbinom(nearbyint(CASE),exp(LOGRHO));
+ REPORTS = rbinom(nearbyint(CASE),RHO);
}
-#undef REPORTS
-
// the process model:
// an SIR model with Euler-multinomial step,
// transmission is seasonal, implemented using B-spline basis functions passed as covariates.
@@ -67,65 +64,38 @@
int nrate = 6;
double rate[nrate]; // transition rates
double trans[nrate]; // transition numbers
- double gamma, mu, iota, beta_sd, beta_var;
double beta; // transmission rate
double dW; // white noise increment
- int nseas = (int) NBASIS; // number of seasonal basis functions
+ int nbasis = (int) NBASIS; // number of seasonal basis functions
int deg = (int) DEGREE; // degree of seasonal basis functions
- double seasonality[nseas];
+ int k;
+ double seasonality[nbasis];
void (*eval_basis)(double,double,int,int,double*);
void (*reulmult)(int,double,double*,double,double*);
- double (*dot)(int,const double *, const double *);
- // untransform the parameters
- gamma = exp(LOGGAMMA);
- mu = exp(LOGMU);
- iota = exp(LOGIOTA);
- beta_sd = exp(LOGBETA_SD);
- beta_var = beta_sd*beta_sd;
-
// to evaluate the basis functions and compute the transmission rate, use some of
// pomp's built-in C-level facilities:
eval_basis = (void (*)(double,double,int,int,double*)) R_GetCCallable("pomp","periodic_bspline_basis_eval");
- dot = (double (*)(int,const double *,const double *)) R_GetCCallable("pomp","dot_product");
// pomp's C-level eulermultinomial simulator
reulmult = (void (*)(int,double,double*,double,double*)) R_GetCCallable("pomp","reulermultinom");
// compute transmission rate from seasonality
- if (nseas <= 0) return;
- (*eval_basis)(t,PERIOD,deg,nseas,&seasonality[0]); // evaluate the periodic B-spline basis
- beta = exp((*dot)(nseas,&seasonality[0],&LOGBETA)); // compute the transmission rate
+ if (nbasis <= 0) return;
+ (*eval_basis)(t,PERIOD,deg,nbasis,&seasonality[0]); // evaluate the periodic B-spline basis
+ for (k = 0, beta = 0; k < nbasis; k++)
+ beta += log(BETA[k])*seasonality[k];
+ beta = exp(beta);
- // test to make sure the parameters and state variable values are sane
- // if not, return without doing any work
- if (!(R_FINITE(beta)) ||
- !(R_FINITE(gamma)) ||
- !(R_FINITE(mu)) ||
- !(R_FINITE(beta_sd)) ||
- !(R_FINITE(iota)) ||
- !(R_FINITE(POPSIZE)) ||
- !(R_FINITE(SUSC)) ||
- !(R_FINITE(INFD)) ||
- !(R_FINITE(RCVD)) ||
- !(R_FINITE(CASE)) ||
- !(R_FINITE(W)))
- return;
-
// compute the environmental stochasticity
- if (beta_sd > 0.0) { // environmental noise is ON
- dW = rgamma(dt/beta_var,beta_var); // gamma noise, mean=dt, variance=(beta_sd^2 dt)
- if (!(R_FINITE(dW))) return;
- } else { // environmental noise is OFF
- dW = dt;
- }
+ dW = rgammawn(BETA_SD,dt);
// compute the transition rates
- rate[0] = mu*POPSIZE; // birth into susceptible class
- rate[1] = (iota+beta*INFD*dW/dt)/POPSIZE; // force of infection
- rate[2] = mu; // death from susceptible class
- rate[3] = gamma; // recovery
- rate[4] = mu; // death from infectious class
- rate[5] = mu; // death from recovered class
+ rate[0] = MU*POPSIZE; // birth into susceptible class
+ rate[1] = (IOTA+beta*INFD*dW/dt)/POPSIZE; // force of infection
+ rate[2] = MU; // death from susceptible class
+ rate[3] = GAMMA; // recovery
+ rate[4] = MU; // death from infectious class
+ rate[5] = MU; // death from recovered class
// compute the transition numbers
trans[0] = rpois(rate[0]*dt); // births are Poisson
@@ -138,17 +108,10 @@
INFD += trans[1]-trans[3]-trans[4];
RCVD += trans[3]-trans[5];
CASE += trans[3]; // cases are cumulative recoveries
- if (beta_sd > 0.0) {
- W += (dW-dt)/beta_sd; // increment mean zero, variance = dt
- }
+ if (BETA_SD > 0.0) W += (dW-dt)/BETA_SD; // increment has mean = 0, variance = dt
}
-#define DSDT (f[stateindex[0]])
-#define DIDT (f[stateindex[1]])
-#define DRDT (f[stateindex[2]])
-#define DCDT (f[stateindex[3]])
-
void sir_ODE (double *f, double *x, const double *p,
const int *stateindex, const int *parindex, const int *covindex,
int covdim, const double *covar, double t)
@@ -156,36 +119,31 @@
int nrate = 6;
double rate[nrate]; // transition rates
double term[nrate]; // terms in the equations
- double gamma, mu, iota;
double beta;
- int nseas = (int) NBASIS; // number of seasonal basis functions
+ int nbasis = (int) NBASIS; // number of seasonal basis functions
int deg = (int) DEGREE; // degree of seasonal basis functions
- double seasonality[nseas];
+ int k;
+ double seasonality[nbasis];
void (*eval_basis)(double,double,int,int,double*);
- double (*dot)(int,const double *, const double *);
- // untransform the parameters
- gamma = exp(LOGGAMMA);
- mu = exp(LOGMU);
- iota = exp(LOGIOTA);
-
// to evaluate the basis functions and compute the transmission rate, use some of
// pomp's built-in C-level facilities:
eval_basis = (void (*)(double,double,int,int,double*)) R_GetCCallable("pomp","periodic_bspline_basis_eval");
- dot = (double (*)(int,const double *,const double *)) R_GetCCallable("pomp","dot_product");
// compute transmission rate from seasonality
- if (nseas <= 0) return;
- (*eval_basis)(t,PERIOD,deg,nseas,&seasonality[0]); // evaluate the periodic B-spline basis
- beta = exp((*dot)(nseas,&seasonality[0],&LOGBETA));
+ if (nbasis <= 0) return;
+ (*eval_basis)(t,PERIOD,deg,nbasis,&seasonality[0]); // evaluate the periodic B-spline basis
+ for (k = 0, beta = 0; k < nbasis; k++)
+ beta += log(BETA[k])*seasonality[k];
+ beta = exp(beta);
// compute the transition rates
- rate[0] = mu*POPSIZE; // birth into susceptible class
- rate[1] = (iota+beta*INFD)/POPSIZE; // force of infection
- rate[2] = mu; // death from susceptible class
- rate[3] = gamma; // recovery
- rate[4] = mu; // death from infectious class
- rate[5] = mu; // death from recovered class
+ rate[0] = MU*POPSIZE; // birth into susceptible class
+ rate[1] = (IOTA+beta*INFD)/POPSIZE; // force of infection
+ rate[2] = MU; // death from susceptible class
+ rate[3] = GAMMA; // recovery
+ rate[4] = MU; // death from infectious class
+ rate[5] = MU; // death from recovered class
// compute the several terms
term[0] = rate[0];
@@ -200,27 +158,6 @@
DIDT = term[1]-term[3]-term[4];
DRDT = term[3]-term[5];
DCDT = term[3]; // accumulate the new I->R transitions
+ DWDT = 0;
}
-
-#undef DSDT
-#undef DIDT
-#undef DRDT
-#undef DCDT
-
-#undef SUSC
-#undef INFD
-#undef RCVD
-#undef CASE
-#undef W
-
-#undef LOGGAMMA
-#undef LOGMU
-#undef LOGIOTA
-#undef LOGBETA
-#undef LOGBETA_SD
-#undef POPSIZE
-#undef LOGRHO
-#undef NBASIS
-#undef DEGREE
-#undef PERIOD
Modified: pkg/inst/include/pomp.h
===================================================================
--- pkg/inst/include/pomp.h 2012-03-27 21:17:56 UTC (rev 634)
+++ pkg/inst/include/pomp.h 2012-03-29 15:30:55 UTC (rev 635)
@@ -7,33 +7,6 @@
#include <Rmath.h>
#include <Rdefines.h>
-// prototypes for C-level access to Euler-multinomial distribution functions
-// NB: 'reulermultinom' does not call GetRNGstate() and PutRNGstate() internally
-void reulermultinom (int ntrans, double size, double *rate, double dt, double *trans);
-double deulermultinom (int ntrans, double size, double *rate, double dt, double *trans, int give_log);
-
-// This function computes r such that if
-// N ~ geometric(prob=1-exp(-r dt)) and T ~ exponential(rate=R),
-// then E[N dt] = E[T]
-// i.e., the rate r for an Euler process that gives the same
-// expected waiting time as the exponential process it approximates.
-// In particular r -> R as dt -> 0.
-inline double exp2geom_rate_correction (double R, double dt) {
- return (dt > 0) ? log(1.0+R*dt)/dt : R;
-}
-
-// This function draws a random increment of a gamma whitenoise process.
-// This will have expectation=dt and variance=(sigma^2*dt)
-// If dW = rgammawn(sigma,dt), then
-// mu dW/dt is a candidate for a random rate process within an
-// Euler-multinomial context, i.e.,
-// E[mu*dW] = mu*dt and Var[mu*dW] = mu*sigma^2*dt
-double rgammawn (double sigma, double dt);
-
-// facility for computing the inner produce of
-// a vector of parameters ('coef') against a vector of basis-function values ('basis')
-double dot_product (int dim, const double *basis, const double *coef);
-
// facility for computing evaluating a basis of periodic bsplines
void periodic_bspline_basis_eval (double x, double period, int degree, int nbasis, double *y);
@@ -188,4 +161,43 @@
// on output:
// lik = pointer to scalar containing (log) likelihood
+// This function computes r such that if
+// N ~ geometric(prob=1-exp(-r dt)) and T ~ exponential(rate=R),
+// then E[N dt] = E[T]
+// i.e., the rate r for an Euler process that gives the same
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 635
More information about the pomp-commits
mailing list