Python计算混合高斯分布微分熵时遇数值错误求助

你的问题出在被积函数的数值稳定性上，而非积分上下限本身。当混合高斯的PDF值极小时，计算1/pdf_gaus_mixture(input)会触发浮点数溢出，再结合对数运算后，会出现0*inf这类无效数值，最终导致积分结果为nan。

具体改进方案 #

1. 重构被积函数，避免溢出

原函数pdf * log(1/pdf)等价于-pdf * log(pdf)，直接用这个形式可以彻底避免计算1/pdf的操作，从根源上消除溢出风险。

2. 处理PDF趋近于0的边界情况

当PDF值小于浮点数能稳定表示的极小值（比如1e-300）时，log(pdf)会变成负无穷，此时-pdf*log(pdf)的极限是0（因为PDF趋近于0的速度远快于log(pdf)趋向负无穷的速度），可以直接返回0，避免无效值。

3. 可选：缩小积分区间（不影响精度）

虽然integrate.quad支持无穷区间，但你可以根据混合高斯的分布范围设置更紧凑的积分区间——比如两个高斯均值各自加减5倍标准差（因为高斯分布99.999%的概率集中在这个范围内），既提升计算速度，又不会损失精度。

修改后的代码 #

import numpy as np
from scipy.stats import norm
from scipy import integrate

def diff_entropy(nu, constant):
    def pdf_gaus_mixture(x):
        return (1 - nu) * norm.pdf(x, loc=0, scale=1) + nu * norm.pdf(x, loc=constant, scale=1)
    
    def func(x):
        pdf = pdf_gaus_mixture(x)
        # 处理PDF极小的情况，直接返回0
        if pdf < 1e-300:
            return 0.0
        # 等价于pdf * log(1/pdf)，但数值更稳定
        return -pdf * np.log(pdf)
    
    # 计算合适的积分区间：覆盖两个高斯分布的99.999%概率范围
    lower_bound = min(0, constant) - 5 * 1
    upper_bound = max(0, constant) + 5 * 1
    return integrate.quad(func, lower_bound, upper_bound)[0]

# 测试计算
nu=0.1
beta=0.01
delta=0.1
sigma=0.01
constant = np.sqrt(1/((beta/delta)+(sigma**2)))
print(diff_entropy(nu, constant))

为什么这样有效？ #

• 重构后的被积函数避免了1/pdf的溢出，直接使用对数的性质转换表达式，数值稳定性大幅提升。
• 边界条件处理消除了log(0)和0*inf这类无效数值的产生。
• 缩小后的积分区间完全覆盖了混合高斯的概率集中区域，积分结果的精度和直接用无穷区间几乎一致，但计算速度更快。

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