如何解决如何为非常小的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)
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。