[Pomp-commits] r887 - in pkg/pomp: . inst src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Mar 2 15:10:33 CET 2014
Author: kingaa
Date: 2014-03-02 15:10:33 +0100 (Sun, 02 Mar 2014)
New Revision: 887
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:
- fix bug in 'dacca' model to do with very occasional state-variable positivity violations. The model is slightly different now, but likelihoods at MLE appear to be much the same. This makes sense since positivity violations are extremely rare near the MLE.
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2014-02-27 02:20:24 UTC (rev 886)
+++ pkg/pomp/DESCRIPTION 2014-03-02 14:10:33 UTC (rev 887)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.47-2
-Date: 2014-02-26
+Version: 0.47-3
+Date: 2014-03-01
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-02-27 02:20:24 UTC (rev 886)
+++ pkg/pomp/inst/NEWS 2014-03-02 14:10:33 UTC (rev 887)
@@ -1,4 +1,10 @@
NEWS
+0.47-3
+ o Fix bug that arises only very occasionally in the 'dacca' cholera example.
+
+ o Modify 'dacca' model so that rare positivity violations are punished in a different way.
+ This modifies the formal model slightly.
+
0.47-2
o By default, 'pompBuilder' now creates all files in the temporary directory.
Modified: pkg/pomp/src/cholmodel.c
===================================================================
--- pkg/pomp/src/cholmodel.c 2014-02-27 02:20:24 UTC (rev 886)
+++ pkg/pomp/src/cholmodel.c 2014-03-02 14:10:33 UTC (rev 887)
@@ -92,31 +92,13 @@
const double *pt;
int j;
- if (!(R_FINITE(SUSCEP))) return;
- if (!(R_FINITE(INFECT))) return;
- if (!(R_FINITE(RSHORT))) return;
- for (pt = &RLONG, j = 0; j < nrstage; j++) {
- if (!(R_FINITE(pt[j]))) return;
- }
+ if (COUNT != 0.0) return;
eps = EPS*NRSTAGE;
- if (!(R_FINITE(GAMMA))) return;
- if (!(R_FINITE(eps))) return;
- if (!(R_FINITE(RHO))) return;
- if (!(R_FINITE(DELTA))) return;
- if (!(R_FINITE(DELTA_I))) return;
- if (!(R_FINITE(CLIN))) return;
- if (!(R_FINITE(SD_BETA))) return;
- if (!(R_FINITE(ALPHA))) return;
- if (!(R_FINITE(BETATREND))) return;
-
beta = exp(dot_product(nbasis,&SEASBASIS,&LOGBETA)+BETATREND*TREND);
omega = exp(dot_product(nbasis,&SEASBASIS,&LOGOMEGA));
- if (!(R_FINITE(beta))) return;
- if (!(R_FINITE(omega))) return;
-
dw = rnorm(0,sqrt(dt)); // white noise
effI = pow(INFECT/POP,ALPHA);
@@ -135,21 +117,6 @@
infections = (omega+(beta+SD_BETA*dw/dt)*effI)*SUSCEP; // infection
sdeaths = DELTA*SUSCEP; // natural S deaths
- if (infections > 0.0) {
- if ((infections+sdeaths)*dt > SUSCEP) { // too many leaving S class
- COUNT += 1.0e10;
- return;
- }
- } else {
- if ((-CLIN*infections+disease+ideaths+passages[0])*dt > INFECT) { // too many leaving I class
- COUNT += 1.0e5;
- return;
- }
- if ((-(1-CLIN)*infections+rsdeaths+wanings)*dt > RSHORT) { // too many leaving Rs class
- COUNT += 1;
- return;
- }
- }
SUSCEP += (births - infections - sdeaths + passages[nrstage] + wanings)*dt;
INFECT += (CLIN*infections - disease - ideaths - passages[0])*dt;
@@ -159,6 +126,15 @@
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;}
+
}
#undef GAMMA
Modified: pkg/pomp/tests/dacca.R
===================================================================
--- pkg/pomp/tests/dacca.R 2014-02-27 02:20:24 UTC (rev 886)
+++ pkg/pomp/tests/dacca.R 2014-03-02 14:10:33 UTC (rev 887)
@@ -14,3 +14,45 @@
pf <- pfilter(dacca,Np=1000,seed=5886855L)
print(round(logLik(pf),digits=1))
+
+## to investigate the rogue crash:
+
+dacca.pars <- c("gamma","eps","deltaI","beta.trend",
+ "log.beta1","log.beta2","log.beta3",
+ "log.beta4","log.beta5","log.beta6",
+ "log.omega1","log.omega2","log.omega3",
+ "log.omega4","log.omega5","log.omega6",
+ "sd.beta","tau")
+dacca.ivps <- c("S.0","I.0","R1.0","R2.0","R3.0")
+dacca.rw.sd <- c(
+ rep(0.1,length(dacca.pars)),
+ rep(0.2,length(dacca.ivps))
+ )
+names(dacca.rw.sd) <- c(dacca.pars,dacca.ivps)
+
+param.tab <- read.csv2(text='
+"";"gamma";"eps";"rho";"delta";"deltaI";"clin";"alpha";"beta.trend";"log.beta1";"log.beta2";"log.beta3";"log.beta4";"log.beta5";"log.beta6";"log.omega1";"log.omega2";"log.omega3";"log.omega4";"log.omega5";"log.omega6";"sd.beta";"tau";"S.0";"I.0";"Rs.0";"R1.0";"R2.0";"R3.0";"nbasis";"nrstage"
+"mle1";20,8;19,1;0;0,02;0,06;1;1;-0,00498;0,747;6,38;-3,44;4,23;3,33;4,55;-1,6928195214;-2,5433835795;-2,8404393891;-4,6918179927;-8,4779724783;-4,3900588064;3,13;0,23;0,621;0,378;0;0,000843;0,000972;1,16e-07;6;3
+"box_min";10;0,2;0;0,02;0,03;1;1;-0,01;-4;0;-4;0;0;0;-10;-10;-10;-10;-10;-10;1;0,1;0;0;0;0;0;0;6;3
+"box_max";40;30;0;0,02;0,6;1;1;0;4;8;4;8;8;8;0;0;0;0;0;0;5;0,5;1;1;0;1;1;1;6;3
+',
+ row.names=1
+ )
+
+dacca.hyperparams <- list(
+ 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)
+ r
+}
+
+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 <- continue(m7)
Modified: pkg/pomp/tests/dacca.Rout.save
===================================================================
--- pkg/pomp/tests/dacca.Rout.save 2014-02-27 02:20:24 UTC (rev 886)
+++ pkg/pomp/tests/dacca.Rout.save 2014-03-02 14:10:33 UTC (rev 887)
@@ -1,6 +1,6 @@
-R version 3.0.2 (2013-09-25) -- "Frisbee Sailing"
-Copyright (C) 2013 The R Foundation for Statistical Computing
+R Under development (unstable) (2014-02-22 r65060) -- "Unsuffered Consequences"
+Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -47,8 +47,52 @@
>
> pf <- pfilter(dacca,Np=1000,seed=5886855L)
> print(round(logLik(pf),digits=1))
-[1] -3747.5
+[1] -3748.6
>
+> ## to investigate the rogue crash:
+>
+> dacca.pars <- c("gamma","eps","deltaI","beta.trend",
++ "log.beta1","log.beta2","log.beta3",
++ "log.beta4","log.beta5","log.beta6",
++ "log.omega1","log.omega2","log.omega3",
++ "log.omega4","log.omega5","log.omega6",
++ "sd.beta","tau")
+> dacca.ivps <- c("S.0","I.0","R1.0","R2.0","R3.0")
+> dacca.rw.sd <- c(
++ rep(0.1,length(dacca.pars)),
++ rep(0.2,length(dacca.ivps))
++ )
+> names(dacca.rw.sd) <- c(dacca.pars,dacca.ivps)
+>
+> param.tab <- read.csv2(text='
++ "";"gamma";"eps";"rho";"delta";"deltaI";"clin";"alpha";"beta.trend";"log.beta1";"log.beta2";"log.beta3";"log.beta4";"log.beta5";"log.beta6";"log.omega1";"log.omega2";"log.omega3";"log.omega4";"log.omega5";"log.omega6";"sd.beta";"tau";"S.0";"I.0";"Rs.0";"R1.0";"R2.0";"R3.0";"nbasis";"nrstage"
++ "mle1";20,8;19,1;0;0,02;0,06;1;1;-0,00498;0,747;6,38;-3,44;4,23;3,33;4,55;-1,6928195214;-2,5433835795;-2,8404393891;-4,6918179927;-8,4779724783;-4,3900588064;3,13;0,23;0,621;0,378;0;0,000843;0,000972;1,16e-07;6;3
++ "box_min";10;0,2;0;0,02;0,03;1;1;-0,01;-4;0;-4;0;0;0;-10;-10;-10;-10;-10;-10;1;0,1;0;0;0;0;0;0;6;3
++ "box_max";40;30;0;0,02;0,6;1;1;0;4;8;4;8;8;8;0;0;0;0;0;0;5;0,5;1;1;0;1;1;1;6;3
++ ',
++ row.names=1
++ )
+>
+> dacca.hyperparams <- list(
++ 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)
++ r
++ }
+>
+> 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 <- continue(m7)
+>
> proc.time()
user system elapsed
- 3.492 0.048 3.555
+ 8.056 0.044 8.137
More information about the pomp-commits
mailing list