如何为非常小的y值计算积分SciPy quad

问题描述

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

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)

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

解决方法

您使用的值对应于具有均值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,mode,args=(mu,sigma))
ans2,err2 = quad(f,np.inf,sigma))

integral = ans1 + ans2
print(integral)