如何计算非常非常小的 y 值的积分（SciPy 四边形）

Question

这是一个对数正态分布的概率密度函数：

from scipy.stats import lognorm
def f(x): return lognorm.pdf(x, s=0.2, loc=0, scale=np.exp(10))

此函数的 y 值非常小（最大值 ~ 1E-5），并且分布在 x 值 ~1E5 上。我们知道PDF的积分应该是1，但是使用下面的代码直接计算积分时，由于计算精度不够，答案是1E-66轮。

from scipy.integrate import quad
import pandas as pd
ans, err = quad(f, -np.inf, np.inf)

您能帮我正确计算这样的积分吗？谢谢。

Answer 1

您使用的值对应于具有平均值mu = 10 和标准差sigma = 0.2 的基础正态分布。使用这些值，分布模式（即 PDF 最大值的位置）位于exp(mu - sigma**2) = 21162.795717500194。函数quad 工作得很好，但它可以被愚弄。在这种情况下，显然quad 只对值非常小的函数进行采样——它永远不会“看到”20000 附近的更高值。

您可以通过计算两个区间的积分来解决此问题，例如 [0, mode] 和 [mode, np.inf]。 （不需要计算负轴上的积分，因为那里的 PDF 为 0。）

例如，这个脚本打印1.0000000000000004

import numpy as np
from scipy.stats import lognorm
from scipy.integrate import quad


def f(x, mu=0, sigma=1):
    return lognorm.pdf(x, s=sigma, loc=0, scale=np.exp(mu))


mu = 10
sigma = 0.2

mode = np.exp(mu - sigma**2)

ans1, err1 = quad(f, 0, mode, args=(mu, sigma))
ans2, err2 = quad(f, mode, np.inf, args=(mu, sigma))

integral = ans1 + ans2
print(integral)