scipy 馈入数组与逐个元素馈入时给出不同的答案

Question

我有这个代码片段：

print(np.sqrt(scipy.stats.chi2.ppf(1-prob,1)))

for i in prob:
    print(np.sqrt(scipy.stats.chi2.ppf(1-i,1)))

我认为这些应该给我相同的答案，但我得到的答案是：

[0. inf 2.27834698 0.22780442 2.21905125]

和

0.0
6.1833132853181185
2.2783471062868474
0.22780441954248226
2.219051610822544

打印概率给了我 [1.0000000e+00 6.2769967e-10 2.2705905e-02 8.1979829e-01 2.6483214e-02]，这很奇怪，因为当我复制这个数组并在控制台中执行 top 命令时 python，我可以得到正确的答案，但是运行脚本总是给我 inf.

我是运行 python 3.9.7，scipy 1.7.3，numpy 1.21.4。 prob.dtype 给我 float32，分别打印数组中每个元素的类型也给我 float32。

有人以前见过这样的东西吗？

Answer 1

我将 prob 转换为 float64 数组并且它起作用了。 float32 的精度似乎有些有趣。

Answer 2

根据问题中显示的值，prob 是：

In [155]: prob = np.array([1.0000000e+00, 6.2769967e-10, 2.2705905e-02, 8.1979829e-01, 2.6483214e-02], dtype=np.float32)

In [156]: prob
Out[156]: 
array([1.0000000e+00, 6.2769967e-10, 2.2705905e-02, 8.1979829e-01,
       2.6483214e-02], dtype=float32)

你给scipy.stats.chi2.ppf的值实际上是1 - prob，这个减法就是问题所在：

In [174]: 1 - prob
Out[174]: 
array([0.        , 1.        , 0.9772941 , 0.18020171, 0.97351676],
      dtype=float32)

请注意，结果中的第二个值为 1.0。这是因为 6.2769967e-10 比 np.float32 的“机器 epsilon”小得多，后者约为 1.19e-7。换句话说，6.27e-10 小于 1.0 周围 np.float32 表示的分辨率。当您将 1 - prob 传递给 ppf() 方法时，输入值 1.0 的结果是 inf（这是正确的）。

正如您在回答中已经指出的那样，您可以通过将 prob 转换为 np.float64 来避免该问题。避免此问题的另一种方法是使用 isf(prob, 1) 而不是 ppf(1 - prob, 1):

In [177]: print(np.sqrt(scipy.stats.chi2.isf(prob, 1)))
[0.         6.18331329 2.27834711 0.22780442 2.21905161]

isf 是 逆生存函数 。该表达式在数学上等同于 ppf(1 - prob, 1)，但它避免了导致精度极度损失的减法。

scipy 馈入数组与逐个元素馈入时给出不同的答案

scipy giving different answer when feeding array vs feeding element by element

python

scipy

scipy.stats