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((1rho^2)*s2^2)) # YX
X[i+1] < rnorm(1,mu1+(s1/s2)*rho*(Y[i]mu2),sqrt((1rho^2)*s1^2)) # XY
}
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(1pow(rho,2))*pow(s2,2)); // YX
x[i+1] = Rf_rnorm(mu1+(s1/s2)*rho*(y[i]mu2),sqrt(1pow(rho,2))*pow(s1,2)); // XY
}
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

The problem here is a problem of parentheses for the
sqrt
in second argument of thernorm
.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*(xmu1),sqrt(rho2*s2*s2)); // YX res(i, 0) = x = R::rnorm(mu1+(s1/s2)*rho*(ymu2),sqrt(rho2*s1*s1)); // XY } 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((1rho^2)*s2^2)) # YX X[i+1] < rnorm(1,mu1+(s1/s2)*rho*(Y[i]mu2),sqrt((1rho^2)*s1^2)) # XY } 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) */

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); // YX x = R::rnorm(mu1 + (s1 / s2) * rho * (y  mu2), sqrt(1.0  rho * rho) * s1); // XY 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((1rho^2)*s2^2)) # YX X[i+1] < rnorm(1,mu1+(s1/s2)*rho*(Y[i]mu2),sqrt((1rho^2)*s1^2)) # XY } 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(1pow(rho,2))*pow(s2,2) ^ ^