在Rcpp中调用cubature库hcubature函数实现四维积分的问题

• 第一步：对齐integrand类型的函数签名要求
  cubature要求传入的被积函数必须符合以下固定签名：

int your_integrand(unsigned dim, const double *x, void *fdata, unsigned fdim, double *fval);

你当前编写的integralFunction参数格式完全不符合要求，这是最核心的问题。

• 第二步：封装固定参数为结构体
  你需要把积分过程中不会变化的固定参数（xA0、xB0、yA0、yB0、t1、sigma）封装为一个自定义结构体，通过hcubature的fdata参数传递给被积函数。
• 第三步：调整数据类型统一为double
  cubature库所有入出参都用double类型，建议把你代码里的float全部替换为double，避免精度损失和类型转换报错。同时不要给临时计算变量加const &引用声明，避免悬空引用风险。

#include <Rcpp.h>
#include <cubature.h>
#include <cmath>

using namespace Rcpp;

// 替换原2dbrownian函数，统一用double，去掉不必要的const引用
double dbrownian2d(double x, double x0, double sigma, double t) {
  double a1 = 1/sqrt(2.0 * M_PI * sigma * t);
  double b1 = exp(-pow(x - x0, 2)/(2.0 * sigma * t));
  return a1 * b1;
}

// 封装固定参数的结构体
struct IntegralParams {
  double xA0;
  double xB0;
  double yA0;
  double yB0;
  double t1;
  double sigma;
};

// 适配hcubature要求的被积函数
int integrand_adapter(unsigned dim, const double *x, void *fdata, unsigned fdim, double *fval) {
  // 转换fdata为我们定义的参数结构体指针
  IntegralParams *params = static_cast<IntegralParams*>(fdata);
  // 拆分x数组为4个积分变量：x[0]=xAt_pos, x[1]=xBt_pos, x[2]=yAt_pos, x[3]=yBt_pos
  double xAt_pos = x[0];
  double xBt_pos = x[1];
  double yAt_pos = x[2];
  double yBt_pos = x[3];
  // 调用原积分计算逻辑
  double temp_pbxA = dbrownian2d(xAt_pos, params->xA0, params->sigma, params->t1);
  double temp_pbxB = dbrownian2d(xBt_pos, params->xB0, params->sigma, params->t1);
  double temp_pbyA = dbrownian2d(yAt_pos, params->yA0, params->sigma, params->t1);
  double temp_pbyB = dbrownian2d(yBt_pos, params->yB0, params->sigma, params->t1);
  // 写入计算结果
  fval[0] = temp_pbxB * temp_pbyB * temp_pbxA * temp_pbyA;
  return 0;
}

// 导出给R调用的积分接口
// [[Rcpp::export]]
NumericVector compute_4d_integral(double xA0, double xB0, double yA0, double yB0, double t1, double sigma,
                                  NumericVector xmin, NumericVector xmax,
                                  double reqAbsError = 1e-8, double reqRelError = 1e-6, size_t maxEval = 0) {
  // 初始化参数结构体
  IntegralParams params = {xA0, xB0, yA0, yB0, t1, sigma};
  // 存储积分结果和误差
  double val, err;
  // 调用hcubature
  int ret = hcubature(1, integrand_adapter, &params,
                      4, xmin.begin(), xmax.begin(),
                      maxEval, reqAbsError, reqRelError,
                      ERROR_INDIVIDUAL, &val, &err);
  // 返回结果给R
  NumericVector res = NumericVector::create(_["value"] = val, _["error"] = err, _["return_code"] = ret);
  return res;
}

# 示例：四个积分变量下限都为-10，上限都为10
xmin <- rep(-10, 4)
xmax <- rep(10, 4)
# 调用积分函数
res <- compute_4d_integral(xA0 = 0, xB0 = 1, yA0 = 0, yB0 = 1, t1 = 0.1, sigma = 1, xmin = xmin, xmax = xmax)
print(res)

hcubature返回码为0表示积分成功收敛，非0表示出现异常。

内容的提问来源于stack exchange，提问作者CafféSospeso