峰度检验#
峰度检验函数 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
np.float64(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 值来量化: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
np.float64(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()

res.pvalue
np.float64(0.0246)
此外,尽管它们具有随机性,但以这种方式计算的 p 值可以精确控制零假设的错误拒绝率 [4]。