1 ## Simple Gibbs Sampler Example
2 ## Adapted from Darren Wilkinson's post at
3 ## http://darrenjw.wordpress.com/2010/04/28/mcmc-programming-in-r-python-java-and-c/
5 ## Sanjog Misra and Dirk Eddelbuettel, June-July 2011
13 ## Actual joint density -- the code which follow implements
14 ## a Gibbs sampler to draw from the following joint density f(x,y)
15 fun <-
function(
x,
y) {
16 x*x * exp(-x*
y*
y -
y*
y + 2*
y - 4*x)
19 ## Note that the full conditionals are propotional to
20 ## f(x|y) = (x^2)*exp(-x*(4+y*y)) : a Gamma density kernel
21 ## f(y|x) = exp(-0.5*2*(x+1)*(y^2 - 2*y/(x+1)) : Normal Kernel
23 ## There is a small typo in Darrens code.
24 ## The full conditional for the normal has the wrong variance
25 ## It should be 1/sqrt(2*(x+1)) not 1/sqrt(1+x)
26 ## This we can verify ...
27 ## The actual conditional (say for x=3) can be computed as follows
28 ## First - Construct the Unnormalized Conditional
29 fy.unnorm <-
function(
y) fun(3,
y)
31 ## Then - Find the appropriate Normalizing Constant
32 K <- integrate(fy.unnorm,-Inf,Inf)
34 ## Finally - Construct Actual Conditional
35 fy <-
function(
y) fy.unnorm(
y)/K$val
37 ## Now - The corresponding Normal should be
38 fy.dnorm <-
function(
y) {
40 dnorm(
y,1/(1+x),sqrt(1/(2*(1+x))))
44 fy.dnorm.wrong <-
function(
y) {
46 dnorm(
y,1/(1+x),sqrt(1/((1+x))))
51 ## Actual (gray thick line)
52 curve(fy,-2,2,
col=
'grey',lwd=5)
54 ## Correct Normal conditional (blue dotted line)
55 curve(fy.dnorm,-2,2,
col=
'blue',add=T,lty=3)
57 ## Wrong Normal (Red line)
58 curve(fy.dnorm.wrong,-2,2,
col=
'red',add=T)
61 ## Here is the actual Gibbs Sampler
62 ## This is Darren Wilkinsons R code (with the corrected variance)
63 ## But we are returning only his columns 2 and 3 as the 1:N sequence
64 ## is never used below
65 Rgibbs <-
function(
N,thin) {
66 mat <- matrix(0,ncol=2,nrow=N)
72 y <-
rnorm(1,1/(x+1),1/sqrt(2*(x+1)))
79 ## We can also try the R compiler on this R function
80 RCgibbs <- cmpfun(Rgibbs)
83 ## mat <- Rgibbs(10000,10); dim(mat)
84 ## would give: [1] 10000 2
87 ## Now for the Rcpp version -- Notice how easy it is to code up!
90 using namespace Rcpp; // inline does that for us already
92 // n and thin are SEXPs which the Rcpp::as function maps to C++ vars
94 int thn = as<int>(thin);
97 NumericMatrix mat(N, 2);
99 RNGScope scope; // Initialize Random number generator
101 // The rest of the code follows the R version
104 for (i=0; i<N; i++) {
105 for (j=0; j<thn; j++) {
106 x = ::Rf_rgamma(3.0,1.0/(y*y+4));
107 y = ::Rf_rnorm(1.0/(x+1),1.0/sqrt(2*x+2));
113 return mat; // Return to R
117 RcppGibbs <- cxxfunction(signature(
n=
"int", thin =
"int"),
118 gibbscode, plugin=
"Rcpp")
122 #include <gsl/gsl_rng.h>
123 #include <gsl/gsl_randist.h>
125 using namespace Rcpp;
130 int thin = as<int>(thns);
132 gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
135 for (i=0; i<
N; i++) {
136 for (
j=0;
j<thin;
j++) {
137 x = gsl_ran_gamma(r,3.0,1.0/(
y*
y+4));
138 y = 1.0/(x+1)+gsl_ran_gaussian(r,1.0/sqrt(2*x+2));
149 GSLGibbs <- cxxfunction(signature(ns="int", thns = "int"),
150 body=gslgibbscode, includes=gslgibbsincl,
153 ## without RcppGSL, using cfunction()
154 #GSLGibbs <- cfunction(signature(ns="int", thns = "int"),
155 # body=gslgibbscode, includes=gslgibbsincl,
157 # cppargs="-I/usr/include",
158 # libargs="-lgsl -lgslcblas")
161 ## Now for some tests
162 ## You can try other values if you like
163 ## Note that the total number of interations are N*thin!
164 Ns <- c(1000,5000,10000,20000)
165 thins <- c(10,50,100,200)
171 for (i in seq_along(Ns)) {
172 tim_R[i] <- system.time(mat <- Rgibbs(Ns[i],thins[i]))[3]
173 tim_RC[i] <- system.time(cmat <- RCgibbs(Ns[i],thins[i]))[3]
174 tim_Rgsl[i] <- system.time(gslmat <- GSLGibbs(Ns[i],thins[i]))[3]
175 tim_Rcpp[i] <- system.time(rcppmat <- RcppGibbs(Ns[i],thins[i]))[3]
176 cat("Replication #", i, "complete \n")
180 speedup <- round(tim_R/tim_Rcpp,2);
181 speedup2 <- round(tim_R/tim_Rgsl,2);
182 speedup3 <- round(tim_R/tim_RC,2);
183 summtab <- round(rbind(tim_R,tim_RC, tim_Rcpp,tim_Rgsl,speedup3,speedup,speedup2),3)
184 colnames(summtab) <- c("N=1000","N=5000","N=10000","N=20000")
185 rownames(summtab) <- c("Elasped Time (R)","Elasped Time (RC)","Elapsed Time (Rcpp)", "Elapsed Time (Rgsl)",
186 "SpeedUp Rcomp.","SpeedUp Rcpp", "SpeedUp GSL")
190 ## Contour Plots -- based on Darren's example
191 if (interactive() &&
require(KernSmooth)) {
192 op <- par(mfrow=c(4,1),mar=c(3,3,3,1))
196 contour(x,y,z,
main="Contours of actual distribution",xlim=c(0,2), ylim=c(-2,4))
197 fit <- bkde2D(
as.matrix(
mat),c(0.1,0.1))
198 contour(drawlabels=T, fit$x1, fit$x2, fit$fhat, xlim=c(0,2), ylim=c(-2,4),
199 main=paste("Contours of empirical distribution:",round(tim_R[4],2)," seconds"))
200 fitc <- bkde2D(
as.matrix(rcppmat),c(0.1,0.1))
201 contour(fitc$x1,fitc$x2,fitc$fhat,xlim=c(0,2), ylim=c(-2,4),
202 main=paste("Contours of Rcpp based empirical distribution:",round(tim_Rcpp[4],2)," seconds"))
203 fitg <- bkde2D(
as.matrix(gslmat),c(0.1,0.1))
204 contour(fitg$x1,fitg$x2,fitg$fhat,xlim=c(0,2), ylim=c(-2,4),
205 main=paste("Contours of GSL based empirical distribution:",round(tim_Rgsl[4],2)," seconds"))
210 ## also use rbenchmark package
213 res <- benchmark(Rgibbs(N,
thn),
217 columns=c(
"test",
"replications",
"elapsed",
218 "relative",
"user.self",
"sys.self"),