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.

enter image description here

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 F. Privé

    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 Ralf Stubner

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