__Last post I set up a for loop for storage for doing MCMC that also prints the iterations every 10,000 iterations. This post will do a simple MCMC with a conjugate prior and a MCMC via a metropolis hastings sampler for the same conjugate prior. This will require us to break into random number generation in c++ . Before jumping into it let’s set up the example, lets use the beta-binomial prior and likelihood. Let__

Then we know that the posterior distribution is

Where s and f are the number success and failures. Don’t believe me?

https://en.wikipedia.org/wiki/Conjugate_prior

Ok so now let’s talk about random number generation because we will need this to sample p from the posterior distribution. C++ has it’s own pseudorandom generator called rand(), but using this to generate random numbers from a distribution has very undesirable properties, even to the point where a uniform distribution will become skewed! Instead we will use the <random> package in c++11 which is described in great detail in this video:

https://channel9.msdn.com/Events/GoingNative/2013/rand-Considered-Harmful

Particularly, they mention that this package has capabilities of randomizing data from the following distributions.

This uses the Mersenne_Twister pseudorandom sampler, which astoundingly has a period of 2^19937 – 1. So there is essentially no chance of repeating the same seed for a random number generator. More info is available on wiki.

https://en.wikipedia.org/wiki/Mersenne_Twister

This is actually the same pseudorandom number generator used in R to generate random numbers.

https://en.wikibooks.org/wiki/R_Programming/Random_Number_Generation

However, it seems to have some trouble working with Rcpp (actually, it doesn’t work at all) and many statisticians recommend using R to source random number generators out to it, we will do that here. One day I’ll look more into random number generation because it seems that things aren’t truly “random”.

For the gibbs sampler, we simply need to use the code we used in the past post and just add

*R::rbeta(a+y,b+n-y); *and store the results in the storage matrix (after burn-in).

For the metropolis hastings, let’s use a uniform proposal (or a beta(1,1) proposal). Then we can generate new observations using R::rbeta() or R::runif(). But first, we need to allocate some memory to variables that we will need for computing the acceptance ratio, add these prior to the If(…) print statements in the for loop.

*double rand;*

* double p;*

* double dnew;*

* double dold;*

* double lknew;*

* double lkold;*

* p=pstart;*

* double cut;*

* double U;*

These are going to be needed in the Metropolis-Hastings sampler. Now for the changes to the for loop.

*rand= R::rbeta(1,1);* This randomly draws a uniform(0,1) a.k.a. a beta(1,1)

*dnew=R::dbeta(rand, a, b, true);* Finds the prior density of the new draw

*dold=R::dbeta(p, a, b, true);* Finds the prior density of the old draw

*lknew=R::dbinom(y,n,rand,true);* Likelihood of new Draw

*lkold=R::dbinom(y,n,p,true);* Likelihood of old Draw

*cut = lknew+dnew-dold-lkold;* Calculates the acceptance ratio

*U=R::runif(0,1);* Draws a random variable for accept/reject

*U=log(U);* takes the log of the random probability

*if(U<cut){p=rand;}* If statement that returns the new value if it’s more likely than the old.

**NOTE: You cant do log(R::unif(0,1)) or even R::log(unif(0,1)). IT DOESN’T WORK! You have to do it separately!!**

Here’s all the code together for those of you following at home:

And here are the results!!

If we just do a for loop in R for this MCMC, we see that the sampling time is

$sampling.time

[1] 9.08

But if we use our c++ version we coded (a.k.a MCMC++ copryright), the sampling time is:

$sampling.time

[1] 1.235

Which is 7.352227 times faster than in R. This is just a simple MCMC++, can you imagine the speed up for a larger Gibbs sampler??

Happy coding!! Check out some of my other posts too if you’re interested in less statistical material!

-Andrew G. Chapple

## 2 thoughts on “Banging My Head Against the Keyboard: MCMC++”