a technical article on how to integrate R and C


http://quantitative-ecology.blogspot.com/2008/11/call-c-from-r-and-r-from-c.html
Several years ago, while a research associate at the University of Chicago, I had the privilege of sitting in on a course taught by Peter Rossi: Bayesian Applications in Marketing and MicroEconometrics. This course — one I recommend to anyone at U Chicago who is interested in statistics — was an incredibly clear treatment of Bayesian statistics, but the aspect I appreciated most was Peter’s careful demonstration of Bayesian theory and methods using R.

One feature of R that I had not made use of up until that point was the ability to call compiled C and Fortran functions from within R (this makes loop-heavy Metropolis-Hastings samplers much, much faster). It turns out that you can also include the R libraries in C source code so that R functions (e.g., random number generators) can be easily accessed. The R-Cran website has an excellent tutorial on how to develop R extensions (here), but I wanted to share an example Peter used in class because it is extremely brief, and for 95% of what I do, this is all I need.

As Peter writes, this is an incredibly inefficient way of simulating from the chisquare distribution, but it demonstrates the point. His more extensive writeup is located here.

Save the following as testfun.c:

#include <R.h>
#include <Rmath.h>
#include <math.h>

/* Function written by Peter Rossi from the University of Chicago GSB */
/*http://faculty.chicagogsb.edu/peter.rossi/teaching/37904/Adding%20Functions%20Written%20in%20C%20to%20R.pdf */

/* include standard C math library and R internal function declarations */
void mychisq(double *vec, double *chisq, int *nu)
/* void means return nothing */
{
int i,iter; /* declare local vars */
/* all statements end in; */
GetRNGstate(); /* set random number seed */
for (i=0 ; i < *nu; ++i)
/* loop over elements of vec */
/*nu “dereferences" the pointer */
{ /* vectors start at 0 location!*/
vec[i] = rnorm(0.0,1.0); /*use R function to draw normals */
Rprintf(“%ith normal draw= %lf \n",(i+1),vec[i]);
/* print out results for “debugging" */
}
*chisq=0.0;
iter=0;
while(iter < *nu) /* “while" version of a loop */
{
if( iter == iter)
{*chisq=*chisq + vec[iter]*vec[iter];}
/* redundant if stmnt */
iter=iter+1; /* note: can’t use ** */
/* if you want to be “cool" use iter += 1 */
}
PutRNGstate(); /* write back ran number seed */
}

To call this function in R, you first need to compile it. To do this you need all the standard compilers and libraries for your operating system. For Debian or Ubuntu, this should do it (if I missed a package, let me know in the comments):

$ sudo aptitude update
$ sudo aptitude install build-essential r-base-dev

Now, you should be able to compile the function:

$ R CMD SHLIB testfun.c

If all goes well, you should see the files testfun.o and testfun.so in the directory. To test the function we will source the following R script into R:

call_mychisq<-function(nu){
##This function is just a wrapper for .C
vector=double(nu); chisq=1
.C(“mychisq",as.double(vector),res=as.double(chisq),as.integer(nu))$res
}

##Load the compiled code (you may need to include
## the explicit file path if it is not local
## NOTE: for Windows machines, you will want to load testfun.dll"
dyn.load(“testfun.so")
result dyn.load(“testfun.so")
> result result
[1] 8.268028

This entry was posted in Stat/Biostat. Bookmark the permalink.

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / 變更 )

Twitter picture

You are commenting using your Twitter account. Log Out / 變更 )

Facebook照片

You are commenting using your Facebook account. Log Out / 變更 )

Google+ photo

You are commenting using your Google+ account. Log Out / 變更 )

連結到 %s