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