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

Gibbs Sampling in R and C++

本文关键字：布斯 采样 C++      更新时间：2023-10-16

我检查了不同编程语言中的吉布斯采样；在R 中

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

您的问题是什么？

R::rgamma()也来自Rcpp。它方便地封装了C级、无名称空间的::Rf_rnorm()。

请注意，您还有矢量化Rcpp::rnorm()，在Darren Wilkinson的最初帖子之后，有很多Gibbs Sampler的例子。我们最好的例子可能是Rcpp画廊的这个页面。

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

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

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>
相关文章：