[Pomp-commits] r727 - in pkg/pomp: . demo inst/include src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jun 1 18:27:00 CEST 2012
Author: kingaa
Date: 2012-06-01 18:27:00 +0200 (Fri, 01 Jun 2012)
New Revision: 727
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/demo/sir.R
pkg/pomp/inst/include/pomp.h
pkg/pomp/src/lookup_table.c
pkg/pomp/src/pomp.h
Log:
- remove superfluous warnings from 'lookup_table'
- fix 'sir' demo to work with covariate table
- add some inline functions to 'pomp.h'
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2012-06-01 13:55:15 UTC (rev 726)
+++ pkg/pomp/DESCRIPTION 2012-06-01 16:27:00 UTC (rev 727)
@@ -1,7 +1,7 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.42-5
+Version: 0.42-6
Date: 2012-06-01
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/pomp/demo/sir.R
===================================================================
--- pkg/pomp/demo/sir.R 2012-06-01 13:55:15 UTC (rev 726)
+++ pkg/pomp/demo/sir.R 2012-06-01 16:27:00 UTC (rev 727)
@@ -1,26 +1,24 @@
require(pomp)
+## negative binomial measurement model
dmeas <- "
- lik = dbinom(cases,nearbyint(incid),rho,give_log);
+ double prob = theta/(theta+rho*incid);
+ lik = dnbinom(cases,theta,prob,give_log);
"
rmeas <- "
- cases = rbinom(nearbyint(incid),rho);
+ double prob = theta/(theta+rho*incid);
+ cases = rnbinom(theta,prob);
"
+## SIR process model with extra-demographic stochasticity
step.fn <- '
int nrate = 6;
double rate[nrate]; // transition rates
double trans[nrate]; // transition numbers
double beta; // transmission rate
double dW; // white noise increment
- double period = 1.0; // period of the seasonality
- int nbasis = 3; // number of seasonality basis functions
- int deg = 3; // degree of the B-spline basis functions
- double seasonality[nbasis];
int k;
- // compute transmission rate from seasonality
- periodic_bspline_basis_eval(t,period,deg,nbasis,&seasonality[0]); // evaluate the periodic B-spline basis
- beta = beta1*seasonality[0]+beta2*seasonality[1]+beta3*seasonality[2];
+ beta = beta1*seas1+beta2*seas2+beta3*seas3;
// compute the environmental stochasticity
dW = rgammawn(beta_sd,dt);
@@ -52,16 +50,9 @@
double term[nrate]; // transition numbers
double beta; // transmission rate
double dW; // white noise increment
- double period = 1.0; // period of the seasonality
- int nbasis = 3;
- int deg = 3;
- double seasonality[nbasis];
int k;
- // compute transmission rate from seasonality
- if (nbasis <= 0) return;
- periodic_bspline_basis_eval(t,period,deg,nbasis,&seasonality[0]); // evaluate the periodic B-spline basis
- beta = beta1*seasonality[0]+beta2*seasonality[1]+beta3*seasonality[2];
+ beta = beta1*seas1+beta2*seas2+beta3*seas3;
// compute the transition rates
rate[0] = mu*popsize; // birth into susceptible class
@@ -86,6 +77,7 @@
Dincid = term[3]; // accumulate the new I->R transitions
DW = 0;
'
+## parameter transformations
partrans <- "
double sum;
Tgamma = exp(gamma);
@@ -96,6 +88,7 @@
Tbeta3 = exp(beta3);
Tbeta_sd = exp(beta_sd);
Trho = expit(rho);
+ Ttheta = exp(theta);
TS_0 = exp(S_0);
TI_0 = exp(I_0);
TR_0 = exp(R_0);
@@ -114,6 +107,7 @@
Tbeta3 = log(beta3);
Tbeta_sd = log(beta_sd);
Trho = logit(rho);
+ Ttheta = log(theta);
sum = S_0+I_0+R_0;
TS_0 = log(S_0/sum);
TI_0 = log(I_0/sum);
@@ -122,6 +116,19 @@
data(LondonYorke)
+cbind(
+ time=seq(from=1928,to=1934,by=0.01),
+ as.data.frame(
+ periodic.bspline.basis(
+ x=seq(from=1928,to=1934,by=0.01),
+ nbasis=3,
+ degree=3,
+ period=1,
+ names="seas%d"
+ )
+ )
+ ) -> covar
+
pompBuilder(
name="SIR",
data=subset(
@@ -137,12 +144,14 @@
step.fn.delta.t=1/52/20,
skeleton.type="vectorfield",
skeleton=skel,
+ covar=covar,
+ tcovar="time",
parameter.transform=partrans,
parameter.inv.transform=paruntrans,
statenames=c("S","I","R","incid","W"),
paramnames=c(
"gamma","mu","iota","beta1","beta2","beta3","beta.sd",
- "popsize","rho","S.0","I.0","R.0"
+ "popsize","rho","theta","S.0","I.0","R.0"
),
zeronames=c("incid","W"),
comp.names=c("S","I","R"),
@@ -159,9 +168,9 @@
coef(po) <- c(
gamma=26,mu=0.02,iota=0.01,
beta1=120,beta2=140,beta3=100,
- beta.sd=1e-3,
+ beta.sd=0.01,
popsize=5e6,
- rho=0.2,
+ rho=0.2,theta=0.001,
S.0=0.22,I.0=0.0018,R.0=0.78
)
@@ -184,12 +193,16 @@
coef(po) <- coef(
traj.match(
pomp(
- window(po,end=1930),
- measurement.model=cases~norm(mean=rho*incid,sd=1000)
+ window(po,end=1930)
+ ## window(po,end=1930),
+ ## measurement.model=cases~norm(mean=rho*incid,sd=100)
),
- est=c("beta1","beta2","beta3","S.0","I.0","R.0"),
+ est=c("S.0","I.0","R.0"),
transform=TRUE
)
)
+pf <- pfilter(po,Np=1000,max.fail=100)
+print(round(logLik(pf),1))
+
dyn.unload("SIR.so")
Modified: pkg/pomp/inst/include/pomp.h
===================================================================
--- pkg/pomp/inst/include/pomp.h 2012-06-01 13:55:15 UTC (rev 726)
+++ pkg/pomp/inst/include/pomp.h 2012-06-01 16:27:00 UTC (rev 727)
@@ -172,7 +172,7 @@
// 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.
-static inline double exp2geom_rate_correction (double R, double dt) {
+static R_INLINE double exp2geom_rate_correction (double R, double dt) {
return (dt > 0) ? log(1.0+R*dt)/dt : R;
}
@@ -182,7 +182,7 @@
// 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
-static inline double rgammawn (double sigma, double dt) {
+static R_INLINE double rgammawn (double sigma, double dt) {
double sigmasq;
sigmasq = sigma*sigma;
return (sigmasq > 0) ? rgamma(dt/sigmasq,sigmasq) : dt;
@@ -205,4 +205,28 @@
// compute probabilities of eulermultinomial transitions
double deulermultinom (int m, double size, double *rate, double dt, double *trans, int give_log);
+static R_INLINE double rbetabinom (double size, double prob, double theta) {
+ return rbinom(size,rbeta(prob*theta,(1.0-prob)*theta));
+}
+
+static R_INLINE double dbetabinom (double x, double size, double prob, double theta, int give_log) {
+ double a = theta*prob;
+ double b = theta*(1.0-prob);
+ double f = lchoose(size,x)-lbeta(a,b)+lbeta(a+x,b+size-x);
+ return (give_log) ? f : exp(f);
+}
+
+static R_INLINE double rbetanbinom (double mu, double size, double theta) {
+ double prob = size/(size+mu);
+ return rnbinom(size,rbeta(prob*theta,(1.0-prob)*theta));
+}
+
+static R_INLINE double dbetanbinom (double x, double mu, double size, double theta, int give_log) {
+ double prob = size/(size+mu);
+ double a = theta*prob;
+ double b = theta*(1.0-prob);
+ double f = lchoose(size+x-1,size-1)-lbeta(a,b)+lbeta(a+size,b+x);
+ return (give_log) ? f : exp(f);
+}
+
#endif
Modified: pkg/pomp/src/lookup_table.c
===================================================================
--- pkg/pomp/src/lookup_table.c 2012-06-01 13:55:15 UTC (rev 726)
+++ pkg/pomp/src/lookup_table.c 2012-06-01 16:27:00 UTC (rev 727)
@@ -58,8 +58,9 @@
double e;
if ((tab == 0) || (tab->length < 1) || (tab->width < 1)) return;
tab->index = findInterval(tab->x,tab->length,x,TRUE,TRUE,tab->index,&flag);
- if (flag != 0) // we are extrapolating
- warning("table_lookup: extrapolating (flag %d) at %le", flag, x);
+ // warn only if we are *outside* the interval
+ if ((x < tab->x[0]) || (x > tab->x[(tab->length)-1]))
+ warning("table_lookup: extrapolating at %le", x);
e = (x - tab->x[tab->index-1]) / (tab->x[tab->index] - tab->x[tab->index-1]);
for (j = 0; j < tab->width; j++) {
k = j*(tab->length)+(tab->index);
Modified: pkg/pomp/src/pomp.h
===================================================================
--- pkg/pomp/src/pomp.h 2012-06-01 13:55:15 UTC (rev 726)
+++ pkg/pomp/src/pomp.h 2012-06-01 16:27:00 UTC (rev 727)
@@ -172,7 +172,7 @@
// 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.
-static inline double exp2geom_rate_correction (double R, double dt) {
+static R_INLINE double exp2geom_rate_correction (double R, double dt) {
return (dt > 0) ? log(1.0+R*dt)/dt : R;
}
@@ -182,7 +182,7 @@
// 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
-static inline double rgammawn (double sigma, double dt) {
+static R_INLINE double rgammawn (double sigma, double dt) {
double sigmasq;
sigmasq = sigma*sigma;
return (sigmasq > 0) ? rgamma(dt/sigmasq,sigmasq) : dt;
@@ -205,4 +205,28 @@
// compute probabilities of eulermultinomial transitions
double deulermultinom (int m, double size, double *rate, double dt, double *trans, int give_log);
+static R_INLINE double rbetabinom (double size, double prob, double theta) {
+ return rbinom(size,rbeta(prob*theta,(1.0-prob)*theta));
+}
+
+static R_INLINE double dbetabinom (double x, double size, double prob, double theta, int give_log) {
+ double a = theta*prob;
+ double b = theta*(1.0-prob);
+ double f = lchoose(size,x)-lbeta(a,b)+lbeta(a+x,b+size-x);
+ return (give_log) ? f : exp(f);
+}
+
+static R_INLINE double rbetanbinom (double mu, double size, double theta) {
+ double prob = size/(size+mu);
+ return rnbinom(size,rbeta(prob*theta,(1.0-prob)*theta));
+}
+
+static R_INLINE double dbetanbinom (double x, double mu, double size, double theta, int give_log) {
+ double prob = size/(size+mu);
+ double a = theta*prob;
+ double b = theta*(1.0-prob);
+ double f = lchoose(size+x-1,size-1)-lbeta(a,b)+lbeta(a+size,b+x);
+ return (give_log) ? f : exp(f);
+}
+
#endif
More information about the pomp-commits
mailing list