Gibbs sampler bivariate normal in Rcpp

I want to generate samples of the bivariate normal distribution (Gibbs Sampler) with fixed parameters in Rcpp. The code in R, is quite simple, but since I am new with Rcpp, it has taken me a lot of work to adapt the code.

My R code

``````gibbsR <- function(n,mu1,mu2,s1,s2,rho){
X <- numeric() ; Y <- numeric()
X[1] <- rnorm(1,mu1,s1) #init value for x_0
for(i in 1:n){
Y[i]   <- rnorm(1,mu2+(s2/s1)*rho*(X[i]-mu1),sqrt((1-rho^2)*s2^2)) # Y|X
X[i+1] <- rnorm(1,mu1+(s1/s2)*rho*(Y[i]-mu2),sqrt((1-rho^2)*s1^2)) # X|Y
}
cbind(x=X[-1],y=Y)}

system.time(resR <- gibbsR(n=1000000,mu1=170,mu2=70,s1=10,s2=5,rho=0.8))
#colMeans(resR);apply(resR, 2, sd);cor(resR)
head(resR)
x        y
[1,] 185.2425 79.27488
[2,] 178.0975 75.53521
[3,] 178.4902 74.29250
[4,] 173.7096 73.37504
[5,] 180.2141 72.89918
[6,] 171.0300 72.66280
``````

My Rcpp code

``````library(Rcpp)
cppFunction("
NumericMatrix gibbsC(int n, double mu1, double mu2, double s1, double s2, double rho) {
NumericVector x; NumericVector y;
x[0] = Rf_rnorm(mu1,s1);
for(int i=0; i<n; i++){
y[i] = Rf_rnorm(mu2+(s2/s1)*rho*(x[i]-mu1),sqrt(1-pow(rho,2))*pow(s2,2)); // Y|X
x[i+1] = Rf_rnorm(mu1+(s1/s2)*rho*(y[i]-mu2),sqrt(1-pow(rho,2))*pow(s1,2)); // X|Y
}
return(cbind(x[-0],y));
}")

system.time(resC <- gibbsC(n=100000,mu1=170,mu2=70,s1=10,s2=5,rho=0.8))
#colMeans(resR);apply(resR, 2, sd);cor(resR)
head(resC)
``````

Any idea how to get a result similar to my gibbsR function

2 answers

• answered 2018-05-16 09:19

The problem here is a problem of parentheses for the `sqrt` in second argument of the `rnorm`.

This code works:

``````#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix gibbsC(int n, double mu1, double mu2, double s1, double s2, double rho) {

NumericMatrix res(n, 2);
double x, y, rho2;
x = R::rnorm(mu1,s1);
for(int i=0; i<n; i++){
rho2 = 1 - rho * rho;
res(i, 1) = y = R::rnorm(mu2+(s2/s1)*rho*(x-mu1),sqrt(rho2*s2*s2)); // Y|X
res(i, 0) = x = R::rnorm(mu1+(s1/s2)*rho*(y-mu2),sqrt(rho2*s1*s1)); // X|Y
}

return res;
}

/*** R
gibbsR <- function(n,mu1,mu2,s1,s2,rho){
X <- numeric() ; Y <- numeric()
X[1] <- rnorm(1,mu1,s1) #init value for x_0
for(i in 1:n){
Y[i]   <- rnorm(1,mu2+(s2/s1)*rho*(X[i]-mu1),sqrt((1-rho^2)*s2^2)) # Y|X
X[i+1] <- rnorm(1,mu1+(s1/s2)*rho*(Y[i]-mu2),sqrt((1-rho^2)*s1^2)) # X|Y
}
cbind(x=X[-1],y=Y)}

N <- 1e5

set.seed(1)
system.time(resR <- gibbsR(n=N,mu1=170,mu2=70,s1=10,s2=5,rho=0.8))
#colMeans(resR);apply(resR, 2, sd);cor(resR)
head(resR)

set.seed(1)
system.time(resC <- gibbsC(n=N,mu1=170,mu2=70,s1=10,s2=5,rho=0.8))
#colMeans(resR);apply(resR, 2, sd);cor(resR)
head(resC)
*/
``````

• answered 2018-05-16 09:22

Here a merge of you code with the example in the Rcpp gallery:

``````#include <Rcpp.h>

using namespace Rcpp;       // shorthand

// [[Rcpp::export]]
NumericMatrix gibbsC(int n, double mu1, double mu2, double s1, double s2, double rho) {

NumericMatrix mat(n, 2);
double x=R::rnorm(mu1, s1); //init value for x_0
double y=0.0;

for (int i=0; i < n; ++i) {
y = R::rnorm(mu2 + (s2 / s1) * rho * (x - mu1),
sqrt(1.0 - rho * rho) * s2); // Y|X
x = R::rnorm(mu1 + (s1 / s2) * rho * (y - mu2),
sqrt(1.0 - rho * rho) * s1); // X|Y
mat(i,0) = x;
mat(i,1) = y;
}
return mat;             // Return to R
}
/*** R
gibbsR <- function(n,mu1,mu2,s1,s2,rho){
X <- numeric() ; Y <- numeric()
X[1] <- rnorm(1,mu1,s1) #init value for x_0
for(i in 1:n){
Y[i]   <- rnorm(1,mu2+(s2/s1)*rho*(X[i]-mu1),sqrt((1-rho^2)*s2^2)) # Y|X
X[i+1] <- rnorm(1,mu1+(s1/s2)*rho*(Y[i]-mu2),sqrt((1-rho^2)*s1^2)) # X|Y
}
cbind(x=X[-1],y=Y)}
set.seed(42)
system.time(resR <- gibbsR(n=1000000,mu1=170,mu2=70,s1=10,s2=5,rho=0.8))
head(resR)

set.seed(42)
system.time(resC <- gibbsC(n=1000000,mu1=170,mu2=70,s1=10,s2=5,rho=0.8))
head(resC)
*/
``````

Note that many computations done at every step here are actually static and could be done before the loop. When I source this file I get as output:

``````> Rcpp::sourceCpp('gibbsC.cpp')

> gibbsR <- function(n,mu1,mu2,s1,s2,rho){
+   X <- numeric() ; Y <- numeric()
+   X[1] <- rnorm(1,mu1,s1) #init value for x_0
+   for(i in 1:n){
+    .... [TRUNCATED]

> set.seed(42)

> system.time(resR <- gibbsR(n=1000000,mu1=170,mu2=70,s1=10,s2=5,rho=0.8))
user  system elapsed
6.221   0.015   6.237

> head(resR)
x        y
[1,] 178.2424 73.78974
[2,] 180.7385 75.19553
[3,] 185.4323 73.97701
[4,] 191.5329 75.88896
[5,] 191.3092 78.42501
[6,] 186.2806 85.38363

> set.seed(42)

> system.time(resC <- gibbsC(n=1000000,mu1=170,mu2=70,s1=10,s2=5,rho=0.8))
user  system elapsed
0.162   0.004   0.166

> head(resC)
[,1]     [,2]
[1,] 178.2424 73.78974
[2,] 180.7385 75.19553
[3,] 185.4323 73.97701
[4,] 191.5329 75.88896
[5,] 191.3092 78.42501
[6,] 186.2806 85.38363
``````

Note that there was a subtle bug in you C++ code concerning parenthesis:

``````sqrt(1-pow(rho,2))*pow(s2,2)
^            ^
``````