[Pomp-commits] r204 - in pkg: R data man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Mar 13 19:35:40 CET 2010
Author: kingaa
Date: 2010-03-13 19:35:40 +0100 (Sat, 13 Mar 2010)
New Revision: 204
Added:
pkg/data/euler.sir.rda
pkg/data/ou2.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/src/euler_sir_density.c
Removed:
pkg/data/euler.sir.R
pkg/data/ou2.R
pkg/data/rw2.R
pkg/data/verhulst.R
Modified:
pkg/R/pomp.R
pkg/man/mif.Rd
pkg/src/eulermultinom.c
pkg/tests/ou2-forecast.R
pkg/tests/sir.R
Log:
- move default initializer to pomp namespace
- go back to binary format data-loadable examples to avoid the ugly problem of tag-along environments
- speed up for degenerate cases in reulermultinom
- put euler-density function for SIR model into source (though it is not used at present)
Modified: pkg/R/pomp.R
===================================================================
--- pkg/R/pomp.R 2010-02-11 02:45:06 UTC (rev 203)
+++ pkg/R/pomp.R 2010-03-13 18:35:40 UTC (rev 204)
@@ -1,4 +1,13 @@
## constructor of the pomp class
+default.initializer <- function (params, t0, ...) {
+ ivpnames <- grep("\\.0$",names(params),val=TRUE)
+ if (length(ivpnames)<1)
+ stop("default initializer error: no parameter names ending in ",sQuote(".0")," found: see ",sQuote("pomp")," documentation")
+ x <- params[ivpnames]
+ names(x) <- sub("\\.0$","",ivpnames)
+ x
+}
+
pomp <- function (data, times, t0, ..., rprocess, dprocess,
rmeasure, dmeasure, measurement.model,
skeleton.map, skeleton.vectorfield, initializer, covar, tcovar,
@@ -114,14 +123,7 @@
}
if (missing(initializer)) {
- initializer <- function (params, t0, ...) {
- ivpnames <- grep("\\.0$",names(params),val=TRUE)
- if (length(ivpnames)<1)
- stop("default initializer error: no parameter names ending in ",sQuote(".0")," found: see ",sQuote("pomp")," documentation")
- x <- params[ivpnames]
- names(x) <- sub("\\.0$","",ivpnames)
- x
- }
+ initializer <- default.initializer
}
if (missing(PACKAGE)) PACKAGE <- character(0)
Deleted: pkg/data/euler.sir.R
===================================================================
--- pkg/data/euler.sir.R 2010-02-11 02:45:06 UTC (rev 203)
+++ pkg/data/euler.sir.R 2010-03-13 18:35:40 UTC (rev 204)
@@ -1,47 +0,0 @@
-require(pomp)
-
-local(
- {
- po <- pomp(
- times=seq(1/52,4,by=1/52),
- data=rbind(measles=numeric(52*4)),
- t0=0,
- delta.t=1/52/20,
- statenames=c("S","I","R","cases","W"),
- paramnames=c(
- "gamma","mu","iota",
- "beta1","beta.sd","pop","rho",
- "nbasis","degree","period"
- ),
- zeronames=c("cases"),
- comp.names=c("S","I","R"),
- step.fun="sir_euler_simulator",
- rprocess=euler.simulate,
- skeleton.vectorfield="sir_ODE",
- rmeasure="binom_rmeasure",
- dmeasure="binom_dmeasure",
- PACKAGE="pomp",
- initializer=function(params, t0, comp.names, ...){
- p <- exp(params)
- snames <- c("S","I","R","cases","W")
- fracs <- p[paste(comp.names,"0",sep=".")]
- x0 <- numeric(length(snames))
- names(x0) <- snames
- x0[comp.names] <- round(p['pop']*fracs/sum(fracs))
- x0
- }
- )
-
- coef(po) <- c(
- gamma=log(26),mu=log(0.02),iota=log(0.01),
- nbasis=3,degree=3,period=1, # NB: all parameters are log-transformed but these
- beta1=log(1200),beta2=log(1800),beta3=log(600),
- beta.sd=log(1e-3),
- pop=log(2.1e6),
- rho=log(0.6),
- S.0=log(26/1200),I.0=log(0.001),R.0=log(1-0.001-26/1200)
- )
-
- simulate(po,nsim=1,seed=329348545L)
- }
- ) -> euler.sir
Added: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Property changes on: pkg/data/euler.sir.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Deleted: pkg/data/ou2.R
===================================================================
--- pkg/data/ou2.R 2010-02-11 02:45:06 UTC (rev 203)
+++ pkg/data/ou2.R 2010-03-13 18:35:40 UTC (rev 204)
@@ -1,76 +0,0 @@
-require(pomp)
-
-local(
- {
- po <- pomp(
- times=seq(1,100),
- data=rbind(
- y1=rep(0,100),
- y2=rep(0,100)
- ),
- t0=0,
- rprocess = function (xstart, times, params, paramnames, ...) {
- nvar <- nrow(xstart)
- npar <- nrow(params)
- nrep <- ncol(xstart)
- ntimes <- length(times)
- ## get indices of the various parameters in the 'params' matrix
- ## C uses zero-based indexing!
- parindex <- match(paramnames,rownames(params))-1
- array(
- .C("ou2_adv",
- X = double(nvar*nrep*ntimes),
- xstart = as.double(xstart),
- par = as.double(params),
- times = as.double(times),
- n = as.integer(c(nvar,npar,nrep,ntimes)),
- parindex = as.integer(parindex),
- DUP = FALSE,
- NAOK = TRUE,
- PACKAGE = "pomp"
- )$X,
- dim=c(nvar,nrep,ntimes),
- dimnames=list(rownames(xstart),NULL,NULL)
- )
- },
- dprocess = function (x, times, params, log, paramnames, ...) {
- nvar <- nrow(x)
- npar <- nrow(params)
- nrep <- ncol(x)
- ntimes <- length(times)
- parindex <- match(paramnames,rownames(params))-1
- array(
- .C("ou2_pdf",
- d = double(nrep*(ntimes-1)),
- X = as.double(x),
- par = as.double(params),
- times = as.double(times),
- n = as.integer(c(nvar,npar,nrep,ntimes)),
- parindex = as.integer(parindex),
- give_log=as.integer(log),
- DUP = FALSE,
- NAOK = TRUE,
- PACKAGE = "pomp"
- )$d,
- dim=c(nrep,ntimes-1)
- )
- },
- dmeasure = "normal_dmeasure",
- rmeasure = "normal_rmeasure",
- paramnames = c(
- "alpha.1","alpha.2","alpha.3","alpha.4",
- "sigma.1","sigma.2","sigma.3",
- "tau"
- ),
- statenames = c("x1","x2")
- )
-
- coef(po) <- c(
- alpha.1=0.9,alpha.2=0,alpha.3=0,alpha.4=0.99,
- sigma.1=1,sigma.2=0,sigma.3=2,
- tau=1,x1.0=50,x2.0=-50
- )
-
- simulate(po,nsim=1,seed=377456545L)
- }
- ) -> ou2
Added: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Property changes on: pkg/data/ou2.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Deleted: pkg/data/rw2.R
===================================================================
--- pkg/data/rw2.R 2010-02-11 02:45:06 UTC (rev 203)
+++ pkg/data/rw2.R 2010-03-13 18:35:40 UTC (rev 204)
@@ -1,39 +0,0 @@
-require(pomp)
-
-local(
- {
- po <- pomp(
- rprocess = onestep.simulate,
- dprocess = onestep.density,
- step.fun = function(x, t, params, delta.t, ...) {
- c(
- x1=rnorm(n=1,mean=x['x1'],sd=params['s1']*delta.t),
- x2=rnorm(n=1,mean=x['x2'],sd=params['s2']*delta.t)
- )
- },
- dens.fun = function (x1, t1, x2, t2, params, ...) {
- sum(
- dnorm(
- x=x2[c('x1','x2')],
- mean=x1[c('x1','x2')],
- sd=params[c('s1','s2')]*(t2-t1),
- log=TRUE
- ),
- na.rm=TRUE
- )
- },
- measurement.model=list(
- y1 ~ norm(mean=x1,sd=tau),
- y2 ~ norm(mean=x2,sd=tau)
- ),
- times=1:100,
- data=rbind(
- y1=rep(0,100),
- y2=rep(0,100)
- ),
- t0=0
- )
-
- simulate(po,params=c(x1.0=0,x2.0=0,s1=1,s2=3,tau=1),nsim=1,seed=738377475L)
- }
- ) -> rw2
Added: pkg/data/rw2.rda
===================================================================
(Binary files differ)
Property changes on: pkg/data/rw2.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Deleted: pkg/data/verhulst.R
===================================================================
--- pkg/data/verhulst.R 2010-02-11 02:45:06 UTC (rev 203)
+++ pkg/data/verhulst.R 2010-03-13 18:35:40 UTC (rev 204)
@@ -1,54 +0,0 @@
-require(pomp)
-
-local(
- {
- po <- pomp(
- data=rbind(obs=rep(0,1000)),
- times=seq(0.1,by=0.1,length=1000),
- t0=0,
- rprocess=euler.simulate,
- step.fun=function(x,t,params,delta.t,...){
- with(
- as.list(c(x,params)),
- rnorm(
- n=1,
- mean=n+r*n*(1-n/K)*delta.t,
- sd=sigma*n*sqrt(delta.t)
- )
- )
- },
- dprocess=onestep.density,
- dens.fun=function(x1,x2,t1,t2,params,log,...){
- delta.t <- t2-t1
- with(
- as.list(c(x1,params)),
- dnorm(
- x=x2['n'],
- mean=n+r*n*(1-n/K)*delta.t,
- sd=sigma*n*sqrt(delta.t),
- log=log
- )
- )
- },
- measurement.model=obs~lnorm(meanlog=log(n),sdlog=log(1+tau)),
- skeleton.vectorfield=function(x,t,params,...){
- with(
- as.list(c(x,params)),
- r*n*(1-n/K)
- )
- },
- delta.t=0.01
- )
- simulate(
- po,
- params=c(
- n.0=10000,
- K=10000,
- r=0.9,
- sigma=0.4,
- tau=0.1
- ),
- seed=73658676L
- )
- }
- ) -> verhulst
Added: pkg/data/verhulst.rda
===================================================================
(Binary files differ)
Property changes on: pkg/data/verhulst.rda
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Modified: pkg/man/mif.Rd
===================================================================
--- pkg/man/mif.Rd 2010-02-11 02:45:06 UTC (rev 203)
+++ pkg/man/mif.Rd 2010-03-13 18:35:40 UTC (rev 204)
@@ -29,8 +29,8 @@
The number of MIF iterations to perform.
}
\item{start}{
- The initial guess of the parameters.
- This must be a named vector.
+ named numerical vector;
+ the starting guess of the parameters.
}
\item{pars}{
optional character vector naming the ordinary parameters to be estimated.
@@ -42,7 +42,7 @@
Every parameter named in \code{ivps} must have a positive random-walk standard deviation specified in \code{rw.sd}.
}
\item{particles}{
- Function of prototype \code{particles(Np,center,sd,...)} which sets up the initial particle matrix by drawing a sample of size \code{Np} from the initial particle distribution centered at \code{center} and of width \code{sd}.
+ Function of prototype \code{particles(Np,center,sd,...)} which sets up the starting particle matrix by drawing a sample of size \code{Np} from the starting particle distribution centered at \code{center} and of width \code{sd}.
If \code{particles} is not supplied by the user, the default behavior is to draw the particles from a multivariate normal distribution with mean \code{center} and standard deviation \code{sd}.
}
\item{rw.sd}{
@@ -66,8 +66,8 @@
}
\item{var.factor}{
a positive number;
- the scaling coefficient relating the width of the initial particle distribution to \code{rw.sd}
- The width of the initial distribution of particles will be \code{random.walk.sd*var.factor}.
+ the scaling coefficient relating the width of the starting particle distribution to \code{rw.sd}.
+ In particular, the width of the distribution of particles at the start of the first MIF iteration will be \code{random.walk.sd*var.factor}.
}
\item{cooling.factor}{
a positive number not greater than 1;
@@ -79,7 +79,7 @@
instead, an unweighed average of the filtering means is used for the update.
}
\item{tol}{
- numeric scalar; particles with log likelihood below \code{tol} are considered to be "lost".
+ numeric scalar; particles with log likelihood below \code{tol} are considered to be \dQuote{lost}.
A filtering failure occurs when, at some time point, all particles are lost.
}
\item{max.fail}{
@@ -120,7 +120,7 @@
The center of the particle distribution returned by \code{particles} should be \code{center}.
The width of the particle distribution should vary monotonically with \code{sd}.
- In particular, when \code{sd=0}, the \code{particles} should return matrices with \code{Np} identical columns, each corresponding to the parameters specified in \code{center}.
+ In particular, when \code{sd=0}, the \code{particles} should return matrices with \code{Np} identical columns, each given by the parameters specified in \code{center}.
}
\references{
E. L. Ionides, C. Bret\'o, & A. A. King,
Added: pkg/src/euler_sir_density.c
===================================================================
--- pkg/src/euler_sir_density.c (rev 0)
+++ pkg/src/euler_sir_density.c 2010-03-13 18:35:40 UTC (rev 204)
@@ -0,0 +1,129 @@
+// dear emacs, please treat this as -*- C++ -*-
+
+#include <Rmath.h>
+
+#include "pomp.h"
+
+static double expit (double x) {
+ return 1.0/(1.0 + exp(-x));
+}
+
+static double logit (double x) {
+ return log(x/(1-x));
+}
+
+static double term_time (double t, double b0, double b1)
+{
+ static double correction = 0.52328767123287671233;
+ double day = 365.0 * (t - floor(t));
+ int tt = (day >= 7.0 && day <= 100.0)
+ || (day >= 115.0 && day <= 200.0)
+ || (day >= 252.0 && day <= 300.0)
+ || (day >= 308.0 && day <= 356.0);
+ return b0 * (1.0 + b1 * (-1.0 + 2.0 * ((double) tt))) / (1.0 + correction * b1);
+}
+
+#define LOGGAMMA (p[parindex[0]]) // recovery rate
+#define LOGMU (p[parindex[1]]) // baseline birth and death rate
+#define LOGIOTA (p[parindex[2]]) // import rate
+#define LOGBETA (p[parindex[3]]) // transmission rate
+#define LOGBETA_SD (p[parindex[4]]) // environmental stochasticity SD in transmission rate
+#define LOGPOPSIZE (p[parindex[5]]) // population size
+#define LOGRHO (p[parindex[6]]) // reporting probability
+#define NBASIS (p[parindex[7]]) // number of periodic B-spline basis functions
+#define DEGREE (p[parindex[8]]) // degree of periodic B-spline basis functions
+#define PERIOD (p[parindex[9]]) // period of B-spline basis functions
+
+#define SUSC (x1[stateindex[0]]) // number of susceptibles
+#define INFD (x1[stateindex[1]]) // number of infectives
+#define RCVD (x1[stateindex[2]]) // number of recovereds
+#define CASE (x1[stateindex[3]]) // number of cases (accumulated per reporting period)
+#define W (x1[stateindex[4]]) // integrated white noise
+#define BIRTHS (x2[stateindex[5]]) // births
+#define dW (x2[stateindex[6]]) // white noise process
+
+// SIR model with Euler multinomial step
+// forced transmission (basis functions passed as covariates)
+// constant population size as a parameter
+// environmental stochasticity on transmission
+void sir_euler_density (double *f, double *x1, double *x2, double t1, double t2, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar)
+{
+ int nrate = 6;
+ double rate[nrate]; // transition rates
+ double *trans; // transition numbers
+ double gamma, mu, iota, beta_sd, popsize;
+ double beta;
+ double dt = t2-t1;
+ int nseas = (int) NBASIS; // number of seasonal basis functions
+ int deg = (int) DEGREE; // degree of seasonal basis functions
+ double seasonality[nseas];
+
+ // untransform the parameters
+ gamma = exp(LOGGAMMA);
+ mu = exp(LOGMU);
+ iota = exp(LOGIOTA);
+ beta_sd = exp(LOGBETA_SD);
+ popsize = exp(LOGPOPSIZE);
+
+ periodic_bspline_basis_eval(t1,PERIOD,deg,nseas,&seasonality[0]);
+ beta = exp(dot_product(nseas,&seasonality[0],&LOGBETA));
+
+ // test to make sure the parameters and state variable values are sane
+ if (!(R_FINITE(beta)) ||
+ !(R_FINITE(gamma)) ||
+ !(R_FINITE(mu)) ||
+ !(R_FINITE(beta_sd)) ||
+ !(R_FINITE(iota)) ||
+ !(R_FINITE(popsize)) ||
+ !(R_FINITE(SUSC)) ||
+ !(R_FINITE(INFD)) ||
+ !(R_FINITE(RCVD)) ||
+ !(R_FINITE(CASE)) ||
+ !(R_FINITE(W)) ||
+ (nseas <= 0)) {
+ *f = R_NaN;
+ return;
+ }
+
+ // compute the transition rates
+ trans = &BIRTHS;
+ if (beta_sd > 0.0) { // environmental noise is ON
+ double beta_var = beta_sd*beta_sd;
+ *f = dgamma(dW,dt/beta_var,beta_var,1);
+ } else { // environmental noise is OFF
+ *f = 0; // THIS ASSUMES THAT dw = dt !!!
+ }
+ rate[0] = mu*popsize; // birth into susceptible class
+ rate[1] = (iota+beta*INFD*dW/dt)/popsize; // force of infection
+ rate[2] = mu; // death from susceptible class
+ rate[3] = gamma; // recovery
+ rate[4] = mu; // death from infectious class
+ rate[5] = mu; // death from recovered class
+
+ // compute the transition numbers
+ *f += dpois(trans[0],rate[0]*dt,1); // births are Poisson
+ *f += deulermultinom(2,SUSC,&rate[1],dt,&trans[1],1);
+ *f += deulermultinom(2,INFD,&rate[3],dt,&trans[3],1);
+ *f += deulermultinom(1,RCVD,&rate[5],dt,&trans[5],1);
+}
+
+#undef SUSC
+#undef INFD
+#undef RCVD
+#undef CASE
+#undef W
+#undef BIRTHS
+#undef dW
+
+#undef LOGGAMMA
+#undef LOGMU
+#undef LOGIOTA
+#undef LOGBETA
+#undef LOGBETA_SD
+#undef LOGPOPSIZE
+#undef LOGRHO
+#undef NBASIS
+#undef DEGREE
+#undef PERIOD
Modified: pkg/src/eulermultinom.c
===================================================================
--- pkg/src/eulermultinom.c 2010-02-11 02:45:06 UTC (rev 203)
+++ pkg/src/eulermultinom.c 2010-03-13 18:35:40 UTC (rev 204)
@@ -9,31 +9,33 @@
double p = 0.0;
int j, k;
if ((size < 0.0) || (dt < 0.0) || (floor(size+0.5) != size)) {
- for (k = 0; k < m; k++)
- trans[k] = R_NaN;
+ for (k = 0; k < m; k++) trans[k] = R_NaN;
return;
}
for (k = 0; k < m; k++) {
if (rate[k] < 0.0) {
- for (j = 0; j < m; j++)
- trans[j] = R_NaN;
+ for (j = 0; j < m; j++) trans[j] = R_NaN;
return;
}
p += rate[k]; // total event rate
}
- size = rbinom(size,1-exp(-p*dt)); // total number of events
- if (!(R_FINITE(size)))
- warning("reulermultinom: result of binomial draw is not finite");
- m -= 1;
- for (k = 0; k < m; k++) {
- if (rate[k] > p) p = rate[k];
- trans[k] = ((size > 0) && (p > 0)) ? rbinom(size,rate[k]/p) : 0;
- if (!(R_FINITE(size)&&R_FINITE(p)&&R_FINITE(rate[k])&&R_FINITE(trans[k])))
+ if (p > 0.0) {
+ size = rbinom(size,1-exp(-p*dt)); // total number of events
+ if (!(R_FINITE(size)))
warning("reulermultinom: result of binomial draw is not finite");
- size -= trans[k];
- p -= rate[k];
+ m -= 1;
+ for (k = 0; k < m; k++) {
+ if (rate[k] > p) p = rate[k];
+ trans[k] = ((size > 0) && (p > 0)) ? rbinom(size,rate[k]/p) : 0;
+ if (!(R_FINITE(size)&&R_FINITE(p)&&R_FINITE(rate[k])&&R_FINITE(trans[k])))
+ warning("reulermultinom: result of binomial draw is not finite");
+ size -= trans[k];
+ p -= rate[k];
+ }
+ trans[m] = size;
+ } else {
+ for (k = 0; k < m; k++) trans[k] = 0.0;
}
- trans[m] = size;
}
// probability density of Euler-multinomial transitions
Modified: pkg/tests/ou2-forecast.R
===================================================================
--- pkg/tests/ou2-forecast.R 2010-02-11 02:45:06 UTC (rev 203)
+++ pkg/tests/ou2-forecast.R 2010-03-13 18:35:40 UTC (rev 204)
@@ -2,6 +2,7 @@
set.seed(921625222L)
+data(ou2)
pf <- pfilter(ou2,Np=1000,save.states=TRUE)
ll <- cumsum(pf$cond.loglik)
pp <- matrix(data=coef(ou2),nrow=length(coef(ou2)),ncol=1000,dimnames=list(names(coef(ou2)),NULL))
Modified: pkg/tests/sir.R
===================================================================
--- pkg/tests/sir.R 2010-02-11 02:45:06 UTC (rev 203)
+++ pkg/tests/sir.R 2010-03-13 18:35:40 UTC (rev 204)
@@ -187,7 +187,7 @@
print(h1[c("S","I","R"),,],digits=4)
## now repeat using the compiled native codes built into the package
-
+data(euler.sir)
po <- euler.sir
set.seed(3049953)
More information about the pomp-commits
mailing list