[Pomp-commits] r889 - in pkg/pomp: . inst/examples src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 10 14:38:56 CET 2014
Author: kingaa
Date: 2014-03-10 14:38:56 +0100 (Mon, 10 Mar 2014)
New Revision: 889
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/inst/examples/dacca.R
pkg/pomp/src/cholmodel.c
pkg/pomp/tests/dacca.R
pkg/pomp/tests/dacca.Rout.save
Log:
- put parameter transforms into C (3X speedup)
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2014-03-09 16:31:45 UTC (rev 888)
+++ pkg/pomp/DESCRIPTION 2014-03-10 13:38:56 UTC (rev 889)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.47-4
-Date: 2014-03-09
+Version: 0.47-5
+Date: 2014-03-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")),
Modified: pkg/pomp/inst/examples/dacca.R
===================================================================
--- pkg/pomp/inst/examples/dacca.R 2014-03-09 16:31:45 UTC (rev 888)
+++ pkg/pomp/inst/examples/dacca.R 2014-03-10 13:38:56 UTC (rev 889)
@@ -122,28 +122,14 @@
statenames = c("S","I","Rs","R1","M","W","count"),
paramnames = c("tau","gamma","eps","delta","deltaI",
"log.omega1","sd.beta","beta.trend","log.beta1",
- "alpha","rho","clin","nbasis","nrstage"),
+ "alpha","rho","clin","nbasis","nrstage",
+ "S.0","I.0","Rs.0","R1.0"),
covarnames = c("pop","dpopdt","seas.1","trend"),
all.state.names=c("S","I","Rs",paste("R",1:nrstage,sep=''),"M","W","count"),
comp.names=c("S","I","Rs",paste("R",1:nrstage,sep='')),
comp.ic.names=c("S.0","I.0","Rs.0",paste("R",1:nrstage,".0",sep='')),
- log.trans=c( # parameters to log transform
- "gamma","eps","rho","delta","deltaI","alpha",
- "tau","sd.beta",
- "S.0","I.0","Rs.0",paste("R",1:nrstage,".0",sep='')
- ),
- logit.trans="clin", # parameters to logit transform
- parameter.transform=function (params, log.trans, logit.trans, comp.ic.names, ...) {
- params[logit.trans] <- plogis(params[logit.trans])
- params[log.trans] <- exp(params[log.trans])
- params[comp.ic.names] <- params[comp.ic.names]/sum(params[comp.ic.names])
- params
- },
- parameter.inv.transform=function (params, log.trans, logit.trans, comp.names, ...) {
- params[logit.trans] <- qlogis(params[logit.trans])
- params[log.trans] <- log(params[log.trans])
- params
- },
+ parameter.transform="_cholmodel_trans",
+ parameter.inv.transform="_cholmodel_untrans",
initializer = function (params, t0, covars, nrstage, comp.ic.names, comp.names, all.state.names, ...) {
states <- numeric(length(all.state.names))
names(states) <- all.state.names
Modified: pkg/pomp/src/cholmodel.c
===================================================================
--- pkg/pomp/src/cholmodel.c 2014-03-09 16:31:45 UTC (rev 888)
+++ pkg/pomp/src/cholmodel.c 2014-03-10 13:38:56 UTC (rev 889)
@@ -18,6 +18,10 @@
#define CLIN (p[parindex[11]])
#define NBASIS (p[parindex[12]])
#define NRSTAGE (p[parindex[13]])
+#define S0 (p[parindex[14]])
+#define I0 (p[parindex[15]])
+#define RS0 (p[parindex[16]])
+#define RL0 (p[parindex[17]])
#define SUSCEP (x[stateindex[0]])
#define INFECT (x[stateindex[1]])
@@ -34,6 +38,52 @@
#define DATADEATHS (y[obsindex[0]])
+void _cholmodel_untrans (double *pt, double *p, int *parindex)
+{
+ int k, nrstage = (int) NRSTAGE;
+ pt[parindex[0]] = log(TAU);
+ pt[parindex[1]] = log(GAMMA);
+ pt[parindex[2]] = log(EPS);
+ pt[parindex[3]] = log(DELTA);
+ pt[parindex[4]] = log(DELTA_I);
+ pt[parindex[6]] = log(SD_BETA);
+ pt[parindex[9]] = log(ALPHA);
+ pt[parindex[10]] = log(RHO);
+ pt[parindex[11]] = logit(CLIN);
+
+ pt[parindex[14]] = log(S0);
+ pt[parindex[15]] = log(I0);
+ pt[parindex[16]] = log(RS0);
+ for (k = 0; k < nrstage; k++)
+ pt[parindex[17]+k] = log((&RL0)[k]);
+}
+
+void _cholmodel_trans (double *pt, double *p, int *parindex)
+{
+ int k, nrstage = (int) NRSTAGE;
+ double sum = 0.0;
+ pt[parindex[0]] = exp(TAU);
+ pt[parindex[1]] = exp(GAMMA);
+ pt[parindex[2]] = exp(EPS);
+ pt[parindex[3]] = exp(DELTA);
+ pt[parindex[4]] = exp(DELTA_I);
+ pt[parindex[6]] = exp(SD_BETA);
+ pt[parindex[9]] = exp(ALPHA);
+ pt[parindex[10]] = exp(RHO);
+ pt[parindex[11]] = expit(CLIN);
+
+ sum += (pt[parindex[14]] = exp(S0));
+ sum += (pt[parindex[15]] = exp(I0));
+ sum += (pt[parindex[16]] = exp(RS0));
+ for (k = 0; k < nrstage; k++)
+ sum += (pt[parindex[17]+k] = exp((&RL0)[k]));
+ pt[parindex[14]] /= sum;
+ pt[parindex[15]] /= sum;
+ pt[parindex[16]] /= sum;
+ for (k = 0; k < nrstage; k++)
+ pt[parindex[17]+k] /= sum;
+}
+
void _cholmodel_norm_rmeasure (double *y, double *x, double *p,
int *obsindex, int *stateindex, int *parindex, int *covindex,
int ncovars, double *covars, double t)
Modified: pkg/pomp/tests/dacca.R
===================================================================
--- pkg/pomp/tests/dacca.R 2014-03-09 16:31:45 UTC (rev 888)
+++ pkg/pomp/tests/dacca.R 2014-03-10 13:38:56 UTC (rev 889)
@@ -15,6 +15,9 @@
pf <- pfilter(dacca,Np=1000,seed=5886855L)
print(round(logLik(pf),digits=1))
+pf1 <- pfilter(simulate(dacca),Np=1000,seed=5886855L)
+print(round(logLik(pf1),digits=1))
+
## to investigate the rogue crash:
dacca.pars <- c("gamma","eps","deltaI","beta.trend",
Modified: pkg/pomp/tests/dacca.Rout.save
===================================================================
--- pkg/pomp/tests/dacca.Rout.save 2014-03-09 16:31:45 UTC (rev 888)
+++ pkg/pomp/tests/dacca.Rout.save 2014-03-10 13:38:56 UTC (rev 889)
@@ -1,5 +1,5 @@
-R version 3.0.2 (2013-09-25) -- "Frisbee Sailing"
+R version 3.0.1 (2013-05-16) -- "Good Sport"
Copyright (C) 2013 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)
@@ -49,6 +49,10 @@
> print(round(logLik(pf),digits=1))
[1] -3748.6
>
+> pf1 <- pfilter(simulate(dacca),Np=1000,seed=5886855L)
+> print(round(logLik(pf1),digits=1))
+[1] -3701.5
+>
> ## to investigate the rogue crash:
>
> dacca.pars <- c("gamma","eps","deltaI","beta.trend",
@@ -131,4 +135,4 @@
>
> proc.time()
user system elapsed
- 28.417 0.076 28.599
+ 36.898 0.108 37.070
More information about the pomp-commits
mailing list