scipy.stats.norm.interval返回inf

如何解决scipy.stats.norm.interval返回inf

我有很多数据点和一个模型，我现在正在寻找可以排除该模型的置信度。

我从卡方值开始，再加上自由度数。

使用scipy，我一直在尝试以下内容：

from scipy.stats import chi2
from scipy.stats import norm

chisq = 74.1
df = 21 #degrees of freedom

cdf = chi2.cdf(chisq,df,scale=1)
sigma = norm.interval(cdf)

这将产生我想要的输出（5.5 sigma置信度），但是对于较高的chi ^ 2值将失败。

在我看来，问题在于scipy使用的数据类型的精度。对于chi ^ 2值334.7（和相同的自由度数），chi2.cdf得出1.0，这又使norm.interval返回-inf / inf。最好返回16.1的西格玛置信度。

我找到了一种使用mpmath提高cdf精度的方法：

import mpmath

mpmath.mp.dps = 200 # decimal digits of precision

def cdf(x,k):
    x,k = mpmath.mpf(x),mpmath.mpf(k)
    return mpmath.gammainc(k/2,x/2,regularized=True)

但是，问题仍然存在，因为scipy.stats.norm.interval将输入舍入为1.0。

是否有一种方法可以通过更改数据类型norm.interval的使用来规避此问题，或者是否有其他方法/软件包来计算具有任意精度输入的正态分布的端点？