Python SciPy stats.quantile_test用法及代码示例

本文简要介绍 python 语言中 scipy.stats.quantile_test 的用法。

用法:

scipy.stats.quantile_test(x, *, q=0, p=0.5, alternative='two-sided')#

执行分位数测试并计算分位数的置信区间。

此函数测试原假设，即 q 是与样本 x 的总体概率 p 相关的分位数值。例如，使用默认参数，它测试 x 背后的总体中位数是否为零。该函数返回一个对象，其中包括检验统计量、p 值以及计算分位数周围置信区间的方法。

参数 ：：

x： array_like

一维样本。

q： 浮点数，默认值：0

分位数的假设值。

p： 浮点数，默认值：0.5

与分位数相关的概率；即小于 q 的人口比例为 p。必须严格介于 0 和 1 之间。

alternative： {‘双面’，‘less’, ‘greater’}，可选

定义备择假设。可以使用以下选项(默认为“双面”)：

• “双面”：与概率 p 相关的分位数不是 q。

• ‘less’：与概率 p 小于 q 相关的分位数。

• ‘greater’：与概率 p 大于 q 相关的分位数。

返回 ：：

result： QuantileTestResult

具有以下属性的对象：

统计 浮点数

可用于分位数检验的两个检验统计量之一。第一个检验统计量，T1, 是样本的比例x小于或等于假设的分位数q。第二个检验统计量，T2, 是样本的比例x严格小于假设的分位数q.

当 alternative = 'greater' 、 T1 用于计算 p 值且 statistic 设置为 T1 时。

当 alternative = 'less' 、 T2 用于计算 p 值且 statistic 设置为 T2 时。

当 alternative = 'two-sided' 时，同时考虑 T1 和 T2，并使用导致最小 p 值的那个。

statistic_type int

任何一个1或者2取决于哪一个T1或者T2用于计算 p 值。

p值 浮点数

与给定替代方案相关的 p 值。

该对象还具有以下方法：

confidence_interval(confidence_level=0.95)

计算与概率相关的总体分位数周围的置信区间p。置信区间以namedtuple带字段低的和高的。值为楠当没有足够的观测值来计算所需置信度的置信区间时。

注意：

该检验及其计算置信区间的方法是非参数的。当且仅当观察值独立同分布时它们才有效。

测试的实施遵循 Conover [1]。考虑两个检验统计量。

T1： 中的观测值数量x小于或等于q.

T1 = (x <= q).sum()

T2： 中的观测值数量x严格小于q.

T2 = (x < q).sum()

有必要使用两个检验统计量来处理 x 是从离散分布或混合分布生成的可能性。

检验的原假设是：

H0: The \(p^{\mathrm{th}}\) population quantile is q.

每个测试统计量的零分布是 \(\mathrm{binom}\left(n, p\right)\) 。当 alternative='less' 时，备择假设为：

H1: The \(p^{\mathrm{th}}\) population quantile is less than q.

p 值是二项式随机变量的概率

\[Y \sim \mathrm{binom}\left(n, p\right)\]

大于或等于观测值T2。

当 alternative='greater' 时，备择假设为：

H1: The \(p^{\mathrm{th}}\) population quantile is greater than q

p 值是二项式随机变量 Y 小于或等于观测值 T1 的概率。

当 alternative='two-sided' 时，备择假设为

H1: q is not the \(p^{\mathrm{th}}\) population quantile.

对于 'less' 和 'greater' 情况，p 值是较小 p 值的两倍。对于相同的数据，这两个 p 值都可能超过 0.5，因此该值被剪裁到区间 \([0, 1]\) 中。

置信区间的方法归功于 Thompson [2]，后来被证明适用于任何 i.i.d 集合。样品[3]。该计算基于以下观察：分位数 \(q\) 大于任何观察值 \(x_m (1\leq m \leq N)\) 的概率可以计算为

\[\mathbb{P}(x_m \leq q) = 1 - \sum_{k=0}^{m-1} \binom{N}{k} q^k(1-q)^{N-k}\]

默认情况下，置信区间是针对 95% 的置信水平计算的。对 95% 置信区间的常见解释是，如果 i.i.d.从同一总体中重复抽取样本，每次都会形成置信区间，置信区间将包含大约 95% 的试验中指定分位数的真实值。

QuantileNPCI R 包 [4] 中提供了类似的函数。基础是相同的，但它通过在样本值之间进行插值来计算置信区间边界，而该函数仅使用样本值作为边界。因此，quantile_test.confidence_interval 返回更保守的间隔(即更大)。

confintr 包 [5] 中包含相同的分位数置信区间计算。

不能保证两侧置信区间是最优的；即，可能存在更紧密的区间，其中可能包含概率大于置信水平的感兴趣分位数。在不对样本进行进一步假设(例如，基础分布的性质)的情况下，单侧间隔是最佳的紧的。

参考：

[1]

23. 10. 科诺弗。实用非参数统计，第三版。 1999 年。

[2]

W. R. Thompson，“关于未知分布形式总体的中位数和其他期望分布的置信范围”，《数理统计年鉴》，卷。 7、没有。 3，第 122-128 页，1936 年，访问时间：2019 年 9 月 18 日。[在线]。可用：https://www.jstor.org/stable/2957563。

[3]

H. A. David 和 H. N. Nagaraja，“订单统计中的非参数推理中的订单统计”，John Wiley & Sons, Ltd，2005 年，第 159-170 页。可用的：https://onlinelibrary.wiley.com/doi/10.1002/0471722162.ch7.

[4]

N. Hutson、A. Hutson、L. Yan，“QuantileNPCI：分位数的非参数置信区间”，R 包，https://cran.r-project.org/package=QuantileNPCI

[5]

M. Mayer，“confintr：置信区间”，R 包，https://cran.r-project.org/package=confintr

例子：

假设我们希望检验总体中位数等于 0.5 的原假设。我们选择99%的置信度；也就是说，如果 p 值小于 0.01，我们将拒绝原假设并支持替代假设。

当测试中位数为 0.5 的标准均匀分布的随机变量时，我们期望数据在大多数情况下与原假设一致。

>>> import numpy as np
>>> from scipy import stats
>>> rng = np.random.default_rng()
>>> rvs = stats.uniform.rvs(size=100, random_state=rng)
>>> stats.quantile_test(rvs, q=0.5, p=0.5)
QuantileTestResult(statistic=45, statistic_type=1, pvalue=0.36820161732669576)

正如预期的那样，p 值不低于我们的阈值 0.01，因此我们不能拒绝原假设。

当测试中位数为 0 的标准正态分布数据时，我们预计零假设会被拒绝。

>>> rvs = stats.norm.rvs(size=100, random_state=rng)
>>> stats.quantile_test(rvs, q=0.5, p=0.5)
QuantileTestResult(statistic=67, statistic_type=2, pvalue=0.0008737198369123724)

事实上，p 值低于我们的阈值 0.01，因此我们拒绝零假设，转而支持默认的 “two-sided” 替代方案：总体中位数不等于 0.5。

然而，假设我们要针对总体中位数大于 0.5 的片面替代方案检验原假设。由于标准正态分布的中位数小于 0.5，因此我们预计原假设不会被拒绝。

>>> stats.quantile_test(rvs, q=0.5, p=0.5, alternative='greater')
QuantileTestResult(statistic=67, statistic_type=1, pvalue=0.9997956114162866)

毫不奇怪，当 p 值大于我们的阈值时，我们不会拒绝零假设而支持所选的替代方案。

分位数检验可用于任何分位数，而不仅仅是中位数。例如，我们可以测试样本分布的第三个四分位数是否大于 0.6。

>>> rvs = stats.uniform.rvs(size=100, random_state=rng)
>>> stats.quantile_test(rvs, q=0.6, p=0.75, alternative='greater')
QuantileTestResult(statistic=64, statistic_type=1, pvalue=0.00940696592998271)

p 值低于阈值。我们拒绝原假设，转而采用另一种假设：样本分布的第三个四分位数大于 0.6。

quantile_test 还可以计算任何分位数的置信区间。

>>> rvs = stats.norm.rvs(size=100, random_state=rng)
>>> res = stats.quantile_test(rvs, q=0.6, p=0.75)
>>> ci = res.confidence_interval(confidence_level=0.95)
>>> ci
ConfidenceInterval(low=0.284491604437432, high=0.8912531024914844)

当测试单边替代方案时，置信区间包含所有观察结果，如果作为 q 传递，测试的 p 值将大于 0.05，因此不会拒绝原假设。例如：

>>> rvs.sort()
>>> q, p, alpha = 0.6, 0.75, 0.95
>>> res = stats.quantile_test(rvs, q=q, p=p, alternative='less')
>>> ci = res.confidence_interval(confidence_level=alpha)
>>> for x in rvs[rvs <= ci.high]:
...     res = stats.quantile_test(rvs, q=x, p=p, alternative='less')
...     assert res.pvalue > 1-alpha
>>> for x in rvs[rvs > ci.high]:
...     res = stats.quantile_test(rvs, q=x, p=p, alternative='less')
...     assert res.pvalue < 1-alpha

此外，如果为随机样本重复生成 95% 置信区间，则置信区间将在大约 95% 的重复中包含真实分位数值。

>>> dist = stats.rayleigh() # our "unknown" distribution
>>> p = 0.2
>>> true_stat = dist.ppf(p) # the true value of the statistic
>>> n_trials = 1000
>>> quantile_ci_contains_true_stat = 0
>>> for i in range(n_trials):
...     data = dist.rvs(size=100, random_state=rng)
...     res = stats.quantile_test(data, p=p)
...     ci = res.confidence_interval(0.95)
...     if ci[0] < true_stat < ci[1]:
...         quantile_ci_contains_true_stat += 1
>>> quantile_ci_contains_true_stat >= 950
True

只要样本是独立同分布的，这适用于任何分布和任何分位数。