Gibbs sampler bivariate normal in Rcpp

Question

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

Answer 1

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)
*/

Answer 2

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)
    ^            ^