[Pomp-commits] r532 - in pkg: data inst/data-R inst/examples man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Aug 8 18:27:24 CEST 2011
Author: kingaa
Date: 2011-08-08 18:27:24 +0200 (Mon, 08 Aug 2011)
New Revision: 532
Modified:
pkg/data/gillespie.sir.rda
pkg/inst/data-R/gillespie.sir.R
pkg/inst/examples/gompertz.R
pkg/inst/examples/sir.R
pkg/inst/examples/sir.c
pkg/man/gompertz.Rd
Log:
- include parameter transformations in 'gillespie.sir'
- bring examples into line with codes in vignettes
Modified: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/data-R/gillespie.sir.R
===================================================================
--- pkg/inst/data-R/gillespie.sir.R 2011-08-08 16:25:54 UTC (rev 531)
+++ pkg/inst/data-R/gillespie.sir.R 2011-08-08 16:27:24 UTC (rev 532)
@@ -1,74 +1,89 @@
library(pomp)
-params <- c(
- nu=log(1/70),
- mu=log(1/70),
- beta1=log(330),
- beta2=log(410),
- beta3=log(490),
- gamma=log(24),
- iota=log(0.1),
- rho=log(0.1),
- S.0=log(0.05),
- I.0=log(1e-4),
- R.0=log(0.95),
- N.0=1000000,
- cases.0=0,
- nbasis=3,
- degree=3,
- period=1
- )
+po <- pomp(
+ data=data.frame(
+ time=seq(from=0,to=10,by=1/52),
+ reports=NA
+ ),
+ times="time",
+ t0=0,
+ rprocess=gillespie.sim(
+ rate.fun="_sir_rates",
+ PACKAGE="pomp",
+ v=cbind(
+ birth=c(1,0,0,1,0),
+ sdeath=c(-1,0,0,-1,0),
+ infection=c(-1,1,0,0,0),
+ ideath=c(0,-1,0,-1,0),
+ recovery=c(0,-1,1,0,1),
+ rdeath=c(0,0,-1,-1,0)
+ ),
+ d=cbind(
+ birth=c(0,0,0,1,0),
+ sdeath=c(1,0,0,0,0),
+ infection=c(1,1,0,1,0),
+ ideath=c(0,1,0,0,0),
+ recovery=c(0,1,0,0,0),
+ rdeath=c(0,0,1,0,0)
+ )
+ ),
+ obsnames="reports",
+ statenames=c("S","I","R","N","cases"),
+ paramnames=c(
+ "gamma","mu","iota",
+ "beta1","nu",
+ "nbasis","degree","period"
+ ),
+ zeronames=c("cases"),
+ measurement.model=reports~binom(size=cases,prob=0.1),
+ to.log.transform=c(
+ "gamma","nu","mu","iota",
+ "beta1","beta2","beta3",
+ "rho",
+ "S.0","I.0","R.0"
+ ),
+ parameter.transform=function (params, to.log.transform, ...) {
+ params[to.log.transform] <- log(params[to.log.transform])
+ params
+ },
+ parameter.inv.transform=function (params, to.log.transform, ...) {
+ params[to.log.transform] <- exp(params[to.log.transform])
+ params
+ },
+ initializer=function(params, t0, ...){
+ comp.names <- c("S","I","R")
+ icnames <- paste(comp.names,"0",sep=".")
+ snames <- c("S","I","R","N","cases")
+ fracs <- exp(params[icnames])
+ x0 <- numeric(length(snames))
+ names(x0) <- snames
+ x0["N"] <- params["pop"]
+ x0[comp.names] <- round(params["pop"]*fracs/sum(fracs))
+ x0
+ }
+ )
+coef(po,transform=TRUE) <- c(
+ gamma=24,
+ nu=1/70,
+ mu=1/70,
+ iota=0.1,
+ beta1=330,
+ beta2=410,
+ beta3=490,
+ rho=0.1,
+ S.0=0.05,
+ I.0=1e-4,
+ R.0=0.95,
+ pop=1000000,
+ nbasis=3,
+ degree=3,
+ period=1
+ )
+
+
simulate(
- pomp(
- data=data.frame(
- time=seq(from=0,to=10,by=1/52),
- reports=NA
- ),
- times="time",
- t0=0,
- rprocess=gillespie.sim(
- rate.fun="_sir_rates",
- PACKAGE="pomp",
- v=cbind(
- birth=c(1,0,0,1,0),
- sdeath=c(-1,0,0,-1,0),
- infection=c(-1,1,0,0,0),
- ideath=c(0,-1,0,-1,0),
- recovery=c(0,-1,1,0,1),
- rdeath=c(0,0,-1,-1,0)
- ),
- d=cbind(
- birth=c(0,0,0,1,0),
- sdeath=c(1,0,0,0,0),
- infection=c(1,1,0,1,0),
- ideath=c(0,1,0,0,0),
- recovery=c(0,1,0,0,0),
- rdeath=c(0,0,1,0,0)
- )
- ),
- obsnames="reports",
- statenames=c("S","I","R","N","cases"),
- paramnames=c(
- "gamma","mu","iota",
- "beta1","nu",
- "nbasis","degree","period"
- ),
- zeronames=c("cases"),
- measurement.model=reports~binom(size=cases,prob=0.1),
- initializer=function(params, t0, ...){
- comp.names <- c("S","I","R")
- icnames <- paste(comp.names,"0",sep=".")
- snames <- c("S","I","R","N","cases")
- fracs <- exp(params[icnames])
- x0 <- numeric(length(snames))
- names(x0) <- snames
- x0["N"] <- params["N.0"]
- x0[comp.names] <- round(params['N.0']*fracs/sum(fracs))
- x0
- }
- ),
- params=params,
+ po,
nsim=1,
seed=1165270654L
) -> gillespie.sir
Modified: pkg/inst/examples/gompertz.R
===================================================================
--- pkg/inst/examples/gompertz.R 2011-08-08 16:25:54 UTC (rev 531)
+++ pkg/inst/examples/gompertz.R 2011-08-08 16:27:24 UTC (rev 532)
@@ -42,10 +42,19 @@
## compute the likelihood of Y|X,tau
f <- dlnorm(x=Y,meanlog=log(X),sdlog=tau,log=log)
return(f)
+ },
+ parameter.transform=function(params,...){
+ params <- c(params["X.0"],log(params[c("r","K","tau","sigma")]))
+ names(params) <- c("X.0","log.r","log.K","log.tau","log.sigma")
+ params
+ },
+ parameter.inv.transform=function(params,...){
+ params <- c(params["X.0"],exp(params[c("log.r","log.K","log.tau","log.sigma")]))
+ names(params) <- c("X.0","r","K","tau","sigma")
+ params
}
) -> gompertz
-
## Now code up the Gompertz example using native routines results in much faster computations.
## The C codes are included in the "examples" directory (file "gompertz.c")
@@ -64,16 +73,21 @@
PACKAGE="gompertz",
paramnames=c("log.r","log.K","log.sigma","log.tau"),
statenames=c("X"),
- obsnames=c("Y")
+ obsnames=c("Y"),
+ parameter.transform=function(params,...){
+ params <- c(params["X.0"],log(params[c("r","K","tau","sigma")]))
+ names(params) <- c("X.0","log.r","log.K","log.tau","log.sigma")
+ params
+ },
+ parameter.inv.transform=function(params,...){
+ params <- c(params["X.0"],exp(params[c("log.r","log.K","log.tau","log.sigma")]))
+ names(params) <- c("X.0","r","K","tau","sigma")
+ params
+ }
)
-params <- c(
- log.K=log(1),
- log.r=log(0.1),
- log.sigma=log(0.1),
- log.tau=log(0.1),
- X.0=1
- )
+## set the parameters
+coef(po,transform=TRUE) <- c(K=1,r=0.1,sigma=0.1,tau=0.1,X.0=1)
if (Sys.info()['sysname']=='Linux') { # only run this under linux
@@ -96,13 +110,13 @@
## compute a trajectory of the deterministic skeleton
tic <- Sys.time()
- X <- trajectory(po,params=params)
+ X <- trajectory(po)
toc <- Sys.time()
print(toc-tic)
## simulate from the model
tic <- Sys.time()
- x <- simulate(po,params=params,nsim=3)
+ x <- simulate(po,nsim=3)
toc <- Sys.time()
print(toc-tic)
Modified: pkg/inst/examples/sir.R
===================================================================
--- pkg/inst/examples/sir.R 2011-08-08 16:25:54 UTC (rev 531)
+++ pkg/inst/examples/sir.R 2011-08-08 16:27:24 UTC (rev 532)
@@ -6,16 +6,6 @@
if (Sys.info()['sysname']=='Linux') { # only run this under linux
- params <- c(
- gamma=26,mu=0.02,iota=0.01,
- beta1=1200,beta2=1800,beta3=600,
- beta.sd=1e-3,
- pop=2.1e6,
- rho=0.6,
- S.0=26/1200,I.0=0.001,R.0=1-0.001-26/1200
- )
-
-
model <- "sir"
pkg <- "pomp"
modelfile <- paste(model,".c",sep="")
@@ -50,37 +40,63 @@
## the order of the state variables assumed in the native routines:
statenames=c("S","I","R","cases","W"),
## the order of the parameters assumed in the native routines:
- paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho","nbasis","degree","period"),
+ paramnames=c(
+ "gamma","mu","iota","beta1","beta.sd","pop","rho",
+ "nbasis","degree","period"
+ ),
## reset cases to zero after each new observation:
zeronames=c("cases"),
- initializer=function(params,t0,...){
- p <- exp(params)
- with(
- as.list(p),
- {
- fracs <- c(S.0,I.0,R.0)
- x0 <- c(
- round(pop*fracs/sum(fracs)), # make sure the three compartments sum to 'pop' initially
- rep(0,2) # zeros for 'cases' and 'W'
- )
- names(x0) <- c("S","I","R","cases","W")
- x0
- }
- )
+ to.log.transform=c(
+ "gamma","mu","iota",
+ "beta1","beta2","beta3","beta.sd",
+ "rho",
+ "S.0","I.0","R.0"
+ ),
+ parameter.transform=function (params, to.log.transform, ...) {
+ params[to.log.transform] <- log(params[to.log.transform])
+ params
+ },
+ parameter.inv.transform=function (params, to.log.transform, ...) {
+ params[to.log.transform] <- exp(params[to.log.transform])
+ params
+ },
+ initializer=function(params, t0, ...) {
+ comp.names <- c("S","I","R")
+ ic.names <- c("S.0","I.0","R.0")
+ snames <- c("S","I","R","cases","W")
+ fracs <- exp(params[ic.names])
+ x0 <- numeric(length(snames))
+ names(x0) <- snames
+ x0[comp.names] <- round(params['pop']*fracs/sum(fracs))
+ ## since 'cases' is in 'zeronames' above, the IC for this variable
+ ## will only matter in trajectory computations
+ ## In trajectory computations, however, 'cases' will be roughly the weekly number of new cases
+ x0["cases"] <- x0["I"]*exp(params["gamma"])/52
+ x0
}
)
+ coef(po,transform=TRUE) <- c(
+ gamma=26,mu=0.02,iota=0.01,
+ beta1=1200,beta2=1800,beta3=600,
+ beta.sd=1e-3,
+ pop=2.1e6,
+ rho=0.6,
+ S.0=26/1200,I.0=0.001,R.0=1-0.001-26/1200,
+ nbasis=3,degree=3,period=1
+ )
+
dyn.load(solib) ## load the shared-object library
## compute a trajectory of the deterministic skeleton
tic <- Sys.time()
- X <- trajectory(po,c(log(params),nbasis=3,degree=3,period=1),hmax=1/52)
+ X <- trajectory(po,hmax=1/52)
toc <- Sys.time()
print(toc-tic)
## simulate from the model
tic <- Sys.time()
- x <- simulate(po,params=c(log(params),nbasis=3,degree=3,period=1),nsim=3)
+ x <- simulate(po,nsim=3)
toc <- Sys.time()
print(toc-tic)
Modified: pkg/inst/examples/sir.c
===================================================================
--- pkg/inst/examples/sir.c 2011-08-08 16:25:54 UTC (rev 531)
+++ pkg/inst/examples/sir.c 2011-08-08 16:27:24 UTC (rev 532)
@@ -23,7 +23,7 @@
#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 POPSIZE (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
@@ -67,7 +67,7 @@
int nrate = 6;
double rate[nrate]; // transition rates
double trans[nrate]; // transition numbers
- double gamma, mu, iota, beta_sd, beta_var, popsize;
+ double gamma, mu, iota, beta_sd, beta_var;
double beta; // transmission rate
double dW; // white noise increment
int nseas = (int) NBASIS; // number of seasonal basis functions
@@ -82,7 +82,6 @@
mu = exp(LOGMU);
iota = exp(LOGIOTA);
beta_sd = exp(LOGBETA_SD);
- popsize = exp(LOGPOPSIZE);
beta_var = beta_sd*beta_sd;
// to evaluate the basis functions and compute the transmission rate, use some of
@@ -104,7 +103,7 @@
!(R_FINITE(mu)) ||
!(R_FINITE(beta_sd)) ||
!(R_FINITE(iota)) ||
- !(R_FINITE(popsize)) ||
+ !(R_FINITE(POPSIZE)) ||
!(R_FINITE(SUSC)) ||
!(R_FINITE(INFD)) ||
!(R_FINITE(RCVD)) ||
@@ -121,8 +120,8 @@
}
// compute the transition rates
- rate[0] = mu*popsize; // birth into susceptible class
- rate[1] = (iota+beta*INFD*dW/dt)/popsize; // force of infection
+ 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
@@ -157,7 +156,7 @@
int nrate = 6;
double rate[nrate]; // transition rates
double term[nrate]; // terms in the equations
- double gamma, mu, iota, popsize;
+ double gamma, mu, iota;
double beta;
int nseas = (int) NBASIS; // number of seasonal basis functions
int deg = (int) DEGREE; // degree of seasonal basis functions
@@ -169,7 +168,6 @@
gamma = exp(LOGGAMMA);
mu = exp(LOGMU);
iota = exp(LOGIOTA);
- popsize = exp(LOGPOPSIZE);
// to evaluate the basis functions and compute the transmission rate, use some of
// pomp's built-in C-level facilities:
@@ -182,8 +180,8 @@
beta = exp((*dot)(nseas,&seasonality[0],&LOGBETA));
// compute the transition rates
- rate[0] = mu*popsize; // birth into susceptible class
- rate[1] = (iota+beta*INFD)/popsize; // force of infection
+ rate[0] = mu*POPSIZE; // birth into susceptible class
+ rate[1] = (iota+beta*INFD)/POPSIZE; // force of infection
rate[2] = mu; // death from susceptible class
rate[3] = gamma; // recovery
rate[4] = mu; // death from infectious class
@@ -224,7 +222,7 @@
#undef LOGIOTA
#undef LOGBETA
#undef LOGBETA_SD
-#undef LOGPOPSIZE
+#undef POPSIZE
#undef LOGRHO
#undef NBASIS
#undef DEGREE
Modified: pkg/man/gompertz.Rd
===================================================================
--- pkg/man/gompertz.Rd 2011-08-08 16:25:54 UTC (rev 531)
+++ pkg/man/gompertz.Rd 2011-08-08 16:27:24 UTC (rev 532)
@@ -9,13 +9,16 @@
\details{
The state process is \eqn{X_{t+1} = K^{(1-S)} X_{t}^S \varepsilon_{t}}, where \eqn{S=e^{-r}} and the \eqn{\varepsilon_t} are i.i.d. lognormal random deviates with variance \eqn{\sigma^2}.
The observed variables \eqn{Y_t} are distributed as \eqn{\mathrm{lognormal}(\log{X_t},\tau)}.
- The model is parameterized by the logarithms of \eqn{r}, \eqn{K}, \eqn{\sigma}, and \eqn{\tau};
+ Parameters include the per-capita growth rate \eqn{r}, the carrying capacity \eqn{K}, the process noise s.d. \eqn{sigma}, the measurement error s.d. \eqn{tau}, and the initial condition \eqn{X_0}.
+ The model is parameterized internally by the logarithms of \eqn{r}, \eqn{K}, \eqn{\sigma}, and \eqn{\tau};
the initial condition is parameterized directly.
+ The \code{pomp} object includes parameter transformations to and from this internal parameterization.
}
\examples{
data(gompertz)
plot(gompertz)
coef(gompertz)
+coef(gompertz,transform=TRUE)
}
\seealso{
\code{\link{pomp-class}} and the introductory vignette \code{vignette("intro_to_pomp")}.
More information about the pomp-commits
mailing list