python - Scipy erfcinv 在 1e-16 附近意外爆炸

Question

我一直在使用 scipy.special.erfcinv 从 pvalues 计算 Z 分数。然而，当 pvalues 变得非常小时，erfcinv 会意外地变大。有什么想法吗？

例子:

In [1]: import numpy as np
In [2]: from scipy.special import erfcinv
In [3]: erfcinv(2e-16) * np.sqrt(2)
Out[3]: 8.2095361516013874
In [4]: erfcinv(1e-16) * np.sqrt(2)
Out[4]: 1.7976931348623155e+308

我正在使用 scipy 0.10.1 运行 python 2.6.6。

Answer 1

简短的回答:你越接近浮点运算精度的极限，就会发生越奇怪的事情:(http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html，http://www.seas.ucla.edu/~vandenbe/103/lectures/flpt.pdf)

有点长。首先让我们看看 erfcinv 函数:

def erfcinv(y):
    return ndtri((2-y)/2.0)/sqrt(2)

如果我们取 y = 2e-16:

In [96]: (2 - 2e-16) / 2
Out[96]: 0.9999999999999999

当我们取 y = 1e-16 时:

In [97]: (2 - 1e-16) / 2
Out[97]: 1.0

现在我们看看 ndtri:

x=ndtri(y) returns the argument x for which the area udnder the
Gaussian probability density function (integrated from minus infinity
to x) is equal to y.

现在一切都应该清楚了，对吗？正如您所怀疑的:

In [99]: ndtri(1)
Out[99]: inf

您的结果可能会有点不同 - 就我而言:

In [101]: erfcinv(1e-16) * np.sqrt(2)
Out[101]: inf