将 R 的集成（）与 R 和 Rcpp 中的函数结合使用

Question

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}'

这个问题与

here所解释的完全一样。

该线程中的最后一篇文章建议使用 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 中的

dnct()

 中使用 

integrate()

。

当我尝试时出现错误。例如：

integrate(function(delta) dnct(1.2, 20, delta), lower = 0, upper = 1) 给出了这个错误：

Error: Expecting a single value: [extent=21]

。

我期待着与

dt

类似的行为：

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

0.2964214 with absolute error < 3.3e-15

我该如何解决这个问题？

预先感谢。

Answer 1

您积分的函数参数必须向量化：

// [[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

 参数进行向量化。