R 和 C++ 中的吉布斯采样

Question

我检查了不同编程语言的 gibbs 采样；在

  x <- rgamma(1,3,y*y+4)
  y <- rnorm(1,1/(x+1),1/sqrt(2*(x+1)))

在 C++ 中

  x = R::rgamma(3.0,1.0/(y*y+4));
  y = R::rnorm(1.0/(x+1),1.0/sqrt(2*x+2));

如果它使用 R 函数，为什么它在 C++ 中有所不同，因为 rgamma 不采用 n= 观察次数，它采用比例而不是速率作为默认输入，并且 rnorm 也没有 n= 观察次数。

对于 Rcpp，它完全不同，例如；

 y = ::Rf_rnorm(1.0/(x+1),1.0/sqrt(2*x+2));

Answer 1

你有什么问题？

R::rgamma() 也来自 Rcpp。它方便地包装了 C 级、非命名空间的 ::Rf_rnorm().

请注意，您还有 vectorised Rcpp::rnorm()——并且在 Darren 最初的 post 之后还有大量 Gibbs 采样器示例威尔金森。我们最好的例子可能是 this page at the Rcpp Gallery.

编辑： 由于您显然对 shape = 1/rate 参数化感到困惑，这里有一个完整且有效的示例：

我们先通过Rcpp编译一个调用C++的方便的R函数：

R> cppFunction("NumericVector callrgamma(int n, double shape, double scale) { 
+                  return(rgamma(n, shape, scale)); }")
R>

然后我们调用 R，确保我们修复了种子：

R> set.seed(42); rgamma(3, 2.0, 2.0)   # calling R
[1] 1.824478 0.444055 0.779610
R>

现在，使用相同的种子，我们调用 C++ 函数并确保遵守“1/over”重新参数化还有：

R> set.seed(42); callrgamma(3, 2.0, 1/2.0)   # calling Rcpp
[1] 1.824478 0.444055 0.779610
R>

Gibbs Sampling in R and C++