R我正在尝试使用我在代码中编写的函数执行 4 维数值积分，C++然后在R使用Rcpp包时获取该函数。

下面是我的代码：

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#define EIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
#include <Eigen/Eigen>
#include <Rcpp.h>
#include <math.h>
#include <iostream>
#include <cstdlib>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/math/special_functions/gamma.hpp>

using namespace std;
using namespace Rcpp;
using namespace Numer;

// [[Rcpp::depends(BH)]]

// [[Rcpp::export]]
double gaussian_free_diffusion(double x,double x0, double sigma, double t) {
  double pi = 2 * acos(0.0);
  double a1 = (1/sqrt(2.0 * pi * sigma * t));
  double b1 = exp(-pow((x - x0), 2.0)/(2.0 * sigma * t));
  double res = a1 * b1;
  return res;
}

// [[Rcpp::export]]
double integral_CppFunctionCUBATURE(NumericVector Zt_pos, const double& xA0, const double& xB0, const double& yA0,
   const double& yB0, const double& t1, const double& sigma){

  double xAt_pos = Zt_pos[0];
  double xBt_pos = Zt_pos[1];
  double yAt_pos = Zt_pos[2];
  double yBt_pos = Zt_pos[3];

  double temp_pbxA = gaussian_free_diffusion(xAt_pos, xA0, sigma,t1);
  double temp_pbxB = gaussian_free_diffusion(xBt_pos, xB0, sigma, t1);
  double temp_pbyA = gaussian_free_diffusion(yAt_pos, yA0, sigma,t1);
  double temp_pbyB = gaussian_free_diffusion(yBt_pos, yB0, sigma, t1);

  return (temp_pbxB * temp_pbyB) * (temp_pbxA * temp_pbyA);

};

integral_CppFunctionCUBATURE是我使用的功能cubintegrate。在 R 中，我会这样做：

sourceCpp('./myCppcode.cpp')

xA0<-3
xB0<-2
yA0<-5
yB0<-10
t<-500
sigma<-1


cubature::cubintegrate(integral_CppFunctionCUBATURE,lower=rep(-1000,4), upper=rep(1000,4),
                       xA0=xA0, xB0=xB0, yA0=yA0,yB0=yB0, t1=t, sigma=sigma,method='cuhre')$integral

结果：0.9999978

多维积分大约需要 4 秒，但我想尽可能加快速度，因为当我对多个 xA0 值运行多维积分时，这可能需要相当长的时间。您对如何提高代码速度有什么建议吗？

C++另外，您会建议一种替代方法来在/中执行快速 4 维积分Rcpp吗？