基于cubature的多元复积分特征函数密度反转问题求助

我在研究中尝试实现多维特征函数反转公式来计算概率密度，该公式出自Sheppard (1991)的论文。调研后发现现有R包存在局限性：

• pracma和elliptic包支持复积分，但最多仅处理2个变量，无法满足多维需求；
• CharFunToolR包实现了反转公式，但仅支持单变量；
• cubature包尝试后存在运行时间长、精度不可靠的问题，且无法妥善处理复数。

我用多元正态特征函数做测试验证（实际研究中特征函数更复杂），测试代码如下：

p <- 6

mu <- c(-0.448, -0.48, -0.661, -0.325, -0.539, -0.297)

Sigma <- matrix(nrow = p, ncol = p, c(
   0.136, -0.013,  0.001, -0.009, -0.005, 0.004,
  -0.013,  0.123,  0.050,  0.002,  0.037, 0.040,
   0.001,  0.050,  0.095, -0.011,  0.056, 0.014,
  -0.009,  0.002, -0.011,  0.099, -0.005, 0.008, 
  -0.005,  0.037,  0.056, -0.005,  0.075, 0.003, 
   0.004,  0.040,  0.014,  0.008,  0.003, 0.102
))

x <- c(0.053, 0.523, -0.041, -0.589, -0.044, -0.066)

使用mvtnorm包的dmvnorm()得到真实结果：

mvtnorm::dmvnorm(x, mu, as.matrix(Sigma))
[1] 0.01410531

但用反转公式实现时，代码如下：

cf <- function(tee, x, mu, Sigma) {
  exp(-1i * crossprod(tee, x)) * exp(1i * crossprod(tee, mu) - 0.5 * t(tee) %*% Sigma %*% tee)
}

cub <- cubature::hcubature(
  f          = cf,
  lowerLimit = rep(-Inf, p),
  upperLimit = rep(Inf, p),
  tol        = 1e-4,
  maxEval    = 1000,
  x          = x,
  mu         = mu,
  Sigma      = Sigma
)

(2 * pi)^(-p) * cub$integral
[1] 0.07572407

得到的结果与真实值偏差极大，增大maxEval后结果完全不稳定；不指定maxEval则运行时间过长，同时出现**"虚部在强制转换中被丢弃"**的警告。即使按包作者指导实现了特征函数的向量化版本，结果仍无改善。

需要解决：

1. 修正代码正确性，使结果接近真实值；
2. 控制运行时间在合理范围；
3. 消除复数相关的警告；
4. 了解其他可用于多维特征函数反转的R包。

解决方案 #

1. 代码正确性修复：分离复积分的实部与虚部 #

cubature包默认仅处理实数输出，直接传入复值函数会导致虚部被丢弃，这是警告和结果错误的核心原因。多维特征函数反转公式的积分是复积分，需分别计算实部和虚部的积分，最终取实部（概率密度为实数）。

修改后的代码：

# 拆分实部和虚部分别定义函数
cf_real <- function(tee, x, mu, Sigma) {
  Re(exp(-1i * crossprod(tee, x)) * exp(1i * crossprod(tee, mu) - 0.5 * t(tee) %*% Sigma %*% tee))
}

cf_imag <- function(tee, x, mu, Sigma) {
  Im(exp(-1i * crossprod(tee, x)) * exp(1i * crossprod(tee, mu) - 0.5 * t(tee) %*% Sigma %*% tee))
}

# 用有限区间近似无穷积分：正态特征函数在|t|>5后衰减极快
cub_real <- cubature::hcubature(
  f          = cf_real,
  lowerLimit = rep(-5, p),
  upperLimit = rep(5, p),
  tol        = 1e-6,
  maxEval    = 5000,
  x          = x,
  mu         = mu,
  Sigma      = Sigma
)

cub_imag <- cubature::hcubature(
  f          = cf_imag,
  lowerLimit = rep(-5, p),
  upperLimit = rep(5, p),
  tol        = 1e-6,
  maxEval    = 5000,
  x          = x,
  mu         = mu,
  Sigma      = Sigma
)

# 计算最终密度，取实部（虚部积分理论上应为0）
final_result <- Re((2 * pi)^(-p) * (cub_real$integral + 1i * cub_imag$integral))
final_result

关键优化点： #

• 拆分实部虚部，避免强制转换丢弃虚部；
• 用[-5,5]^p近似无穷区间：多元正态特征函数的指数项为负二次型，t分量绝对值超过5后函数值可忽略，足够满足精度；
• 调整tol和maxEval平衡精度与速度。

2. 运行时间优化 #

• 使用cubature包的pcubature（并行积分）替代hcubature，利用多核CPU加速：

cub_real <- cubature::pcubature(
  f          = cf_real,
  lowerLimit = rep(-5, p),
  upperLimit = rep(5, p),
  tol        = 1e-6,
  maxEval    = 5000,
  x          = x,
  mu         = mu,
  Sigma      = Sigma
)

• 高维场景下，可使用稀疏网格积分（如sparseGrid包），比普通自适应积分效率更高。

3. 消除警告 #

通过拆分实部和虚部的积分，避免了复值到实值的强制转换，自然消除"虚部被丢弃"的警告。

4. 其他可用R包 #

• sparseGrid：支持高维数值积分，适配特征函数反转的高维积分场景；
• mvQuad：提供多种高维积分方法，包括自适应和稀疏网格，支持自定义函数；
• RcppNumerical：基于C的数值积分库，可实现高效复积分，适合极致性能需求，需编写少量C代码。

内容的提问来源于stack exchange，提问作者little_sky