Rcpp 中的 Gibbs 采样器双变量正态分布

Question

我想在 Rcpp 中生成具有固定参数的二元正态分布样本 (Gibbs Sampler)。 R 中的代码非常简单，但由于我是 Rcpp 新手，因此我花了很多工作来调整代码。

我的 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)}

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

我的 Rcpp 代码

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)

知道如何获得类似于我的 gibbsR 函数的结果

Answer 1

这里的问题是sqrt 在rnorm 的第二个参数中的括号问题。

此代码有效：

#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

这里是你的代码与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)
*/

请注意，这里每一步完成的许多计算实际上都是静态的，可以在循环之前完成。当我获取这个文件时，我得到了输出：

> 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

请注意，您的 C++ 代码中有一个关于括号的细微错误：

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