RcppNumerical 中的二重积分 [英] Double integral in RcppNumerical

问题描述

我必须计算函数的二重积分:

I have to calculate the double integral of the function:

 > DIntegral <- function(x,y){res <- pnorm(x,1,0.1) * dexp(y-2,1.2)
                              return(res)
                              }

x 和 y 的上限分别为:10 和 Infinity.x 和 y 的下限分别为:1 和 2.

The upper limit for x and y are: 10 and Infinity respectively. The lower limit for x and y are: 1 and 2 respectively.

如何在 RcppNumerical 中进行这种双重积分?

How can I do this double integration in RcppNumerical?

对于一维集成，我的 C++ 文件如下所示:

For one dimensional integration my C++ file looks like:

    // [[Rcpp::depends(RcppEigen)]]
   // [[Rcpp::depends(RcppNumerical)]]
   #include <RcppNumerical.h>
   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);

}

有限限制的二维积分代码

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



class PDF: public MFunc
{
private:
  double mu;
  double sigma;
  double beta;
  double M0;

 public:
  PDF( double mu_, double sigma_, double beta_, double M0_): mu(mu_),     sigma(sigma_), beta(beta_), M0(M0_) {};

  double operator()(Constvec& x)
  { 
    return  R::pnorm(x[0],mu,sigma,1,0) * R::dexp(x[1]-M0,beta,0);

  }
};

// [[Rcpp::export]]
double integrate_test2( double mu, double sigma, double beta, double M0)
{
  Eigen::VectorXd lower(2);
  lower << 1, 2;
  Eigen::VectorXd upper(2);
  upper << 10, 50;

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

}

推荐答案

我已更新您的代码，以将积分限制作为参数并返回错误代码和估计值:

I have updated your code to take the integration limits as arguments and to return the error code and estimate:

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

class PDF: public MFunc
{
private:
  double mu;
  double sigma;
  double beta;
  double M0;

public:
  PDF( double mu_, double sigma_, double beta_, double M0_): mu(mu_),     sigma(sigma_), beta(beta_), M0(M0_) {};

  double operator()(Constvec& x)
  { 
    return  R::pnorm(x[0],mu,sigma,1,0) * R::dexp(x[1]-M0,beta,0);

  }
};

// [[Rcpp::export]]
Rcpp::List integrate_test2( double mu, double sigma, double beta, double M0, Eigen::VectorXd lower, Eigen::VectorXd upper)
{
  PDF f( mu, sigma, beta, M0);
  double err_est;
  int err_code;
  double res = integrate ( f, lower, upper,err_est,err_code);
  return Rcpp::List::create(
    Rcpp::Named("result") = res,
    Rcpp::Named("error_estimate") = err_est,
    Rcpp::Named("error_code") = err_code
  );
}
/*** R
integrate_test2(1, 0.1, 1.2, 2, c(1, 2), c(10, 50))
integrate_test2(1, 0.1, 1.2, 2, c(1, 2), c(10, 1e4))
integrate_test2(1, 0.1, 1.2, 2, c(1, 2), c(10, 1e6))
integrate_test2(1, 0.1, 1.2, 2, c(1, 2), c(10, 1e8))
*/

结果:

> Rcpp::sourceCpp('2d_int.cpp')

> integrate_test2(1, 0.1, 1.2, 2, c(1, 2), c(10, 50))
$result
[1] 8.950068

$error_estimate
[1] 0.3570577

$error_code
[1] 1


> integrate_test2(1, 0.1, 1.2, 2, c(1, 2), c(10, 1e4))
$result
[1] 4.787999

$error_estimate
[1] 16.12484

$error_code
[1] 1


> integrate_test2(1, 0.1, 1.2, 2, c(1, 2), c(10, 1e6))
$result
[1] 1.605216e-314

$error_estimate
[1] 4.320299e-313

$error_code
[1] 0


> integrate_test2(1, 0.1, 1.2, 2, c(1, 2), c(10, 1e8))
$result
[1] 0

$error_estimate
[1] 0

$error_code
[1] 0

因此对于小值 y 的上限，积分不会收敛.当它收敛时，结果(几乎)为零.当将上限增加到大约 1e307 时，这不会改变，即几乎 .Machine$double.xmax.之后我得到 NaN.

So for small values the upper limit of y, the integration does not converge. And when it does converge, the result is (nearly) zero. This does not change when one increases the upper limit up to about 1e307, i.e. almost .Machine$double.xmax. After that I get NaN.

然而，如果我使用 cubature 包，结果就大不相同了:

However, if I use the cubature package the result is quite different:

library(cubature)
DIntegral <- function(x){
  res <- pnorm(x[1],1,0.1) * dexp(x[2]-2,1.2)
  return(res)
}
cubintegrate(f = DIntegral, lower = c(1, 2), upper = c(10, 50), method = "hcubature")
#> $integral
#> [1] 8.961023
#> 
#> $error
#> [1] 4.888071e-05
#> 
#> $neval
#> [1] 983
#> 
#> $returnCode
#> [1] 0
cubintegrate(f = DIntegral, lower = c(1, 2), upper = c(10, 1000), method = "hcubature")
#> $integral
#> [1] 8.960415
#> 
#> $error
#> [1] 7.145898e-05
#> 
#> $neval
#> [1] 1701611595
#> 
#> $returnCode
#> [1] 0
cubintegrate(f = DIntegral, lower = c(1, 2), upper = c(10, Inf), method = "hcubature")
#> $integral
#> [1] 8.960105
#> 
#> $error
#> [1] 8.515124e-05
#> 
#> $neval
#> [1] 1706522167
#> 
#> $returnCode
#> [1] 0

我不确定这里发生了什么.

I am not sure what is going on here.

这篇关于RcppNumerical 中的二重积分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！