scipy中非中心卡方分布实现的准确性答案

【问题标题】：Accuracy of non-central chi-squared distribution implementation in scipyscipy中非中心卡方分布实现的准确性
【发布时间】：2018-04-11 08:16:53
【问题描述】：

在 scipy (scipy.stats.ncx2, scipy.special.chndtr) 中的非中心卡方分布的实现对于非中心参数的大值是不准确的。例如。考虑

from scipy.stats import ncx2 

print(ncx2.cdf(1200,2,1000))
print(ncx2.cdf(1500,2,1000))
print(ncx2.cdf(2000,2,1000))
print(ncx2.cdf(5000,2,1000))

产生

0.998604279948
0.999933004449
0.999933004449
0.999933004449

scipy 中的实现基于 Abramowitz 和 Stegun 的公式 26.4.25，Handbook of Mathematical Functions (1966) （见https://github.com/scipy/scipy/blob/master/scipy/special/cdflib/cdfchn.f）。

有更快、更准确的算法来实现这个功能，尤其是对于大的非中心性。例如，R 显然使用了 Ding (1992)（参见 http://stat.ethz.ch/R-manual/R-devel/library/stats/html/Chisquare.html）。我想知道是否有人知道替代实现是否可以作为标准 python 库的一部分？

【问题讨论】：

也许在 SciPy 跟踪器中打开一个问题是在 Python 科学堆栈中实现它的最佳机会。
好的，就这么做了。

标签： python scipy

【解决方案1】：

回答我自己的问题，以防它对某人有用：可以使用 rpy2 从 python 访问非中心卡方分布的 R 实现，如下所示（改编上面的示例）：

import rpy2.robjects as ro
from rpy2.robjects.packages import importr

stats = importr('stats')

print(stats.pchisq(1200,2,1000)[0])
print(stats.pchisq(1500,2,1000)[0])
print(stats.pchisq(2000,2,1000)[0])
print(stats.pchisq(5000,2,1000)[0])

输出是

0.998663933426865
1.0
1.0
1.0

【讨论】：