[Pomp-commits] r897 - in pkg/pomp: R inst inst/examples src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 20 19:31:26 CET 2014
Author: kingaa
Date: 2014-03-20 19:31:26 +0100 (Thu, 20 Mar 2014)
New Revision: 897
Modified:
pkg/pomp/R/minim.R
pkg/pomp/inst/NEWS
pkg/pomp/inst/examples/bbs.R
pkg/pomp/src/cholmodel.c
pkg/pomp/src/pomp.h
pkg/pomp/src/sir.c
pkg/pomp/tests/dacca.R
pkg/pomp/tests/dacca.Rout.save
pkg/pomp/tests/ou2-trajmatch.Rout.save
Log:
- simplify dacca test
- updates tests
- add barycentric transformation functions to pomp.h
- rewrite sir.c and cholmodel.c to use the barycentric transformations
- use negative binomial reporting model in bbs
Modified: pkg/pomp/R/minim.R
===================================================================
--- pkg/pomp/R/minim.R 2014-03-20 17:57:01 UTC (rev 896)
+++ pkg/pomp/R/minim.R 2014-03-20 18:31:26 UTC (rev 897)
@@ -22,7 +22,6 @@
if (length(est)==0) {
- params <- c(A=3)
val <- objfun(guess)
conv <- NA
evals <- as.integer(c(1,0))
Modified: pkg/pomp/inst/NEWS
===================================================================
--- pkg/pomp/inst/NEWS 2014-03-20 17:57:01 UTC (rev 896)
+++ pkg/pomp/inst/NEWS 2014-03-20 18:31:26 UTC (rev 897)
@@ -5,6 +5,9 @@
o 'pomp' now depends on 'nloptr', which provides a suite of optimization algorithms.
This package can now be used in various methods for optimization of an objective function.
+ o New inline C functions 'to_log_barycentric' and 'from_log_barycentric' are provided in 'pomp.h' to facilitate log-barycentric transformations.
+ These have proven very useful in dealing with parameters constrained to sum to one (e.g., initial conditions of compartmental models).
+
0.48-3
o Correct a bug in 'abc' to do with parameter transformation.
Modified: pkg/pomp/inst/examples/bbs.R
===================================================================
--- pkg/pomp/inst/examples/bbs.R 2014-03-20 17:57:01 UTC (rev 896)
+++ pkg/pomp/inst/examples/bbs.R 2014-03-20 18:31:26 UTC (rev 897)
@@ -37,7 +37,7 @@
),
skeleton.type="vectorfield",
skeleton="_sir_ODE",
- measurement.model=reports~norm(mean=rho*cases,sd=1+sigma*cases),
+ measurement.model=reports~nbinom(mu=rho*cases,size=1/sigma),
PACKAGE="pomp",
obsnames = c("reports"),
statenames=c("S","I","R","cases","W"),
Modified: pkg/pomp/src/cholmodel.c
===================================================================
--- pkg/pomp/src/cholmodel.c 2014-03-20 17:57:01 UTC (rev 896)
+++ pkg/pomp/src/cholmodel.c 2014-03-20 18:31:26 UTC (rev 897)
@@ -51,17 +51,12 @@
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]);
+ to_log_barycentric(&pt[parindex[14]],&S0,3+nrstage);
}
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);
@@ -72,20 +67,12 @@
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;
+ from_log_barycentric(&pt[parindex[14]],&S0,3+nrstage);
}
void _cholmodel_norm_rmeasure (double *y, double *x, double *p,
- int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int *obsindex, int *stateindex,
+ int *parindex, int *covindex,
int ncovars, double *covars, double t)
{
double v, tol = 1.0e-18;
@@ -97,8 +84,10 @@
}
}
-void _cholmodel_norm_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
- int *obsindex, int *stateindex, int *parindex, int *covindex,
+void _cholmodel_norm_dmeasure (double *lik, double *y, double *x,
+ double *p, int give_log,
+ int *obsindex, int *stateindex,
+ int *parindex, int *covindex,
int ncovars, double *covars, double t)
{
double v, tol = 1.0e-18;
@@ -119,7 +108,8 @@
// truncation is not used
// instead, particles with negative states are killed
void _cholmodel_one (double *x, const double *p,
- const int *stateindex, const int *parindex, const int *covindex,
+ const int *stateindex, const int *parindex,
+ const int *covindex,
int covdim, const double *covar,
double t, double dt)
{ // implementation of the SIRS cholera model
Modified: pkg/pomp/src/pomp.h
===================================================================
--- pkg/pomp/src/pomp.h 2014-03-20 17:57:01 UTC (rev 896)
+++ pkg/pomp/src/pomp.h 2014-03-20 18:31:26 UTC (rev 897)
@@ -227,6 +227,8 @@
return(trans);
}
+
+// LOGIT AND INVERSE LOGIT TRANSFORMATIONS
static R_INLINE double logit (double p) {
return log(p/(1.0-p));
}
@@ -242,7 +244,8 @@
// this must be done by the calling program
// But note that when reulermultinom is called inside a pomp 'rprocess', there is no need to call
// {Get,Put}RNGState() as this is handled by pomp
-static void reulermultinom (int m, double size, double *rate, double dt, double *trans) {
+static void reulermultinom (int m, double size, double *rate,
+ double dt, double *trans) {
double p = 0.0;
int j, k;
if ((size < 0.0) || (dt < 0.0) || (floor(size+0.5) != size)) {
@@ -276,7 +279,8 @@
}
// COMPUTE PROBABILITIES OF EULER-MULTINOMIAL TRANSITIONS
-static double deulermultinom (int m, double size, double *rate, double dt, double *trans, int give_log) {
+static double deulermultinom (int m, double size, double *rate, double dt,
+ double *trans, int give_log) {
double p = 0.0;
double n = 0.0;
double ff = 0.0;
@@ -318,11 +322,42 @@
return ff;
}
+// C-LEVEL DEFINITIONS OF LOG-BARYCENTRIC TRANSFORMATION.
+// USEFUL FOR WORKING WITH PARAMETERS CONSTRAINED TO SUM TO 1
+
+// TRANSFORMS TO LOG BARYCENTRIC COORDINATES
+// on input:
+// x = pointer to vector of parameters to be tranformed to
+// log barycentric coordinates,
+// n = length of vector.
+// on output:
+// xt = pointer to vector of log barycentric coordinates
+static R_INLINE void to_log_barycentric (double *xt, const double *x, int n) {
+ double sum;
+ int i;
+ for (i = 0, sum = 0.0; i < n; i++) sum += x[i];
+ for (i = 0; i < n; i++) xt[i] = log(x[i]/sum);
+}
+
+// TRANSFORMS FROM LOG BARYCENTRIC COORDINATES
+// on input:
+// x = pointer to vector of parameters in log barycentric coordinates,
+// n = length of vector.
+// on output:
+// xt = pointer to vector of coordinates on unit simplex
+static R_INLINE void from_log_barycentric (double *xt, const double *x, int n) {
+ double sum;
+ int i;
+ for (i = 0, sum = 0.0; i < n; i++) sum += (xt[i] = exp(x[i]));
+ for (i = 0; i < n; i++) xt[i] /= sum;
+}
+
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) {
+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);
@@ -334,7 +369,8 @@
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) {
+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);
Modified: pkg/pomp/src/sir.c
===================================================================
--- pkg/pomp/src/sir.c 2014-03-20 17:57:01 UTC (rev 896)
+++ pkg/pomp/src/sir.c 2014-03-20 18:31:26 UTC (rev 897)
@@ -52,15 +52,13 @@
pt[parindex[3]+k] = log(BETA[k]);
pt[parindex[4]] = log(BETA_SD);
pt[parindex[6]] = logit(RHO);
- pt[parindex[7]] = log(S0);
- pt[parindex[8]] = log(I0);
- pt[parindex[9]] = log(R0);
+ to_log_barycentric(pt+parindex[7],&S0,3);
+
}
void _sir_par_trans (double *pt, double *p, int *parindex)
{
int nbasis = *(get_pomp_userdata_int("nbasis"));
- double sum = 0.0;
int k;
pt[parindex[0]] = exp(GAMMA);
pt[parindex[1]] = exp(MU);
@@ -69,12 +67,7 @@
pt[parindex[3]+k] = exp(BETA[k]);
pt[parindex[4]] = exp(BETA_SD);
pt[parindex[6]] = expit(RHO);
- sum += (pt[parindex[7]] = exp(S0));
- sum += (pt[parindex[8]] = exp(I0));
- sum += (pt[parindex[9]] = exp(R0));
- pt[parindex[7]] /= sum;
- pt[parindex[8]] /= sum;
- pt[parindex[9]] /= sum;
+ from_log_barycentric(pt+parindex[7],&S0,3);
}
void _sir_binom_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
Modified: pkg/pomp/tests/dacca.R
===================================================================
--- pkg/pomp/tests/dacca.R 2014-03-20 17:57:01 UTC (rev 896)
+++ pkg/pomp/tests/dacca.R 2014-03-20 18:31:26 UTC (rev 897)
@@ -7,18 +7,10 @@
pompExample(dacca)
x <- as.data.frame(dacca)
- print(names(x))
- print(dim(x))
-
x <- simulate(dacca,nsim=3,as.data.frame=TRUE)
- print(names(x))
- print(dim(x))
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:
@@ -55,6 +47,8 @@
r
}
+ op <- options(warn=-1)
+
set.seed(7777+7)
params.tricky <- dacca.rprior(dacca.hyperparams)
m7 <- mif(
@@ -90,4 +84,6 @@
transform=TRUE
)
+ options(op)
+
}
Modified: pkg/pomp/tests/dacca.Rout.save
===================================================================
--- pkg/pomp/tests/dacca.Rout.save 2014-03-20 17:57:01 UTC (rev 896)
+++ pkg/pomp/tests/dacca.Rout.save 2014-03-20 18:31:26 UTC (rev 897)
@@ -24,18 +24,10 @@
+ pompExample(dacca)
+
+ x <- as.data.frame(dacca)
-+ print(names(x))
-+ print(dim(x))
-+
+ x <- simulate(dacca,nsim=3,as.data.frame=TRUE)
-+ print(names(x))
-+ print(dim(x))
+
+ 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:
+
@@ -110,22 +102,10 @@
+ }
Loading required package: mvtnorm
Loading required package: subplex
+Loading required package: nloptr
Loading required package: deSolve
newly created pomp object(s):
dacca
- [1] "time" "cholera.deaths" "seas.1" "seas.2"
- [5] "seas.3" "seas.4" "seas.5" "seas.6"
- [9] "pop" "dpopdt" "trend"
-[1] 600 11
- [1] "time" "cholera.deaths" "S" "I"
- [5] "Rs" "R1" "R2" "R3"
- [9] "M" "W" "count" "seas.1"
-[13] "seas.2" "seas.3" "seas.4" "seas.5"
-[17] "seas.6" "pop" "dpopdt" "trend"
-[21] "sim"
-[1] 1800 21
-[1] -3748.6
-[1] -3701.5
Warning messages:
1: 11 filtering failures occurred in 'pfilter'
2: 4 filtering failures occurred in 'pfilter'
@@ -138,4 +118,4 @@
>
> proc.time()
user system elapsed
- 11.056 0.060 11.164
+ 11.092 0.056 11.206
Modified: pkg/pomp/tests/ou2-trajmatch.Rout.save
===================================================================
--- pkg/pomp/tests/ou2-trajmatch.Rout.save 2014-03-20 17:57:01 UTC (rev 896)
+++ pkg/pomp/tests/ou2-trajmatch.Rout.save 2014-03-20 18:31:26 UTC (rev 897)
@@ -133,7 +133,7 @@
Error in tmof.internal(object = object, params = params, est = est, transform = transform, :
parameter(s): 'bob''harry' not found in 'params'
> try(traj.match(ou2,est=c('x1.0','x2.0','bob','alpha.4','tau')))
-Error in maxim.internal(objfun = traj.match.objfun(object = object, params = start, :
+Error in minim.internal(objfun = traj.match.objfun(object = object, params = start, :
'est' must refer to parameters named in 'start'
>
> fit <- optim(par=c(0,0,1,0.5,0.5),fn=ofun,method="Nelder-Mead",control=list(maxit=400))
@@ -177,4 +177,4 @@
>
> proc.time()
user system elapsed
- 1.588 0.084 1.700
+ 1.672 0.028 1.723
More information about the pomp-commits
mailing list