[Pomp-commits] r725 - in pkg/pomp: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue May 29 20:12:41 CEST 2012
Author: kingaa
Date: 2012-05-29 20:12:40 +0200 (Tue, 29 May 2012)
New Revision: 725
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/R/builder.R
pkg/pomp/R/pomp.R
pkg/pomp/demo/gompertz.R
pkg/pomp/demo/sir.R
pkg/pomp/man/builder.Rd
pkg/pomp/man/pomp.Rd
Log:
- add support for parameter transformations to 'pompBuilder'
- more error-trapping in 'pomp'
- update demos to use compiled parameter transformations
- fix minor omission in 'pomp' man page
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2012-05-15 16:28:37 UTC (rev 724)
+++ pkg/pomp/DESCRIPTION 2012-05-29 18:12:40 UTC (rev 725)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.42-4
-Date: 2012-05-14
+Version: 0.42-5
+Date: 2012-05-29
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/pomp/R/builder.R
===================================================================
--- pkg/pomp/R/builder.R 2012-05-15 16:28:37 UTC (rev 724)
+++ pkg/pomp/R/builder.R 2012-05-29 18:12:40 UTC (rev 725)
@@ -22,8 +22,9 @@
pompBuilder <- function (data, times, t0, name,
statenames, paramnames,
rmeasure, dmeasure, step.fn, step.fn.delta.t,
- skeleton, skeleton.type, skelmap.delta.t = 1, ...,
- link = TRUE) {
+ skeleton, skeleton.type, skelmap.delta.t = 1,
+ parameter.transform, parameter.inv.transform,
+ ..., link = TRUE) {
obsnames <- names(data)
obsnames <- setdiff(obsnames,times)
solib <- pompCBuilder(
@@ -34,7 +35,9 @@
rmeasure=rmeasure,
dmeasure=dmeasure,
step.fn=step.fn,
- skeleton=skeleton
+ skeleton=skeleton,
+ parameter.transform=parameter.transform,
+ parameter.inv.transform=parameter.inv.transform
)
if (link) pompLink(name)
pomp(
@@ -49,6 +52,8 @@
skeleton=render("{%name%}_skelfn",name=name),
skeleton.type=skeleton.type,
skelmap.delta.t=skelmap.delta.t,
+ parameter.transform=render("{%name%}_par_trans",name=name),
+ parameter.inv.transform=render("{%name%}_par_untrans",name=name),
PACKAGE=name,
obsnames=obsnames,
statenames=statenames,
@@ -77,34 +82,52 @@
)
header <- list(
- file="/* pomp model file: {%name%} */\n\n#include <pomp.h>\n#include <R_ext/Rdynload.h>\n",
+ file="/* pomp model file: {%name%} */\n\n#include <pomp.h>\n#include <R_ext/Rdynload.h>\n\n",
rmeasure="\nvoid {%name%}_rmeasure (double *__y, double *__x, double *__p, int *__obsindex, int *__stateindex, int *__parindex, int *__covindex, int __ncovars, double *__covars, double t)\n{\n",
dmeasure= "\nvoid {%name%}_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)\n{\n",
step.fn="\nvoid {%name%}_stepfn (double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, int __covdim, const double *__covar, double t, double dt)\n{\n",
- skeleton="\nvoid {%name%}_skelfn (double *__f, double *__x, double *__p, int *__stateindex, int *__parindex, int *__covindex, int __ncovars, double *__covars, double t)\n{\n"
+ skeleton="\nvoid {%name%}_skelfn (double *__f, double *__x, double *__p, int *__stateindex, int *__parindex, int *__covindex, int __ncovars, double *__covars, double t)\n{\n",
+ parameter.transform="\nvoid {%name%}_par_trans (double *__pt, double *__p, int *__parindex)\n{\n",
+ parameter.inv.transform="\nvoid {%name%}_par_untrans (double *__pt, double *__p, int *__parindex)\n{\n"
)
decl <- list(
- periodic_bspline_basis_eval="void (*periodic_bspline_basis_eval)(double,double,int,int,double*);\nperiodic_bspline_basis_eval = (void (*)(double,double,int,int,double*)) R_GetCCallable(\"pomp\",\"periodic_bspline_basis_eval\");\n",
- reulermultinom="void (*reulermultinom)(int,double,double*,double,double*);\nreulermultinom = (void (*)(int,double,double*,double,double*)) R_GetCCallable(\"pomp\",\"reulermultinom\");\n",
- get_pomp_userdata="const SEXP (*get_pomp_userdata)(const char *);\npomp_get_userdata = (const SEXP (*)(const char*)) R_GetCCallable(\"pomp\",\"get_pomp_userdata\");\n",
- get_pomp_userdata_int="const int * (*get_pomp_userdata_int)(const char *);\npomp_get_userdata_int = (const int *(*)(const char*)) R_GetCCallable(\"pomp\",\"get_pomp_userdata_int\");\n",
- get_pomp_userdata_double="const double * (*get_pomp_userdata_double)(const char *);\npomp_get_userdata_double = (const double *(*)(const char*)) R_GetCCallable(\"pomp\",\"get_pomp_userdata_double\");\n"
- )
+ periodic_bspline_basis_eval="\tvoid (*periodic_bspline_basis_eval)(double,double,int,int,double*);\nperiodic_bspline_basis_eval = (void (*)(double,double,int,int,double*)) R_GetCCallable(\"pomp\",\"periodic_bspline_basis_eval\");\n",
+ reulermultinom="\tvoid (*reulermultinom)(int,double,double*,double,double*);\nreulermultinom = (void (*)(int,double,double*,double,double*)) R_GetCCallable(\"pomp\",\"reulermultinom\");\n",
+ get_pomp_userdata="\tconst SEXP (*get_pomp_userdata)(const char *);\npomp_get_userdata = (const SEXP (*)(const char*)) R_GetCCallable(\"pomp\",\"get_pomp_userdata\");\n",
+ get_pomp_userdata_int="\tconst int * (*get_pomp_userdata_int)(const char *);\npomp_get_userdata_int = (const int *(*)(const char*)) R_GetCCallable(\"pomp\",\"get_pomp_userdata_int\");\n",
+ get_pomp_userdata_double="\tconst double * (*get_pomp_userdata_double)(const char *);\npomp_get_userdata_double = (const double *(*)(const char*)) R_GetCCallable(\"pomp\",\"get_pomp_userdata_double\");\n",
+ logit="\tdouble logit(double p){return log(p/(1.0-p));}\n\n",
+ expit="\tdouble expit(double x){return 1.0/(1.0+exp(-x));}\n\n"
+ )
footer <- list(
rmeasure="\n}\n\n",
dmeasure="\n}\n\n",
step.fn="\n}\n\n",
- skeleton="\n}\n\n"
+ skeleton="\n}\n\n",
+ parameter.transform="\n}\n\n",
+ parameter.inv.transform="\n}\n\n"
)
-pompCBuilder <- function (name, statenames, paramnames, obsnames, rmeasure, dmeasure, step.fn, skeleton)
+utility.fns <- list(
+ )
+
+
+pompCBuilder <- function (name, statenames, paramnames, obsnames, rmeasure, dmeasure,
+ step.fn, skeleton, parameter.transform, parameter.inv.transform)
{
if (missing(name)) stop(sQuote("name")," must be supplied");
if (missing(statenames)) stop(sQuote("name")," must be supplied");
if (missing(paramnames)) stop(sQuote("name")," must be supplied");
if (missing(obsnames)) stop(sQuote("name")," must be supplied");
+
+ mpt <- missing(parameter.transform)
+ mpit <- missing(parameter.inv.transform)
+ if (xor(mpt,mpit))
+ stop("if you supply one transformation function, you must supply its inverse")
+ has.trans <- !mpt
+
name <- cleanForC(name)
statenames <- cleanForC(statenames)
paramnames <- cleanForC(paramnames)
@@ -116,6 +139,11 @@
out <- file(description=modelfile,open="w")
cat(file=out,render(header$file,name=name))
+
+ for (f in utility.fns) {
+ cat(file=out,f)
+ }
+
## variable/parameter/observations definitions
for (v in seq_along(paramnames)) {
cat(file=out,render(define$var,variable=paramnames[v],ptr='__p',ilist='__parindex',index=v-1))
@@ -129,8 +157,28 @@
for (v in seq_along(statenames)) {
cat(file=out,render(define$var,variable=paste0("D",statenames[v]),ptr='__f',ilist='__stateindex',index=v-1))
}
+ for (v in seq_along(paramnames)) {
+ cat(file=out,render(define$var,variable=paste0("T",paramnames[v]),ptr='__pt',ilist='__parindex',index=v-1))
+ }
cat(file=out,render(define$var.alt,variable="lik",ptr='__lik',index=0))
+ if (has.trans) {
+ ## parameter transformation function
+ cat(file=out,render(header$parameter.transform,name=name))
+ for (fn in names(decl)) {
+ if (grepl(fn,parameter.transform))
+ cat(file=out,decl[[fn]])
+ }
+ cat(file=out,parameter.transform,footer$parameter.transform)
+ ## inverse parameter transformation function
+ cat(file=out,render(header$parameter.inv.transform,name=name))
+ for (fn in names(decl)) {
+ if (grepl(fn,parameter.inv.transform))
+ cat(file=out,decl[[fn]])
+ }
+ cat(file=out,parameter.inv.transform,footer$parameter.inv.transform)
+ }
+
## rmeasure function
cat(file=out,render(header$rmeasure,name=name),rmeasure,footer$rmeasure)
@@ -166,6 +214,9 @@
for (v in seq_along(statenames)) {
cat(file=out,render(undefine$var,variable=paste0("D",statenames[v])))
}
+ for (v in seq_along(paramnames)) {
+ cat(file=out,render(undefine$var,variable=paste0("T",paramnames[v])))
+ }
close(out)
cflags <- paste0("PKG_CFLAGS=\"",
Modified: pkg/pomp/R/pomp.R
===================================================================
--- pkg/pomp/R/pomp.R 2012-05-15 16:28:37 UTC (rev 724)
+++ pkg/pomp/R/pomp.R 2012-05-29 18:12:40 UTC (rev 725)
@@ -54,6 +54,10 @@
obsnames, statenames, paramnames, covarnames, zeronames,
PACKAGE, parameter.transform, parameter.inv.transform) {
+ if (missing(data)) stop(sQuote("data")," is a required argument")
+ if (missing(times)) stop(sQuote("times")," is a required argument")
+ if (missing(t0)) stop(sQuote("t0")," is a required argument")
+
## check the data
if (is.data.frame(data)) {
if (!is.character(times) || length(times)!=1 || !(times%in%names(data)))
@@ -231,21 +235,14 @@
" 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
- }
+ mpt <- missing(parameter.transform)
+ mpit <- missing(parameter.inv.transform)
+ if (xor(mpt,mpit)) {
+ stop("if one of ",sQuote("parameter.transform"),", ",
+ sQuote("parameter.inv.transform"),
+ " is supplied, then so must the other")
}
+ has.trans <- !mpt
if (has.trans) {
par.trans <- pomp.fun(
f=parameter.transform,
@@ -261,12 +258,7 @@
par.trans <- pomp.fun()
par.untrans <- pomp.fun()
}
- if (
- has.trans &&
- par.trans at mode<0 &&
- par.untrans at mode<0
- )
- has.trans <- FALSE
+ if (has.trans && par.trans at mode<0 && par.untrans at mode<0) has.trans <- FALSE
new(
'pomp',
Modified: pkg/pomp/demo/gompertz.R
===================================================================
--- pkg/pomp/demo/gompertz.R 2012-05-15 16:28:37 UTC (rev 724)
+++ pkg/pomp/demo/gompertz.R 2012-05-29 18:12:40 UTC (rev 725)
@@ -58,7 +58,6 @@
"
rmeas <- "
Y = rlnorm(log(X),tau);
-
"
step.fn <- "
double S = exp(-r*dt);
@@ -72,6 +71,20 @@
/* note that X is not over-written in the skeleton function */
DX = pow(K,1-S)*pow(X,S);
"
+partrans <- "
+ Tr = exp(r);
+ TK = exp(K);
+ Tsigma = exp(sigma);
+ TX_0 = exp(X_0);
+ Ttau = exp(tau);
+"
+paruntrans <- "
+ Tr = log(r);
+ TK = log(K);
+ Tsigma = log(sigma);
+ TX_0 = log(X_0);
+ Ttau = log(tau);
+"
pompBuilder(
name="Gompertz",
@@ -87,12 +100,8 @@
skeleton=skel,
skeleton.type="map",
skelmap.delta.t=1,
- parameter.inv.transform=function(params,...){
- log(params)
- },
- parameter.transform=function(params,...){
- exp(params)
- }
+ parameter.inv.transform=partrans,
+ parameter.transform=paruntrans
) -> Gompertz
## simulate some data
Modified: pkg/pomp/demo/sir.R
===================================================================
--- pkg/pomp/demo/sir.R 2012-05-15 16:28:37 UTC (rev 724)
+++ pkg/pomp/demo/sir.R 2012-05-29 18:12:40 UTC (rev 725)
@@ -1,10 +1,10 @@
require(pomp)
dmeas <- "
- lik = dbinom(reports,nearbyint(cases),rho,give_log);
+ lik = dbinom(cases,nearbyint(incid),rho,give_log);
"
rmeas <- "
- reports = rbinom(nearbyint(cases),rho);
+ cases = rbinom(nearbyint(incid),rho);
"
step.fn <- '
int nrate = 6;
@@ -20,8 +20,7 @@
// compute transmission rate from seasonality
periodic_bspline_basis_eval(t,period,deg,nbasis,&seasonality[0]); // evaluate the periodic B-spline basis
- for (k = 0, beta = 0; k < nbasis; k++)
- beta += (&beta1)[k]*seasonality[k];
+ beta = beta1*seasonality[0]+beta2*seasonality[1]+beta3*seasonality[2];
// compute the environmental stochasticity
dW = rgammawn(beta_sd,dt);
@@ -44,7 +43,7 @@
S += trans[0]-trans[1]-trans[2];
I += trans[1]-trans[3]-trans[4];
R += trans[3]-trans[5];
- cases += trans[3]; // cases are cumulative recoveries
+ incid += trans[3]; // cases are cumulative recoveries
if (beta_sd > 0.0) W += (dW-dt)/beta_sd; // increment has mean = 0, variance = dt
'
skel <- '
@@ -62,8 +61,7 @@
// compute transmission rate from seasonality
if (nbasis <= 0) return;
periodic_bspline_basis_eval(t,period,deg,nbasis,&seasonality[0]); // evaluate the periodic B-spline basis
- for (k = 0, beta = 0; k < nbasis; k++)
- beta += (&beta1)[k]*seasonality[k];
+ beta = beta1*seasonality[0]+beta2*seasonality[1]+beta3*seasonality[2];
// compute the transition rates
rate[0] = mu*popsize; // birth into susceptible class
@@ -85,52 +83,73 @@
DS = term[0]-term[1]-term[2];
DI = term[1]-term[3]-term[4];
DR = term[3]-term[5];
- Dcases = term[3]; // accumulate the new I->R transitions
+ Dincid = term[3]; // accumulate the new I->R transitions
DW = 0;
'
+partrans <- "
+ double sum;
+ Tgamma = exp(gamma);
+ Tmu = exp(mu);
+ Tiota = exp(iota);
+ Tbeta1 = exp(beta1);
+ Tbeta2 = exp(beta2);
+ Tbeta3 = exp(beta3);
+ Tbeta_sd = exp(beta_sd);
+ Trho = expit(rho);
+ TS_0 = exp(S_0);
+ TI_0 = exp(I_0);
+ TR_0 = exp(R_0);
+ sum = TS_0+TI_0+TR_0;
+ TS_0 /= sum;
+ TI_0 /= sum;
+ TR_0 /= sum;
+"
+paruntrans <- "
+ double sum;
+ Tgamma = log(gamma);
+ Tmu = log(mu);
+ Tiota = log(iota);
+ Tbeta1 = log(beta1);
+ Tbeta2 = log(beta2);
+ Tbeta3 = log(beta3);
+ Tbeta_sd = log(beta_sd);
+ Trho = logit(rho);
+ sum = S_0+I_0+R_0;
+ TS_0 = log(S_0/sum);
+ TI_0 = log(I_0/sum);
+ TR_0 = log(R_0/sum);
+"
+data(LondonYorke)
+
pompBuilder(
name="SIR",
- data=data.frame(
- time=seq(from=1/52,to=4,by=1/52),
- reports=NA
- ),
+ data=subset(
+ LondonYorke,
+ subset=town=="New York"&disease=="measles"&year>=1928&year<=1933,
+ select=c(time,cases)
+ ),
times="time",
- t0=0,
+ t0=1928,
dmeasure=dmeas,
rmeasure=rmeas,
step.fn=step.fn,
step.fn.delta.t=1/52/20,
skeleton.type="vectorfield",
skeleton=skel,
- statenames=c("S","I","R","cases","W"),
+ parameter.transform=partrans,
+ parameter.inv.transform=paruntrans,
+ statenames=c("S","I","R","incid","W"),
paramnames=c(
- "gamma","mu","iota","beta1","beta.sd","popsize","rho"
+ "gamma","mu","iota","beta1","beta2","beta3","beta.sd",
+ "popsize","rho","S.0","I.0","R.0"
),
- zeronames=c("cases","W"),
- logvar=c(
- "gamma","mu","iota",
- "beta1","beta2","beta3","beta.sd",
- "S.0","I.0","R.0"
- ),
- logitvar="rho",
+ zeronames=c("incid","W"),
comp.names=c("S","I","R"),
ic.names=c("S.0","I.0","R.0"),
- parameter.transform=function (params, logvar, logitvar, ic.names, ...) {
- params[logvar] <- exp(params[logvar])
- params[logitvar] <- plogis(params[logitvar])
- params[ic.names] <- params[ic.names]/sum(params[ic.names])
- params
- },
- parameter.inv.transform=function (params, logvar, logitvar, ic.names, ...) {
- params[ic.names] <- params[ic.names]/sum(params[ic.names])
- params[logvar] <- log(params[logvar])
- params[logitvar] <- qlogis(params[logitvar])
- params
- },
initializer=function(params, t0, comp.names, ic.names, ...) {
x0 <- numeric(5)
- names(x0) <- c("S","I","R","cases","W")
+ names(x0) <- c("S","I","R","incid","W")
fracs <- params[ic.names]
x0[comp.names] <- round(params['popsize']*fracs/sum(fracs))
x0
@@ -141,9 +160,9 @@
gamma=26,mu=0.02,iota=0.01,
beta1=120,beta2=140,beta3=100,
beta.sd=1e-3,
- popsize=5e5,
- rho=0.6,
- S.0=26/120,I.0=0.001,R.0=1-26/120
+ popsize=5e6,
+ rho=0.2,
+ S.0=0.22,I.0=0.0018,R.0=0.78
)
## compute a trajectory of the deterministic skeleton
@@ -152,7 +171,7 @@
toc <- Sys.time()
print(toc-tic)
-plot(cases~time,data=X,type='l')
+plot(incid~time,data=X,type='l')
## simulate from the model
tic <- Sys.time()
@@ -160,6 +179,17 @@
toc <- Sys.time()
print(toc-tic)
-plot(cases~time,data=x,col=as.factor(x$sim),pch=16)
+plot(incid~time,data=x,col=as.factor(x$sim),pch=16)
+coef(po) <- coef(
+ traj.match(
+ pomp(
+ window(po,end=1930),
+ measurement.model=cases~norm(mean=rho*incid,sd=1000)
+ ),
+ est=c("beta1","beta2","beta3","S.0","I.0","R.0"),
+ transform=TRUE
+ )
+ )
+
dyn.unload("SIR.so")
Modified: pkg/pomp/man/builder.Rd
===================================================================
--- pkg/pomp/man/builder.Rd 2012-05-15 16:28:37 UTC (rev 724)
+++ pkg/pomp/man/builder.Rd 2012-05-29 18:12:40 UTC (rev 725)
@@ -8,6 +8,7 @@
pompBuilder(data, times, t0, name, statenames, paramnames,
rmeasure, dmeasure, step.fn, step.fn.delta.t,
skeleton, skeleton.type, skelmap.delta.t = 1,
+ parameter.transform, parameter.inv.transform,
\dots, link = TRUE)
}
\arguments{
@@ -34,6 +35,11 @@
As in \code{pomp}, \code{skeleton.type} indicates whether the skeleton is a map (discrete-time) or vectorfield (continuous-time).
If the former, \code{skelmap.delta.t} is the time-step of the map.
}
+ \item{parameter.transform, parameter.inv.transform}{
+ optional C codes that implement parameter transformations.
+ \code{parameter.transform} maps parameters from the estimation scale to the natural scale;
+ \code{parameter.inv.transformation} maps them from the natural scale to the estimation scale.
+ }
\item{\dots}{
additional arguments are passed to \code{\link{pomp}}
}
Modified: pkg/pomp/man/pomp.Rd
===================================================================
--- pkg/pomp/man/pomp.Rd 2012-05-15 16:28:37 UTC (rev 724)
+++ pkg/pomp/man/pomp.Rd 2012-05-29 18:12:40 UTC (rev 725)
@@ -113,7 +113,7 @@
\code{covar} can be specified as either a matrix or a data frame.
In either case the columns are taken to be distinct covariates.
If \code{covar} is a data frame, \code{tcovar} can be either the name or the index of the time variable.
- If a covariate table is supplied, then the value of each of the covariates is interpolated as needed, i.e., whenever \code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}, or \code{init.state} is evaluated.
+ If a covariate table is supplied, then the value of each of the covariates is interpolated as needed, i.e., whenever \code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}, \code{skeleton}, or \code{init.state} is evaluated.
The resulting interpolated values are passed to the corresponding functions as a numeric vector named \code{covars}.
}
\item{obsnames, statenames, paramnames, covarnames}{
More information about the pomp-commits
mailing list