Adaptive tau-leaping in SSA
Just a brief update on my activities 11 days before the R package GillespieSSA is to be released. I have been hacking away on the code, debugging, and running test simulations like there is no tomorrow. The last SSA implementation I added is the optimized tau-leaping algorithm by Cao et al. (2006) (see the reference list on the GillespieSSA page). Initially it was a bit daunting trying to wrap ones head the underlying, especially after implementing all the “straightforward” methods. For a while I was asking myself “Can I build this?” and then my two year old replied in god old Bob style “Yes we can!”. Actually, once one figures it out it is a rather ingenious implementation of the SSA.
The optimized tau-leap method (or OTL for short) choses the time increment 
Below is the trajectory of Lotka’s predator-prey prey model illustrating this point. The red point is the neutrally stable-point for the deterministic ODE and also the starting point for the simulation. The blue points is the trajectory of the populations plotted every 5 time steps. There are several interesting patterns one can see in this plot. First, while the finite population starts on the neutrally-stable equilibrium it rapidly starts to oscillate with ever increasing amplitude, inevitably leading to extinction. This nicely illustrates the important concept of deterministic models not always being able to accurately predict the dynamics of finite (and hence stochastic) populations . The second interesting feature in the plot is that can see the adaptive switches between the Poisson tau-leap method and the direct method that the optimized tau-leap method relies on, these are the regions with the dense point trajectories (lower-left corner).
By the way, Gillespie uses this model in his 1976 paper to illustrate the direct method and this is his view of the trajectory,
Thus, one way of viewing the behavior of the Lotka model is to say that the microscopic random fluctuations inherent in the reactions cause the system to execute a “drunkard’s walk” over the continuum of concentric, neutrally stable solution orbits, staggering sometimes outward and sometimes inward.
Variance-Covariance Matrix in glm
vcov.glm<-function(obj){
#return the variance-covariance matrix of a glm object
#from p. 188 in Venables and Ripley. 2002.
#Modern Applied Statistics With S. Springer. New York.
so <- summary(obj,corr=F)
so$dispersion * so$cov.unscaled
}
Visualizing Data
files within R libraries. Of course you can paste them directly
into R, but this web page allows you to view the images within the
context of the help files:
R Graphical Manuals
Also, check out this link for the R Graph Gallery.
Testing for missing values of numeric (double) vectors in R
I just came across this bit of R trickiness when trying to improve the error messages for one of our packages. I tried adding a check at the C level for NA_REAL, but to no avail. The test failed even when I was certain that I was passing down a vector containing missing values. Consulting the Writing R Extensions Manual (WREM) led me quickly to ISNA. The WREM says that ISNA only applies to numeric values of type double and that other types can be checked by equality comparison to NA_INTEGER, NA_LOGICAL, etc. What could perhaps be better emphasized is that ISNA is the only appropriate way to test for missingness of REALSXP elements and that equality testing with NA_REAL does not work.
It is an easy mistake to make since one might be lulled into complacency by the repeating patterns of a switch statement. Here’s how NOT to do it:
SEXP hello_hasNA(SEXP v)
{
int i, found_na = 0;
for (i = 0; i < length(v); i++) {
switch (TYPEOF(v)) {
case INTSXP:
if (INTEGER(v)[i] == NA_INTEGER)
found_na = 1;
break;
case LGLSXP:
if (LOGICAL(v)[i] == NA_LOGICAL)
found_na = 1;
break;
case REALSXP:
if (REAL(v)[i] == NA_REAL) /* WRONG, must use ISNA() */
found_na = 1;
break;
case STRSXP:
if (STRING_ELT(v, i) == NA_STRING)
found_na = 1;
break;
default:
error("no support for type");
}
if (found_na)
break;
}
return ScalarLogical(found_na);
}
To fix things, replace the REALSXP case like this:
case REALSXP:
if (ISNA(REAL(v)[i]))
found_na = 1;
break;
Technorati Tags: programming, R
Testing for missing values of numeric (double) vectors in R
I just came across this bit of R trickiness when trying to improve the error messages for one of our packages. I tried adding a check at the C level for NA_REAL, but to no avail. The test failed even when I was certain that I was passing down a vector containing missing values. Consulting the Writing R Extensions Manual (WREM) led me quickly to ISNA. The WREM says that ISNA only applies to numeric values of type double and that other types can be checked by equality comparison to NA_INTEGER, NA_LOGICAL, etc. What could perhaps be better emphasized is that ISNA is the only appropriate way to test for missingness of REALSXP elements and that equality testing with NA_REAL does not work.
It is an easy mistake to make since one might be lulled into complacency by the repeating patterns of a switch statement. Here’s how NOT to do it:
SEXP hello_hasNA(SEXP v)
{
int i, found_na = 0;
for (i = 0; i < length(v); i++) {
switch (TYPEOF(v)) {
case INTSXP:
if (INTEGER(v)[i] == NA_INTEGER)
found_na = 1;
break;
case LGLSXP:
if (LOGICAL(v)[i] == NA_LOGICAL)
found_na = 1;
break;
case REALSXP:
if (REAL(v)[i] == NA_REAL) /* WRONG, must use ISNA() */
found_na = 1;
break;
case STRSXP:
if (STRING_ELT(v, i) == NA_STRING)
found_na = 1;
break;
default:
error("no support for type");
}
if (found_na)
break;
}
return ScalarLogical(found_na);
}
To fix things, replace the REALSXP case like this:
case REALSXP:
if (ISNA(REAL(v)[i]))
found_na = 1;
break;
Technorati Tags: programming, R
E-mail updates
