[Pomp-commits] r594 - pkg/src

[noreply at r-forge.r-project.org]

Author: kingaa
Date: 2012-01-15 20:30:18 +0100 (Sun, 15 Jan 2012)
New Revision: 594

Modified:
   pkg/src/pomp.h
Log:
- add two inline helper functions: one for doing the exponential to geometric rate correction (see He et al. 2010),
the other for computing the increment of an integrated gamma white-noise process (see Breto et al.)


Modified: pkg/src/pomp.h
===================================================================
--- pkg/src/pomp.h	2012-01-12 14:56:13 UTC (rev 593)
+++ pkg/src/pomp.h	2012-01-15 19:30:18 UTC (rev 594)
@@ -8,11 +8,34 @@
 #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
+inline double rgammawn (double sigma, double dt) {
+  double sigmasq;
+  sigmasq = sigma*sigma;
+  return (sigmasq > 0) ? rgamma(dt/sigmasq,sigmasq) : 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