Below is an example of writing the loop for your Metropolis Hastings code
in c and accessing it from R/S.  Note that all random numbers are
generated in R/S and passed to C so that we can take advantage of the fact
that R/S does a much better/easier/more valid job of generating random
numbers than C does. Feel free to just put a copy of this e-mail on the
web page and direct specific questions to me.



First we need to creat the c code to do the looping, put the following
lines in a file called mh2.c (not including the lines of -------):

---------------mh2.c---------------------

#include <math.h>

double g( double x ){
  return( exp(-0.5*pow(x,2)) + 0.5 * exp(-0.5*pow((x-3),2) ) );
}

void mh2(int *m, double *norms, double *unifs, double *out, double *canda,
double *accepta)
{
  int i;
  double cand,accept;

  for (i=1 ; i < *m ; i++){
    out[i] = out[i-1];
    cand = out[i-1] + norms[i];
    canda[i] = cand;
    accept = g(cand)/g(out[i-1]);
    accepta[i] = accept;
    if (unifs[i] < accept){
      out[i] = cand;
    }
  }

}

-----------end mh2.c--------------------------

This looks just like Dr. Reese's code except that the normal and uniform
random numbers are passed in as vectors rather than generated on the fly
(because R/S does a much better job of generating random numbers).  canda
and accepta are empty vectors provided by R/S and filled with values by
the C code.  Note also that the vectors go from 0 to m-1 in C but from 1
to m in R/S.

Now we need to compile the C code into a .dll file using gcc (gcc is a
free C compiler that should be installed on all lab computers).

from a dos window in the directory where mh2.c is saved, run the following
2 commands:

gcc -c mh2.c
gcc -shared -o mh2.dll mh2.o

you should now have a file in that directory named mh2.dll

Now start up R or S-PLUS and define the following function to call the
compiled C code:

mh2 <- function(m,start){
out <- numeric(m)
out[1] <- start
.C("mh2",
    m=as.integer(m), norms=as.double(rnorm(m)),
    unifs=as.double(runif(m)), out=as.double(out),
    cand=as.double(numeric(m)), accept=as.double(numeric(m)))
}


The first argument is M, the number of steps to take and the second
argument is the starting value. We first create the out vector and put our
starting value in as the first element, then the .C function passes all of
our information (including creating vectors of random numbers and passing
them) to the compiled "mh2" function.  The return value is a list with a
component for each thing we passed in (with the vectors updated).


Next we need to load the compiled function into R/S by using one of the
following 2 commands (change the path to reflect where you stored your
files):

dyn.load("c:/classes/s595/mh2.dll")  # R
dll.load("c:/classes/s595/mh2.dll")  # S-PLUS


Now (if everything worked) you can run the function like:

temp <- mh2(3000, 3)
plot(density(temp$out))


I also created (copied) the following function for time comparisons:

mh1 <- function(m,start){
out <- numeric(m)
out[1] <- start
for(i in 2:m){
   out[i]<-out[i-1]
   cand <- rnorm(1,out[i-1],1)
   accept <- g(cand)/g(out[i-1])
   u <- runif(1,0,1)
   if (u < accept) {out[i] <- cand}
}
out
}

here are some timed comparisons:

> system.time(temp <- mh1(3000,3))
[1] 0.97 0.00 1.23   NA   NA
> system.time(temp2 <- mh2(3000,3))
[1] 0.02 0.00 0.15   NA   NA
> system.time(temp <- mh1(30000,3))
[1] 11.14  0.20 27.86    NA    NA
> system.time(temp2 <- mh2(30000,3))
[1] 0.29 0.01 0.92   NA   NA
> system.time(temp <- mh1(300000,3))
[1] 107.96   1.02 337.60     NA     NA
> system.time(temp2 <- mh2(300000,3))
[1]  4.88  0.18 28.29    NA    NA


mh2 (the C compiled version) is a lot faster than mh1 (R only).  Times for
S-PLUS are similar for the compiled code, I did not try the S only version
in S-PLUS.


Hope someone finds this useful.