[Pomp-commits] r888 - in pkg/pomp: . inst src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Mar 9 17:31:46 CET 2014
Author: kingaa
Date: 2014-03-09 17:31:45 +0100 (Sun, 09 Mar 2014)
New Revision: 888
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/inst/NEWS
pkg/pomp/src/cholmodel.c
pkg/pomp/tests/dacca.R
pkg/pomp/tests/dacca.Rout.save
Log:
- in 'dacca' model, positivity violations are now treated a bit differently
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2014-03-02 14:10:33 UTC (rev 887)
+++ pkg/pomp/DESCRIPTION 2014-03-09 16:31:45 UTC (rev 888)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.47-3
-Date: 2014-03-01
+Version: 0.47-4
+Date: 2014-03-09
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/NEWS
===================================================================
--- pkg/pomp/inst/NEWS 2014-03-02 14:10:33 UTC (rev 887)
+++ pkg/pomp/inst/NEWS 2014-03-09 16:31:45 UTC (rev 888)
@@ -1,4 +1,7 @@
NEWS
+0.47-4
+ o Revisit 'dacca' bug. Set negative compartments to zero along with compartments immediately downstream.
+
0.47-3
o Fix bug that arises only very occasionally in the 'dacca' cholera example.
Modified: pkg/pomp/src/cholmodel.c
===================================================================
--- pkg/pomp/src/cholmodel.c 2014-03-02 14:10:33 UTC (rev 887)
+++ pkg/pomp/src/cholmodel.c 2014-03-09 16:31:45 UTC (rev 888)
@@ -89,7 +89,7 @@
double beta;
double omega;
double dw;
- const double *pt;
+ double *pt;
int j;
if (COUNT != 0.0) return;
@@ -110,9 +110,9 @@
rsdeaths = DELTA*RSHORT; // natural Rs deaths
wanings = RHO*RSHORT; // loss of immunity
- for (pt = &RLONG, j = 0; j < nrstage; j++) {
- rldeaths[j] = DELTA*pt[j]; // natural R deaths
- passages[j+1] = eps*pt[j]; // passage to the next immunity class
+ for (pt = &RLONG, j = 0; j < nrstage; j++, pt++) {
+ rldeaths[j] = *pt*DELTA; // natural R deaths
+ passages[j+1] = *pt*eps; // passage to the next immunity class
}
infections = (omega+(beta+SD_BETA*dw/dt)*effI)*SUSCEP; // infection
@@ -121,20 +121,39 @@
SUSCEP += (births - infections - sdeaths + passages[nrstage] + wanings)*dt;
INFECT += (CLIN*infections - disease - ideaths - passages[0])*dt;
RSHORT += ((1-CLIN)*infections - rsdeaths - wanings)*dt;
- for (j = 0; j < nrstage; j++)
- (&RLONG)[j] += (passages[j] - passages[j+1] - rldeaths[j])*dt;
+ for (pt = &RLONG, j = 0; j < nrstage; j++, pt++)
+ *pt += (passages[j] - passages[j+1] - rldeaths[j])*dt;
DEATHS += disease*dt; // cumulative deaths due to disease
NOISE += dw;
// check for violations of positivity constraints
// nonzero COUNT variable signals violation
- if (SUSCEP < 0.0) { COUNT += 1; return;}
- if (INFECT < 0.0) { COUNT += 1e3; return;}
- if (RSHORT < 0.0) { COUNT += 1e6; return;}
- if (DEATHS < 0.0) { COUNT += 1e9; return;}
- for (j = 0; j < nrstage; j++)
- if ((&RLONG)[j] < 0.0) { COUNT += 1e12; return;}
-
+ if (SUSCEP < 0.0) {
+ SUSCEP = 0.0; INFECT = 0.0; RSHORT = 0.0;
+ COUNT += 1;
+ }
+ if (INFECT < 0.0) {
+ INFECT = 0.0; SUSCEP = 0.0;
+ COUNT += 1e3;
+ }
+ if (RSHORT < 0.0) {
+ RSHORT = 0.0; SUSCEP = 0.0;
+ COUNT += 1e6;
+ }
+ if (DEATHS < 0.0) {
+ DEATHS = 0.0;
+ COUNT += 1e9;
+ }
+ for (pt = &RLONG, j = 0; j < nrstage-1; j++, pt++) {
+ if (*pt < 0.0) {
+ *pt = 0.0; *(pt+1) = 0.0;
+ COUNT += 1e12;
+ }
+ }
+ if (*pt < 0.0) {
+ *pt = 0.0; SUSCEP = 0.0;
+ COUNT += 1e12;
+ }
}
#undef GAMMA
Modified: pkg/pomp/tests/dacca.R
===================================================================
--- pkg/pomp/tests/dacca.R 2014-03-02 14:10:33 UTC (rev 887)
+++ pkg/pomp/tests/dacca.R 2014-03-09 16:31:45 UTC (rev 888)
@@ -43,6 +43,7 @@
min=unlist(param.tab["box_min",]),
max=unlist(param.tab["box_max",])
)
+
dacca.rprior <- function (hyperparams, ...) {
r <- runif(length(hyperparams$min),min=hyperparams$min,max=hyperparams$max)
names(r) <- names(hyperparams$min)
@@ -51,8 +52,35 @@
set.seed(7777+7)
params.tricky <- dacca.rprior(dacca.hyperparams)
-m7 <- mif(dacca,Nmif=2,start=params.tricky,pars=dacca.pars,ivps=dacca.ivps,
- Np=100,ic.lag=600,method="mif2",rw.sd=dacca.rw.sd,
- cooling.type="geometric",cooling.fraction=sqrt(0.1),var.factor=2,
- transform=TRUE)
+m7 <- mif(
+ dacca,
+ Nmif=2,
+ start=params.tricky,
+ pars=dacca.pars,
+ ivps=dacca.ivps,
+ Np=100,
+ method="mif2",
+ rw.sd=dacca.rw.sd,
+ cooling.type="geometric",
+ cooling.fraction=sqrt(0.1),
+ var.factor=2,
+ transform=TRUE
+ )
m7 <- continue(m7)
+
+set.seed(12350)
+th.draw <- dacca.rprior(dacca.hyperparams)
+m1 <- mif(
+ dacca,
+ Nmif=10,
+ Np=100,
+ start=th.draw,
+ pars=dacca.pars,
+ ivps=dacca.ivps,
+ method="mif2",
+ rw.sd=dacca.rw.sd,
+ cooling.type="geometric",
+ cooling.fraction=sqrt(0.1),
+ var.factor=2,
+ transform=TRUE
+ )
Modified: pkg/pomp/tests/dacca.Rout.save
===================================================================
--- pkg/pomp/tests/dacca.Rout.save 2014-03-02 14:10:33 UTC (rev 887)
+++ pkg/pomp/tests/dacca.Rout.save 2014-03-09 16:31:45 UTC (rev 888)
@@ -1,6 +1,6 @@
-R Under development (unstable) (2014-02-22 r65060) -- "Unsuffered Consequences"
-Copyright (C) 2014 The R Foundation for Statistical Computing
+R version 3.0.2 (2013-09-25) -- "Frisbee Sailing"
+Copyright (C) 2013 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -77,6 +77,7 @@
+ min=unlist(param.tab["box_min",]),
+ max=unlist(param.tab["box_max",])
+ )
+>
> dacca.rprior <- function (hyperparams, ...) {
+ r <- runif(length(hyperparams$min),min=hyperparams$min,max=hyperparams$max)
+ names(r) <- names(hyperparams$min)
@@ -85,14 +86,49 @@
>
> set.seed(7777+7)
> params.tricky <- dacca.rprior(dacca.hyperparams)
-> m7 <- mif(dacca,Nmif=2,start=params.tricky,pars=dacca.pars,ivps=dacca.ivps,
-+ Np=100,ic.lag=600,method="mif2",rw.sd=dacca.rw.sd,
-+ cooling.type="geometric",cooling.fraction=sqrt(0.1),var.factor=2,
-+ transform=TRUE)
-Warning message:
-1 filtering failure occurred in 'pfilter'
+> m7 <- mif(
++ dacca,
++ Nmif=2,
++ start=params.tricky,
++ pars=dacca.pars,
++ ivps=dacca.ivps,
++ Np=100,
++ method="mif2",
++ rw.sd=dacca.rw.sd,
++ cooling.type="geometric",
++ cooling.fraction=sqrt(0.1),
++ var.factor=2,
++ transform=TRUE
++ )
+Warning messages:
+1: 11 filtering failures occurred in 'pfilter'
+2: 4 filtering failures occurred in 'pfilter'
> m7 <- continue(m7)
>
+> set.seed(12350)
+> th.draw <- dacca.rprior(dacca.hyperparams)
+> m1 <- mif(
++ dacca,
++ Nmif=10,
++ Np=100,
++ start=th.draw,
++ pars=dacca.pars,
++ ivps=dacca.ivps,
++ method="mif2",
++ rw.sd=dacca.rw.sd,
++ cooling.type="geometric",
++ cooling.fraction=sqrt(0.1),
++ var.factor=2,
++ transform=TRUE
++ )
+Warning messages:
+1: 1 filtering failure occurred in 'pfilter'
+2: 16 filtering failures occurred in 'pfilter'
+3: 2 filtering failures occurred in 'pfilter'
+4: 2 filtering failures occurred in 'pfilter'
+5: 1 filtering failure occurred in 'pfilter'
+6: 7 filtering failures occurred in 'pfilter'
+>
> proc.time()
user system elapsed
- 8.056 0.044 8.137
+ 28.417 0.076 28.599
More information about the pomp-commits
mailing list