47

我已经计算了一个测试统计量，该统计量分布为一个自由度为 1 的卡方，并且想找出使用​​ python 对应的 P 值。

我是 python 和数学/统计新手，所以我想我想要的是 SciPy 中 chi2 分布的概率密度函数。但是，当我这样使用它时：

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

然而，一些谷歌搜索和一些懂数学但不懂python的同事说它应该是0.05。

有任何想法吗？干杯，戴维

4

7 回答 7

60

在这里快速复习：

概率密度函数：将其视为一个点值；给定点的概率有多密集？

累积分布函数：这是函数到给定点的概率质量；分布的百分比在这一点的一侧？

在您的情况下，您选择了 PDF，并得到了正确答案。如果您尝试 1 - CDF：

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

PDF CDF

于 2012-07-30T19:14:30.527 回答
27

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

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

注意：

chisqprob已弃用！stats.chisqprob 在 scipy 0.17.0 中已弃用；改用stats.distributions.chi2.sf _

于 2013-11-23T17:25:02.260 回答
25

更新：如前所述， 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() 对于小的卡方值看起来相当，但对于大的卡方值，前者更可取。后者无法提供小于机器 epsilon 的 p 值，并且会给出接近机器 epsilon 的非常不准确的答案。如其他人所示，使用两种方法可得出较小卡方值的可比较值：

>>>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
于 2015-05-22T16:15:33.940 回答
7

你的意思是：

>>> 1 - stats.chi2.cdf(3.84, 1)
0.050043521248705147
于 2012-07-30T19:19:13.027 回答
5

其他一些解决方案已被弃用。使用scipy.stats.chi2生存功能。这与1 - cdf(chi_statistic, df)

例子：

from scipy.stats import chi2
p_value = chi2.sf(chi_statistic, df)
于 2016-01-23T16:03:20.027 回答
3

如果您想了解数学，样本的 p 值 x（固定）为

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

在哪里，

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

因此，如果您要计算固定观测值 x 的 p 值，则计算 m(x)（广义 z 分数）和 1-G(m(x)^2)。

例如，众所周知，如果 x 是从单变量 (k = 1) 正态分布中采样的，并且 z 分数 = 2（与平均值相差 2 个标准差），那么 p 值约为 0.046（参见 a z 分数表） 

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
于 2016-10-11T20:01:01.730 回答
2

对于超高精度，当 scipy'schi2.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 。

于 2020-08-19T06:24:06.497 回答