1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
|
## Simple Gibbs Sampler Example
## Adapted from Darren Wilkinson's post at
## http://darrenjw.wordpress.com/2010/04/28/mcmc-programming-in-r-python-java-and-c/
##
## Sanjog Misra and Dirk Eddelbuettel, June-July 2011
## Updated by Dirk Eddelbuettel, Mar 2020
suppressMessages({
library(Rcpp)
library(rbenchmark)
})
## Actual joint density -- the code which follow implements
## a Gibbs sampler to draw from the following joint density f(x,y)
fun <- function(x,y) {
x*x * exp(-x*y*y - y*y + 2*y - 4*x)
}
## Note that the full conditionals are propotional to
## f(x|y) = (x^2)*exp(-x*(4+y*y)) : a Gamma density kernel
## f(y|x) = exp(-0.5*2*(x+1)*(y^2 - 2*y/(x+1)) : Normal Kernel
## There is a small typo in Darrens code.
## The full conditional for the normal has the wrong variance
## It should be 1/sqrt(2*(x+1)) not 1/sqrt(1+x)
## This we can verify ...
## The actual conditional (say for x=3) can be computed as follows
## First - Construct the Unnormalized Conditional
fy.unnorm <- function(y) fun(3,y)
## Then - Find the appropriate Normalizing Constant
K <- integrate(fy.unnorm,-Inf,Inf)
## Finally - Construct Actual Conditional
fy <- function(y) fy.unnorm(y)/K$val
## Now - The corresponding Normal should be
fy.dnorm <- function(y) {
x <- 3
dnorm(y,1/(1+x),sqrt(1/(2*(1+x))))
}
## and not ...
fy.dnorm.wrong <- function(y) {
x <- 3
dnorm(y,1/(1+x),sqrt(1/((1+x))))
}
if (interactive()) {
## Graphical check
## Actual (gray thick line)
curve(fy,-2,2,col='grey',lwd=5)
## Correct Normal conditional (blue dotted line)
curve(fy.dnorm,-2,2,col='blue',add=T,lty=3)
## Wrong Normal (Red line)
curve(fy.dnorm.wrong,-2,2,col='red',add=T)
}
## Here is the actual Gibbs Sampler
## This is Darren Wilkinsons R code (with the corrected variance)
## But we are returning only his columns 2 and 3 as the 1:N sequence
## is never used below
Rgibbs <- function(N,thin) {
mat <- matrix(0,ncol=2,nrow=N)
x <- 0
y <- 0
for (i in 1:N) {
for (j in 1:thin) {
x <- rgamma(1,3,y*y+4)
y <- rnorm(1,1/(x+1),1/sqrt(2*(x+1)))
}
mat[i,] <- c(x,y)
}
mat
}
## Now for the Rcpp version -- Notice how easy it is to code up!
cppFunction("NumericMatrix RcppGibbs(int N, int thn){
NumericMatrix mat(N, 2); // Setup storage
double x = 0, y = 0; // The rest follows the R version
for (int i = 0; i < N; i++) {
for (int j = 0; j < thn; j++) {
x = R::rgamma(3.0,1.0/(y*y+4));
y = R::rnorm(1.0/(x+1),1.0/sqrt(2*x+2));
}
mat(i,0) = x;
mat(i,1) = y;
}
return mat; // Return to R
}")
## Use of the sourceCpp() is preferred for users who wish to source external
## files or specify their headers and Rcpp attributes within their code.
## Code here is able to easily be extracted and placed into its own C++ file.
## Compile and Load
sourceCpp(code="
#include <RcppGSL.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
using namespace Rcpp; // just to be explicit
// [[Rcpp::depends(RcppGSL)]]
// [[Rcpp::export]]
NumericMatrix GSLGibbs(int N, int thin){
gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
double x = 0, y = 0;
NumericMatrix mat(N, 2);
for (int i = 0; i < N; i++) {
for (int j = 0; j < thin; j++) {
x = gsl_ran_gamma(r,3.0,1.0/(y*y+4));
y = 1.0/(x+1)+gsl_ran_gaussian(r,1.0/sqrt(2*x+2));
}
mat(i,0) = x;
mat(i,1) = y;
}
gsl_rng_free(r);
return mat; // Return to R
}")
## Now for some tests
## You can try other values if you like
## Note that the total number of interations are N*thin!
Ns <- c(1000,5000,10000,20000)
thins <- c(10,50,100,200)
tim_R <- rep(0,4)
tim_Rgsl <- rep(0,4)
tim_Rcpp <- rep(0,4)
for (i in seq_along(Ns)) {
tim_R[i] <- system.time(mat <- Rgibbs(Ns[i],thins[i]))[3]
tim_Rgsl[i] <- system.time(gslmat <- GSLGibbs(Ns[i],thins[i]))[3]
tim_Rcpp[i] <- system.time(rcppmat <- RcppGibbs(Ns[i],thins[i]))[3]
cat("Replication #", i, "complete \n")
}
## Comparison
speedup <- round(tim_R/tim_Rcpp,2);
speedup2 <- round(tim_R/tim_Rgsl,2);
summtab <- round(rbind(tim_R, tim_Rcpp,tim_Rgsl,speedup,speedup2),3)
colnames(summtab) <- c("N=1000","N=5000","N=10000","N=20000")
rownames(summtab) <- c("Elasped Time (R)","Elapsed Time (Rcpp)", "Elapsed Time (Rgsl)",
"SpeedUp Rcpp", "SpeedUp GSL")
print(summtab)
## Contour Plots -- based on Darren's example
if (interactive() && require(KernSmooth)) {
op <- par(mfrow=c(4,1),mar=c(3,3,3,1))
x <- seq(0,4,0.01)
y <- seq(-2,4,0.01)
z <- outer(x,y,fun)
contour(x,y,z,main="Contours of actual distribution",xlim=c(0,2), ylim=c(-2,4))
fit <- bkde2D(as.matrix(mat),c(0.1,0.1))
contour(drawlabels=T, fit$x1, fit$x2, fit$fhat, xlim=c(0,2), ylim=c(-2,4),
main=paste("Contours of empirical distribution:",round(tim_R[4],2)," seconds"))
fitc <- bkde2D(as.matrix(rcppmat),c(0.1,0.1))
contour(fitc$x1,fitc$x2,fitc$fhat,xlim=c(0,2), ylim=c(-2,4),
main=paste("Contours of Rcpp based empirical distribution:",round(tim_Rcpp[4],2)," seconds"))
fitg <- bkde2D(as.matrix(gslmat),c(0.1,0.1))
contour(fitg$x1,fitg$x2,fitg$fhat,xlim=c(0,2), ylim=c(-2,4),
main=paste("Contours of GSL based empirical distribution:",round(tim_Rgsl[4],2)," seconds"))
par(op)
}
## also use rbenchmark package
N <- 20000
thn <- 200
res <- benchmark(Rgibbs(N, thn),
RcppGibbs(N, thn),
GSLGibbs(N, thn),
columns=c("test", "replications", "elapsed",
"relative", "user.self", "sys.self"),
order="relative",
replications=10)
print(res)
## And we are done
|
