1

I have a Rcpp file where I am trying to integrate the function-

f<- function(x){dexp(x-a,beta)}

.

Where a=2 and beta=1.2. The file looks like-

    // [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#include <RcppNumerical.h>
#include <limits>
using namespace Numer;



class PDF: public Func
{
private:

  double beta;
  double M0;

public:
  PDF( double beta_, double M0_): beta(beta_), M0(M0_) {};

  double operator()(const double& x) const
  { 
    return   R::dexp(x-M0,beta,0);

  }
};

// [[Rcpp::export]]
double integrate_test2( double beta, double M0, double upper, double lower)
{

  PDF f( beta, M0);
  double err_est;
  int err_code;
  double res =  integrate ( f, lower,upper,err_est, err_code);
  return(res);

}

Now if I want to integrate the function from lower limit 2 to upper limit 5 using the code :

>  integrate_test2(  beta=1.2, M0=2, upper=5,  lower=2)
    [1] 0.917915

it gives me value "0.917915". But when I use R code to integrate using the code=

> integrate(function(x){dexp(x-2,1.2,0)},2,5)
  0.9726763 with absolute error < 1.1e-14

It gives me result " 0.9726763". Why this happens?

Geotas
  • 35
  • 6
  • Why don't you use the distribution function? `pexp(5-2,1.2) - pexp(2-2,1.2)` – Roland Dec 14 '18 at 12:40
  • @Roland I need to use RcppNumerical for my work – Geotas Dec 14 '18 at 12:44
  • 1
    This reminds me f https://stackoverflow.com/questions/53648355/double-integral-in-rcppnumerical. Anyway, @Roland's comment tells us that RcppNumerical is wrong here. It is probably best to [open an issue](https://github.com/yixuan/RcppNumerical/issues). – Ralf Stubner Dec 14 '18 at 13:15
  • 3
    Change to `return R::dexp(x-M0,1/beta,0);` – Roland Dec 14 '18 at 13:23
  • 2
    @RalfStubner Definitely not a bug, just the wrong parameterization. – Roland Dec 14 '18 at 13:36

0 Answers0