Seminar 9: Проверка статистических гипотез

Author

Carlos Buitrago

Show Python code

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.special import gammaln
rng = np.random.default_rng(2026)
plt.rcParams.update({'figure.figsize': (9, 4.8), 'axes.grid': True, 'grid.alpha': .22, 'font.size': 11})

1 1. Постановка задачи

Пусть $\mathcal P$ — семейство распределений наблюдения $X$. Основная гипотеза $H_0:P\in\mathcal P_0$ проверяется против альтернативы $H_1:P\in\mathcal P_1$, где множества распределений не пересекаются. Для параметрического семейства $\{P_\theta:\theta\in\Theta\}$ записываем $H_0:\theta\in\Theta_0$, $H_1:\theta\in\Theta_1$.

Критерий задаётся критической областью $S$ в выборочном пространстве: если $X\in S$, отвергаем $H_0$; иначе не отвергаем. Неотвержение не доказывает истинность $H_0$.

Две ошибки

Ошибка I рода: отвергнуть верную $H_0$.
Ошибка II рода: не отвергнуть $H_0$, когда верна $H_1$.
При фиксированном ограничении ошибки I рода стараемся максимизировать вероятность обнаружения альтернативы.

2 2. Функция мощности, размер и уровень

Для критической области $S$ функция мощности определяется формулой \[\beta(\theta,S)=P_\theta(X\in S).\] При $\theta\in\Theta_0$ это вероятность ошибки I рода, при $\theta\in\Theta_1$ — мощность. Вероятность ошибки II рода при альтернативе $\theta$ равна $1-\beta(\theta,S)$.

Размер критерия: $\alpha_0(S)=\sup_{\theta\in\Theta_0}\beta(\theta,S)$. Критерий имеет уровень значимости $\alpha$, если $\alpha_0(S)\le\alpha$ (равенство не обязательно).

Равномерно наиболее мощный (р.н.м., UMP) критерий уровня $\alpha$: для любого другого критерия $R$ уровня $\alpha$ и каждого $\theta\in\Theta_1$ выполнено $\beta(\theta,S)\ge\beta(\theta,R)$.

Стратегия решения

Выписать $H_0,H_1$ и плотность/вероятности выборки.
Найти статистику, монотонную по отношению правдоподобия.
Определить направление критической области.
Выбрать порог из распределения при $H_0$.
Вычислить мощность при произвольном параметре альтернативы.
Обосновать оптимальность леммой Неймана—Пирсона или теоремой о монотонном отношении правдоподобия.

2.1 2.1. Почему уровень и мощность — разные вещи

Для нормальной модели $X_i\sim N(\mu,1)$ проверяем $H_0:\mu\le0$ против $H_1:\mu>0$. Правосторонний критерий имеет вид $\sqrt n\,\overline X>z_{1-\alpha}$, где $z_q=\Phi^{-1}(q)$. Его функция мощности: \[\beta(\mu)=1-\Phi(z_{1-\alpha}-\sqrt n\,\mu).\] При $\mu=0$ она равна $\alpha$; при $\mu>0$ растёт к единице.

Show Python code

alpha=.05
mu=np.linspace(-.8,1.3,400)
fig,ax=plt.subplots()
for nn in [5,20,80]:
    ax.plot(mu, stats.norm.sf(stats.norm.isf(alpha)-np.sqrt(nn)*mu),label=f'n = {nn}',lw=2)
ax.axvline(0,color='black',ls='--',label='граница H₀')
ax.axhline(alpha,color='gray',ls=':',label='α = 0.05')
ax.set(xlabel='Истинное μ',ylabel='Вероятность отклонить H₀',title='Функция мощности и влияние размера выборки',ylim=(0,1.02))
ax.legend();plt.show()

3 3. Лемма Неймана—Пирсона

Лемма (простая гипотеза против простой альтернативы)

Пусть $H_0:\theta=\theta_0$, $H_1:\theta=\theta_1$ и существуют плотности $p_0,p_1$ относительно общей меры. Критерий, отвергающий $H_0$ при большом отношении $p_1(x)/p_0(x)$, является наиболее мощным среди критериев данного уровня $\alpha$. В общем случае на границе $p_1=cp_0$ допускается рандомизация, чтобы получить точный уровень.

\[\varphi(x)=\begin{cases}1,&p_1(x)>cp_0(x),\\ \gamma,&p_1(x)=cp_0(x),\\0,&p_1(x)<cp_0(x).\end{cases}\]

Идея доказательства. Для произвольного конкурирующего критерия $\psi$ уровня не выше $\alpha$ имеем $(\varphi-\psi)(p_1-cp_0)\ge0$ поточечно. Интегрируя, получаем \[E_1\varphi-E_1\psi\ge c(E_0\varphi-E_0\psi)\ge0.\] Последнее неравенство использует $E_0\varphi=\alpha$ (при точной калибровке). Для нерандомизированных дискретных критериев может получиться уровень строго ниже $\alpha$.

Важная оговорка: если $p_0(x)=0$ и $p_1(x)>0$, такие точки всегда выгодно включить в критическую область: это даёт мощность без ошибки I рода.

3.1 3.1. Бернулли: почему статистика $\sum X_i$?

Для $X_i\stackrel{iid}{\sim}\operatorname{Bern}(p)$ положим $K=\sum_{i=1}^n X_i$. Тогда \[L(p;x)=p^K(1-p)^{n-K},\] \[\frac{L(p_1;x)}{L(p_0;x)}=\left(\frac{p_1(1-p_0)}{p_0(1-p_1)}\right)^K\left(\frac{1-p_1}{1-p_0}\right)^n.\] Если $p_1>p_0$, отношение растёт по $K$: отвергаем при большом $K$. Если $p_1<p_0$, при маленьком $K$.

При дискретном распределении порог без рандомизации обычно даёт консервативный уровень. Для точного уровня можно рандомизировать на одном граничном значении $K$.

Show Python code

n=25;p0=.35;p1=.60;alpha=.05
k=np.arange(n+1)
llr=stats.binom.logpmf(k,n,p1)-stats.binom.logpmf(k,n,p0)
cut=next(j for j in range(n+1) if stats.binom.sf(j-1,n,p0)<=alpha)
fig,axs=plt.subplots(1,2,figsize=(12,4.4))
axs[0].plot(k,llr,'o-',lw=2);axs[0].axvline(cut,color='crimson',ls='--');axs[0].set(xlabel='K',ylabel='log(L₁/L₀)',title='Отношение правдоподобия монотонно')
axs[1].bar(k,stats.binom.pmf(k,n,p0),alpha=.6,label='H₀')
axs[1].bar(k,stats.binom.pmf(k,n,p1),alpha=.5,label='H₁')
axs[1].axvline(cut-.5,color='crimson',ls='--',label=f'Критическая область K ≥ {cut}')
axs[1].legend();axs[1].set(xlabel='K',ylabel='Вероятность',title='Ошибка I рода и мощность')
plt.tight_layout();plt.show()
print('Фактический уровень:',stats.binom.sf(cut-1,n,p0))
print('Мощность при p₁:',stats.binom.sf(cut-1,n,p1))

Фактический уровень: 0.02546060824372137
Мощность при p₁: 0.732282173361202

4 4. Монотонное отношение правдоподобия (MLR)

Семейство $\{P_\theta\}$ имеет неубывающее отношение правдоподобия по статистике $T(X)$, если при любых $\theta_1>\theta_0$ отношение $p_{\theta_1}(x)/p_{\theta_0}(x)$ является неубывающей функцией $T(x)$ (на общей области определения плотностей). Если отношение невозрастает, направление неравенств меняется.

Теорема о р.н.м. критерии (Карлин—Рубин)

При неубывающем отношении правдоподобия по $T$ критерий $H_0:\theta\le\theta_0$ против $H_1:\theta>\theta_0$ отвергает при $T>c$ (при необходимости с рандомизацией на $T=c$). Порог выбирается по граничному распределению при $\theta_0$.

Для альтернативы $\theta<\theta_0$ критерий отвергает при малых $T$. Это можно получить заменой параметра $\eta=-\theta$ и статистики $-T$.

Почему оптимален для всех альтернатив сразу? Лемма Неймана—Пирсона для каждой фиксированной $\theta_1>\theta_0$ приводит к одной и той же форме критической области $\{T>c\}$; критический порог определяется только $\theta_0$ и $\alpha$. Монотонность мощности по параметру гарантирует контроль ошибки на всём составном $H_0$.

Пример: в $N(\mu,\sigma^2)$ с известной $\sigma^2$ статистика $T=\sum X_i$ имеет MLR по $\mu$. Для $\operatorname{Exp}(\text{rate }\lambda)$ отношение правдоподобия по $T=\sum X_i$ невозрастает по $T$ при росте $\lambda$.

5 5. Равномерное распределение и зависимость носителя от параметра

Если $X_i\sim U[0,\theta]$, то совместная плотность \[p_\theta(x)=\theta^{-n}\mathbf1\{0\le X_{(1)},\ X_{(n)}\le\theta\},\qquad X_{(n)}=\max_iX_i.\] Для $0\le t\le\theta$ имеем \[P_\theta(X_{(n)}\le t)=\left(\frac t\theta\right)^n.\] Для проверки $H_0:\theta=\theta_0$ против $H_1:\theta<\theta_0$ естественна область $X_{(n)}\le c\theta_0$, где $c=\alpha^{1/n}$.

Поддержка зависит от параметра

При $\theta>\theta_0$ наблюдение $X_{(n)}>\theta_0$ невозможно при $H_0$ и возможно при альтернативе. Такие точки следует включать в критическую область без увеличения уровня. Это особенно важно для двусторонней альтернативы в задаче 6.

Show Python code

n=12;theta0=1.;alpha=.05;B=40000
mx=rng.random((B,n)).max(axis=1)*theta0
x=np.linspace(0,1,500)
fig,ax=plt.subplots()
ax.hist(mx,bins=65,density=True,alpha=.45,label='Монте-Карло')
ax.plot(x,n*x**(n-1),lw=2,label='Теоретическая плотность максимума')
c=alpha**(1/n);ax.axvline(c,color='crimson',ls='--',label=f'c = α^(1/n) = {c:.3f}')
ax.set(xlabel='Максимум выборки',ylabel='Плотность',title='Распределение максимума из U[0,1]')
ax.legend();plt.show()
print('Эмпирический уровень:',np.mean(mx<=c),'Теоретический:',alpha)

Эмпирический уровень: 0.049825 Теоретический: 0.05

6 6. Экспоненциальное семейство и точные критические значения

Если $X_i\sim\operatorname{Exp}(\lambda)$, где $\lambda$ — интенсивность (rate), то \[L(\lambda;x)=\lambda^n\exp\left(-\lambda\sum_i x_i\right),\qquad T=\sum_iX_i\sim\operatorname{Gamma}(n,\text{scale}=1/\lambda).\] Эквивалентно, $2\lambda T\sim\chi^2_{2n}$.

Для $H_1:\lambda>\lambda_0$ критическая область $T\le\chi^2_{2n;\alpha}/(2\lambda_0)$.
Для $H_1:\lambda<\lambda_0$ критическая область $T\ge\chi^2_{2n;1-\alpha}/(2\lambda_0)$.

Для $N(\mu,1)$ статистика $\sqrt n(\overline X-\mu_0)\sim N(0,1)$ при границе $H_0$; для правосторонней альтернативы сравниваем с $z_{1-\alpha}$, для левосторонней — с $z_\alpha$.

Show Python code

n=8;lam0=1.;alpha=.05
xx=np.linspace(0,20,600)
left=stats.gamma.ppf(alpha,a=n,scale=1/lam0)
right=stats.gamma.isf(alpha,a=n,scale=1/lam0)
fig,ax=plt.subplots()
for lam in [.65,1.,1.6]:
    ax.plot(xx,stats.gamma.pdf(xx,a=n,scale=1/lam),label=f'λ = {lam}')
ax.axvline(left,color='crimson',ls='--',label='Левый порог')
ax.axvline(right,color='darkgreen',ls='--',label='Правый порог')
ax.set(title='Сумма экспоненциальных наблюдений',xlabel='T = ΣXᵢ',ylabel='Плотность')
ax.legend();plt.show()

7 7. Двусторонние альтернативы и отсутствие UMP

Для $H_1:\theta\ne\theta_0$ наиболее мощные критерии против альтернатив справа и слева могут требовать противоположных критических областей. Поэтому UMP-критерий для двусторонней альтернативы часто не существует. Это нужно доказывать, а не просто утверждать.

Для Бернулли можно сравнить единственные (или вынужденно совпадающие с ними) наиболее мощные области для $p_1>p_0$ и $p_2<p_0$ по лемме Неймана—Пирсона. При дискретности нужно аккуратно учитывать граничные уровни и возможность рандомизации.

Рандомизированный критерий — функция $\varphi(x)\in[0,1]$, задающая вероятность отвергнуть $H_0$ при наблюдении $x$. Его мощность $\beta(\theta,\varphi)=E_\theta\varphi(X)$. Это расширение полезно для точного уровня в дискретных моделях.

Show Python code

n=20;p0=.5
k=np.arange(n+1);a=.05
cl=next(j for j in range(n+1) if stats.binom.cdf(j,n,p0)>a)-1
cr=next(j for j in range(n+1) if stats.binom.sf(j-1,n,p0)<=a)
ps=np.linspace(.02,.98,300)
fig,ax=plt.subplots()
ax.plot(ps,stats.binom.cdf(cl,n,ps),label=f'Левый критерий: K ≤ {cl}')
ax.plot(ps,stats.binom.sf(cr-1,n,ps),label=f'Правый критерий: K ≥ {cr}')
ax.axvline(p0,color='black',ls='--');ax.set(xlabel='p',ylabel='Мощность',title='Противоположные направления оптимальности')
ax.legend();plt.show()

8 8. Дисперсия нормального распределения: простая против простой

В задаче 5 рассматривается $X_i\sim N(0,\theta)$, где $\theta>0$ — дисперсия, а не стандартное отклонение. Плотность выборки \[L(\theta;x)=(2\pi\theta)^{-n/2}\exp\left(-\frac{\sum_i x_i^2}{2\theta}\right).\] Для $H_0:\theta=\theta_0$ против $H_1:\theta=\theta_1$ логарифм отношения правдоподобия: \[\log\frac{L(\theta_1;x)}{L(\theta_0;x)}=-\frac n2\log\frac{\theta_1}{\theta_0}+\frac12\left(\frac1{\theta_0}-\frac1{\theta_1}\right)\sum_i x_i^2.\] Следовательно, при $\theta_1>\theta_0$ отвергаем для больших $T=\sum X_i^2$, при $\theta_1<\theta_0$ — для малых. При $H_0$ $T/\theta_0\sim\chi_n^2$, поэтому порог вычисляется точно. Мощность находится из распределения $T/\theta_1\sim\chi_n^2$.

Show Python code

n=15;theta0=1.;theta1=1.8;alpha=.05
threshold=theta0*stats.chi2.isf(alpha,n)
thetas=np.linspace(.45,2.6,350)
power=stats.chi2.sf(threshold/thetas,n)
fig,ax=plt.subplots()
ax.plot(thetas,power,lw=2)
ax.axvline(theta0,color='black',ls='--',label='H₀')
ax.axhline(alpha,color='gray',ls=':',label='α')
ax.scatter([theta1],[stats.chi2.sf(threshold/theta1,n)],s=75,zorder=5,label='θ₁')
ax.set(title='Мощность критерия дисперсии при θ₁ > θ₀',xlabel='Истинная дисперсия θ',ylabel='Мощность',ylim=(0,1))
ax.legend();plt.show()

9 9. Особая задача 1: равномерное против показательного

Для одного наблюдения $X$ гипотезы $H_0:U[0,1]$ и $H_1:\operatorname{Exp}(1)$ — простые. На $[0,1]$ отношение плотностей равно $e^{-x}$, а при $x>1$ плотность при $H_0$ равна нулю, при $H_1$ положительна.

Поэтому при построении критерия Неймана—Пирсона обязательно рассмотрите две части критической области: $x>1$ (бесплатная мощность) и участок внутри $[0,1]$, выбираемый по отношению правдоподобия. Мощность вычисляется интегрированием экспоненциальной плотности по обеим частям.

Show Python code

x=np.linspace(0,3,700)
fig,ax=plt.subplots()
ax.plot(x,np.where(x<=1,1.,0.),lw=2,label='H₀: U[0,1]')
ax.plot(x,np.exp(-x),lw=2,label='H₁: Exp(1)')
ax.axvspan(1,3,color='orange',alpha=.13,label='x > 1: невозможно при H₀')
ax.set(title='Разные носители и «бесплатная» мощность',xlabel='x',ylabel='Плотность',ylim=(0,1.12))
ax.legend();plt.show()

10 10. Допустимость и несмещённость критериев

Критерий $S$ называется несмещённым (в смысле мощности), если \[\inf_{\theta\in\Theta_1}\beta(\theta,S)\ge\sup_{\theta\in\Theta_0}\beta(\theta,S).\] Это понятие не связано с несмещённостью оценок вида $E\widehat\theta=\theta$.

Критерий $S$ допустим, если не существует другого критерия $R$ с $\beta(\theta,R)\ge\beta(\theta,S)$ для всех параметров альтернативы (и при сопоставимом контроле ошибки I рода), причём строгое неравенство выполняется хотя бы при одной альтернативе.

Подсказка к задаче 8

Предположите, что другой критерий $R$ уровня $\alpha$ не менее мощен на всей альтернативе и строго лучше где-то. Сопоставьте нижнюю грань мощности $S$ на альтернативе с верхней гранью мощности $S$ на нулевой гипотезе. Для строгого доказательства в формулировке задачи важно учитывать равномерную наиболее мощность $S$: она сразу запрещает строгое улучшение среди критериев уровня $\alpha$.

11 11. Практикум: проверяем теорию методом Монте-Карло

11.1 11.1. Эмпирическая ошибка I рода и мощность

Смоделируем правосторонний критерий нормального среднего и сравним с точной формулой.

Show Python code

B=30000;n=30;alpha=.05;mu0=0.;mus=[0.,.15,.3,.5]
crit=stats.norm.isf(alpha)
for mu in mus:
    z=np.sqrt(n)*(rng.normal(mu,1,size=(B,n)).mean(axis=1)-mu0)
    empirical=np.mean(z>crit)
    theory=stats.norm.sf(crit-np.sqrt(n)*(mu-mu0))
    print(f'μ={mu:.2f} | эмпирическая мощность={empirical:.4f} | теория={theory:.4f}')

μ=0.00 | эмпирическая мощность=0.0495 | теория=0.0500
μ=0.15 | эмпирическая мощность=0.2079 | теория=0.2052
μ=0.30 | эмпирическая мощность=0.4971 | теория=0.4993
μ=0.50 | эмпирическая мощность=0.8653 | теория=0.8630

11.2 11.2. Размер выборки и вероятность ошибки II рода

При фиксированной альтернативе $\mu_1>\mu_0$ вероятность ошибки II рода равна $\Phi(z_{1-\alpha}-\sqrt n(\mu_1-\mu_0))$. Это позволяет планировать объём выборки для заданной мощности.

Show Python code

ns=np.arange(2,201);delta=.35;alpha=.05
power=stats.norm.sf(stats.norm.isf(alpha)-np.sqrt(ns)*delta)
fig,ax=plt.subplots()
ax.plot(ns,power,lw=2);ax.axhline(.8,color='crimson',ls='--',label='Мощность 80%')
ax.set(title='Планирование эксперимента',xlabel='Размер выборки n',ylabel='Мощность',ylim=(0,1))
ax.legend();plt.show()
print('Минимальное n для мощности ≥ 80%:',ns[np.flatnonzero(power>=.8)[0]])

Минимальное n для мощности ≥ 80%: 51

12 12. Задачи для самостоятельного решения (стр. 75–76)

Ниже приведены все восемь задач из пособия. Полные решения намеренно не включены. Под каждой задачей — раскрываемая подсказка с нужным инструментом.

Задача 1

Имеется $X_1$ — выборка объёма 1. Основная гипотеза $H_0$ состоит в том, что $X_1$ имеет равномерное распределение на отрезке $[0,1]$, альтернатива — в том, что $X_1$ имеет показательное распределение с параметром 1. Постройте наиболее мощный критерий уровня значимости $\alpha$ для различения этих гипотез и вычислите его мощность.