Scipy lognormal CFD和SF大部分时间都是0

Question

根据Scipy reference的说法，给定一个lognormal distribution，CDF给出X等于或低于某个特定值的概率，而SF给出X高于某个特定值的概率。

因为我的发行版以x=195为中心，所以我在调用这些函数时使用了loc参数。我没有使用scale参数。问题可能在这里，因为对这两个参数如何工作的解释不好。

因此，我编写了一小段代码来确定x在两个值x1和x2之间的概率，其中x1<x2

from scipy.stats import lognorm

sigma = 0.149
mu = 195
x1 = 188.78
x2 = 201.68
p_below = lognorm.cdf(x1, sigma, loc=mu)
p_above = lognorm.sf(x2, sigma, loc=mu)
p_between = 1 - p_above - p_below

我要说的是，在方差如此低的x=195，x低于x1或高于x2的概率将很高。尽管如此，cdf和sf返回的值大多数时候是零。我只在将sigma降低到0.068时才成功地获得了大于零的sf结果，这很奇怪(方差越大，越容易越界)。此外，我永远不会期望cdf或sf返回零，而是返回一个非常低的正float。

你知道我做错了什么吗？

Answer 1

让我们看看：

In [41]: st.lognorm.sf(x2-mu, s=sigma)
Out[41]: 0.0

In [42]: a = np.log(x2-mu) / sigma

In [43]: st.norm.sf(a)
Out[43]: 1.6458720667149256e-37

In [44]: st.norm.logsf(a)
Out[44]: -84.697378065290948

但

In [45]: st.norm.cdf(a)
Out[45]: 1.0

查看source code，可以清楚地看到scipy.lognorm的实现可以通过定义显式的_sf、_logcdf和_logsf方法来改进，以转发到正态分布的相应方法，就像_cdf一样。

目前，解决方法是手动获取日志并使用scipy.stats.norm方法。

编辑:这是scipy票证：https://github.com/scipy/scipy/issues/5940 (欢迎使用补丁！)

EDIT2:已在scipy主分支中修复此问题：

In [19]: scipy.__version__
Out[19]: '0.18.0.dev0+f63d0a6'

In [20]: lognorm.sf(x2-mu, s=sigma)
Out[20]: 1.6458720667149256e-37

EDIT3: OP中的代码片段中仍然存在问题:浮点中的1 - 1.e-37等于1。

Answer 2

我在scipy中遇到了同样的问题，即sf和cdf函数不是很准确，在变为0之前只允许~10^-30的概率。作为另一种选择，您可以尝试执行数值积分，可能是使用scipy/numpy或mpmath (http://mpmath.org/doc/current/calculus/integration.html)。

或者，您可以尝试重新表述您的问题，以便您可以使用另一个统计数据进行测试。这就是我解决问题的方法。然而，这可能并不总是一种选择。