如何计算非常小的 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 值 (max ~ 1E-5) 并分布在 ~1E5 的 x 值上。我们知道一个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)

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

How to calculate integral for very very small y values (SciPy quad)

python

integral

scipy

numerical-methods

numerical-integration