[Pomp-commits] r654 - in pkg: . data inst inst/data-R inst/doc man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Apr 12 04:10:04 CEST 2012
Author: kingaa
Date: 2012-04-12 04:10:04 +0200 (Thu, 12 Apr 2012)
New Revision: 654
Modified:
pkg/DESCRIPTION
pkg/data/bbs.rda
pkg/data/blowflies.rda
pkg/data/dacca.rda
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/data/gompertz.rda
pkg/data/ou2.rda
pkg/data/ricker.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/TODO
pkg/inst/data-R/blowflies.R
pkg/inst/data-R/ricker.R
pkg/inst/data-R/sir.R
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/bsmc-ricker-flat-prior.rda
pkg/inst/doc/intro_to_pomp.Rnw
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/doc/ricker-mif.rda
pkg/inst/doc/ricker-probe-match.rda
pkg/man/blowflies.Rd
pkg/man/gompertz.Rd
pkg/man/parmat.Rd
pkg/man/ricker.Rd
pkg/man/sir.Rd
pkg/man/traj-match.Rd
pkg/src/blowfly.c
pkg/src/ricker.c
pkg/src/sir.c
pkg/tests/bbs.R
pkg/tests/bbs.Rout.save
pkg/tests/dacca.Rout.save
pkg/tests/dimchecks.Rout.save
pkg/tests/fhn.Rout.save
pkg/tests/filtfail.Rout.save
pkg/tests/gillespie.R
pkg/tests/gillespie.Rout.save
pkg/tests/gompertz.Rout.save
pkg/tests/logistic.Rout.save
pkg/tests/ou2-bsmc.Rout.save
pkg/tests/ou2-forecast.Rout.save
pkg/tests/ou2-icfit.Rout.save
pkg/tests/ou2-kalman.Rout.save
pkg/tests/ou2-mif-fp.Rout.save
pkg/tests/ou2-mif.Rout.save
pkg/tests/ou2-nlf.Rout.save
pkg/tests/ou2-pmcmc.Rout.save
pkg/tests/ou2-probe.Rout.save
pkg/tests/ou2-procmeas.Rout.save
pkg/tests/ou2-simulate.Rout.save
pkg/tests/ou2-trajmatch.Rout.save
pkg/tests/partrans.Rout.save
pkg/tests/pfilter.Rout.save
pkg/tests/pomppomp.Rout.save
pkg/tests/ricker-bsmc.R
pkg/tests/ricker-bsmc.Rout.save
pkg/tests/ricker-probe.R
pkg/tests/ricker-probe.Rout.save
pkg/tests/ricker-spect.R
pkg/tests/ricker-spect.Rout.save
pkg/tests/ricker.R
pkg/tests/ricker.Rout.save
pkg/tests/rw2.Rout.save
pkg/tests/sir.Rout.save
pkg/tests/skeleton.R
pkg/tests/skeleton.Rout.save
pkg/tests/steps.Rout.save
pkg/tests/synlik.Rout.save
pkg/tests/verhulst.Rout.save
Log:
- update blowflies example: better parameterization, better documentation
- improve parameterization of ricker, euler.sir, gillespie.sir, and bbs examples
- these changes are not backward compatible!
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/DESCRIPTION 2012-04-12 02:10:04 UTC (rev 654)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.41-3
-Date: 2012-04-10
+Version: 0.41-4
+Date: 2012-04-12
Revision: $Rev$
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing, Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
Maintainer: Aaron A. King <kingaa at umich.edu>
Modified: pkg/data/bbs.rda
===================================================================
(Binary files differ)
Modified: pkg/data/blowflies.rda
===================================================================
(Binary files differ)
Modified: pkg/data/dacca.rda
===================================================================
(Binary files differ)
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gompertz.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ricker.rda
===================================================================
(Binary files differ)
Modified: pkg/data/rw2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/verhulst.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/TODO
===================================================================
--- pkg/inst/TODO 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/inst/TODO 2012-04-12 02:10:04 UTC (rev 654)
@@ -1,5 +1,9 @@
-* put BSMC example into the intro vignette
+* take away confusing name-changes in introductory examples (ricker, sir)
+* fix bbs parameter transformations to include sigma
+
+* pompBuilder function for automatic pomp construction
+
* effects of replacing 'sample.int' with systematic resampling inside 'bsmc'
* sorting vs no-sorting systematic resampling algorithm?
@@ -29,8 +33,6 @@
* Add LPA model examples.
-* Add BBS example.
-
* SDE examples.
* Extended Kalman filter.
Modified: pkg/inst/data-R/blowflies.R
===================================================================
--- pkg/inst/data-R/blowflies.R 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/inst/data-R/blowflies.R 2012-04-12 02:10:04 UTC (rev 654)
@@ -21,20 +21,25 @@
step.fun="_blowfly_model_simulator",
delta.t=1,
),
- paramnames=c("log.P","log.N0","log.delta","log.sigma.P","log.sigma.d","tau","log.sigma.y"),
+ paramnames=c("P","N0","delta","sigma.P","sigma.d","tau","sigma.y"),
statenames=c("N1","R","S","e","eps"),
obsnames=c("y"),
- measurement.model=y~nbinom(mu=N1,size=exp(-2*log.sigma.y)),
+ measurement.model=y~nbinom(mu=N1,size=1/(sigma.y^2)),
y.init=with( ## initial data
raw.data,
- approx(
- x=day,
- y=y,
- xout=seq(from=0,to=14,by=1),
- rule=2
- )$y
+ approx(x=day,y=y,xout=seq(from=0,to=14,by=1),rule=2)$y
),
# y.init=c(948, 948, 942, 930, 911, 885, 858, 833.7, 801, 748.3, 676, 589.8, 504, 434.9, 397),
+ parameter.inv.transform=function(params,...) {
+ p <- c(log(params[c("delta","sigma.P","sigma.d","sigma.y","P","N0")]),params["tau"])
+ names(p) <- c(paste("log",c("delta","sigma.P","sigma.d","sigma.y","P","N0"),sep="."),"tau")
+ p
+ },
+ parameter.transform=function(params,...) {
+ p <- c(exp(params[paste("log",c("delta","sigma.P","sigma.d","sigma.y","P","N0"),sep=".")]),params["tau"])
+ names(p) <- c(c("delta","sigma.P","sigma.d","sigma.y","P","N0"),"tau")
+ p
+ },
initializer=function (params, t0, y.init, ...) {
ntau <- length(y.init)
n <- y.init[ntau:1]
@@ -44,25 +49,14 @@
) -> blowflies1
pomp(
- data=subset(raw.data[c("day","y")],day>14&day<400),
- times="day",
- t0=14,
+ blowflies1,
rprocess=discrete.time.sim(
step.fun="_blowfly_model_simulator",
delta.t=2,
),
- paramnames=c("log.P","log.N0","log.delta","log.sigma.P","log.sigma.d","tau","log.sigma.y"),
- statenames=c("N1","R","S","e","eps"),
- obsnames=c("y"),
- measurement.model=y~nbinom(mu=N1,size=exp(-2*log.sigma.y)),
y.init=with( ## initial data
raw.data,
- approx(
- x=day,
- y=y,
- xout=seq(from=0,to=14,by=2),
- rule=2
- )$y
+ approx(x=day,y=y,xout=seq(from=0,to=14,by=2),rule=2)$y
),
#y.init=c(948, 942, 911, 858, 801, 676, 504, 397),
initializer=function (params, t0, y.init, ...) {
@@ -74,26 +68,26 @@
) -> blowflies2
## mle from search to date
-coef(blowflies1) <- c(
- log.P = 1.189 ,
- log.delta = -1.828 ,
- log.N0 = 6.522 ,
- log.sigma.P = 0.301 ,
- log.sigma.d = -0.292 ,
- log.sigma.y = -3.625 ,
- tau = 14
- )
+coef(blowflies1,transform=TRUE) <- c(
+ log.P = 1.189,
+ log.delta = -1.828,
+ log.N0 = 6.522,
+ log.sigma.P = 0.301,
+ log.sigma.d = -0.292,
+ log.sigma.y = -3.625,
+ tau = 14
+ )
## mle from search to date
-coef(blowflies2) <- c(
- log.P = 1.005 ,
- log.delta = -1.75 ,
- log.N0 = 6.685 ,
- log.sigma.P = 0.366 ,
- log.sigma.d = -0.274 ,
- log.sigma.y = -4.524 ,
- tau = 7
- )
+coef(blowflies2,transform=TRUE) <- c(
+ log.P = 1.005,
+ log.delta = -1.75,
+ log.N0 = 6.685,
+ log.sigma.P = 0.366,
+ log.sigma.d = -0.274,
+ log.sigma.y = -4.524,
+ tau = 7
+ )
test <- FALSE
if(test){
@@ -109,9 +103,9 @@
## check that it matches the deterministic skeleton when noise is small
params.1.skel <- coef(blowflies1)
- params.1.skel["log.sigma.P"] <- log(0.00001)
- params.1.skel["log.sigma.d"] <- log(0.00001)
- params.1.skel["log.sigma.y"] <- log(0.00001)
+ params.1.skel["sigma.P"] <- 0.00001
+ params.1.skel["sigma.d"] <- 0.00001
+ params.1.skel["sigma.y"] <- 0.00001
simulate(blowflies1,params=params.1.skel,nsim=1,seed=73691676L) -> b1.skel
plot(obs(blowflies1)['y',],ty='l',lty="dashed")
lines(obs(b1.skel)['y',],ty='l')
Modified: pkg/inst/data-R/ricker.R
===================================================================
--- pkg/inst/data-R/ricker.R 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/inst/data-R/ricker.R 2012-04-12 02:10:04 UTC (rev 654)
@@ -12,23 +12,21 @@
dmeasure="_ricker_poisson_dmeasure",
skeleton.type="map",
skeleton="_ricker_skeleton",
- paramnames=c("log.r","sigma","phi"),
+ paramnames=c("r","sigma","phi"),
statenames=c("N","e"),
obsnames=c("y"),
parameter.inv.transform=function(params,...) {
- params <- c(params["log.r"],log(params[c("sigma","phi","N.0")]),params["e.0"])
- names(params) <- c("log.r","log.sigma","log.phi","log.N.0","e.0")
+ params[c("r","sigma","phi","N.0")] <- log(params[c("r","sigma","phi","N.0")])
params
},
parameter.transform=function(params,...) {
- params <- c(params["log.r"],exp(params[c("log.sigma","log.phi","log.N.0")]),params["e.0"])
- names(params) <- c("log.r","sigma","phi","N.0","e.0")
+ params[c("r","sigma","phi","N.0")] <- exp(params[c("r","sigma","phi","N.0")])
params
},
PACKAGE="pomp"
),
params=c(
- log.r=3.8,
+ r=exp(3.8),
sigma=0.3,
phi=10,
N.0=7,
Modified: pkg/inst/data-R/sir.R
===================================================================
--- pkg/inst/data-R/sir.R 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/inst/data-R/sir.R 2012-04-12 02:10:04 UTC (rev 654)
@@ -21,7 +21,7 @@
statenames=c("S","I","R","cases","W"),
paramnames=c(
"gamma","mu","iota",
- "log.beta1","beta.sd","pop","rho",
+ "beta1","beta.sd","pop","rho",
"nbasis","degree","period",
"S.0","I.0","R.0"
),
@@ -41,14 +41,14 @@
)
coef(po) <- c(
- gamma=26,mu=0.02,iota=0.01,
- nbasis=3,degree=3,period=1,
- log.beta1=log(1200),log.beta2=log(1800),log.beta3=log(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
- )
+ gamma=26,mu=0.02,iota=0.01,
+ nbasis=3,degree=3,period=1,
+ 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
+ )
simulate(po,nsim=1,seed=329348545L) -> euler.sir
@@ -94,7 +94,7 @@
coef(po) <- c(
gamma=24,mu=1/70,iota=0.1,
- log.beta1=log(330),log.beta2=log(410),log.beta3=log(490),
+ beta1=330,beta2=410,beta3=490,
rho=0.1,
S.0=0.05,I.0=1e-4,R.0=0.95,
pop=1000000,
@@ -111,7 +111,7 @@
save(gillespie.sir,file="gillespie.sir.rda",compress="xz")
tc <- textConnection("
-day;flu
+day;reports
1;3
2;8
3;28
@@ -142,21 +142,34 @@
),
skeleton.type="vectorfield",
skeleton="_sir_ODE",
- measurement.model=flu~norm(mean=rho*cases,sd=1+sigma*cases),
+ measurement.model=reports~norm(mean=rho*cases,sd=1+sigma*cases),
PACKAGE="pomp",
- obsnames = c("flu"),
+ obsnames = c("reports"),
statenames=c("S","I","R","cases","W"),
paramnames=c(
"gamma","mu","iota",
- "log.beta1","beta.sd","pop","rho",
+ "beta","beta.sd","pop","rho",
"nbasis","degree","period",
"S.0","I.0","R.0"
),
zeronames=c("cases"),
comp.names=c("S","I","R"),
ic.names=c("S.0","I.0","R.0"),
- parameter.transform="_sir_par_trans",
- parameter.inv.transform="_sir_par_untrans",
+ logvar=c(
+ "beta","gamma","mu","iota","sigma","beta.sd",
+ "S.0","I.0","R.0"
+ ),
+ logitvar="rho",
+ parameter.inv.transform=function (params, logvar, logitvar, ...) {
+ params[logvar] <- log(params[logvar])
+ params[logitvar] <- qlogis(params[logitvar])
+ params
+ },
+ parameter.transform=function (params, logvar, logitvar, ...) {
+ params[logvar] <- exp(params[logvar])
+ params[logitvar] <- plogis(params[logitvar])
+ params
+ },
initializer=function(params, t0, comp.names, ic.names, ...) {
snames <- c("S","I","R","cases","W")
fracs <- params[ic.names]
@@ -170,7 +183,7 @@
coef(po) <- c(
gamma=1/3,mu=0.0,iota=0.0,
nbasis=1,degree=0,period=1,
- log.beta1=0.33,
+ beta=1.4,
beta.sd=0,
pop=1400,
rho=0.9,sigma=3.6,
@@ -179,4 +192,3 @@
bbs <- po
save(bbs,file="bbs.rda",compress="xz")
-
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/bsmc-ricker-flat-prior.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.Rnw
===================================================================
--- pkg/inst/doc/intro_to_pomp.Rnw 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/inst/doc/intro_to_pomp.Rnw 2012-04-12 02:10:04 UTC (rev 654)
@@ -838,7 +838,7 @@
ricker.sim <- function (x, t, params, delta.t, ...) {
e <- rnorm(n=1,mean=0,sd=params["sigma"])
xnew <- c(
- exp(params["log.r"]+log(x["N"])-x["N"]+e),
+ params["r"]*x["N"]*exp(-x["N"]+e),
e
)
names(xnew) <- c("N","e")
@@ -860,7 +860,7 @@
measurement.model=y~pois(lambda=N*phi)
)
coef(ricker) <- c(
- log.r=3.8,
+ r=exp(3.8),
sigma=0.3,
phi=10,
N.0=7,
@@ -900,7 +900,7 @@
To see this, we'll apply probe to the Ricker model at the true parameters and at a wild guess.
<<first-probe>>=
pb.truth <- probe(ricker,probes=plist,nsim=1000,seed=1066L)
-guess <- c(log.r=log(20),sigma=1,phi=20,N.0=7,e.0=0)
+guess <- c(r=20,sigma=1,phi=20,N.0=7,e.0=0)
pb.guess <- probe(ricker,params=guess,probes=plist,nsim=1000,seed=1066L)
@
Results summaries and diagnostic plots showing the model-data agreement and correlations among the probes can be obtained by
@@ -942,7 +942,7 @@
<<ricker-probe-match-calc,eval=F>>=
pm <- probe.match(
pb.guess,
- est=c("log.r","log.sigma","log.phi"),
+ est=c("r","sigma","phi"),
transform=TRUE,
method="Nelder-Mead",
maxit=2000,
@@ -980,7 +980,7 @@
var.factor=2,
ic.lag=3,
max.fail=50,
- rw.sd=c(log.r=0.1,log.sigma=0.1,log.phi=0.1)
+ rw.sd=c(r=0.1,sigma=0.1,phi=0.1)
)
mf <- continue(mf,Nmif=500,max.fail=20)
@
@@ -1076,7 +1076,7 @@
start=starts[[j]],
transform=log,
transform.params=TRUE,
- est=c("log.K","log.r"),
+ est=c("K","r"),
lags=c(1,2),
seed=7639873L,
method="Nelder-Mead",
@@ -1252,7 +1252,7 @@
true.fit <- nlf(
gompertz,
transform.params=TRUE,
- est=c("log.K","log.r"),
+ est=c("K","r"),
lags=2,
seed=7639873,
method="Nelder-Mead",
@@ -1270,8 +1270,8 @@
}
@
From \verb+true.fit$params+ and \verb+true.fit$se+ we get the estimates ($\pm$ 1 standard error)
-${r}=$~\Sexpr{signif(true.fit$params["r"],2)}~$\pm$~\Sexpr{signif(true.fit$params["r"]*true.fit$se["log.r"],2)}
-and ${K}=$~\Sexpr{signif(true.fit$params["K"],2)}~$\pm$~\Sexpr{signif(true.fit$params["K"]*true.fit$se["log.K"],2)}.
+${r}=$~\Sexpr{signif(true.fit$params["r"],2)}~$\pm$~\Sexpr{signif(true.fit$params["r"]*true.fit$se["r"],2)}
+and ${K}=$~\Sexpr{signif(true.fit$params["K"],2)}~$\pm$~\Sexpr{signif(true.fit$params["K"]*true.fit$se["K"],2)}.
%%$\log K = 0.081 \pm 0.064$, $\log r = -2.2 \pm 0.51$
The standard errors provided by \code{nlf} are based on a Newey-West estimate of the variance-covariance matrix that is generally
@@ -1298,7 +1298,7 @@
gompertz,
start=coef(gompertz),
transform.params=TRUE,
- est=c("log.K","log.r"),
+ est=c("K","r"),
lags=lags,
seed=7639873,
bootstrap=TRUE,
@@ -1351,7 +1351,7 @@
fit <- nlf(
gompertz,
transform.params=TRUE,
- est=c("log.K","log.r"),
+ est=c("K","r"),
lags=lags,
seed=7639873L,
bootstrap=TRUE,
@@ -1395,7 +1395,7 @@
set.seed(136872209L)
Np <- 10000
prior1 <- parmat(coef(ricker),nrep=Np)
-prior1["log.r",] <- runif(n=Np,min=2,max=5)
+prior1["r",] <- exp(runif(n=Np,min=2,max=5))
prior1["sigma",] <- runif(n=Np,min=0.1,max=1)
@
Now \code{params} holds these 10,000 samples.
@@ -1403,7 +1403,7 @@
Note that, by specifying \code{transform=TRUE}, we cause the estimation to proceed on the transformed scale.
<<bsmc-example-flat-prior-3,eval=F>>=
fit1 <- bsmc(ricker,params=prior1,transform=TRUE,
- est=c("log.r","log.sigma"),smooth=0.2)
+ est=c("r","sigma"),smooth=0.2)
@
<<bsmc-example-flat-prior-eval,eval=T,echo=F>>=
binary.file <- "bsmc-ricker-flat-prior.rda"
@@ -1425,7 +1425,7 @@
\begin{figure}
<<bsmc-example-flat-prior-plot,fig=T,pdf=F,png=T,height=6,width=6,echo=F>>=
- plot(fit1,pars=c("log.r","sigma"),thin=5000)
+ plot(fit1,pars=c("r","sigma"),thin=5000)
@
\caption{
Results of \code{bsmc} on the Ricker model.
@@ -1442,12 +1442,12 @@
set.seed(90348704L)
Np <- 10000
prior2 <- parmat(coef(ricker),nrep=Np)
-## normal prior on log.r
-prior2["log.r",] <- rnorm(n=Np,mean=4,sd=3)
+## log-normal prior on r
+prior2["r",] <- rlnorm(n=Np,meanlog=4,sdlog=3)
## log-normal prior on sigma
prior2["sigma",] <- rlnorm(n=Np,meanlog=log(0.5),sdlog=5)
fit2 <- bsmc(ricker,params=prior2,transform=TRUE,
- est=c("log.r","log.sigma"),smooth=0.2)
+ est=c("r","sigma"),smooth=0.2)
signif(coef(fit2),digits=2)
@
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/ricker-mif.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/ricker-probe-match.rda
===================================================================
(Binary files differ)
Modified: pkg/man/blowflies.Rd
===================================================================
--- pkg/man/blowflies.Rd 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/man/blowflies.Rd 2012-04-12 02:10:04 UTC (rev 654)
@@ -19,6 +19,24 @@
Unlimited quantities of larval food were provided;
the adult food supply (ground liver) was constant at 0.4g per day.
The data were taken from the table provided by Brillinger et al. (1980).
+
+ The models are discrete delay equations:
+ \deqn{R(t+1) \sim \mathrm{Poisson}(P N(t-\tau) \exp{(-N(t-\tau)/N_{0})} e(t+1) {\Delta}t)}{%
+ R[t+1] ~ Poisson(P N[t-tau] exp(-N[t-tau]/N0) e[t+1] dt)}
+ \deqn{S(t+1) \sim \mathrm{binomial}(N(t),\exp{(-\delta \epsilon(t+1) {\Delta}t)})}{%
+ S[t+1] ~ binomial(N[t],exp(-delta eps[t+1] dt))}
+ \deqn{N(t) = R(t)+S(t)}{N[t]=R[t]+S[t]}
+ where \eqn{e(t)}{e[t]} and \eqn{\epsilon(t)}{eps[t]} are Gamma-distributed i.i.d. random variables with mean 1 and variances \eqn{{\sigma_p^2}/{{\Delta}t}}{sigma.p^2/dt}, \eqn{{\sigma_d^2}/{{\Delta}t}}{sigma.d^2/dt}, respectively.
+ \code{blowflies1} has a timestep (\eqn{{\Delta}t}{dt}) of 1 day, and \code{blowflies2} has a timestep of 2 days.
+ The process model in \code{blowflies1} thus corresponds exactly to that studied by Wood (2010).
+ The measurement model in both cases is taken to be
+ \deqn{y(t) \sim \mathrm{negbin}(N(t),1/\sigma_y^2)}{y[t] ~ negbin(N[t],1/sigma.y^2)}, i.e., the observations are assumed to be negative-binomially distributed with mean \eqn{N(t)}{N[t]} and variance \eqn{N(t)+(\sigma_y N(t))^2}{N[t]+(sigma.y N[t])^2}.
+
+ Do
+
+ \code{file.show(system.file("data-R","blowflies.R",package="pomp"))}
+
+ to view the code that constructs these pomp objects.
}
\seealso{\code{\link{pomp-class}} and the vignettes}
\references{
Modified: pkg/man/gompertz.Rd
===================================================================
--- pkg/man/gompertz.Rd 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/man/gompertz.Rd 2012-04-12 02:10:04 UTC (rev 654)
@@ -7,12 +7,10 @@
}
\usage{data(gompertz)}
\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)}.
- 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.
+ The state process is \eqn{X_{t+1} = K^{1-S} X_{t}^S \epsilon_{t}}{X[t+1]=K^(1-S) X[t]^S eps[t]}, where \eqn{S=e^{-r}}{S=e^{-r}} and the \eqn{\epsilon_t}{eps[t]} are i.i.d. lognormal random deviates with variance \eqn{\sigma^2}{sigma^2}.
+ The observed variables \eqn{Y_t} are distributed as \eqn{\mathrm{lognormal}(\log{X_t},\tau)}{lognormal(log(X[t]),tau)}.
+ Parameters include the per-capita growth rate \eqn{r}, the carrying capacity \eqn{K}, the process noise s.d. \eqn{\sigma}{sigma}, the measurement error s.d. \eqn{\tau}{tau}, and the initial condition \eqn{X_0}{X[0]}.
+ The \code{pomp} object includes parameter transformations that log-transform the parameters for estimation purposes.
}
\examples{
data(gompertz)
Modified: pkg/man/parmat.Rd
===================================================================
--- pkg/man/parmat.Rd 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/man/parmat.Rd 2012-04-12 02:10:04 UTC (rev 654)
@@ -23,7 +23,7 @@
## generate a bifurcation diagram for the Ricker map
data(ricker)
p <- parmat(coef(ricker),nrep=500)
- p["log.r",] <- seq(from=1.5,to=4,length=500)
+ p["r",] <- exp(seq(from=1.5,to=4,length=500))
x <- trajectory(ricker,times=seq(from=1000,to=2000,by=1),params=p)
- matplot(p["log.r",],x["N",,],pch='.',col='black',xlab="log(r)",ylab="N")
+ matplot(p["r",],x["N",,],pch='.',col='black',xlab="log(r)",ylab="N",log='x')
}
Modified: pkg/man/ricker.Rd
===================================================================
--- pkg/man/ricker.Rd 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/man/ricker.Rd 2012-04-12 02:10:04 UTC (rev 654)
@@ -7,7 +7,7 @@
}
\usage{data(ricker)}
\details{
- The state process is \eqn{N_{t+1} = r N_{t} \exp(-N_{t}+e_{t})}{N[t+1] = r N[t] exp(-N[t]+e[t])}, where the \eqn{e_t}{e[t]} are i.i.d. normal random deviates with variance \eqn{\sigma^2}{sigma^2}.
+ The state process is \eqn{N_{t+1} = r N_{t} \exp(-N_{t}+e_{t})}{N[t+1] = r N[t] exp(-N[t]+e[t])}, where the \eqn{e_t}{e[t]} are i.i.d. normal random deviates with zero mean and variance \eqn{\sigma^2}{sigma^2}.
The observed variables \eqn{y_t}{y[t]} are distributed as \eqn{\mathrm{Poisson}(\phi N_t)}{Poisson(phi N[t])}.
}
\examples{
Modified: pkg/man/sir.Rd
===================================================================
--- pkg/man/sir.Rd 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/man/sir.Rd 2012-04-12 02:10:04 UTC (rev 654)
@@ -25,13 +25,16 @@
data(euler.sir)
plot(euler.sir)
plot(simulate(euler.sir,seed=20348585))
+coef(euler.sir)
data(gillespie.sir)
plot(gillespie.sir)
plot(simulate(gillespie.sir,seed=20348585))
+coef(gillespie.sir)
data(bbs)
plot(bbs)
+coef(bbs)
}
\references{
Anonymous (1978).
Modified: pkg/man/traj-match.Rd
===================================================================
--- pkg/man/traj-match.Rd 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/man/traj-match.Rd 2012-04-12 02:10:04 UTC (rev 654)
@@ -135,7 +135,7 @@
lines(x2~time,data=as(res,"data.frame"),col='red')
data(ricker)
- ofun <- traj.match.objfun(ricker,est=c("log.r","log.phi"),transform=TRUE)
+ ofun <- traj.match.objfun(ricker,est=c("r","phi"),transform=TRUE)
optim(fn=ofun,par=c(2,0),method="BFGS")
}
\seealso{
Modified: pkg/src/blowfly.c
===================================================================
--- pkg/src/blowfly.c 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/src/blowfly.c 2012-04-12 02:10:04 UTC (rev 654)
@@ -4,35 +4,31 @@
#include "pomp.h"
-#define LOG_P (p[parindex[0]]) // growth rate
-#define LOG_NZERO (p[parindex[1]]) // density-dependence parameter
-#define LOG_DELTA (p[parindex[2]]) // survival parameter
-#define LOG_SIGMAP (p[parindex[3]]) // recruitment noise SD
-#define LOG_SIGMAD (p[parindex[4]]) // survivorship noise SD
-#define TAU (p[parindex[5]]) // delay
+#define P (p[parindex[0]]) // growth rate
+#define NZERO (p[parindex[1]]) // density-dependence parameter
+#define DELTA (p[parindex[2]]) // survival parameter
+#define SIGMAP (p[parindex[3]]) // recruitment noise SD
+#define SIGMAD (p[parindex[4]]) // survivorship noise SD
+#define TAU (p[parindex[5]]) // delay
-#define N (&x[stateindex[0]]) // total population
-#define R (x[stateindex[1]]) // recruits
-#define S (x[stateindex[2]]) // survivors
-#define E (x[stateindex[3]]) // recruitment noise
-#define EPS (x[stateindex[4]]) // survival noise
+#define N (&x[stateindex[0]]) // total population
+#define R (x[stateindex[1]]) // recruits
+#define S (x[stateindex[2]]) // survivors
+#define E (x[stateindex[3]]) // recruitment noise
+#define EPS (x[stateindex[4]]) // survival noise
// Ricker model with log-normal process noise
void _blowfly_model_simulator (double *x, const double *p,
- const int *stateindex, const int *parindex, const int *covindex,
- int covdim, const double *covar,
- double t, double dt)
+ const int *stateindex, const int *parindex, const int *covindex,
+ int covdim, const double *covar,
+ double t, double dt)
{
- double var_p = exp(2*LOG_SIGMAP)/dt;
- double var_d = exp(2*LOG_SIGMAD)/dt;
- double nzero = exp(LOG_NZERO);
- double e = (var_p > 0.0) ? rgamma(1.0/var_p,var_p) : 1.0;
- double eps = (var_d > 0.0) ? rgamma(1.0/var_d,var_d) : 1.0;
- double P = exp(LOG_P);
- int tau = (int) TAU;
+ double e = rgammawn(SIGMAP,dt)/dt;
+ double eps = rgammawn(SIGMAD,dt)/dt;
+ int tau = nearbyint(TAU);
int k;
- R = rpois(P*N[tau]*exp(-N[tau]/nzero)*dt*e);
- S = rbinom(N[0],exp(-exp(LOG_DELTA)*dt*eps));
+ R = rpois(P*N[tau]*exp(-N[tau]/NZERO)*dt*e);
+ S = rbinom(N[0],exp(-DELTA*dt*eps));
E = e;
EPS = eps;
for (k = tau; k > 0; k--) N[k] = N[k-1];
Modified: pkg/src/ricker.c
===================================================================
--- pkg/src/ricker.c 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/src/ricker.c 2012-04-12 02:10:04 UTC (rev 654)
@@ -4,7 +4,7 @@
#include "pomp.h"
-#define LOG_R (p[parindex[0]]) // log growth rate
+#define R (p[parindex[0]]) // growth rate
#define SIGMA (p[parindex[1]]) // process noise level
#define PHI (p[parindex[2]]) // measurement scale parameter
@@ -34,7 +34,7 @@
double t, double dt)
{
double e = (SIGMA > 0.0) ? rnorm(0,SIGMA) : 0.0;
- N = exp(LOG_R+log(N)-N+e);
+ N = exp(log(R)+log(N)-N+e);
E = e;
}
@@ -42,13 +42,13 @@
const int *stateindex, const int *parindex, const int *covindex,
int covdim, const double *covar, double t)
{
- f[0] = exp(LOG_R+log(N)-N);
+ f[0] = exp(log(R)+log(N)-N);
f[1] = 0.0;
}
#undef N
#undef E
-#undef LOG_R
+#undef R
#undef SIGMA
#undef PHI
Modified: pkg/src/sir.c
===================================================================
--- pkg/src/sir.c 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/src/sir.c 2012-04-12 02:10:04 UTC (rev 654)
@@ -27,7 +27,7 @@
#define GAMMA (p[parindex[0]]) // recovery rate
#define MU (p[parindex[1]]) // baseline birth and death rate
#define IOTA (p[parindex[2]]) // import rate
-#define LOGBETA (p[parindex[3]]) // transmission rate
+#define BETA (&p[parindex[3]]) // transmission rates
#define BETA_SD (p[parindex[4]]) // environmental stochasticity SD in transmission rate
#define POPSIZE (p[parindex[5]]) // population size
#define RHO (p[parindex[6]]) // reporting probability
@@ -54,9 +54,13 @@
void _sir_par_untrans (double *pt, double *p, int *parindex)
{
+ int nbasis = (int) NBASIS;
+ int k;
pt[parindex[0]] = log(GAMMA);
pt[parindex[1]] = log(MU);
pt[parindex[2]] = log(IOTA);
+ for (k = 0; k < nbasis; k++)
+ pt[parindex[3]+k] = log(BETA[k]);
pt[parindex[4]] = log(BETA_SD);
pt[parindex[6]] = logit(RHO);
pt[parindex[10]] = log(S0);
@@ -66,9 +70,13 @@
void _sir_par_trans (double *pt, double *p, int *parindex)
{
+ int nbasis = (int) NBASIS;
+ int k;
pt[parindex[0]] = exp(GAMMA);
pt[parindex[1]] = exp(MU);
pt[parindex[2]] = exp(IOTA);
+ for (k = 0; k < nbasis; k++)
+ pt[parindex[3]+k] = exp(BETA[k]);
pt[parindex[4]] = exp(BETA_SD);
pt[parindex[6]] = expit(RHO);
pt[parindex[10]] = exp(S0);
@@ -116,13 +124,16 @@
double trans[nrate]; // transition numbers
double beta;
double dW;
- int nseas = (int) NBASIS; // number of seasonal basis functions
+ int nbasis = (int) NBASIS; // number of seasonal basis functions
int deg = (int) DEGREE; // degree of seasonal basis functions
- double seasonality[nseas];
+ double seasonality[nbasis];
+ int k;
- if (nseas <= 0) return;
- periodic_bspline_basis_eval(t,PERIOD,deg,nseas,&seasonality[0]);
- beta = exp(dot_product(nseas,&seasonality[0],&LOGBETA));
+ if (nbasis <= 0) return;
+ periodic_bspline_basis_eval(t,PERIOD,deg,nbasis,&seasonality[0]);
+ for (k = 0, beta = 0; k < nbasis; k++)
+ beta += seasonality[k]*log(BETA[k]);
+ beta = exp(beta);
// test to make sure the parameters and state variable values are sane
if (!(R_FINITE(beta)) ||
@@ -171,14 +182,17 @@
double rate[nrate]; // transition rates
double term[nrate]; // terms in the equations
double beta;
- int nseas = (int) NBASIS; // number of seasonal basis functions
+ int nbasis = (int) NBASIS; // number of seasonal basis functions
int deg = (int) DEGREE; // degree of seasonal basis functions
- double seasonality[nseas];
-
- if (nseas <= 0) return;
- periodic_bspline_basis_eval(t,PERIOD,deg,nseas,&seasonality[0]);
- beta = exp(dot_product(nseas,&seasonality[0],&LOGBETA));
+ double seasonality[nbasis];
+ int k;
+ if (nbasis <= 0) return;
+ periodic_bspline_basis_eval(t,PERIOD,deg,nbasis,&seasonality[0]);
+ for (k = 0, beta = 0; k < nbasis; k++)
+ beta += seasonality[k]*log(BETA[k]);
+ beta = exp(beta);
+
// compute the transition rates
rate[0] = MU*POPSIZE; // birth into susceptible class
rate[1] = (IOTA+beta*INFD)/POPSIZE; // force of infection
@@ -221,9 +235,10 @@
int ncovar, double *covar) {
double beta;
double rate = 0.0;
- int nseas = (int) NBASIS; // number of seasonal basis functions
+ int nbasis = (int) NBASIS; // number of seasonal basis functions
int deg = (int) DEGREE; // degree of seasonal basis functions
- double seasonality[nseas];
+ double seasonality[nbasis];
+ int k;
switch (j) {
case 1: // birth
@@ -233,8 +248,10 @@
rate = MU*SUSC;
break;
case 3: // infection
- periodic_bspline_basis_eval(t,PERIOD,deg,nseas,&seasonality[0]);
- beta = exp(dot_product(nseas,&seasonality[0],&LOGBETA));
+ periodic_bspline_basis_eval(t,PERIOD,deg,nbasis,&seasonality[0]);
+ for (k = 0, beta = 0; k < nbasis; k++)
+ beta += seasonality[k]*log(BETA[k]);
+ beta = exp(beta);
rate = (beta*INFD+IOTA)*SUSC/POPSIZE;
break;
case 4: // infected death
@@ -252,4 +269,3 @@
}
return rate;
}
-
Modified: pkg/tests/bbs.R
===================================================================
--- pkg/tests/bbs.R 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/tests/bbs.R 2012-04-12 02:10:04 UTC (rev 654)
@@ -8,10 +8,10 @@
coef(bbs,transform=TRUE)
prior <- parmat(coef(bbs),nrep=1000)
-prior["log.beta1",] <- runif(n=1000,min=1,max=2)
+prior["beta",] <- exp(runif(n=1000,min=1,max=2))
prior["sigma",] <- runif(n=1000,min=2,max=4)
-fit1 <- bsmc(bbs,params=prior,transform=TRUE,est=c("log.beta1","sigma"),smooth=0.2)
+fit1 <- bsmc(bbs,params=prior,transform=TRUE,est=c("beta","sigma"),smooth=0.2)
signif(coef(fit1),3)
-fit2 <- traj.match(bbs,est=c("log.beta1","sigma"),transform=TRUE)
+fit2 <- traj.match(bbs,est=c("beta","sigma"),transform=TRUE)
signif(coef(fit2),3)
Modified: pkg/tests/bbs.Rout.save
===================================================================
--- pkg/tests/bbs.Rout.save 2012-04-10 19:48:02 UTC (rev 653)
+++ pkg/tests/bbs.Rout.save 2012-04-12 02:10:04 UTC (rev 654)
@@ -1,5 +1,5 @@
-R version 2.14.2 (2012-02-29)
+R version 2.15.0 (2012-03-30)
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-unknown-linux-gnu (64-bit)
@@ -28,36 +28,38 @@
> coef(bbs)
gamma mu iota nbasis degree period
0.3333333 0.0000000 0.0000000 1.0000000 0.0000000 1.0000000
- log.beta1 beta.sd pop rho sigma S.0
- 0.3300000 0.0000000 1400.0000000 0.9000000 3.6000000 0.9990000
+ beta beta.sd pop rho sigma S.0
+ 1.4000000 0.0000000 1400.0000000 0.9000000 3.6000000 0.9990000
I.0 R.0
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 654
More information about the pomp-commits
mailing list