[Pomp-commits] r598 - in pkg: R inst/include man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jan 19 13:03:04 CET 2012
Author: kingaa
Date: 2012-01-19 13:03:03 +0100 (Thu, 19 Jan 2012)
New Revision: 598
Modified:
pkg/R/eulermultinom.R
pkg/inst/include/pomp.h
pkg/man/eulermultinom.Rd
pkg/src/eulermultinom.c
pkg/src/pomp.h
Log:
- add an R interface to the 'rgammawn' function
Modified: pkg/R/eulermultinom.R
===================================================================
--- pkg/R/eulermultinom.R 2012-01-19 11:44:49 UTC (rev 597)
+++ pkg/R/eulermultinom.R 2012-01-19 12:03:03 UTC (rev 598)
@@ -3,3 +3,6 @@
deulermultinom <- function (x, size, rate, dt, log = FALSE)
.Call(D_Euler_Multinom,as.matrix(x),size,rate,dt,log)
+
+rgammawn <- function (n = 1, sigma, dt)
+ .Call(R_GammaWN,n,sigma,dt)
Modified: pkg/inst/include/pomp.h
===================================================================
--- pkg/inst/include/pomp.h 2012-01-19 11:44:49 UTC (rev 597)
+++ pkg/inst/include/pomp.h 2012-01-19 12:03:03 UTC (rev 598)
@@ -8,11 +8,30 @@
#include <Rdefines.h>
// prototypes for C-level access to Euler-multinomial distribution functions
-// NB: 'reulermultinom' does not call GetRNGstate() and PutRNGstate() internally, so the user must do so
+// NB: 'reulermultinom' does not call GetRNGstate() and PutRNGstate() internally
void reulermultinom (int ntrans, double size, double *rate, double dt, double *trans);
double deulermultinom (int ntrans, double size, double *rate, double dt, double *trans, int give_log);
-// facility for dotting a vector of parameters ('coef') against a vector of basis-function values ('basis')
+// This function computes r such that if
+// N ~ geometric(prob=1-exp(-r dt)) and T ~ exponential(rate=R),
+// then E[N dt] = E[T]
+// i.e., the rate r for an Euler process that gives the same
+// expected waiting time as the exponential process it approximates.
+// In particular r -> R as dt -> 0.
+inline double exp2geom_rate_correction (double R, double dt) {
+ return (dt > 0) ? log(1.0+R*dt)/dt : R;
+}
+
+// This function draws a random increment of a gamma whitenoise process.
+// This will have expectation=dt and variance=(sigma^2*dt)
+// If dW = rgammawn(sigma,dt), then
+// mu dW/dt is a candidate for a random rate process within an
+// Euler-multinomial context, i.e.,
+// E[mu*dW] = mu*dt and Var[mu*dW] = mu*sigma^2*dt
+double rgammawn (double sigma, double dt);
+
+// facility for computing the inner produce of
+// a vector of parameters ('coef') against a vector of basis-function values ('basis')
double dot_product (int dim, const double *basis, const double *coef);
// facility for computing evaluating a basis of periodic bsplines
Modified: pkg/man/eulermultinom.Rd
===================================================================
--- pkg/man/eulermultinom.Rd 2012-01-19 11:44:49 UTC (rev 597)
+++ pkg/man/eulermultinom.Rd 2012-01-19 12:03:03 UTC (rev 598)
@@ -2,6 +2,7 @@
\alias{eulermultinom}
\alias{reulermultinom}
\alias{deulermultinom}
+\alias{rgammawn}
\title{Euler-multinomial death process}
\description{
Density and random-deviate generation for the Euler-multinomial death process with parameters \code{size}, \code{rate}, and \code{dt}.
@@ -9,11 +10,13 @@
\usage{
reulermultinom(n = 1, size, rate, dt)
deulermultinom(x, size, rate, dt, log = FALSE)
+rgammawn(n = 1, sigma, dt)
}
\arguments{
\item{n}{integer; number of random variates to generate.}
\item{size}{scalar integer; number of individuals at risk.}
\item{rate}{numeric vector of hazard rates.}
+ \item{sigma}{numeric scalar; intensity of the Gamma white noise process.}
\item{dt}{numeric scalar; duration of Euler step.}
\item{x}{matrix or vector containing number of individuals that have succumbed to each death process.}
\item{log}{logical; if TRUE, return logarithm(s) of probabilities.}
@@ -28,14 +31,37 @@
where \eqn{r=(r_1,\dots,r_k)}{r=(r1,\dots,rk)}.
Draw \eqn{m} random samples from this distribution by doing
- \code{reulermultinom(n=m,size=N,rate=r,dt=dt)},
+ \code{dn <- reulermultinom(n=m,size=N,rate=r,dt=dt)},
where \code{r} is the vector of rates.
Evaluate the probability that \eqn{x=(x_1,\dots,x_k)}{x=(x1,\dots,xk)} are the numbers of individuals who have died in each of the \eqn{k} ways over the interval \eqn{\Delta t=}{}\code{dt}, by doing
\code{deulermultinom(x=x,size=N,rate=r,dt=dt)}.
- Direct access to the underlying C routines is available: see the header file \dQuote{pomp.h}, included with the package.
+ Breto & Ionides discuss how an infinitesimally overdispersed death process can be constructed by compounding a binomial process with a Gamma white noise process.
+ The Euler approximation of the resulting process can be obtained as follows.
+ Let the increments of the equidispersed process be given by
+
+ \code{reulermultinom(size=N,rate=r,dt=dt)}.
+
+ In this expression, replace the rate \eqn{r} with \eqn{r {\Delta W}/{\Delta t}},
+ where \eqn{\Delta W ~ Gamma(dt/\sigma^2,\sigma^2)} is the increment of an integrated Gamma white noise process with intensity \eqn{\sigma}.
+ That is, \eqn{\Delta W} has mean \eqn{\Delta t} and variance \eqn{\sigma^2 \Delta t}.
+ The resulting process is overdispersed and converges (as \eqn{\Delta t} goes to zero) to a well-defined process.
+ The following lines of \R code accomplish this:
+
+ \code{dW <- rgammawn(sigma=sigma,dt=dt)}
+
+ \code{dn <- reulermultinom(size=N,rate=r,dt=dW)}
+
+ or
+
+ \code{dn <- reulermultinom(size=N,rate=r*dW/dt,dt=dt)}.
+
+ He et al. use such overdispersed death processes in modeling measles.
+
+ For all of the functions described here, direct access to the underlying C routines is available:
+ see the header file \dQuote{pomp.h}, included with the package.
}
\value{
\item{reulermultinom}{
@@ -46,10 +72,26 @@
\item{deulermultinom}{
Returns a vector (of length equal to the number of columns of \code{x}) containing the probabilities of observing each column of \code{x} given the specified parameters (\code{size}, \code{rate}, \code{dt}).
}
+ \item{rgammawn}{
+ Returns a vector of length \code{n} containing random increments of the integrated Gamma white noise process with intensity \code{sigma}.
+ }
}
\author{Aaron A. King \email{kingaa at umich dot edu}}
\examples{
-print(x <- reulermultinom(5,size=100,rate=c(a=1,b=2,c=3),dt=0.1))
-deulermultinom(x,size=100,rate=c(1,2,3),dt=0.1)
+print(dn <- reulermultinom(5,size=100,rate=c(a=1,b=2,c=3),dt=0.1))
+deulermultinom(x=dn,size=100,rate=c(1,2,3),dt=0.1)
+## an Euler-multinomial with overdispersed transitions:
+dt <- 0.1
+dW <- rgammawn(sigma=0.1,dt=dt)
+print(dn <- reulermultinom(5,size=100,rate=c(a=1,b=2,c=3),dt=dW))
}
+\references{
+ C. Bret\\'o & E. L. Ionides,
+ Compound Markov counting processes and their applications to modeling infinitesimally over-dispersed systems.
+ Stoch. Proc. Appl., 121:2571--2591, 2011.
+
+ D. He, E. L. Ionides, & A. A. King,
+ Plug-and-play inference for disease dynamics: measles in large and small populations as a case study.
+ J. R. Soc. Interface, 7:271--283, 2010.
+ }
\keyword{distribution}
Modified: pkg/src/eulermultinom.c
===================================================================
--- pkg/src/eulermultinom.c 2012-01-19 11:44:49 UTC (rev 597)
+++ pkg/src/eulermultinom.c 2012-01-19 12:03:03 UTC (rev 598)
@@ -134,3 +134,37 @@
UNPROTECT(nprotect);
return f;
}
+
+// This function draws a random increment of a gamma whitenoise process.
+// This will have expectation=dt and variance=(sigma^2*dt)
+// If dW = rgammawn(sigma,dt), then
+// mu dW/dt is a candidate for a random rate process within an
+// Euler-multinomial context, i.e.,
+// E[mu*dW] = mu*dt and Var[mu*dW] = mu*sigma^2*dt
+double rgammawn (double sigma, double dt) {
+ double sigmasq;
+ sigmasq = sigma*sigma;
+ return (sigmasq > 0) ? rgamma(dt/sigmasq,sigmasq) : dt;
+}
+
+SEXP R_GammaWN (SEXP n, SEXP sigma, SEXP deltat) {
+ int nprotect = 0;
+ int k, nval, nsig, ndt;
+ double *x, *sig, *dt;
+ SEXP ans;
+ PROTECT(n = AS_INTEGER(n)); nprotect++;
+ nval = INTEGER(n)[0];
+ PROTECT(sigma = AS_NUMERIC(sigma)); nprotect++;
+ nsig = LENGTH(sigma);
+ sig = REAL(sigma);
+ PROTECT(deltat = AS_NUMERIC(deltat)); nprotect++;
+ ndt = LENGTH(deltat);
+ dt = REAL(deltat);
+ PROTECT(ans = NEW_NUMERIC(nval)); nprotect++;
+ x = REAL(ans);
+ for (k = 0; k < nval; k++)
+ x[k] = rgammawn(sig[k%nsig],dt[k%ndt]);
+ UNPROTECT(nprotect);
+ return ans;
+}
+
Modified: pkg/src/pomp.h
===================================================================
--- pkg/src/pomp.h 2012-01-19 11:44:49 UTC (rev 597)
+++ pkg/src/pomp.h 2012-01-19 12:03:03 UTC (rev 598)
@@ -28,11 +28,7 @@
// mu dW/dt is a candidate for a random rate process within an
// Euler-multinomial context, i.e.,
// E[mu*dW] = mu*dt and Var[mu*dW] = mu*sigma^2*dt
-inline double rgammawn (double sigma, double dt) {
- double sigmasq;
- sigmasq = sigma*sigma;
- return (sigmasq > 0) ? rgamma(dt/sigmasq,sigmasq) : dt;
-}
+double rgammawn (double sigma, double dt);
// facility for computing the inner produce of
// a vector of parameters ('coef') against a vector of basis-function values ('basis')
More information about the pomp-commits
mailing list