[Pomp-commits] r528 - in pkg: . R 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 Jul 28 23:19:28 CEST 2011
Author: kingaa
Date: 2011-07-28 23:19:27 +0200 (Thu, 28 Jul 2011)
New Revision: 528
Modified:
pkg/DESCRIPTION
pkg/R/pomp-methods.R
pkg/R/pomp.R
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/ChangeLog
pkg/inst/NEWS
pkg/inst/data-R/dacca.R
pkg/inst/data-R/euler.sir.R
pkg/inst/data-R/gompertz.R
pkg/inst/doc/advanced_topics_in_pomp.Rnw
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/gompertz-multi-mif.rda
pkg/inst/doc/gompertz-trajmatch.rda
pkg/inst/doc/intro_to_pomp.Rnw
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/doc/nlf-block-boot.rda
pkg/inst/doc/nlf-boot.rda
pkg/inst/doc/nlf-fit-from-truth.rda
pkg/inst/doc/nlf-fits.rda
pkg/inst/doc/nlf-lag-tests.rda
pkg/inst/doc/nlf-multi-short.rda
pkg/inst/doc/ricker-mif.rda
pkg/inst/doc/ricker-probe-match.rda
pkg/man/dacca.Rd
pkg/man/pomp-class.Rd
pkg/man/pomp-methods.Rd
pkg/man/pomp.Rd
pkg/src/sir.c
pkg/tests/ricker.Rout.save
pkg/tests/rw2.Rout.save
pkg/tests/sir-icfit.Rout.save
pkg/tests/sir.R
pkg/tests/sir.Rout.save
Log:
- add facility for parameter transformations to pomp:
- 'pomp' takes two new optional arguments: 'parameter.transform', 'parameter.inv.transform'
- 'coef' and 'coef<-' methods take new optional argument 'transform' which applies the appropriate transform when TRUE
- all data()-loadable objects are changed accordingly: 'euler.sir', 'dacca', and 'gompertz' make use of this facility
- the section on parameter transformation in the 'intro_to_pomp' vignette has been changed accordingly
- add a section on byte-compiled code to the advanced topics vignette
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2011-07-24 00:12:25 UTC (rev 527)
+++ pkg/DESCRIPTION 2011-07-28 21:19:27 UTC (rev 528)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.38-5
-Date: 2011-07-23
+Version: 0.39-1
+Date: 2011-07-28
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>
URL: http://pomp.r-forge.r-project.org
Modified: pkg/R/pomp-methods.R
===================================================================
--- pkg/R/pomp-methods.R 2011-07-24 00:12:25 UTC (rev 527)
+++ pkg/R/pomp-methods.R 2011-07-28 21:19:27 UTC (rev 528)
@@ -161,24 +161,56 @@
}
)
+pomp.transform <- function (object, params, dir = c("forward","inverse")) {
+ dir <- match.arg(dir)
+ r <- length(dim(params))
+ nm <- if (r>0) rownames(params) else names(params)
+ tfunc <- switch(
+ dir,
+ forward=function (x) do.call(object at par.trans,c(list(x),object at userdata)),
+ inverse=function (x) do.call(object at par.untrans,c(list(x),object at userdata))
+ )
+ if (r > 1)
+ retval <- apply(params,2:r,tfunc)
+ else
+ retval <- tfunc(params)
+ if (is.null(names(retval)))
+ switch(
+ dir,
+ forward=stop(
+ "invalid ",sQuote("pomp")," object: ",
+ sQuote("parameter.transform")," must return a named numeric vector"
+ ),
+ inverse=stop(
+ "invalid ",sQuote("pomp")," object: ",
+ sQuote("parameter.inv.transform")," must return a named numeric vector"
+ )
+ )
+ retval
+}
+
## extract the coefficients
setMethod(
"coef",
"pomp",
- function (object, pars, ...) {
- if (missing(pars)) {
- pars <- names(object at params)
- } else {
- excl <- !(pars%in%names(object at params))
- if (any(excl)) {
+ function (object, pars, transform = FALSE, ...) {
+ if (transform)
+ params <- pomp.transform(object,params=object at params,dir="inverse")
+ else
+ params <- object at params
+ if (missing(pars))
+ pars <- names(params)
+ else {
+ excl <- setdiff(pars,names(params))
+ if (length(excl)>0) {
stop(
"in ",sQuote("coef"),": name(s) ",
- paste(sapply(pars[excl],sQuote),collapse=","),
+ paste(sQuote(excl),collapse=","),
" correspond to no parameter(s)"
)
}
}
- object at params[pars]
+ params[pars]
}
)
@@ -186,36 +218,54 @@
setMethod(
"coef<-",
"pomp",
- function (object, pars, ..., value) {
+ function (object, pars, transform = FALSE, ..., value) {
if (missing(pars)) { ## replace the whole params slot with 'value'
+ if (transform)
+ value <- pomp.transform(object,params=value,dir="forward")
pars <- names(value)
- if (is.null(pars))
- stop("in ",sQuote("coef<-"),": ",sQuote("value")," must be a named vector")
- object at params <- numeric(length(pars))
- names(object at params) <- pars
- object at params[] <- as.numeric(value)
+ if (is.null(pars)) {
+ if (transform)
+ stop(sQuote("parameter.transform(value)")," must be a named vector")
+ else
+ stop(sQuote("value")," must be a named vector")
+ }
+ object at params <- value
} else { ## replace or append only the parameters named in 'pars'
if (!is.null(names(value))) ## we ignore the names of 'value'
- warning("in ",sQuote("coef<-"),": names of ",sQuote("value")," are being discarded",call.=FALSE)
+ warning(
+ "in ",sQuote("coef<-"),
+ " names of ",sQuote("value")," are being discarded",
+ call.=FALSE
+ )
+## if (length(pars)!=length(value))
+## stop(sQuote("pars")," and ",sQuote("value")," must be of equal length")
if (length(object at params)==0) { ## no pre-existing 'params' slot
- object at params <- numeric(length(pars))
- names(object at params) <- pars
- object at params[] <- as.numeric(value)
+ val <- numeric(length(pars))
+ names(val) <- pars
+ val[] <- value
+ if (transform)
+ value <- pomp.transform(object,params=val,dir="forward")
+ object at params <- value
} else { ## pre-existing params slot
- excl <- !(pars%in%names(object at params)) ## new parameters
- if (any(excl)) { ## append parameters
+ params <- coef(object,transform=transform)
+ val <- numeric(length(pars))
+ names(val) <- pars
+ val[] <- value
+ excl <- !(pars%in%names(params)) ## new parameter names
+ if (any(excl)) { ## append parameters
warning(
"in ",sQuote("coef<-"),": name(s) ",
paste(sQuote(pars[excl]),collapse=","),
- " are not existing parameter(s);",
+ " do not refer to existing parameter(s);",
" they are being concatenated",
call.=FALSE
)
- x <- c(object at params,numeric(length(excl)))
- names(x) <- c(names(object at params),pars[excl])
- object at params <- x
+ params <- c(params,val[excl])
}
- object at params[pars] <- as.numeric(value)
+ params[pars] <- val
+ if (transform)
+ params <- pomp.transform(object,params=params,dir="forward")
+ object at params <- params
}
}
object
@@ -258,6 +308,10 @@
}
cat("initializer = \n")
show(object at initializer)
+ cat("parameter transform function = \n")
+ show(object at par.trans)
+ cat("parameter inverse transform function = \n")
+ show(object at par.untrans)
if (length(object at userdata)>0) {
cat("userdata = \n")
show(object at userdata)
Modified: pkg/R/pomp.R
===================================================================
--- pkg/R/pomp.R 2011-07-24 00:12:25 UTC (rev 527)
+++ pkg/R/pomp.R 2011-07-28 21:19:27 UTC (rev 528)
@@ -20,6 +20,8 @@
statenames = 'character',
paramnames = 'character',
covarnames = 'character',
+ par.trans = 'function',
+ par.untrans = 'function',
PACKAGE = 'character',
userdata = 'list'
)
@@ -46,7 +48,7 @@
skeleton = NULL, skeleton.type = c("map","vectorfield"),
initializer, covar, tcovar,
obsnames, statenames, paramnames, covarnames,
- PACKAGE) {
+ PACKAGE, parameter.transform, parameter.inv.transform) {
## check the data
if (is.data.frame(data)) {
@@ -219,6 +221,41 @@
" do not embrace the data times: covariates may be extrapolated"
)
+ if (missing(parameter.transform)) {
+ if (missing(parameter.inv.transform)) {
+ has.trans <- FALSE
+ } else {
+ stop("pomp error: if ",sQuote("parameter.inv.transform")," is supplied, then " ,
+ sQuote("parameter.transform")," must also be supplied")
+ }
+ } else {
+ if (missing(parameter.inv.transform)) {
+ stop("pomp error: if ",sQuote("parameter.transform")," is supplied, then " ,
+ sQuote("parameter.inv.transform")," must also be supplied")
+ } else {
+ has.trans <- TRUE
+ }
+ }
+ if (has.trans) {
+ par.trans <- match.fun(parameter.transform)
+ par.untrans <- match.fun(parameter.inv.transform)
+ if (!all(c('params','...')%in%names(formals(par.trans))))
+ stop(
+ "pomp error: ",sQuote("parameter.transform")," must be a function of prototype ",
+ sQuote("parameter.transform(params,...)"),
+ call.=TRUE
+ )
+ if (!all(c('params','...')%in%names(formals(par.untrans))))
+ stop(
+ "pomp error: ",sQuote("parameter.inv.transform")," must be a function of prototype ",
+ sQuote("parameter.inv.transform(params,...)"),
+ call.=TRUE
+ )
+ } else {
+ par.untrans <- par.trans <- function(params, ...) params
+ }
+
+
new(
'pomp',
rprocess = rprocess,
@@ -237,6 +274,8 @@
statenames = statenames,
paramnames = paramnames,
covarnames = covarnames,
+ par.trans = par.trans,
+ par.untrans = par.untrans,
PACKAGE = PACKAGE,
userdata = list(...)
)
@@ -366,7 +405,7 @@
skeleton.map = NULL, skeleton.vectorfield = NULL,
initializer, covar, tcovar,
obsnames, statenames, paramnames, covarnames,
- PACKAGE) {
+ PACKAGE, parameter.transform, parameter.inv.transform) {
skel <- skeleton.jigger(
skeleton=skeleton,
skeleton.type=skeleton.type,
@@ -392,6 +431,8 @@
paramnames=paramnames,
covarnames=covarnames,
PACKAGE=PACKAGE,
+ parameter.transform=parameter.transform,
+ parameter.inv.transform=parameter.inv.transform,
...
)
}
@@ -406,7 +447,7 @@
skeleton.map = NULL, skeleton.vectorfield = NULL,
initializer, covar, tcovar,
obsnames, statenames, paramnames, covarnames,
- PACKAGE) {
+ PACKAGE, parameter.transform, parameter.inv.transform) {
skel <- skeleton.jigger(
skeleton=skeleton,
skeleton.type=skeleton.type,
@@ -432,6 +473,8 @@
paramnames=paramnames,
covarnames=covarnames,
PACKAGE=PACKAGE,
+ parameter.transform=parameter.transform,
+ parameter.inv.transform=parameter.inv.transform,
...
)
}
@@ -447,7 +490,7 @@
skeleton.map = NULL, skeleton.vectorfield = NULL,
initializer, covar, tcovar,
obsnames, statenames, paramnames, covarnames,
- PACKAGE) {
+ PACKAGE, parameter.transform, parameter.inv.transform) {
skel <- skeleton.jigger(
skeleton=skeleton,
skeleton.type=skeleton.type,
@@ -473,6 +516,8 @@
paramnames=paramnames,
covarnames=covarnames,
PACKAGE=PACKAGE,
+ parameter.transform=parameter.transform,
+ parameter.inv.transform=parameter.inv.transform,
...
)
}
@@ -486,7 +531,7 @@
skeleton, skeleton.type,
initializer, covar, tcovar,
obsnames, statenames, paramnames, covarnames,
- PACKAGE) {
+ PACKAGE, parameter.transform, parameter.inv.transform) {
mmg <- !missing(measurement.model)
dmg <- !missing(dmeasure)
rmg <- !missing(rmeasure)
@@ -517,6 +562,25 @@
if (missing(PACKAGE)) PACKAGE <- data at PACKAGE
if (missing(skeleton.type)) skeleton.type <- data at skeleton.type
if (missing(skeleton)) skeleton <- data at skeleton
+
+ if (missing(parameter.transform)) {
+ if (missing(parameter.inv.transform)) {
+ par.trans <- data at par.trans
+ par.untrans <- data at par.untrans
+ } else {
+ stop("pomp error: if ",sQuote("parameter.inv.transform")," is supplied, then " ,
+ sQuote("parameter.transform")," must also be supplied")
+ }
+ } else {
+ if (missing(parameter.inv.transform)) {
+ stop("pomp error: if ",sQuote("parameter.transform")," is supplied, then " ,
+ sQuote("parameter.inv.transform")," must also be supplied")
+ } else {
+ par.trans <- match.fun(parameter.transform)
+ par.untrans <- match.fun(parameter.inv.transform)
+ }
+ }
+
userdata <- data at userdata
added.userdata <- list(...)
userdata[names(added.userdata)] <- added.userdata
@@ -540,7 +604,9 @@
statenames=statenames,
paramnames=paramnames,
covarnames=covarnames,
- PACKAGE=PACKAGE
+ PACKAGE=PACKAGE,
+ parameter.transform=par.trans,
+ parameter.inv.transform=par.untrans
),
userdata
)
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/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2011-07-24 00:12:25 UTC (rev 527)
+++ pkg/inst/ChangeLog 2011-07-28 21:19:27 UTC (rev 528)
@@ -1,5 +1,6 @@
2011-07-24 kingaa
+ * [r527] inst/ChangeLog: - update ChangeLog
* [r526] DESCRIPTION, src/SSA.f, tests/fhn.Rout.save,
tests/filtfail.Rout.save, tests/gillespie.Rout.save,
tests/logistic.Rout.save, tests/ou2-bsmc.Rout.save,
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2011-07-24 00:12:25 UTC (rev 527)
+++ pkg/inst/NEWS 2011-07-28 21:19:27 UTC (rev 528)
@@ -1,4 +1,10 @@
NEWS
+0.39-1
+ o New facilities for parameter transformation are provided.
+ New optional arguments 'parameter.transform' and 'parameter.inv.transform' to 'pomp' allow the user to specify transformations between natural and internal parameter scalings.
+ A new option 'transfom' to 'coef' and 'coef<-' allows the user to get and set parameters on the natural scale.
+ See 'pomp?coef' and the "intro_to_pomp" vignette for details and examples.
+
0.38-5
o 'pfilter' will now give an error if ever a non-finite likelihood (dmeasure value) is returned.
Before, errors were only generated when dmeasure returned NA.
Modified: pkg/inst/data-R/dacca.R
===================================================================
--- pkg/inst/data-R/dacca.R 2011-07-24 00:12:25 UTC (rev 527)
+++ pkg/inst/data-R/dacca.R 2011-07-28 21:19:27 UTC (rev 528)
@@ -123,25 +123,23 @@
"omega1","sd.beta","beta.trend","log.beta1",
"alpha","rho","clin","nbasis","nrstage"),
covarnames = c("pop","dpopdt","seas.1","trend"),
- ## log.trans=c( # parameters to log transform
- ## "gamma","eps","rho","delta","deltaI","alpha",
- ## "tau","sd.beta",
- ## paste("omega",seq(length=nbasis),sep=''),
- ## "S.0","I.0","Rs.0","R1.0","R2.0","R3.0"
- ## ),
- ## logit.trans="clin", # parameters to logit transform
- ## transform=function (params, log.trans, logit.trans, ...) {
- ## logit <- function(p){log(p/(1-p))} # (0,1) -> (-inf,inf)
- ## params[logit.trans] <- logit(params[logit.trans])
- ## params[log.trans] <- log(params[log.trans])
- ## params
- ## },
- ## inv.transform=function (params, log.trans, logit.trans, ...) {
- ## expit <- function(r){1/(1+exp(-r))} # (-inf,inf) -> (0,1)
- ## params[logit.trans] <- expit(params[logit.trans])
- ## params[log.trans] <- exp(params[log.trans])
- ## params
- ## },
+ log.trans=c( # parameters to log transform
+ "gamma","eps","rho","delta","deltaI","alpha",
+ "tau","sd.beta",
+ paste("omega",seq(length=nbasis),sep=''),
+ "S.0","I.0","Rs.0","R1.0","R2.0","R3.0"
+ ),
+ logit.trans="clin", # parameters to logit transform
+ parameter.transform=function (params, log.trans, logit.trans, ...) {
+ params[logit.trans] <- qlogis(params[logit.trans])
+ params[log.trans] <- log(params[log.trans])
+ params
+ },
+ parameter.inv.transform=function (params, log.trans, logit.trans, ...) {
+ params[logit.trans] <- plogis(params[logit.trans])
+ params[log.trans] <- exp(params[log.trans])
+ params
+ },
initializer = function (params, t0, covars, all.state.names, comp.names, nrstage, ...) {
all.state.names <- c("S","I","Rs","R1","R2","R3","M","W","count")
comp.names <- c("S","I","Rs",paste("R",1:nrstage,sep=''))
@@ -190,4 +188,4 @@
}
coef(dacca) <- dacca.transform(mle,dir="forward")
-save(dacca,dacca.transform,file="dacca.rda")
+save(dacca,file="dacca.rda")
Modified: pkg/inst/data-R/euler.sir.R
===================================================================
--- pkg/inst/data-R/euler.sir.R 2011-07-24 00:12:25 UTC (rev 527)
+++ pkg/inst/data-R/euler.sir.R 2011-07-28 21:19:27 UTC (rev 528)
@@ -1,57 +1,68 @@
require(pomp)
-simulate(
- pomp(
- data=data.frame(
- time=seq(from=1/52,to=4,by=1/52),
- reports=NA
- ),
- times="time",
- t0=0,
- rprocess=euler.sim(
- step.fun="_sir_euler_simulator",
- delta.t=1/52/20,
- PACKAGE="pomp"
- ),
- skeleton.type="vectorfield",
- skeleton="_sir_ODE",
- rmeasure="_sir_binom_rmeasure",
- dmeasure="_sir_binom_dmeasure",
- PACKAGE="pomp",
- obsnames = c("reports"),
- 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"),
- 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))
- ## 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
- }
- ),
- params=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)
- ),
- nsim=1,
- seed=329348545L
- ) -> euler.sir
+po <- pomp(
+ data=data.frame(
+ time=seq(from=1/52,to=4,by=1/52),
+ reports=NA
+ ),
+ times="time",
+ t0=0,
+ rprocess=euler.sim(
+ step.fun="_sir_euler_simulator",
+ delta.t=1/52/20,
+ PACKAGE="pomp"
+ ),
+ skeleton.type="vectorfield",
+ skeleton="_sir_ODE",
+ rmeasure="_sir_binom_rmeasure",
+ dmeasure="_sir_binom_dmeasure",
+ PACKAGE="pomp",
+ obsnames = c("reports"),
+ 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"),
+ 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, ...) {
+ ic.names <- paste(comp.names,".0",sep="")
+ 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,
+ 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
+
save(euler.sir,file="euler.sir.rda")
Modified: pkg/inst/data-R/gompertz.R
===================================================================
--- pkg/inst/data-R/gompertz.R 2011-07-24 00:12:25 UTC (rev 527)
+++ pkg/inst/data-R/gompertz.R 2011-07-28 21:19:27 UTC (rev 528)
@@ -1,30 +1,33 @@
require(pomp)
-simulate(
- pomp(
- data=data.frame(time=seq(0,100,by=1),Y=NA),
- times="time",
- t0=0,
- rprocess=discrete.time.sim(
- step.fun="_gompertz_simulator"
- ),
- rmeasure="_gompertz_normal_rmeasure",
- dmeasure="_gompertz_normal_dmeasure",
- skeleton.type="map",
- skeleton="_gompertz_skeleton",
- paramnames=c("log.r","log.K","log.sigma","log.tau"),
- statenames=c("X"),
- obsnames=c("Y")
- ),
- params=c(
- log.K=log(1),
- log.r=log(0.1),
- log.sigma=log(0.1),
- log.tau=log(0.1),
- X.0=1
- ),
- nsim=1,
- seed=299438676L
- ) -> gompertz
+po <- pomp(
+ data=data.frame(time=seq(0,100,by=1),Y=NA),
+ times="time",
+ t0=0,
+ rprocess=discrete.time.sim(
+ step.fun="_gompertz_simulator"
+ ),
+ rmeasure="_gompertz_normal_rmeasure",
+ dmeasure="_gompertz_normal_dmeasure",
+ skeleton.type="map",
+ skeleton="_gompertz_skeleton",
+ paramnames=c("log.r","log.K","log.sigma","log.tau"),
+ statenames=c("X"),
+ 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
+ }
+ )
+coef(po,transform=TRUE) <- c(K=1,r=0.1,sigma=0.1,tau=0.1,X.0=1)
+
+simulate(po,nsim=1,seed=299438676L) -> gompertz
+
save(gompertz,file="gompertz.rda")
Modified: pkg/inst/doc/advanced_topics_in_pomp.Rnw
===================================================================
--- pkg/inst/doc/advanced_topics_in_pomp.Rnw 2011-07-24 00:12:25 UTC (rev 527)
+++ pkg/inst/doc/advanced_topics_in_pomp.Rnw 2011-07-28 21:19:27 UTC (rev 528)
@@ -205,6 +205,57 @@
Doing \Sexpr{dim(simdat.Rvect)[2]} simulations of \code{ou2.Rvect} took \Sexpr{round(etime.Rvect,2)}~\Sexpr{units(etime.Rvect)}.
Compared to the \Sexpr{round(etime.Rplug,2)}~\Sexpr{units(etime.Rplug)} it took to run \Sexpr{dim(simdat.Rplug)[2]} simulations of \code{ou2.Rplug}, this is a \Sexpr{round(as.numeric(etime.Rplug)/as.numeric(etime.Rvect))}-fold speed-up.
+\subsection{Using \R's byte compiler}
+
+From version {2.13}, \R\ has provided byte-compilation facilities, in the base package \code{compiler}.
+Let's see to what extent we can speed up our codes by byte-compiling the components of our \code{pomp} object.
+<<plugin-byte-code,echo=T>>=
+require(compiler)
+pomp(
+ ou2.Rplug,
+ rprocess=discrete.time.sim(
+ step.fun=cmpfun(
+ function (x, t, params, ...) {
+ eps <- rnorm(n=2,mean=0,sd=1) # noise terms
+ xnew <- c(
+ x1=params["alpha.1"]*x["x1"]+params["alpha.3"]*x["x2"]+
+ params["sigma.1"]*eps[1],
+ x2=params["alpha.2"]*x["x1"]+params["alpha.4"]*x["x2"]+
+ params["sigma.2"]*eps[1]+params["sigma.3"]*eps[2]
+ )
+ names(xnew) <- c("x1","x2")
+ xnew
+ },
+ options=list(optimize=3)
+ )
+ )
+ ) -> ou2.Bplug
+@
+<<plugin-byte-code-sim,echo=F>>=
+tic <- Sys.time()
+simdat.Bplug <- simulate(ou2.Bplug,params=coef(ou2),nsim=1000,states=T)
+toc <- Sys.time()
+etime.Bplug <- toc-tic
+@
+Doing these \Sexpr{dim(simdat.Bplug)[2]} simulations of \code{ou2.Bplug} took \Sexpr{round(etime.Bplug,2)}~\Sexpr{units(etime.Bplug)}.
+This is a \Sexpr{round(as.numeric(etime.Rplug)/as.numeric(etime.Bplug),1)}-fold speed-up relative to the plug-in code written in \R.
+
+We can byte-compile the vectorized \R\ code, too, and compare its performance:
+<<vectorized-byte-code>>=
+ou2.Bvect <- pomp(ou2.Rplug,rprocess=cmpfun(ou2.Rvect.rprocess,options=list(optimize=3)))
+@
+<<vectorized-byte-code-sim>>=
+tic <- Sys.time()
+simdat.Bvect <- simulate(ou2.Bvect,params=theta,states=T,nsim=1000)
+toc <- Sys.time()
+etime.Bvect <- toc-tic
+@
+<<echo=F>>=
+units(etime.Bvect) <- units(etime.Rplug)
+@
+This code shows a \Sexpr{round(as.numeric(etime.Rplug)/as.numeric(etime.Bvect))}-fold speed-up relative to the plug-in code written in \R.
+
+
\subsection{Accelerating the code using C: a plug-in based implementation}
As we've seen, we can usually acheive big accelerations using compiled native code.
@@ -342,17 +393,33 @@
),
## reset cases to zero at 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 <- round(c(pop*fracs/sum(fracs),0,0))
- names(x0) <- c("S","I","R","cases","W")
- x0
- }
- )
+ comp.names=c("S","I","R"),
+ 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, ...) {
+ ic.names <- paste(comp.names,".0",sep="")
+ 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
}
) -> sir
@
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/gompertz-multi-mif.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/gompertz-trajmatch.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.Rnw
===================================================================
--- pkg/inst/doc/intro_to_pomp.Rnw 2011-07-24 00:12:25 UTC (rev 527)
+++ pkg/inst/doc/intro_to_pomp.Rnw 2011-07-28 21:19:27 UTC (rev 528)
@@ -494,6 +494,7 @@
coef(gompertz)
coef(gompertz,c("sigma","tau")) <- c(1,0)
@
+See below for more information on the \code{coef} and \code{coef<-} methods for getting and setting parameters.
Finally, one has access to the unobserved states via, e.g.,
<<eval=F>>=
states(gompertz)
@@ -508,6 +509,7 @@
The parameters in the Gompertz model above are constrained to be positive.
When we estimate these parameters using numerical search algorithms, we must find some way to ensure that these constraints will be honored.
A straightforward way to accomplish this is to transform the parameters so that they become unconstrained.
+To introduce some terminology, we want to transform the parameters from the \emph{natural scale} ($r$, $K$, $\sigma$, $\tau$), to another, \emph{internal scale}, on which they will have some desirable property, e.g., they will be unconstrained.
The following codes re-implement the Gompertz model using transformed parameters.
<<loggomp-proc-sim-def>>=
@@ -553,26 +555,71 @@
}
@
-Note that, in each of the above functions, we \emph{untransform} the parameters before we do any computations.
+Note that, in each of the above functions, we \emph{untransform} the parameters back onto the natural scale before we do any computations.
-<<loggomp-pomp-construction,eval=F>>=
+\pkg{pomp} provides a facility to make it easier to work with parameter transformations.
+Specifically, we can specify the optional functions \code{parameter.transform} and \code{parameter.inv.transform} when we create the \code{pomp} object.
+The first function will take transform our parameters to the internal scale;
+the second one will invert that operation, transforming the parameters from the internal back to the natural scale.
+
+<<loggomp-pomp-construction,eval=T>>=
dat <- as(gompertz,"data.frame")
-theta <- c(
- log(coef(gompertz,c("r","K","tau","sigma"))),
- coef(gompertz,"X.0")
- )
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 528
More information about the pomp-commits
mailing list