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

Question

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.

Wow, I see now that the answer you refer to was written by me. I don't remember that at all. :) — Roland, Jun 14 '23 at 06:09

Roland · Accepted Answer · 2023-06-16T05:15:36.987

5

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.

edited Jun 16 '23 at 05:15

answered Jun 14 '23 at 05:43

Roland

127,288
10
191
288

It works, thanks! I would not be able to code this as I do not know C nor C++. – Bruce James Jun 14 '23 at 06:13
It's really not that difficult. I'm not good in C++ myself but you can learn the basics in a few hours. – Roland Jun 14 '23 at 06:17
Nice Boost integration example. Only nit I'd have is the dual `using namespace` which I find dicey. – Dirk Eddelbuettel Jun 14 '23 at 12:52
@DirkEddelbuettel Thanks for the feedback. You mean the risk of name collisions, right? – Roland Jun 15 '23 at 05:00
1

Right. Even now in your edited example: is `pdf()` really from `Rcpp` or does it need the Boost namespave? – Dirk Eddelbuettel Jun 15 '23 at 12:10

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

1 Answers1