使用R的integrate()函数与来自R和Rcpp的函数一起使用

Question

使用R的integrate()函数与来自R和Rcpp的函数一起使用

4

在R中，函数dt有时会根据计算中使用的值发出警告信息。例如，dt(1.424781, 1486, -5)显示

Warning messages:
1: In dt(1.424781, 1486, -5) :
  full precision may not have been achieved in 'pnt{final}'
2: In dt(1.424781, 1486, -5) :
  full precision may not have been achieved in 'pnt{final}'

这个问题与这里所解释的完全相同。

那个帖子中的最后一篇建议使用rcpp来获得更好的精度。我尝试了一下，效果很好。
我将rcpp代码保存在文件boost_t.cpp中：

// [[Rcpp::depends(BH)]]

#include <Rcpp.h>
#include <boost/math/distributions/non_central_t.hpp> 

using namespace boost::math;

// [[Rcpp::export]]
double dnct(const double x, const double df, const double ncp) {
  non_central_t dist(df, ncp);
  return pdf(dist, x);
}

然后在R中只需执行以下操作

library(Rcpp)
sourceCpp("boost_t.cpp")
dnct(1.424781, 1486, -5)
[1] 4.393078e-10

我想做的是在R中在integrate()内部使用dnct()。当我尝试时，出现错误。例如：

integrate(function(delta) dnct(1.2, 20, delta), lower = 0, upper = 1)

出现了以下错误：错误：预期只有一个值：[范围=21]。

我原本希望能得到与dt相似的行为：

integrate(function(delta) dt(1.2, 20, delta), lower = 0, upper = 1)

0.2964214 with absolute error < 3.3e-15

我该如何解决这个问题？
提前感谢。

- Bruce James

哇，我现在看到你所提到的答案是我写的。我完全不记得了。 :) - Roland

哇，我现在明白了，你所提到的答案是我写的。我完全不记得了。 :) - Roland

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Roland · Accepted Answer

函数参数在进行积分时必须是矢量化的。

// [[Rcpp::depends(BH)]]

#include <Rcpp.h>
#include <boost/math/distributions/non_central_t.hpp> 

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector dnct(const double x, const double df, const NumericVector ncp) {

  R_xlen_t n = ncp.length();
  NumericVector y(n);

  for (R_xlen_t i = 0; i < n; ++i) {
    boost::math::non_central_t dist(df, ncp[i]);
    y[i] = boost::math::pdf(dist, x);
  }

  return y;
}

/*** R
integrate(function(delta) dnct(1.2, 20, delta), lower = 0, upper = 1)
*/

留给读者作为练习的是将 x 和 df 参数也进行向量化处理。