Jarque-Bera 拟合优度检验#
假设我们希望从测量数据推断一项医学研究中成年男性的人体体重是否不服从正态分布[1]。体重（磅）记录在下面的数组 x 中。
import numpy as np
x = np.array([148, 154, 158, 160, 161, 162, 166, 170, 182, 195, 236])
Jarque-Bera 检验scipy.stats.jarque_bera 首先计算一个基于样本偏度和峰度的统计量。
from scipy import stats
res = stats.jarque_bera(x)
res.statistic
6.982848237344646
由于正态分布的偏度为零，并且（“超额”或“Fisher”）峰度为零，因此对于从正态分布中抽取的样本，此统计量的值往往较低。
该检验是通过将观测到的统计量的值与零分布进行比较来进行的：零分布是在零假设下导出的统计量值分布，该零假设是体重是从正态分布中抽取的。
对于 Jarque-Bera 检验，对于非常大的样本，零分布是具有两个自由度的卡方分布。
import matplotlib.pyplot as plt
dist = stats.chi2(df=2)
jb_val = np.linspace(0, 11, 100)
pdf = dist.pdf(jb_val)
fig, ax = plt.subplots(figsize=(8, 5))
def jb_plot(ax): # we'll reuse this
ax.plot(jb_val, pdf)
ax.set_title("Jarque-Bera Null Distribution")
ax.set_xlabel("statistic")
ax.set_ylabel("probability density")
jb_plot(ax)
plt.show();
通过 p 值量化比较：零分布中大于或等于统计量观测值的值的比例。
fig, ax = plt.subplots(figsize=(8, 5))
jb_plot(ax)
pvalue = dist.sf(res.statistic)
annotation = (f'p-value={pvalue:.6f}\n(shaded area)')
props = dict(facecolor='black', width=1, headwidth=5, headlength=8)
_ = ax.annotate(annotation, (7.5, 0.01), (8, 0.05), arrowprops=props)
i = jb_val >= res.statistic # indices of more extreme statistic values
ax.fill_between(jb_val[i], y1=0, y2=pdf[i])
ax.set_xlim(0, 11)
ax.set_ylim(0, 0.3)
plt.show()
res.pvalue
0.03045746622458189
如果 p 值“小” - 也就是说，如果从正态分布的总体中抽样产生如此极端的统计量值的概率很低 - 这可以作为反对零假设的证据，从而支持备择假设：体重不是从正态分布中抽取的。请注意
反之则不然；也就是说，该检验不用来提供支持零假设的证据。
被认为是“小”的值的阈值是一个选择，应该在分析数据之前做出[2]，同时考虑假阳性（错误地拒绝零假设）和假阴性（未能拒绝错误的零假设）的风险。
请注意，卡方分布提供了零分布的渐近近似；它仅对具有许多观测值的样本准确。对于像我们这样的小样本，scipy.stats.monte_carlo_test 可能会提供更准确（尽管是随机的）的精确 p 值近似。
def statistic(x, axis):
# underlying calculation of the Jarque Bera statistic
s = stats.skew(x, axis=axis)
k = stats.kurtosis(x, axis=axis)
return x.shape[axis]/6 * (s**2 + k**2/4)
res = stats.monte_carlo_test(x, stats.norm.rvs, statistic,
alternative='greater')
fig, ax = plt.subplots(figsize=(8, 5))
jb_plot(ax)
ax.hist(res.null_distribution, np.linspace(0, 10, 50),
density=True)
ax.legend(['asymptotic approximation (many observations)',
'Monte Carlo approximation (11 observations)'])
plt.show()
res.pvalue
0.0095
此外，尽管它们的随机性，以这种方式计算的 p 值可以用来精确控制错误拒绝零假设的速率[3]。