rrcppnumerical-integration

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.


Solution

  • 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.