峰度检验#
峰度检验 scipy.stats.kurtosistest
函数检验零假设:样本所来自的总体峰度与正态分布的峰度相同。
假设我们希望从测量值推断一项医学研究中成年男性人类的体重是否不是正态分布的 [1]。体重(磅)记录在下面的数组 x
中。
import numpy as np
x = np.array([148, 154, 158, 160, 161, 162, 166, 170, 182, 195, 236])
来自 [2] 的峰度检验首先计算一个基于样本(超额/Fisher)峰度的统计量。
from scipy import stats
res = stats.kurtosistest(x)
res.statistic
/home/treddy/github_projects/scipy/doc/build/env/lib/python3.10/site-packages/scipy/stats/_axis_nan_policy.py:586: UserWarning: `kurtosistest` p-value may be inaccurate with fewer than 20 observations; only n=11 observations were given.
res = hypotest_fun_out(*samples, **kwds)
2.3048235214240873
(该测试警告说我们的样本观测值太少,无法执行测试。我们将在示例末尾回到这一点。)由于正态分布的超额峰度为零(根据定义),因此对于从正态分布中抽取的样本,此统计量的大小往往较低。
通过将观察到的统计量值与零分布进行比较来执行测试:零分布是在零假设(权重是从正态分布中提取的)下导出的统计量值的分布。
对于此测试,非常大样本的统计量零分布是标准正态分布。
import matplotlib.pyplot as plt
dist = stats.norm()
kt_val = np.linspace(-5, 5, 100)
pdf = dist.pdf(kt_val)
fig, ax = plt.subplots(figsize=(8, 5))
def kt_plot(ax): # we'll reuse this
ax.plot(kt_val, pdf)
ax.set_title("Kurtosis Test Null Distribution")
ax.set_xlabel("statistic")
ax.set_ylabel("probability density")
kt_plot(ax)
plt.show()
比较由 p 值量化:零分布中与观察到的统计量值一样极端或更极端的值的比例。在统计量为正的双尾测试中,大于观察到的统计量的零分布元素和小于观察到的统计量的负值的零分布元素都被认为是“更极端”。
fig, ax = plt.subplots(figsize=(8, 5))
kt_plot(ax)
pvalue = dist.cdf(-res.statistic) + dist.sf(res.statistic)
annotation = (f'p-value={pvalue:.3f}\n(shaded area)')
props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
_ = ax.annotate(annotation, (3, 0.005), (3.25, 0.02), arrowprops=props)
i = kt_val >= res.statistic
ax.fill_between(kt_val[i], y1=0, y2=pdf[i], color='C0')
i = kt_val <= -res.statistic
ax.fill_between(kt_val[i], y1=0, y2=pdf[i], color='C0')
ax.set_xlim(-5, 5)
ax.set_ylim(0, 0.1)
plt.show()
res.pvalue
0.0211764592113868
如果 p 值“很小”,即,如果从正态分布总体中抽样数据并产生如此极端的统计量值的概率很低,则这可以作为反对零假设而支持备择假设的证据:权重不是从正态分布中提取的。请注意
反之则不然;也就是说,该测试不用于为零假设提供证据。
将被视为“小”的值的阈值是一个应该在分析数据之前做出的选择 [3],同时考虑到假阳性(错误地拒绝零假设)和假阴性(未能拒绝错误的零假设)的风险。
请注意,标准正态分布提供了零分布的渐近近似;它仅对具有许多观测值的样本准确。这就是我们在示例开头收到警告的原因;我们的样本很小。在这种情况下,scipy.stats.monte_carlo_test
可以提供更准确的,尽管是随机的,对精确 p 值的近似。
def statistic(x, axis):
# get just the skewtest statistic; ignore the p-value
return stats.kurtosistest(x, axis=axis).statistic
res = stats.monte_carlo_test(x, stats.norm.rvs, statistic)
fig, ax = plt.subplots(figsize=(8, 5))
kt_plot(ax)
ax.hist(res.null_distribution, np.linspace(-5, 5, 50),
density=True)
ax.legend(['asymptotic approximation\n(many observations)',
'Monte Carlo approximation\n(11 observations)'])
plt.show()
/home/treddy/github_projects/scipy/doc/build/env/lib/python3.10/site-packages/scipy/stats/_axis_nan_policy.py:586: UserWarning: `kurtosistest` p-value may be inaccurate with fewer than 20 observations; only n=11 observations were given.
res = hypotest_fun_out(*samples, **kwds)
/home/treddy/github_projects/scipy/doc/build/env/lib/python3.10/site-packages/scipy/stats/_axis_nan_policy.py:618: UserWarning: `kurtosistest` p-value may be inaccurate with fewer than 20 observations; only n=11 observations were given.
res = hypotest_fun_out(*samples, axis=axis, **kwds)
res.pvalue
0.0228
此外,尽管具有随机性质,但以这种方式计算的 p 值可用于精确控制零假设的错误拒绝率 [4]。