我检查了不同编程语言中的吉布斯采样;在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>