CODEWITHSUNDEEP
rrcppnumerical-integration
Bruce James
Bruce James

Reputation: 131

Using R's integrate() with functions from R and Rcpp

Function dt in R can sometimes throw warning messages depending on the values used in the computation.
For example, dt(1.424781, 1486, -5) shows

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

This issue is exactly the same as explained here.

The last post in that thread suggests to use rcpp to get better precision. I did, and it works well.
I saved the rcpp code in file 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);
}

and then in R simply do

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

What I would like to do is to use dnct() inside integrate() in R.
When I try I get an error. For example:

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

gives this error: Error: Expecting a single value: [extent=21].

I was expecting a similar behavior that I have with dt:

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

0.2964214 with absolute error < 3.3e-15

How can I fix this?
Thanks in advance.

Upvotes: 4

Views: 105

Answers (1)

Roland
Roland

Reputation: 132969

The function argument over which you integrate must be vectorized:

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

It is left as an exercise for the reader to also vectorize the x and df arguments.

Upvotes: 5

Related Questions