Критерий согласия Харке-Бера#
Предположим, мы хотим сделать вывод из измерений, отличаются ли веса взрослых мужчин
в медицинском исследовании от нормального распределения [1]. Веса (фунты)
записываются в массив x ниже.
import numpy as np
x = np.array([148, 154, 158, 160, 161, 162, 166, 170, 182, 195, 236])
Тест Жарке-Бера scipy.stats.jarque_bera начинается с вычисления статистики на основе выборочной асимметрии и эксцесса.
from scipy import stats
res = stats.jarque_bera(x)
res.statistic
np.float64(6.982848237344646)
Поскольку нормальное распределение имеет нулевую асимметрию и нулевой («избыточный» или «Фишеровский») эксцесс, значение этой статистики обычно низкое для выборок, взятых из нормального распределения.
Тест выполняется путем сравнения наблюдаемого значения статистики с нулевым распределением: распределением значений статистики, полученным при нулевой гипотезе о том, что веса были взяты из нормального распределения.
Для теста Жарке-Бера нулевое распределение для очень больших выборок — это
chi-squared distribution с двумя степенями свободы.
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
np.float64(0.03045746622458189)
Если p-значение «маленькое» - то есть если существует низкая вероятность выборки данных из нормально распределенной совокупности, которая дает такое экстремальное значение статистики - это может быть принято как свидетельство против нулевой гипотезы в пользу альтернативы: веса не были взяты из нормального распределения. Обратите внимание, что:
Обратное неверно; то есть тест не используется для подтверждения нулевой гипотезы.
Порог для значений, которые будут считаться «малыми», — это выбор, который следует сделать до анализа данных [2] с учётом рисков как ложных с positives (ошибочное отклонение нулевой гипотезы), так и ложных negatives (неспособность отклонить ложную нулевую гипотезу).
Обратите внимание, что распределение хи-квадрат даёт асимптотическое приближение
нулевого распределения; оно точно только для выборок с большим
количеством наблюдений. Для малых выборок, как наша, 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
np.float64(0.0068)
Кроме того, несмотря на их стохастическую природу, p-значения, вычисленные таким образом, могут использоваться для точного контроля уровня ложных отклонений нулевой гипотезы [3].