两种数值积分方法结果差异：Gauss-Hermite是否更可靠？

阿华AIGC实验室

2026-5-12

关于高斯密度与函数乘积积分的数值方法准确性疑问解答

咱们来聊聊你遇到的这个积分计算问题——用rmutil::int2()和Gauss-Hermite求积法计算同一个高斯密度与函数乘积的积分，结果存在差异，想知道是不是GH方法更准确，而int2的数值积分只是近似？

首先得明确：两种方法本质上都是数值近似，但Gauss-Hermite求积法在处理这类与高斯密度相关的积分时，理论上有天生的精度优势，下面具体拆解：

两种方法的核心差异

rmutil::int2()：这是一个通用的自适应数值积分工具，针对无穷区间的积分，它会先做变量变换把无穷区间转成有限区间，然后通过逐步细化区间、调整采样点来逼近积分值。它的精度严重依赖于你设置的参数（比如eps、max、d），如果参数不够严格，或者被积函数在某些区域（比如高斯分布的尾部）变化特殊，就容易产生较大的近似误差。
Gauss-Hermite求积法：这是专门为计算「函数与高斯密度乘积的积分」设计的正交多项式求积方法。对于单变量的$\int_{-\infty}^{\infty} f(x) e^{-x2} dx$，它能通过精心选择的节点和权重，对多项式函数的积分做到精确计算；扩展到多元情况时，通过张量积节点和变量变换（就像你代码里做的那样，把多元正态转换成标准正态形式），只要你的被积函数能被多项式很好地近似，用足够多的节点就能得到极高的精度。

针对你的示例分析

你的被积函数正好是gfun(B0,B1)乘以多元正态密度，完全匹配GH方法的适用场景——GH的权重天生就对应了高斯密度的结构，所以计算效率和精度都比通用的int2高很多。而int2作为通用工具，面对这种带高斯尾的积分，需要更严格的参数设置（比如把eps调到1e-8，增大max和d）才能逼近GH的结果，否则很容易因为区间采样不够细致，导致结果偏差。

验证建议

如果想进一步确认，可以试试这两个操作：

调优int2的参数：把eps设得更小（比如1e-8），同时增大max和d的值，看看结果是否更接近GH的输出。如果是的话，说明之前的int2参数设置不够严格，导致了较大的近似误差。
增加GH的节点数：比如把mgauss.hermite()里的n改成15或20，观察result_GH是否稳定。如果结果几乎不变，说明GH已经收敛到了非常准确的积分值。

你的示例代码

1. rmutil::int2() 实现

library(rmutil)
Sig <- matrix (c(0.2^2, 0, 0, 0.8^2), ncol=2)
Mu<- c(2, 0)
to.integrate <- function(B0, B1) {
  first.int= 1/0.8 * (1.2 * exp(B0 + B1 * 0.5))^(-1/0.8) * gamma(1/0.8)
  B=matrix(c(B0, B1), ncol=1)
  multi.norm=1 / (2 * pi * det(Sig)^(1/2)) * exp (- 0.5 * t( B - Mu ) %*% solve(Sig) %*%( B - Mu ) )
  return (first.int %*% multi.norm)
}
result_int2 <- int2(to.integrate, a=c(-Inf, -Inf), b=c(Inf, Inf), eps=1.0e-6, max=16, d=5)

2. 多元高斯求积点实现

library(statmod)
mgauss.hermite <- function(n, mu, sigma) {
  dm <- length(mu)
  gh <- gauss.quad(n, 'hermite')
  gh <- cbind(gh$nodes, gh$weights)
  idx <- as.matrix(expand.grid(rep(list(1:n), dm)))
  pts <- matrix(gh[idx, 1], nrow(idx), dm)
  wts <- apply(matrix(gh[idx, 2], nrow(idx), dm), 1, prod)
  eig <- eigen(sigma)
  rot <- eig$vectors %*% diag(sqrt(eig$values))
  pts <- t(rot %*% t(pts) + mu)
  return(list(points=pts, weights=wts))
}
nod_wei <- mgauss.hermite(10, mu=Mu, sigma=Sig)
gfun <- function(B0, B1) {
  first.int <- 1/0.8 *(1.2 * exp(B0 + B1 * 0.5))^(-1/0.8)* gamma(1/0.8)
  return(first.int)
}
result_GH <- sum(gfun(nod_wei$points[, 1], nod_wei$points[, 2]) * nod_wei$weights)/pi