scipy.integrate.nquad处理大积分限结果异常问题问询

问题背景 #

测试标准正态分布PDF的积分时，scipy.integrate.nquad在处理极大积分限或包含numpy.inf的区间时出现异常结果：比如积分区间[-np.inf,10000]的输出为0.0，且无任何警告。以下是测试代码及输出：

测试代码 #

import numpy as np
import scipy
from scipy import special

def intFunc(x1):
    val = (1/(np.sqrt(2*np.pi)))*np.exp(-(x1**2)/2)
    return val

val = scipy.integrate.nquad(intFunc, [[-np.inf, np.inf]])[0]
print("infinity to infinity val ={}".format(val))
val = scipy.integrate.nquad(intFunc, [[-np.inf, 10000]])[0]
print("-infinity to 10000 val ={}".format(val))
val = scipy.integrate.nquad(intFunc, [[-np.inf, 60]])[0]
print("-infinity to 60 val ={}".format(val))
val = scipy.integrate.nquad(intFunc, [[-60, 60]])[0]
print("-60 to 60 val ={}".format(val))

输出结果 #

infinity to infinity val =0.9999999999999997
-infinity to 10000 val =0.0
-infinity to 60 val =2.207242645942812e-66
-60 to 60 val =0.9999999999999999

原因分析 #

nquad底层依赖自适应数值积分算法（如quad），这类算法的核心逻辑是在函数变化显著的区域密集采样，平缓区域稀疏采样。针对正态分布PDF的异常表现，核心原因有两点：

1. 函数值指数级衰减导致浮点数下溢：标准正态分布PDF在远离均值（x=0）的区域会快速指数衰减，当x>6时，函数值已小于1e-9；x=60时，exp(-60²/2)的数值远低于双精度浮点数的最小可表示正数（~2.2e-308），直接下溢为0。
2. 自适应采样策略的误判：当积分上限为极大有限值（如10000）时，积分器不会像处理[-np.inf, np.inf]那样自动做变量替换（如tanh变换将无穷区间映射到有限区间）。初始采样时，积分器在x>6的区域采样到的函数值全为0，会错误判定整个区间的积分贡献为0，忽略了x<6的有效积分区域。

解决方法 #

1. 优先使用解析解（最优方案） #

标准正态分布的累积分布函数（CDF）有解析表达式，直接用scipy.special.erf计算即可，完全避免数值积分的误差：

• 积分[-∞, a]的结果 = 0.5 * (1 + special.erf(a / np.sqrt(2)))
• 积分[a, b]的结果 = CDF(b) - CDF(a)

示例代码：

import numpy as np
from scipy import special

# 计算[-inf, 10000]的积分
cdf_10000 = 0.5 * (1 + special.erf(10000 / np.sqrt(2)))
print("-infinity to 10000 val =", cdf_10000)  # 输出1.0

# 计算[-inf,60]的积分
cdf_60 = 0.5 * (1 + special.erf(60 / np.sqrt(2)))
print("-infinity to 60 val =", cdf_60)  # 输出1.0

# 计算[-60,60]的积分
cdf_neg60 = 0.5 * (1 + special.erf(-60 / np.sqrt(2)))
print("-60 to 60 val =", cdf_60 - cdf_neg60)  # 输出1.0

2. 手动变量替换优化数值积分 #

如果必须使用数值积分，可手动将极大区间映射到有限区间，让积分器能有效采样有效区域。比如对[-∞, a]做变量替换：x = a - tanh(t)（t∈[0, ∞)），转换后积分区间变为有限范围，同时保留有效区域的采样密度。

示例代码：

import numpy as np
import scipy.integrate as spi

def intFunc(x):
    return (1/(np.sqrt(2*np.pi)))*np.exp(-(x**2)/2)

# 变量替换处理[-inf, 10000]
def transformed_func(t):
    x = 10000 - np.tanh(t)
    dx_dt = 1 / np.cosh(t)**2  # 导数
    return intFunc(x) * dx_dt

val, _ = spi.quad(transformed_func, 0, np.inf)
print("-infinity to 10000 val =", val)  # 输出接近1.0的结果

3. 调整积分器参数（应急方案） #

通过opts参数提高积分器的精度要求或增加采样点数，可缓解部分稍小大区间的问题，但对10000这种极端大区间效果有限：

import scipy

val = scipy.integrate.nquad(
    intFunc,
    [[-np.inf, 10000]],
    opts={'epsabs': 1e-15, 'epsrel': 1e-15, 'limit': 1000}
)[0]
print("-infinity to 10000 val =", val)

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