2

我正在评估具有不同支持的伽马分布的密度。这是我的 Rcpp 代码。

// [[Rcpp::export]] 
NumericVector fdensv(NumericVector w, NumericMatrix pard){
nj = w.size();
NumericVector out(nj);
for (int j=0;j<nj;j++){
    out[j] = R::dgamma(log(w[j]),pard(0,j),pard(1,j),0);
}
return out;
}

sourceCpp("test2.cpp")

现在测试代码

nj = 200
dr = exp(rgamma(nj,2,3))
pr = matrix(runif(400*2,2,4),2,200)
gg = fdensv(dr,pr)

gg2 = NULL
for (i in 1:nj)
 {
 gg2[i] = dgamma(log(dr[i]),pr[1,i],pr[2,i])
  }

cbind(gg,gg2)
all.equal(gg,gg2)

我得到“平均相对差异：32.77”......知道这种差异来自哪里吗？谢谢！

4

2 回答 2

2

R 对函数在 R 级别和 C 级别使用不同的参数化这一事实让您感到困扰gamma。如果您将dgamma调用更改为

R::dgamma(log(w[j]), pard(0,j), 1/pard(1,j), 0);

这在R-exts 6.7.1中有简要讨论——注意它需要scale而不是rate.

于 2014-05-03T09:25:56.537 回答
0

你把它弄得太复杂了。开始更简单：

// [[Rcpp::export]] 
double mydgamma(double a, double b, double c, int d) {
  return R::dgamma(a, b, c, d);
}

并查看有关形状和速率参数的标题。

然后

R> mydgamma(0.5, 1, 1, 0)
[1] 0.606531
R> dgamma(0.5,1)
[1] 0.606531
R>

或者

R> mydgamma(0.25, 0.5, 1, 0)
[1] 0.878783
R> dgamma(0.25,0.5)
[1] 0.878783
R>
于 2014-05-03T02:41:40.590 回答