[Pomp-commits] r717 - in pkg/pompExamples: inst inst/data-R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat May 12 17:35:25 CEST 2012
Author: kingaa
Date: 2012-05-12 17:35:24 +0200 (Sat, 12 May 2012)
New Revision: 717
Modified:
pkg/pompExamples/inst/ChangeLog
pkg/pompExamples/inst/data-R/pertussis.sim.R
pkg/pompExamples/src/pertussis.c
Log:
- move parameter transformations into C for the pertussis example
Modified: pkg/pompExamples/inst/ChangeLog
===================================================================
--- pkg/pompExamples/inst/ChangeLog 2012-05-12 15:09:24 UTC (rev 716)
+++ pkg/pompExamples/inst/ChangeLog 2012-05-12 15:35:24 UTC (rev 717)
@@ -1,5 +1,16 @@
+2012-05-08 kingaa
+
+ * [r712] DESCRIPTION, tests/pertussis.Rout.save: - bring unit tests
+ up to date
+ * [r710] DESCRIPTION: - update version number
+ * [r709] inst/doc/budmoth-model.pdf, inst/doc/pertussis-model.pdf:
+ - update vignettes
+ * [r708] data/budmoth.sim.rda, data/pertussis.sim.rda: - update for
+ compatibility with pomp version 0.42-2
+
2012-04-21 kingaa
+ * [r672] inst/ChangeLog: - add ChangeLog
* [r671] DESCRIPTION, data/budmoth.sim.rda, data/pertussis.sim.rda,
inst/NEWS, inst/data-R/budmoth.sim.R,
inst/doc/budmoth-model-slices.rda,
Modified: pkg/pompExamples/inst/data-R/pertussis.sim.R
===================================================================
--- pkg/pompExamples/inst/data-R/pertussis.sim.R 2012-05-12 15:09:24 UTC (rev 716)
+++ pkg/pompExamples/inst/data-R/pertussis.sim.R 2012-05-12 15:35:24 UTC (rev 717)
@@ -43,7 +43,8 @@
"birthrate","deathrate","mean.beta","ampl.beta",
"imports","sigma","gamma","alpha","alpha.ratio",
"report.prob","boost.prob","polar.prob",
- "foi.mod","noise.sigma","popsize","tau"
+ "foi.mod","noise.sigma","popsize","tau",
+ "S.0","E.0","I.0","R1.0","R2.0"
),
statenames=c("S","E","I","R1","R2","cases","W","err","simpop"),
zeronames=c("cases","err"),
@@ -51,12 +52,8 @@
comp.names=c("S","E","I","R1","R2"),
rmeasure = "negbin_rmeasure",
dmeasure = "negbin_dmeasure",
- logitvar=c("report.prob","boost.prob","polar.prob"),
- logvar=c(
- "birthrate","deathrate","mean.beta","ampl.beta","imports","sigma","gamma",
- "alpha","alpha.ratio","foi.mod","noise.sigma","tau",
- "S.0","E.0","I.0","R1.0","R2.0"
- ),
+ parameter.inv.transform="pertussis_par_untrans",
+ parameter.transform="pertussis_par_trans",
initializer = function (params, t0, statenames, comp.names, ivps, ...) {
states <- numeric(length(statenames))
names(states) <- statenames
@@ -65,18 +62,6 @@
states[comp.names] <- round(params['popsize']*frac/sum(frac))
states["simpop"] <- params["popsize"]
states
- },
- parameter.inv.transform=function (params, logitvar, logvar, ivps, ...) {
- params[ivps] <- params[ivps]/sum(params[ivps],na.rm=TRUE)
- params[logitvar] <- qlogis(params[logitvar])
- params[logvar] <- log(params[logvar])
- params
- },
- parameter.transform=function (params, logitvar, logvar, ivps, ...) {
- params[logitvar] <- plogis(params[logitvar])
- params[logvar] <- exp(params[logvar])
- params[ivps] <- params[ivps]/sum(params[ivps],na.rm=TRUE)
- params
}
)
coef(po) <- unlist(params[d,])
Modified: pkg/pompExamples/src/pertussis.c
===================================================================
--- pkg/pompExamples/src/pertussis.c 2012-05-12 15:09:24 UTC (rev 716)
+++ pkg/pompExamples/src/pertussis.c 2012-05-12 15:35:24 UTC (rev 717)
@@ -5,6 +5,10 @@
#include <Rdefines.h>
#include <R_ext/Rdynload.h>
+inline double expit (double x) {return 1.0/(1.0+exp(-x));}
+
+inline double logit (double p) {return log(p/(1.0-p));}
+
static inline double rgammawn (double sigma, double dt) {
double sigmasq;
sigmasq = sigma*sigma;
@@ -36,6 +40,11 @@
#define NOISE_SIGMA (p[parindex[13]]) // white noise intensity
#define POP (p[parindex[14]]) // population size
#define TAU (p[parindex[15]]) // reporting SD
+#define S0 (p[parindex[16]]) // initial S
+#define E0 (p[parindex[17]]) // initial E
+#define I0 (p[parindex[18]]) // initial I
+#define R10 (p[parindex[19]]) // initial R1
+#define R20 (p[parindex[20]]) // initial R2
#define SUS (x[stateindex[0]]) // susceptibles
#define EXP (x[stateindex[1]]) // exposed
@@ -59,6 +68,69 @@
#define REPORTS (y[obsindex[0]])
+void pertussis_par_untrans (double *pt, double *p, int *parindex)
+{
+ double sum;
+ pt[parindex[0]] = log(BIRTHRATE);
+ pt[parindex[1]] = log(DEATHRATE);
+ pt[parindex[2]] = log(MEANBETA);
+ pt[parindex[3]] = log(AMPLBETA);
+ pt[parindex[4]] = log(IMPORTS);
+ pt[parindex[5]] = log(SIGMA);
+ pt[parindex[6]] = log(GAMMA);
+ pt[parindex[7]] = log(ALPHA);
+ pt[parindex[8]] = log(ALPHA_RATIO);
+ pt[parindex[9]] = logit(REPORT_PROB);
+ pt[parindex[10]] = logit(BOOST_PROB);
+ pt[parindex[11]] = logit(POLAR_PHI);
+ pt[parindex[12]] = log(FOI_MOD);
+ pt[parindex[13]] = log(NOISE_SIGMA);
+ pt[parindex[15]] = log(TAU);
+
+ sum = S0+E0+I0+R10+R20;
+ pt[parindex[16]] = log(S0/sum);
+ pt[parindex[17]] = log(E0/sum);
+ pt[parindex[18]] = log(I0/sum);
+ pt[parindex[19]] = log(R10/sum);
+ pt[parindex[20]] = log(R20/sum);
+
+}
+
+void pertussis_par_trans (double *pt, double *p, int *parindex)
+{
+ double sum;
+ double s, e, i, r1, r2;
+
+ pt[parindex[0]] = exp(BIRTHRATE);
+ pt[parindex[1]] = exp(DEATHRATE);
+ pt[parindex[2]] = exp(MEANBETA);
+ pt[parindex[3]] = exp(AMPLBETA);
+ pt[parindex[4]] = exp(IMPORTS);
+ pt[parindex[5]] = exp(SIGMA);
+ pt[parindex[6]] = exp(GAMMA);
+ pt[parindex[7]] = exp(ALPHA);
+ pt[parindex[8]] = exp(ALPHA_RATIO);
+ pt[parindex[9]] = expit(REPORT_PROB);
+ pt[parindex[10]] = expit(BOOST_PROB);
+ pt[parindex[11]] = expit(POLAR_PHI);
+ pt[parindex[12]] = exp(FOI_MOD);
+ pt[parindex[13]] = exp(NOISE_SIGMA);
+ pt[parindex[15]] = exp(TAU);
+
+ s = exp(S0);
+ e = exp(E0);
+ i = exp(I0);
+ r1 = exp(R10);
+ r2 = exp(R20);
+ sum = s+e+i+r1+r2;
+ pt[parindex[16]] = s/sum;
+ pt[parindex[17]] = e/sum;
+ pt[parindex[18]] = i/sum;
+ pt[parindex[19]] = r1/sum;
+ pt[parindex[20]] = r2/sum;
+
+}
+
void pertussis_sveirr_EM (double *x, double *p,
int *stateindex, int *parindex, int *covindex,
int covdim, double *covar, double t, double dt)
More information about the pomp-commits
mailing list