Статистические тесты позволяют принимать обоснованные решения на основе данных. В финансовом анализе, исследованиях и разработке моделей машинного обучения проверка гипотез определяет валидность предположений о данных, выявляет значимые различия между группами и подтверждает применимость математических методов.
Python предоставляет инструменты для реализации статистических тестов через библиотеки scipy, statsmodels и numpy. Правильное применение этих методов снижает риск ошибочных выводов и повышает надежность аналитических результатов.
Основы статистических гипотез
Статистическая гипотеза — утверждение о параметрах или свойствах генеральной совокупности, которое можно проверить на данных выборки. Тестирование гипотез формализует процесс принятия решений в условиях неопределенности.
- Нулевая гипотеза (H₀) представляет утверждение об отсутствии эффекта или различий.
- Альтернативная гипотеза (H₁) противоположна нулевой и отражает наличие эффекта. Например, при проверке нормальности распределения H₀ утверждает, что данные следуют нормальному распределению, а H₁ — что не следуют.
- P-value (вероятностное значение) показывает вероятность получить наблюдаемый или более экстремальный результат при условии истинности нулевой гипотезы. Значение p-value < 0.05 традиционно считается основанием для отклонения H₀. Уровень значимости α (обычно 0.05 или 0.01) определяет порог, ниже которого результат считается статистически значимым.
- Ошибка первого рода (Type I error) возникает при отклонении истинной нулевой гипотезы. Вероятность такой ошибки равна уровню значимости α.
- Ошибка второго рода (Type II error) происходит при принятии ложной нулевой гипотезы. Мощность теста (1 — β) определяет вероятность правильного отклонения ложной H₀.
Выбор уровня значимости зависит от задачи. В медицинских исследованиях используют α = 0.01 для минимизации ошибок первого рода. В разведочном анализе допустим α = 0.10. Баланс между ошибками первого и второго рода определяется контекстом применения.
Тесты нормальности распределения
Проверка нормальности распределения важна для применения параметрических статистических методов. Многие тесты и модели предполагают нормальность данных: t-тест, ANOVA, линейная регрессия. Нарушение этого предположения приводит к искаженным результатам.
Критерий Шапиро-Уилка
Тест Шапиро-Уилка оценивает, насколько выборка согласуется с нормальным распределением. Метод основан на корреляции между данными и соответствующими квантилями нормального распределения. Тест обладает высокой мощностью для выборок размером от 3 до 5000 наблюдений.
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
# Генерация финансовых временных рядов
np.random.seed(42)
n = 500
# Ряд 1: доходности с нормальным распределением
returns_normal = np.random.normal(0.0005, 0.02, n)
# Ряд 2: доходности с тяжелыми хвостами (t-распределение)
returns_heavy_tail = stats.t.rvs(df=3, loc=0.0005, scale=0.02, size=n)
# Ряд 3: цены с трендом (преобразуем доходности в цены)
prices = 100 * np.exp(np.cumsum(returns_normal))
# Тест Шапиро-Уилка для доходностей
stat_normal, p_normal = stats.shapiro(returns_normal)
stat_heavy, p_heavy = stats.shapiro(returns_heavy_tail)
stat_prices, p_prices = stats.shapiro(prices)
print("Тест Шапиро-Уилка:")
print(f"Нормальные доходности: статистика={stat_normal:.4f}, p-value={p_normal:.4f}")
print(f"Доходности с тяжелыми хвостами: статистика={stat_heavy:.4f}, p-value={p_heavy:.4f}")
print(f"Ценовой ряд: статистика={stat_prices:.4f}, p-value={p_prices:.4f}")
# Визуализация
fig, axes = plt.subplots(2, 3, figsize=(14, 8))
# Гистограммы
axes[0, 0].hist(returns_normal, bins=30, color='darkgray', alpha=0.7, edgecolor='black')
axes[0, 0].set_title('Нормальные доходности')
axes[0, 0].set_xlabel('Доходность')
axes[0, 0].set_ylabel('Частота')
axes[0, 1].hist(returns_heavy_tail, bins=30, color='darkgray', alpha=0.7, edgecolor='black')
axes[0, 1].set_title('Доходности с тяжелыми хвостами')
axes[0, 1].set_xlabel('Доходность')
axes[0, 2].hist(prices, bins=30, color='darkgray', alpha=0.7, edgecolor='black')
axes[0, 2].set_title('Ценовой ряд')
axes[0, 2].set_xlabel('Цена')
# Q-Q графики
stats.probplot(returns_normal, dist="norm", plot=axes[1, 0])
axes[1, 0].set_title(f'Q-Q график: нормальные (p={p_normal:.3f})')
axes[1, 0].get_lines()[0].set_color('black')
axes[1, 0].get_lines()[1].set_color('red')
stats.probplot(returns_heavy_tail, dist="norm", plot=axes[1, 1])
axes[1, 1].set_title(f'Q-Q график: тяжелые хвосты (p={p_heavy:.3f})')
axes[1, 1].get_lines()[0].set_color('black')
axes[1, 1].get_lines()[1].set_color('red')
stats.probplot(prices, dist="norm", plot=axes[1, 2])
axes[1, 2].set_title(f'Q-Q график: цены (p={p_prices:.3f})')
axes[1, 2].get_lines()[0].set_color('black')
axes[1, 2].get_lines()[1].set_color('red')
plt.tight_layout()
plt.show()Тест Шапиро-Уилка:
Нормальные доходности: статистика=0.9967, p-value=0.4013
Доходности с тяжелыми хвостами: статистика=0.9028, p-value=0.0000
Ценовой ряд: статистика=0.9327, p-value=0.0000
Рис. 1: Гистограммы и Q-Q графики для трех типов данных. Верхний ряд показывает распределения через гистограммы. Нижний ряд содержит Q-Q графики с p-value из теста Шапиро-Уилка. Нормальные доходности демонстрируют точки вдоль красной диагонали, тяжелые хвосты отклоняются на концах, ценовой ряд показывает систематическое искривление из-за тренда
Код генерирует три типа данных для демонстрации различных свойств распределений. Первый ряд представляет доходности с нормальным распределением, второй — доходности с тяжелыми хвостами через t-распределение с 3 степенями свободы, третий — ценовой ряд с трендом.
Тест Шапиро-Уилка возвращает тестовую статистику и p-value. Для нормально распределенных доходностей p-value превышает 0.05, что не дает оснований отклонить нулевую гипотезу о нормальности. Доходности с тяжелыми хвостами показывают низкий p-value (обычно < 0.01), указывая на отклонение от нормальности. Ценовой ряд также демонстрирует значимое отклонение из-за наличия тренда.
Q-Q графики визуализируют соответствие квантилей данных квантилям нормального распределения. Точки, лежащие на диагональной линии, указывают на нормальность. Отклонения в хвостах распределения (концах графика) свидетельствуют о тяжелых или легких хвостах относительно нормального распределения.
Критерий Харке-Бера
Тест Харке-Бера (Jarque-Bera) проверяет нормальность через асимметрию (skewness) и эксцесс (kurtosis) распределения. Метод особенно эффективен для больших выборок (n > 2000) и чувствителен к отклонениям в хвостах распределения.
Формула теста:
JB = (n/6) × (S² + (K — 3)²/4)
где:
- n — размер выборки;
- S — коэффициент асимметрии;
- K — коэффициент эксцесса.
Для нормального распределения S = 0 и K = 3. Статистика JB следует распределению χ² с двумя степенями свободы.
from scipy.stats import jarque_bera, skew, kurtosis
# Тест Харке-Бера
jb_normal, p_jb_normal = jarque_bera(returns_normal)
jb_heavy, p_jb_heavy = jarque_bera(returns_heavy_tail)
# Вычисление асимметрии и эксцесса
skew_normal = skew(returns_normal)
kurt_normal = kurtosis(returns_normal, fisher=True) # fisher=True дает excess kurtosis
skew_heavy = skew(returns_heavy_tail)
kurt_heavy = kurtosis(returns_heavy_tail, fisher=True)
print("\nТест Харке-Бера:")
print(f"Нормальные доходности: JB={jb_normal:.4f}, p-value={p_jb_normal:.4f}")
print(f" Асимметрия: {skew_normal:.4f}, Эксцесс: {kurt_normal:.4f}")
print(f"Доходности с тяжелыми хвостами: JB={jb_heavy:.4f}, p-value={p_jb_heavy:.4f}")
print(f" Асимметрия: {skew_heavy:.4f}, Эксцесс: {kurt_heavy:.4f}")Тест Харке-Бера:
Нормальные доходности: JB=4.0581, p-value=0.1315
Асимметрия: 0.1796, Эксцесс: 0.2564
Доходности с тяжелыми хвостами: JB=1146.3943, p-value=0.0000
Асимметрия: -0.5854, Эксцесс: 7.3250Тест Харке-Бера дополняет метод Шапиро-Уилка, предоставляя информацию о конкретных характеристиках распределения. Высокий эксцесс (> 3) указывает на тяжелые хвосты — частое свойство финансовых доходностей. Асимметрия показывает смещение распределения влево (отрицательные значения) или вправо (положительные значения).
Параметр fisher=True в функции kurtosis возвращает excess kurtosis (эксцесс минус 3), где значение 0 соответствует нормальному распределению. Положительный excess kurtosis означает более тяжелые хвосты по сравнению с нормальным распределением, что характерно для финансовых рядов.
Выбор между тестами Шапиро-Уилка и Харке-Бера зависит от размера выборки и цели анализа. Для малых выборок (n < 50) предпочтителен метод Шапиро-Уилка. Тест Харке-Бера эффективен на больших выборках и дает понимание механизма отклонения от нормальности через асимметрию и эксцесс.
Тесты на стационарность временных рядов
Стационарность временного ряда подразумевает постоянство статистических свойств (среднего, дисперсии, автокорреляции) во времени. Нестационарные ряды усложняют прогнозирование и моделирование, приводя к ложной корреляции между независимыми рядами.
Расширенный тест Дики-Фуллера (ADF)
Тест Дики-Фуллера проверяет наличие единичного корня в авторегрессионной модели. Единичный корень указывает на нестационарность ряда. Расширенная версия (Augmented Dickey-Fuller, ADF) учитывает более сложную структуру автокорреляции через добавление лагированных разностей.
ADF тестирует уравнение:
Δy_t = α + βt + γy_{t-1} + δ₁Δy_{t-1} + … + δ_pΔy_{t-p} + ε_t
где:
- Δy_t — первая разность ряда (y_t — y_{t-1});
- α — константа;
- βt — детерминированный тренд;
- γ — коэффициент при лагированном значении;
- δ_i — коэффициенты при лагированных разностях;
- ε_t — случайная ошибка.
Нулевая гипотеза: γ = 0 (наличие единичного корня, ряд нестационарен). Альтернативная гипотеза: γ < 0 (ряд стационарен).
from statsmodels.tsa.stattools import adfuller, kpss
# Генерация данных для тестов стационарности
np.random.seed(123)
n = 300
# Стационарный ряд (доходности)
stationary_series = np.random.normal(0, 1, n)
# Нестационарный ряд с трендом (цены)
trend = np.linspace(100, 150, n)
noise = np.random.normal(0, 2, n)
non_stationary_trend = trend + noise
# Нестационарный ряд со случайным блужданием
random_walk = 100 + np.cumsum(np.random.normal(0, 1, n))
# Тест ADF
def run_adf_test(series, name):
result = adfuller(series, autolag='AIC')
print(f"\n{name}:")
print(f" ADF статистика: {result[0]:.4f}")
print(f" P-value: {result[1]:.4f}")
print(f" Количество лагов: {result[2]}")
print(f" Критические значения:")
for key, value in result[4].items():
print(f" {key}: {value:.4f}")
return result[1]
p_stationary = run_adf_test(stationary_series, "Стационарный ряд")
p_trend = run_adf_test(non_stationary_trend, "Ряд с трендом")
p_random_walk = run_adf_test(random_walk, "Случайное блуждание")
# Визуализация
fig, axes = plt.subplots(3, 1, figsize=(12, 10))
axes[0].plot(stationary_series, color='black', linewidth=1)
axes[0].set_title(f'Стационарный ряд (ADF p-value: {p_stationary:.4f})')
axes[0].set_ylabel('Значение')
axes[0].axhline(y=0, color='red', linestyle='--', alpha=0.5)
axes[1].plot(non_stationary_trend, color='black', linewidth=1)
axes[1].set_title(f'Ряд с трендом (ADF p-value: {p_trend:.4f})')
axes[1].set_ylabel('Значение')
axes[2].plot(random_walk, color='black', linewidth=1)
axes[2].set_title(f'Случайное блуждание (ADF p-value: {p_random_walk:.4f})')
axes[2].set_xlabel('Время')
axes[2].set_ylabel('Значение')
plt.tight_layout()
plt.show()Стационарный ряд:
ADF статистика: -16.5125
P-value: 0.0000
Количество лагов: 0
Критические значения:
1%: -3.4524
5%: -2.8713
10%: -2.5719
Ряд с трендом:
ADF статистика: 0.3017
P-value: 0.9774
Количество лагов: 16
Критические значения:
1%: -3.4537
5%: -2.8718
10%: -2.5722
Случайное блуждание:
ADF статистика: -0.5244
P-value: 0.8872
Количество лагов: 0
Критические значения:
1%: -3.4524
5%: -2.8713
10%: -2.5719
Рис. 2: Три временных ряда с различной стационарностью. Стационарный ряд колеблется вокруг нуля без тренда. Ряд с трендом показывает систематический рост. Случайное блуждание демонстрирует накопление случайных изменений без возврата к среднему
Код генерирует три типа временных рядов с различными свойствами стационарности. Функция adfuller принимает параметр autolag=’AIC’, который автоматически выбирает оптимальное количество лагов через информационный критерий Акаике. Это снижает риск неправильной спецификации модели.
Стационарный ряд показывает низкую ADF статистику (сильно отрицательную) и p-value < 0.05, что позволяет отклонить нулевую гипотезу о наличии единичного корня. Ряд с детерминированным трендом и случайное блуждание демонстрируют высокий p-value (> 0.05), подтверждая нестационарность.
Критические значения на уровнях 1%, 5% и 10% позволяют оценить силу отклонения нулевой гипотезы. Если ADF статистика меньше критического значения, гипотеза о единичном корне отклоняется. Сравнение с несколькими уровнями значимости дает более полную картину.
Тест Квятковского-Филлипса-Шмидта-Шина (KPSS)
Тест KPSS использует противоположную логику относительно ADF. Нулевая гипотеза утверждает стационарность ряда, альтернативная — наличие единичного корня. Комбинация ADF и KPSS дает более надежную оценку стационарности.
# Тест KPSS
def run_kpss_test(series, name):
result = kpss(series, regression='c', nlags='auto')
print(f"\n{name}:")
print(f" KPSS статистика: {result[0]:.4f}")
print(f" P-value: {result[1]:.4f}")
print(f" Количество лагов: {result[2]}")
print(f" Критические значения:")
for key, value in result[3].items():
print(f" {key}: {value:.4f}")
return result[1]
print("\nТест KPSS:")
kpss_stationary = run_kpss_test(stationary_series, "Стационарный ряд")
kpss_trend = run_kpss_test(non_stationary_trend, "Ряд с трендом")
kpss_random_walk = run_kpss_test(random_walk, "Случайное блуждание")
# Сводная таблица результатов
print("\n" + "="*60)
print("Сводная таблица тестов стационарности:")
print("="*60)
print(f"{'Ряд':<25} {'ADF p-value':<15} {'KPSS p-value':<15} {'Вывод'}")
print("-"*60)
print(f"{'Стационарный':<25} {p_stationary:<15.4f} {kpss_stationary:<15.4f} {'Стационарен'}")
print(f"{'С трендом':<25} {p_trend:<15.4f} {kpss_trend:<15.4f} {'Нестационарен'}")
print(f"{'Случайное блуждание':<25} {p_random_walk:<15.4f} {kpss_random_walk:<15.4f} {'Нестационарен'}")Тест KPSS:
Стационарный ряд:
KPSS статистика: 0.1566
P-value: 0.1000
Количество лагов: 2
Критические значения:
10%: 0.3470
5%: 0.4630
2.5%: 0.5740
1%: 0.7390
Ряд с трендом:
KPSS статистика: 2.8245
P-value: 0.0100
Количество лагов: 10
Критические значения:
10%: 0.3470
5%: 0.4630
2.5%: 0.5740
1%: 0.7390
Случайное блуждание:
KPSS статистика: 1.6830
P-value: 0.0100
Количество лагов: 10
Критические значения:
10%: 0.3470
5%: 0.4630
2.5%: 0.5740
1%: 0.7390
============================================================
Сводная таблица тестов стационарности:
============================================================
Ряд ADF p-value KPSS p-value Вывод
------------------------------------------------------------
Стационарный 0.0000 0.1000 Стационарен
С трендом 0.9774 0.0100 Нестационарен
Случайное блуждание 0.8872 0.0100 НестационаренПараметр regression=’c’ в функции kpss указывает на модель с константой без тренда. Для рядов с явным трендом используют regression=’ct’ (константа и тренд). Параметр nlags=’auto’ автоматически определяет количество лагов по формуле Ньюи-Уэста.
Интерпретация комбинации тестов:
- ADF p-value < 0.05, KPSS p-value > 0.05: ряд стационарен;
- ADF p-value > 0.05, KPSS p-value < 0.05: ряд нестационарен;
- Оба p-value > 0.05: результаты неоднозначны, требуется дополнительный анализ;
- Оба p-value < 0.05: возможна разностно-стационарная модель (difference-stationary model).
Стационарный ряд проходит оба теста согласованно. Ряд с трендом и случайное блуждание показывают нестационарность в обоих тестах. Использование двух тестов с противоположными нулевыми гипотезами снижает риск ошибочных выводов.
Сравнение выборок
Сравнение двух или более выборок определяет, принадлежат ли они одной генеральной совокупности или имеют значимые различия. Задача возникает при A/B тестировании, сравнении результатов экспериментов, оценке эффективности различных методов.
T-тест для независимых выборок
T-тест сравнивает средние значения двух выборок при предположении нормальности распределений. Метод подходит для выборок с приблизительно равными дисперсиями (гомоскедастичность). Для разных дисперсий используют Т-тест Уэлча.
from scipy.stats import ttest_ind, mannwhitneyu, levene
# Генерация данных для сравнения
np.random.seed(456)
n_samples = 100
# Выборка A: доходности стратегии A
strategy_a = np.random.normal(0.001, 0.015, n_samples)
# Выборка B: доходности стратегии B (немного выше среднего)
strategy_b = np.random.normal(0.0015, 0.015, n_samples)
# Выборка C: стратегия с другой дисперсией
strategy_c = np.random.normal(0.001, 0.025, n_samples)
# Тест Левене для проверки равенства дисперсий
levene_ab = levene(strategy_a, strategy_b)
levene_ac = levene(strategy_a, strategy_c)
print("Тест Левене (равенство дисперсий):")
print(f"Стратегии A vs B: статистика={levene_ab.statistic:.4f}, p-value={levene_ab.pvalue:.4f}")
print(f"Стратегии A vs C: статистика={levene_ac.statistic:.4f}, p-value={levene_ac.pvalue:.4f}")
# T-тест для независимых выборок
# equal_var=True для обычного t-теста, False для Welch's t-test
t_stat_ab, p_value_ab = ttest_ind(strategy_a, strategy_b, equal_var=True)
t_stat_ac, p_value_ac = ttest_ind(strategy_a, strategy_c, equal_var=False)
print("\nT-тест для независимых выборок:")
print(f"Стратегии A vs B: t-статистика={t_stat_ab:.4f}, p-value={p_value_ab:.4f}")
print(f"Стратегии A vs C (Welch): t-статистика={t_stat_ac:.4f}, p-value={p_value_ac:.4f}")
# Описательная статистика
print("\nОписательная статистика:")
print(f"Стратегия A: среднее={np.mean(strategy_a):.6f}, std={np.std(strategy_a, ddof=1):.6f}")
print(f"Стратегия B: среднее={np.mean(strategy_b):.6f}, std={np.std(strategy_b, ddof=1):.6f}")
print(f"Стратегия C: среднее={np.mean(strategy_c):.6f}, std={np.std(strategy_c, ddof=1):.6f}")
# Визуализация
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# Boxplot
axes[0].boxplot([strategy_a, strategy_b, strategy_c], labels=['Стратегия A', 'Стратегия B', 'Стратегия C'])
axes[0].set_ylabel('Доходность')
axes[0].set_title('Распределение доходностей стратегий')
axes[0].grid(True, alpha=0.3)
# Гистограммы
axes[1].hist(strategy_a, bins=20, alpha=0.5, label='Стратегия A', color='gray', edgecolor='black')
axes[1].hist(strategy_b, bins=20, alpha=0.5, label='Стратегия B', color='blue', edgecolor='black')
axes[1].hist(strategy_c, bins=20, alpha=0.5, label='Стратегия C', color='red', edgecolor='black')
axes[1].set_xlabel('Доходность')
axes[1].set_ylabel('Частота')
axes[1].set_title('Наложенные распределения')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()Тест Левене (равенство дисперсий):
Стратегии A vs B: статистика=2.6367, p-value=0.1060
Стратегии A vs C: статистика=32.5772, p-value=0.0000
T-тест для независимых выборок:
Стратегии A vs B: t-статистика=1.5353, p-value=0.1263
Стратегии A vs C (Welch): t-статистика=-0.4556, p-value=0.6493
Описательная статистика:
Стратегия A: среднее=0.003171, std=0.013876
Стратегия B: среднее=-0.000132, std=0.016441
Стратегия C: среднее=0.004468, std=0.024867
Рис. 3: Сравнение трех стратегий. Боксплот слева показывает медианы, квартили и выбросы. Наложенные гистограммы справа демонстрируют формы распределений. Стратегия C имеет более широкое распределение при схожем среднем со стратегией A
Тест Левене проверяет гомоскедастичность — равенство дисперсий между выборками. Результат определяет выбор между стандартным t-тестом (equal_var=True) и t-тестом Уэлча (equal_var=False). Метод Уэлча не требует равенства дисперсий и более устойчив в случае их нарушения.
Параметр ddof=1 в функции np.std обеспечивает несмещенную оценку стандартного отклонения через деление на (n-1) вместо n. Это корректирует систематическое занижение дисперсии в выборках.
P-value < 0.05 указывает на статистически значимое различие средних. Стратегии A и B могут не показать значимых различий несмотря на разные средние из-за высокой дисперсии доходностей. Стратегия C с большей дисперсией требует проведения t-теста Уэлча для корректного сравнения.
Тест Манна-Уитни
Тест Манна-Уитни (Mann-Whitney U test) — непараметрическая альтернатива t-тесту, не требующая нормальности распределений. Метод сравнивает ранги наблюдений вместо их значений, что делает его устойчивым к выбросам и применимым для ординальных данных.
# Генерация данных с ненормальным распределением
np.random.seed(789)
# Экспоненциальное распределение (сильно асимметричное)
sample_exp1 = np.random.exponential(scale=2.0, size=80)
sample_exp2 = np.random.exponential(scale=2.5, size=80)
# Тест Манна-Уитни
u_stat, p_value_mw = mannwhitneyu(sample_exp1, sample_exp2, alternative='two-sided')
print("\nТест Манна-Уитни (непараметрический):")
print(f"U-статистика: {u_stat:.4f}")
print(f"P-value: {p_value_mw:.4f}")
# Сравнение с t-тестом на тех же данных
t_stat_exp, p_value_t = ttest_ind(sample_exp1, sample_exp2)
print("\nДля сравнения, t-тест на тех же данных:")
print(f"T-статистика: {t_stat_exp:.4f}")
print(f"P-value: {p_value_t:.4f}")
print("\nОписательная статистика экспоненциальных выборок:")
print(f"Выборка 1: среднее={np.mean(sample_exp1):.4f}, медиана={np.median(sample_exp1):.4f}")
print(f"Выборка 2: среднее={np.mean(sample_exp2):.4f}, медиана={np.median(sample_exp2):.4f}")Тест Манна-Уитни (непараметрический):
U-статистика: 2902.0000
P-value: 0.3100
Для сравнения, t-тест на тех же данных:
T-статистика: -1.5443
P-value: 0.1245
Описательная статистика экспоненциальных выборок:
Выборка 1: среднее=1.8784, медиана=1.5342
Выборка 2: среднее=2.3214, медиана=1.7470Параметр alternative в mannwhitneyu определяет тип альтернативной гипотезы: ‘two-sided’ для двустороннего теста (распределения различаются), ‘less’ или ‘greater’ для односторонних тестов (одна выборка систематически больше другой).
Тест Манна-Уитни основан на ранжировании объединенных данных обеих выборок. U-статистика подсчитывает количество пар, где значение из первой выборки меньше значения из второй. При равных распределениях ожидается U ≈ n₁ × n₂ / 2.
Для экспоненциально распределенных данных тест Манна-Уитни дает более надежные результаты, чем t-тест. T-тест предполагает нормальность и чувствителен к асимметрии распределения. Непараметрический подход избегает этих ограничений, сохраняя статистическую мощность.
Выбор между t-тестом и Манна-Уитни зависит от свойств данных. При нормальности распределений t-тест обладает большей мощностью. Для ненормальных данных, выбросов или малых выборок предпочтителен Манна-Уитни. Проверка нормальности через тесты Шапиро-Уилка или Харке-Бера помогает принять обоснованное решение.
Тесты на автокорреляцию
Автокорреляция временного ряда показывает связь между текущими и предыдущими значениями. Наличие автокорреляции нарушает предположение о независимости наблюдений в регрессионных моделях и влияет на точность стандартных ошибок коэффициентов.
Критерий Льюнга-Бокса
Тест Льюнга-Бокса проверяет, значимо ли отличаются от нуля первые k автокорреляций. Метод применяется для диагностики остатков моделей временных рядов и проверки гипотезы белого шума.
Статистика теста рассчитывается по формуле:
Q = n(n + 2) × Σ(ρ²ₖ / (n — k))
где:
- n — размер выборки;
- ρₖ — автокорреляция на лаге k;
- суммирование по k от 1 до m.
Статистика Q следует распределению χ² с m степенями свободы при нулевой гипотезе об отсутствии автокорреляции.
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
# Генерация рядов с различной автокорреляцией
np.random.seed(101)
n = 200
# Белый шум (нет автокорреляции)
white_noise = np.random.normal(0, 1, n)
# AR(1) процесс с автокорреляцией
ar1_series = np.zeros(n)
ar1_series[0] = np.random.normal(0, 1)
phi = 0.7 # коэффициент автокорреляции
for t in range(1, n):
ar1_series[t] = phi * ar1_series[t-1] + np.random.normal(0, 1)
# MA(1) процесс
ma1_series = np.zeros(n)
errors = np.random.normal(0, 1, n)
theta = 0.6
for t in range(1, n):
ma1_series[t] = errors[t] + theta * errors[t-1]
# Тест Льюнга-Бокса
def ljung_box_test(series, name, lags=10):
result = acorr_ljungbox(series, lags=lags, return_df=True)
print(f"\n{name} (первые {lags} лагов):")
print(result[['lb_stat', 'lb_pvalue']].head(lags))
return result
lb_white = ljung_box_test(white_noise, "Белый шум")
lb_ar1 = ljung_box_test(ar1_series, "AR(1) процесс")
lb_ma1 = ljung_box_test(ma1_series, "MA(1) процесс")
# Визуализация автокорреляционных функций
fig, axes = plt.subplots(3, 3, figsize=(14, 10))
# Временные ряды
axes[0, 0].plot(white_noise, color='black', linewidth=0.8)
axes[0, 0].set_title('Белый шум')
axes[0, 0].set_ylabel('Значение')
axes[1, 0].plot(ar1_series, color='black', linewidth=0.8)
axes[1, 0].set_title('AR(1) процесс')
axes[1, 0].set_ylabel('Значение')
axes[2, 0].plot(ma1_series, color='black', linewidth=0.8)
axes[2, 0].set_title('MA(1) процесс')
axes[2, 0].set_xlabel('Время')
axes[2, 0].set_ylabel('Значение')
# ACF графики
plot_acf(white_noise, lags=20, ax=axes[0, 1], color='black')
axes[0, 1].set_title('ACF: Белый шум')
plot_acf(ar1_series, lags=20, ax=axes[1, 1], color='black')
axes[1, 1].set_title('ACF: AR(1)')
plot_acf(ma1_series, lags=20, ax=axes[2, 1], color='black')
axes[2, 1].set_title('ACF: MA(1)')
# PACF графики
plot_pacf(white_noise, lags=20, ax=axes[0, 2], color='black', method='ywm')
axes[0, 2].set_title('PACF: Белый шум')
plot_pacf(ar1_series, lags=20, ax=axes[1, 2], color='black', method='ywm')
axes[1, 2].set_title('PACF: AR(1)')
plot_pacf(ma1_series, lags=20, ax=axes[2, 2], color='black', method='ywm')
axes[2, 2].set_title('PACF: MA(1)')
plt.tight_layout()
plt.show()Белый шум (первые 10 лагов):
lb_stat lb_pvalue
1 1.447212 0.228976
2 1.783037 0.410033
3 2.578395 0.461290
4 3.271393 0.513478
5 6.742075 0.240540
6 6.861447 0.333851
7 6.905907 0.438742
8 11.719403 0.164169
9 14.848603 0.095179
10 14.946305 0.134032
AR(1) процесс (первые 10 лагов):
lb_stat lb_pvalue
1 89.998819 2.383022e-21
2 118.799245 1.596143e-26
3 130.142268 5.039811e-28
4 132.752973 1.003614e-27
5 133.113759 5.195390e-27
6 133.195003 2.728971e-26
7 133.256440 1.311156e-25
8 133.586675 5.101737e-25
9 136.281179 6.035908e-25
10 138.158927 1.004298e-24
MA(1) процесс (первые 10 лагов):
lb_stat lb_pvalue
1 41.584263 1.128984e-10
2 41.631316 9.117477e-10
3 44.641782 1.102534e-09
4 47.149906 1.419125e-09
5 49.365843 1.868175e-09
6 49.668824 5.477335e-09
7 51.036384 9.037416e-09
8 54.476757 5.579043e-09
9 56.507844 6.277245e-09
10 56.820572 1.440450e-08
Рис. 4: Три типа временных рядов с ACF и PACF графиками. Левый столбец показывает исходные ряды. Средний столбец содержит ACF с доверительными интервалами (синие области). Правый столбец — PACF графики. Белый шум имеет незначимые автокорреляции. AR(1) показывает экспоненциальное затухание ACF и один значимый лаг PACF. MA(1) демонстрирует один значимый лаг ACF и затухание PACF
Функция acorr_ljungbox вычисляет статистику Льюнга-Бокса для указанного количества лагов. Параметр return_df=True возвращает результаты в формате датафрейма с колонками lb_stat (статистика теста) и lb_pvalue (p-value).
Белый шум показывает высокие p-value для всех лагов, подтверждая отсутствие автокорреляции. AR(1) процесс демонстрирует значимую автокорреляцию на первых лагах с низкими p-value. MA(1) процесс показывает автокорреляцию только на первом лаге, что соответствует теоретическим свойствам модели скользящего среднего.
ACF (Autocorrelation Function) отображает корреляцию ряда с его лагированными значениями. PACF (Partial Autocorrelation Function) показывает корреляцию после удаления влияния промежуточных лагов. AR процессы имеют экспоненциально затухающую ACF и резкий обрыв PACF на лаге p. MA процессы демонстрируют обратную картину: обрыв ACF на лаге q и затухающую PACF.
Тест Дарбина-Уотсона
Тест Дарбина-Уотсона специально разработан для обнаружения автокорреляции первого порядка в остатках регрессионной модели. Метод широко используется в эконометрике для диагностики линейных регрессий.
Статистика теста рассчитывается по формуле:
DW = Σ(eₜ — eₜ₋₁)² / Σeₜ²
где eₜ — остатки модели в момент времени t.
Значение DW находится в диапазоне [0, 4]. DW ≈ 2 указывает на отсутствие автокорреляции. DW < 2 свидетельствует о положительной автокорреляции, DW > 2 — об отрицательной.
from statsmodels.regression.linear_model import OLS
from statsmodels.tools import add_constant
from statsmodels.stats.stattools import durbin_watson
# Генерация данных для регрессии
np.random.seed(202)
n = 150
# Независимая переменная
x = np.linspace(0, 10, n)
# Зависимая переменная с автокоррелированными остатками
true_beta = 2.5
true_alpha = 10
# Модель 1: остатки без автокорреляции
errors_no_ac = np.random.normal(0, 2, n)
y_no_ac = true_alpha + true_beta * x + errors_no_ac
# Модель 2: остатки с положительной автокорреляцией
errors_ac = np.zeros(n)
errors_ac[0] = np.random.normal(0, 2)
rho = 0.8 # коэффициент автокорреляции остатков
for t in range(1, n):
errors_ac[t] = rho * errors_ac[t-1] + np.random.normal(0, 2)
y_with_ac = true_alpha + true_beta * x + errors_ac
# Регрессионные модели
X = add_constant(x)
model_no_ac = OLS(y_no_ac, X).fit()
model_with_ac = OLS(y_with_ac, X).fit()
# Тест Дарбина-Уотсона
dw_no_ac = durbin_watson(model_no_ac.resid)
dw_with_ac = durbin_watson(model_with_ac.resid)
print("Тест Дарбина-Уотсона:")
print(f"Модель без автокорреляции: DW = {dw_no_ac:.4f}")
print(f"Модель с автокорреляцией: DW = {dw_with_ac:.4f}")
print("\nИнтерпретация DW статистики:")
print("DW ≈ 2.0: нет автокорреляции")
print("DW < 1.5: положительная автокорреляция") print("DW > 2.5: отрицательная автокорреляция")
# Визуализация остатков
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Регрессии
axes[0, 0].scatter(x, y_no_ac, alpha=0.5, color='gray', s=20)
axes[0, 0].plot(x, model_no_ac.predict(X), color='red', linewidth=2)
axes[0, 0].set_title(f'Регрессия без автокорреляции (DW={dw_no_ac:.2f})')
axes[0, 0].set_xlabel('X')
axes[0, 0].set_ylabel('Y')
axes[0, 1].scatter(x, y_with_ac, alpha=0.5, color='gray', s=20)
axes[0, 1].plot(x, model_with_ac.predict(X), color='red', linewidth=2)
axes[0, 1].set_title(f'Регрессия с автокорреляцией (DW={dw_with_ac:.2f})')
axes[0, 1].set_xlabel('X')
axes[0, 1].set_ylabel('Y')
# Остатки во времени
axes[1, 0].plot(model_no_ac.resid, color='black', linewidth=1)
axes[1, 0].axhline(y=0, color='red', linestyle='--', alpha=0.5)
axes[1, 0].set_title('Остатки без автокорреляции')
axes[1, 0].set_xlabel('Наблюдение')
axes[1, 0].set_ylabel('Остаток')
axes[1, 1].plot(model_with_ac.resid, color='black', linewidth=1)
axes[1, 1].axhline(y=0, color='red', linestyle='--', alpha=0.5)
axes[1, 1].set_title('Остатки с автокорреляцией')
axes[1, 1].set_xlabel('Наблюдение')
axes[1, 1].set_ylabel('Остаток')
plt.tight_layout()
plt.show()Тест Дарбина-Уотсона:
Модель без автокорреляции: DW = 1.8046
Модель с автокорреляцией: DW = 0.4821
Интерпретация DW статистики:
DW ≈ 2.0: нет автокорреляции
DW < 1.5: положительная автокорреляция DW > 2.5: отрицательная автокорреляция
Рис. 5: Сравнение регрессий с различными свойствами остатков. Верхний ряд показывает диаграммы рассеяния с линиями регрессии. Нижний ряд отображает остатки во времени. Модель без автокорреляции имеет хаотично распределенные остатки. Модель с автокорреляцией демонстрирует плавные волны остатков, подтверждая зависимость между соседними значениями
Функция add_constant добавляет столбец единиц к матрице регрессоров для оценки свободного члена. Класс OLS реализует метод наименьших квадратов, метод fit() возвращает объект с результатами оценки модели.
Модель без автокорреляции показывает DW ≈ 2.0, что соответствует независимым остаткам. Модель с автокоррелированными остатками демонстрирует DW < 1.5, указывая на положительную автокорреляцию. Визуальный паттерн остатков во времени подтверждает результаты: независимые остатки хаотично колеблются вокруг нуля, автокоррелированные показывают плавные волны с периодами устойчивых положительных и отрицательных значений.
Автокорреляция остатков нарушает эффективность оценок метода наименьших квадратов. Стандартные ошибки коэффициентов становятся заниженными, что приводит к завышенной статистической значимости предикторов. Обнаружение автокорреляции требует применения альтернативных методов оценки: обобщенный метод наименьших квадратов (GLS), модели ARIMA для остатков, робастные стандартные ошибки Ньюи-Уэста.
Множественное тестирование и корректировка p-value
При одновременном проведении множества статистических тестов возрастает вероятность ошибки первого рода. Если проводить 20 независимых тестов на уровне значимости α = 0.05, ожидаемое количество ложноположительных результатов составит 1 тест. Для 100 тестов — уже 5 ложных открытий.
Проблема множественного тестирования чаще всего возникает в разведочном анализе данных и отборе признаков для машинного обучения. Без корректировки p-value результаты теряют статистическую валидность.
Метод Бонферрони
Корректировка Бонферрони — простейший метод контроля групповой вероятности ошибки (Family-wise error rate, FWER). FWER определяет вероятность сделать хотя бы одну ошибку первого рода среди всех тестов. Метод Бонферрони делит уровень значимости на количество тестов: α_adj = α / m.
from statsmodels.stats.multitest import multipletests
# Генерация данных для множественного тестирования
np.random.seed(303)
n_tests = 50
n_samples = 100
# Создание p-values для разных сценариев
# 45 тестов без эффекта (нулевая гипотеза верна)
p_values_null = []
for i in range(45):
sample1 = np.random.normal(0, 1, n_samples)
sample2 = np.random.normal(0, 1, n_samples)
_, p = ttest_ind(sample1, sample2)
p_values_null.append(p)
# 5 тестов с реальным эффектом
p_values_effect = []
for i in range(5):
sample1 = np.random.normal(0, 1, n_samples)
sample2 = np.random.normal(0.5, 1, n_samples) # сдвиг среднего
_, p = ttest_ind(sample1, sample2)
p_values_effect.append(p)
# Объединение всех p-values
p_values = np.array(p_values_null + p_values_effect)
# Без корректировки
significant_uncorrected = np.sum(p_values < 0.05)
# Корректировка Бонферрони
reject_bonf, pvals_bonf, _, _ = multipletests(p_values, alpha=0.05, method='bonferroni')
significant_bonf = np.sum(reject_bonf)
# Корректировка FDR (Benjamini-Hochberg)
reject_fdr, pvals_fdr, _, _ = multipletests(p_values, alpha=0.05, method='fdr_bh')
significant_fdr = np.sum(reject_fdr)
print("Результаты множественного тестирования:")
print(f"Всего тестов: {n_tests}")
print(f"Тесты с реальным эффектом: 5")
print(f"Тесты без эффекта: 45")
print(f"\nБез корректировки (α=0.05):")
print(f" Значимых результатов: {significant_uncorrected}")
print(f"\nМетод Бонферрони:")
print(f" Скорректированный α: {0.05/n_tests:.4f}")
print(f" Значимых результатов: {significant_bonf}")
print(f"\nМетод FDR (Benjamini-Hochberg):")
print(f" Значимых результатов: {significant_fdr}")
# Визуализация p-values
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# Сортировка p-values для визуализации
sorted_idx = np.argsort(p_values)
sorted_p = p_values[sorted_idx]
# График 1: Все p-values без корректировки
axes[0].scatter(range(n_tests), sorted_p, color='black', s=30)
axes[0].axhline(y=0.05, color='red', linestyle='--', label='α=0.05')
axes[0].set_xlabel('Ранг теста')
axes[0].set_ylabel('P-value')
axes[0].set_title(f'Без корректировки ({significant_uncorrected} значимых)')
axes[0].legend()
axes[0].set_yscale('log')
# График 2: Корректировка Бонферрони
axes[1].scatter(range(n_tests), sorted_p, color='black', s=30)
axes[1].axhline(y=0.05/n_tests, color='red', linestyle='--',
label=f'α_Bonf={0.05/n_tests:.4f}')
axes[1].set_xlabel('Ранг теста')
axes[1].set_ylabel('P-value')
axes[1].set_title(f'Бонферрони ({significant_bonf} значимых)')
axes[1].legend()
axes[1].set_yscale('log')
# График 3: FDR
axes[2].scatter(range(n_tests), sorted_p, color='black', s=30)
# Линия FDR
fdr_line = 0.05 * np.arange(1, n_tests + 1) / n_tests
axes[2].plot(range(n_tests), fdr_line, color='red', linestyle='--',
label='FDR порог')
axes[2].set_xlabel('Ранг теста')
axes[2].set_ylabel('P-value')
axes[2].set_title(f'FDR ({significant_fdr} значимых)')
axes[2].legend()
axes[2].set_yscale('log')
plt.tight_layout()
plt.show()Результаты множественного тестирования:
Всего тестов: 50
Тесты с реальным эффектом: 5
Тесты без эффекта: 45
Без корректировки (α=0.05):
Значимых результатов: 7
Метод Бонферрони:
Скорректированный α: 0.0010
Значимых результатов: 4
Метод FDR (Benjamini-Hochberg):
Значимых результатов: 4
Рис. 6: Сравнение методов корректировки множественного тестирования. Все графики используют логарифмическую шкалу для p-значений. Слева — необработанные p-values с порогом α = 0.05. В центре — жесткий порог Бонферрони. Справа — адаптивный порог FDR, растущий с рангом теста. FDR находит больше эффектов, чем Бонферрони, контролируя долю ложных открытий вместо вероятности хотя бы одной ошибки
Код выше демонстрирует различия методов на логарифмической шкале. FDR использует адаптивный порог, который растет с рангом теста: pₖ ≤ (k/m) × α. Это позволяет обнаружить больше эффектов по сравнению с Бонферрони, сохраняя контроль над долей ложных открытий.
Без корректировки обнаруживается около 5-7 значимых результатов: 5 истинно положительных (реальный эффект) плюс 2-3 ложноположительных из 45 тестов без эффекта при α = 0.05. Метод Бонферрони снижает порог до α / 50 = 0.001, что резко уменьшает количество обнаруженных эффектов. Метод находит только самые сильные эффекты, минимизируя ложноположительные результаты, но увеличивая ошибки второго рода.
Бонферрони излишне консервативен при большом количестве тестов. Для 1000 тестов скорректированный уровень значимости составляет 0.00005, что практически исключает обнаружение умеренных эффектов. Метод подходит для критических применений, где ошибки первого рода недопустимы.
False Discovery Rate (FDR)
FDR контролирует ожидаемую долю ложных открытий среди всех отклоненных нулевых гипотез. В отличие от FWER, который контролирует вероятность хотя бы одной ошибки, FDR допускает некоторое количество ложноположительных результатов, повышая мощность теста.
Метод Benjamini-Hochberg реализует FDR через следующую процедуру:
- Отсортировать p-значения по возрастанию: p₁ ≤ p₂ ≤ … ≤ pₘ;
- Найти максимальный индекс k, для которого pₖ ≤ (k/m) × α;
- Отклонить нулевые гипотезы для всех тестов с индексами от 1 до k.
FDR = 0.05 означает, что среди всех отклоненных гипотез ожидается не более 5% ложных открытий. Если обнаружено 20 значимых результатов при FDR = 0.05, ожидаемое количество ложноположительных составляет 1 результат.
Выбор между Бонферрони и FDR определяется контекстом:
- Бонферрони: медицинские исследования, регуляторные решения, ситуации с высокой ценой ошибок первого рода;
- FDR: разведочный анализ, скрининг переменных для моделей, геномика, ситуации с приемлемой долей ложных открытий.
Для разведочного анализа с последующей валидацией FDR предпочтительнее. Метод обеспечивает баланс между обнаружением эффектов и контролем ошибок. Для подтверждающих исследований используют Бонферрони или более мощные методы контроля FWER (Holm, Hochberg).
Заключение
Статистические тесты превращают данные в обоснованные выводы, снижая субъективность аналитических решений. Проверка нормальности через Шапиро-Уилка или Харке-Бера определяет применимость параметрических методов. Тесты стационарности ADF и KPSS выявляют структурные свойства временных рядов, критичные для моделирования. Сравнение выборок через t-тесты или Манна-Уитни количественно оценивает различия между группами. Диагностика автокорреляции обеспечивает валидность регрессионных моделей.
Python через scipy и statsmodels предоставляет инструменты для всех перечисленных задач с минимальным кодом. Корректное применение тестов требует понимания их предположений и ограничений. Множественное тестирование усложняет интерпретацию, но методы Бонферрони и FDR контролируют ошибки. Комбинация нескольких тестов дает более надежную картину: ADF с KPSS для стационарности, Шапиро-Уилка с Q-Q графиками для нормальности. Статистическая строгость на этапе анализа данных закладывает фундамент для качественных моделей и обоснованных бизнес-решений.