偏度检验#

此函数检验样本所来自的总体偏度是否与相应的正态分布相同的零假设。

假设我们希望从测量结果推断医学研究中成年男性的人体重量是否不呈正态分布 [1]。 重量(磅)记录在下面的数组 x 中。

import numpy as np
x = np.array([148, 154, 158, 160, 161, 162, 166, 170, 182, 195, 236])

来自 [2] 的偏度检验 scipy.stats.skewtest 首先基于样本偏度计算一个统计量。

from scipy import stats
res = stats.skewtest(x)
res.statistic
2.7788579769903414

由于正态分布的偏度为零,因此对于从正态分布中抽取的样本,此统计量的大小往往较小。

通过将观察到的统计量值与零分布进行比较来执行测试:零分布是在零假设(即权重是从正态分布中抽取的)下得出的统计量值的分布。

对于此检验,非常大样本的统计量的零分布是标准正态分布。

import matplotlib.pyplot as plt
dist = stats.norm()
st_val = np.linspace(-5, 5, 100)
pdf = dist.pdf(st_val)
fig, ax = plt.subplots(figsize=(8, 5))

def st_plot(ax):  # we'll reuse this
    ax.plot(st_val, pdf)
    ax.set_title("Skew Test Null Distribution")
    ax.set_xlabel("statistic")
    ax.set_ylabel("probability density")

st_plot(ax)
plt.show()
../../_images/5787b267e279ab2972d6063d4dcf22fbe18f8e8712b30e447381527ff54b42cd.png

比较由 p 值量化:零分布中与观察到的统计量值一样极端或更极端的值的比例。在双侧检验中,大于观察到的统计量的零分布的元素和小于观察到的统计量的负数的零分布的元素都被认为是“更极端”的。

fig, ax = plt.subplots(figsize=(8, 5))
st_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 = st_val >= res.statistic
ax.fill_between(st_val[i], y1=0, y2=pdf[i], color='C0')
i = st_val <= -res.statistic
ax.fill_between(st_val[i], y1=0, y2=pdf[i], color='C0')
ax.set_xlim(-5, 5)
ax.set_ylim(0, 0.1)
plt.show()
../../_images/539c9a141bd83ae52193ba114a45bbfd7c2f0bb11e12606ec20f867f927a98d9.png
res.pvalue
0.005455036974740185

如果 p 值“小”——也就是说,如果从正态分布总体中采样产生如此极端的统计量值的概率很低——这可以被视为不利于零假设而支持替代假设的证据:权重不是从正态分布中抽取的。请注意,

  • 反之则不然;也就是说,该检验不用于为零假设提供证据。

  • 将被视为“小”的值的阈值是一个选择,应在分析数据之前做出 [3],并考虑假阳性(错误地拒绝零假设)和假阴性(未能拒绝错误的零假设)的风险。

请注意,标准正态分布提供了零分布的渐近近似;它仅对具有许多观察值的样本准确。对于我们这样的小样本,scipy.stats.monte_carlo_test 可以提供更准确的,尽管是随机的,对精确 p 值的近似。

def statistic(x, axis):
    # get just the skewtest statistic; ignore the p-value
    return stats.skewtest(x, axis=axis).statistic

res = stats.monte_carlo_test(x, stats.norm.rvs, statistic)
fig, ax = plt.subplots(figsize=(8, 5))
st_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()
../../_images/f7a5652423cde0d7c4134a9ff2f8e9c3db5d4f1ceb2c767d202af661ede39230.png
res.pvalue
0.0042

在这种情况下,即使对于我们的小样本,渐近近似和蒙特卡罗近似也相当接近。

参考文献#