Python中卡方检验统计量的P值

Question

我计算了一个检验统计量，其服从自由度为1的卡方分布，并且想要使用Python找出对应的P值。作为Python和数学/统计方面的新手，我认为我需要使用SciPy中卡方分布的概率密度函数。但是，当我像这样使用它时：

from scipy import stats
stats.chi2.pdf(3.84 , 1)
0.029846

然而，一些谷歌搜索和与一些了解数学但不了解Python的同事交流后，他们说应该是0.05。

有什么想法吗？ 谢谢， Davy

- Davy Kavanagh

---

如果您使用scipy.stats.chisquare运行测试，是否能得到所需的结果？ - Fred Foo

顺便提一句，当我按照维基百科计算概率密度函数时，我得到的结果与 SciPy 相同：x = 3.84; reciprocal(2**.5 * gamma(.5)) * x ** (.5 - 1) * exp(- x / 2)。 - Fred Foo

我认为你正在使用错误的函数……正如@larsmans所提到的，你应该使用chisquare函数，但要确保将实际值和期望值的数组传递给它，它会返回你所寻找的3.84和p-value。 - ernie

我一开始没有做这个的原因是我没有实际和预期数组。这是一个同事开发的新型分析方法。它遵循卡方分布，但不是经典的卡方列联表检验。所以我认为我不能使用它。我需要能够仅从测试统计量中推导出p值。 - Davy Kavanagh

那么，pval = 1 - stats.chi2.cdf(3.84, 1) 这一行代码怎么样？（我在这个帖子中看到的）。 - ernie

甚至更好的是 >>> stats.chi2.sf(3.84, 1) 0.050043521248705106，以增加尾部精度为代价 - Josef

Answer 1

这里是一个简要回顾：

概率密度函数：可以将其视为点值；在给定点处的概率有多密集？

累积分布函数：这是函数在给定点之前的概率质量；分布百分之几位于该点的一侧？

在您的情况下，您使用了概率密度函数（PDF），并得出了正确答案。如果您尝试 1 - CDF：

>>> 1 - stats.chi2.cdf(3.84, 1)
0.050043521248705147

概率密度函数 累积分布函数

- learner

---

这并没有真正回答问题，因为你假设op对于p值的右侧配置感兴趣。 - Tobbey

Answer 2

更新：请注意，从scipy版本0.17.0开始，chisqprob()已被弃用。现在可以通过scipy.stats.distributions.chi2.sf()来获得高精度的卡方值，例如：

>>>from scipy.stats.distributions import chi2
>>>chi2.sf(3.84,1)
0.050043521248705189
>>>chi2.sf(1424,1)
1.2799986253099803e-311

---

虽然对于小的卡方值，stats.chisqprob()和1-stats.chi2.cdf()看起来相似，但是对于大的卡方值，前者更可取。后者无法提供比机器精度更小的p值，并且在接近机器精度时会给出非常不准确的答案。正如其他人所示，对于小的卡方值，两种方法得到的结果是可比较的：

>>>from scipy.stats import chisqprob, chi2
>>>chisqprob(3.84,1)
0.050043521248705189
>>>1 - chi2.cdf(3.84,1)
0.050043521248705147

使用1-chi2.cdf()在这里出现问题：

>>>1 - chi2.cdf(67,1)
2.2204460492503131e-16
>>>1 - chi2.cdf(68,1)
1.1102230246251565e-16
>>>1 - chi2.cdf(69,1)
1.1102230246251565e-16
>>>1 - chi2.cdf(70,1)
0.0

虽然 chisqprob() 可以为更大范围的卡方值提供准确结果，产生的 p 值几乎与大于零的最小浮点数一样小，但它也会在下溢之前停止计算：

>>>chisqprob(67,1)
2.7150713219425247e-16
>>>chisqprob(68,1)
1.6349553217245471e-16
>>>chisqprob(69,1)
9.8463440314253303e-17    
>>>chisqprob(70,1)
5.9304458500824782e-17
>>>chisqprob(500,1)
9.505397766554137e-111
>>>chisqprob(1000,1)
1.7958327848007363e-219
>>>chisqprob(1424,1)
1.2799986253099803e-311
>>>chisqprob(1425,1)
0.0

- Christophe Lambert

---

一个非常好的回答 :) - jb.

一个快速的提示：chisqprob已被弃用，http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chisqprob.html#scipy.stats.chisqprob - ZK Zhao

谢谢 @cqcn1991 -- 我已经更新了解决方案。 - Christophe Lambert

Answer 3

要计算给定卡方和和自由度时的零假设概率，您还可以调用chisqprob：

>>> from scipy.stats import chisqprob
>>> chisqprob(3.84, 1)
0.050043521248705189

---

注意：

chisqprob已经被弃用！自scipy 0.17.0版本开始，使用stats.distributions.chi2.sf代替stats.chisqprob。

- jb.

---

感谢您关于chi2.sf()的说明。 - Ivan

Answer 4

您想要的是：

您想执行的操作是：

>>> 1 - stats.chi2.cdf(3.84, 1)
0.050043521248705147

- David Robinson

Answer 5

其他一些解决方案已被废弃，请使用scipy.stats.chi2生存函数，其与1 - cdf(chi_statistic, df)相同。

示例：

from scipy.stats import chi2
p_value = chi2.sf(chi_statistic, df)

- SARose

Answer 6

如果你想了解数学知识，样本x的p值（固定）是：

P[P(X) <= P(x)] = P[m(X) >= m(x)] = 1 - G(m(x)^2)

其中，

P是已知协方差（cov）和均值的（例如k元）正态分布的概率， X是来自该正态分布的随机变量， m(x)是马氏距离= sqrt(.请注意，在1-d中，这只是z-score的绝对值。 G是具有k个自由度的chi^2分布的CDF。

因此，如果你正在计算固定观察值x的p值，那么你需要计算m(x)（广义z-score），并且是1-G(m(x)^2)。

例如，众所周知，如果x是从单变量（k = 1）正态分布中抽样的，并且具有z-score = 2（它距平均值2个标准偏差），则p值约为.046（请参见z-score表）。

In [7]: from scipy.stats import chi2

In [8]: k = 1

In [9]: z = 2

In [10]: 1-chi2.cdf(z**2, k)
Out[10]: 0.045500263896358528

- travelingbones

Answer 7

如果需要超高精度，当scipy的chi2.sf()无法满足需求时，就要动用大杀器了：

>>> import numpy as np
>>> from rpy2.robjects import r
>>> np.exp(np.longdouble(r.pchisq(19000, 2, lower_tail=False, log_p=True)[0]))
1.5937563168532229629e-4126

由另一位用户（WestCoastProjects）更新使用OP中的值时，我们得到：

np.exp(np.longdouble(r.pchisq(3.84,1, lower_tail=False, log_p=True)[0]))
Out[5]: 0.050043521248705198928

所以这就是你要找的那个0.05。

- 1''

---

是的，R 是统计学中的大杀器。 - WestCoastProjects

但是p值不是P(chisquared >= x)而是P(chisquared > x)吗？ - user3494047