§2. Сложная гипотеза: когда параметры неизвестны

Метод группировки с последующим применением критерия хи-квадрат годится и для проверки сложной гипотезы $H_0'$ о принадлежности неизвестной функции распределения элементов выборки некоторому заданному классу $\mathcal F = \{F(x, \boldsymbol\theta),\ \boldsymbol\theta \in \Theta \subseteq \mathbb R^k\}$.

В этом случае общая (при всевозможных $\boldsymbol\theta$) область значений $\xi_1$ тоже разбивается на $N$ промежутков $\Delta_j = (a_j, b_j]$, и $\nu_j$ — число элементов выборки в $\Delta_j$. Но теперь вероятности $\mathbf P(\xi_1 \in \Delta_j)$ при $H_0'$ не заданы однозначно, а являются функциями от $\boldsymbol\theta$: $p_j(\boldsymbol\theta) = F(b_j, \boldsymbol\theta) - F(a_j, \boldsymbol\theta)$. Из-за этого нельзя просто подставить $p_j(\boldsymbol\theta)$ вместо $p_j^0$ в (1).

Р. Фишер (1924) доказал: если подставить $p_j(\tilde{\boldsymbol\theta})$, где $\tilde{\boldsymbol\theta}$ — оценка максимального правдоподобия, основанная на частотах, то при некоторых условиях на класс $\mathcal F$ статистика

снова имеет предельным законом хи-квадрат, только уже с $(N - 1 - k)$ степенями свободы, где $k$ — размерность вектора $\boldsymbol\theta$.

Смысл по-пацански: каждый оценённый по данным параметр «подгоняет» теорию под данные и съедает одну степень свободы. Подогнали два параметра — минус две.

Определение. Значением оценки максимального правдоподобия, основанной на частотах, $\tilde{\boldsymbol\theta}$ служит вектор $\boldsymbol\theta = (\theta_1, \dots, \theta_k)$, на котором достигается максимум вероятности

$$\mathbf P(\nu_1 = l_1, \dots, \nu_N = l_N) = \frac{n!}{l_1!\dots l_N!}[p_1(\boldsymbol\theta)]^{l_1}\dots[p_N(\boldsymbol\theta)]^{l_N}.$$

Это равносильно максимизации по $\boldsymbol\theta$ функции

$$\sum_{j=1}^Nl_j\ln p_j(\boldsymbol\theta) \tag{7}$$

или (для гладких моделей) решению системы, вообще говоря, нелинейных уравнений

$$\sum_{j=1}^Nl_j\frac{\partial\ln p_j(\boldsymbol\theta)}{\partial\theta_m} = 0, \qquad m = 1, \dots, k. \tag{8}$$

Пример 2. Пуассоновская модель

Положим $\pi_m(\theta) = e^{-\theta}\theta^m/m!$, $m \ge 0$. Возьмём промежутки $\Delta_j = [j - 1, j)$, $j = 1, \dots, N - 1$; $\Delta_N = [N - 1, \infty)$. Тогда $p_j(\theta) = \pi_{j-1}(\theta)$ при $j = 1, \dots, N - 1$ и $p_N(\theta) = \sum_{m=N-1}^\infty\pi_m(\theta)$.

Так как $\theta$ — скалярный параметр и $(d/d\theta)\ln\pi_m(\theta) = m/\theta - 1$, система (8) сводится к одному уравнению. С учётом $l_1 + \dots + l_N = n$ получаем

Первый член в скобках равен сумме всех значений $\xi_i$, меньших или равных $N - 2$. Второй — это $l_N\mathbf M(\xi_1 \mid \xi_1 \ge N - 1)$ (П7), он приближённо равен сумме всех значений $\xi_i \ge N - 1$. Поэтому решение $\tilde\theta$ уравнения (9) близко к среднему арифметическому $\bar\xi$ — оценке максимального правдоподобия по всей выборке. (Почему $\bar\xi$ — ОМП для пуассоновской модели — вопрос 3.)

По кирпичикам: уравнение (9) — это «среднее», в котором для хвоста, где мы знаем только «попало в последнюю ячейку», подставлено ожидаемое значение внутри этой ячейки. Решается простой итерацией: подставил $\theta$ справа — получил новое $\theta$ слева, и так до сходимости.

Самолёты-снаряды над Лондоном

Применим критерий хи-квадрат к данным о падениях самолётов-снарядов на юг Лондона во время Второй мировой войны (у Лагутина — по книге [81]). Опасность попадания в жилые дома вместо военных объектов велика при низкой точности стрельбы (при так называемой стрельбе по площадной цели).

Карта южной части Лондона была разбита на $n = 24 \times 24 = 576$ небольших участков площадью по $1/4$ кв. км, и на ней отмечены места падения самолётов-снарядов. В таблице — количества участков $l_{j+1}$ ровно с $j$ падениями. Так как участков много, а вероятность попадания самолёта-снаряда в отдельный участок мала, то при справедливости гипотезы о низкой точности стрельбы можно воспользоваться законом редких событий (§1 главы 5): число попаданий на любой участок — приближённо пуассоновская величина с общим для всех участков параметром $\theta$. Предположим также, что попадания на разные участки независимы.

Общее число падений $M = \sum jl_{j+1} = 537$. Начальная оценка — среднее число падений на участок $\hat\theta = M/n \approx 0{,}932$. Тогда ожидаемые количества участков ровно с $j$ падениями примерно равны $n\pi_j(\hat\theta)$:

Прежде чем вычислять статистику, надо объединить последние 4 столбца, чтобы ожидаемое количество оказалось не меньше 5: $l_4 + \dots + l_7 = 8$ и $n(\pi_4 + \dots + \pi_7) = 8{,}71$.

Теперь заменим начальную оценку $\hat\theta$ на $\tilde\theta$, максимизируя по $\theta$ функцию (7) на компьютере (удобно вычислять $p_5(\theta) = 1 - p_1(\theta) - \dots - p_4(\theta)$). Ответ: $\tilde\theta = 0{,}930$, отличие от $\hat\theta$ всего на 0,002. Значение статистики (6) равно $\tilde X_n^2 = 1{,}02$. Поскольку $N = 5$ и $k = 1$, предельный закон имеет $N - k - 1 = 3$ степени свободы. Фактический уровень значимости 0,80, так что гипотеза о низкой точности стрельбы принимается.

(В книге $\tilde\theta = 0{,}924$, ожидаемые количества 228,6; 211,3; 97,6; 30,1; 8,46 и $\tilde X_n^2 = 1{,}05$. Это не максимум функции (7): по уравнению (9) и прямым перебором с шагом 0,001 максимум достигается при 0,930. Вывод от этого не меняется — уровень 0,80 против 0,79.)

Феллер по поводу этого примера отмечал (перескажем своими словами): большинство жителей верило, что снаряды ложатся кучно, — тогда было бы больше участков без попаданий и больше участков с многими попаданиями. А данные оказались совершенно случайными. Это поучительная иллюстрация: неискушённому человеку случайность часто кажется регулярностью или стремлением к скоплению.

Смысл по-пацански: точности у этого оружия не хватало, чтобы выбрать конкретный квартал, — снаряды ложились как попало, по закону редких событий. А жителям казалось, что их район выбрали нарочно: честная случайность на глаз всегда выглядит кучной.

Объединять маловероятные промежутки необходимо. Если оставить $N = 8$, то $\tilde\theta \approx \hat\theta = 0{,}932$ и $\tilde X_n^2 = 32{,}6$. Это значимо велико для $\chi^2_6$ даже на уровне $10^{-4}$ (фактический уровень около $1{,}2\cdot10^{-5}$). Причина резкого роста — малая величина $np_8(\tilde\theta) \approx 0{,}03$: она придаёт слишком большой вес квадрату отклонения наблюдаемого количества $l_8 = 1$ от ожидаемого. Один участок с 7 попаданиями «стоит» в статистике $(1 - 0{,}031)^2/0{,}031 \approx 30$!

Пример 3. Проверка нормальности по сгруппированным данным

Если данные предварительно группируются, то $\boldsymbol\theta$ можно оценить и до группировки, например обычным методом максимального правдоподобия (§2 и §4 главы 9). Однако, как показывает следующий пример, в этом случае $\tilde X_n^2$ будет сходиться, вообще говоря, к другому предельному закону.

Пусть $\xi_i \sim \mathcal N(\mu, \sigma^2)$, $i = 1, \dots, n$, и оба параметра неизвестны. Для разбиения прямой на промежутки оценим неизвестную функцию распределения $\Phi((x - \mu)/\sigma)$ через $\Phi((x - \bar\xi)/S)$, где $\bar\xi = \frac1n\sum\xi_i$ и $S^2 = \frac1n\sum(\xi_i - \bar\xi)^2$. Чтобы вероятности попадания в промежутки были примерно одинаковы, возьмём в качестве $b_j$ решения уравнений

(таблица T2 или приближение Хамакера для $\Phi^{-1}$ из §5 главы 4). Далее подсчитаем $\nu_j$ — количества попаданий в построенные промежутки. Затем найдём ОМП по частотам $\tilde{\boldsymbol\theta} = (\tilde\theta_1, \tilde\theta_2)$ численным поиском максимума функции (7), исходя из точки $(\bar\xi, S)$.

Для вычисления $p_j(\boldsymbol\theta)$ понадобится приближённо вычислять $y = \Phi(x)$, например алгоритмом Морана: $s = 0$, $t = x\sqrt2/3$; для $i = 0, \dots, 12$ положить $z = i + 0{,}5$ и прибавить к $s$ величину $\sin(zt)e^{-z^2/9}/z$; ответ $y = 0{,}5 + s/\pi$. Он обеспечивает 9 точных десятичных цифр у $\Phi(x)$ при $|x| \le 7$.

Важно отметить, что сами оценки $\bar\xi$ и $S$ использовать в формуле (6) нельзя. При нарушении этого запрета статистика $\tilde X_n^2$ не будет (асимптотически) следовать хи-квадрат с $N - 3$ степенями свободы: график её функции распределения пройдёт несколько ниже графика функции распределения $\chi^2_{N-3}$. Не будет она следовать и $\chi^2_{N-1}$ (как было бы при точно известных параметрах) — график пройдёт несколько выше.

Как показали Г. Чернов и Э. Леман (1954), статистика $\tilde X_n^2$ асимптотически распределена как сумма $\xi_1^2 + \dots + \xi_{N-3}^2 + \gamma_1\xi_{N-2}^2 + \gamma_2\xi_{N-1}^2$, где $\xi_i$ — независимые $\mathcal N(0, 1)$, а числа $\gamma_1$ и $\gamma_2$ лежат между 0 и 1 и зависят от проверяемого закона и способа разбиения.

Предельный закон $\tilde X_n^2$ при $N = 10$ (если подставлять $\bar\xi$ и $S$) зажат в полосе между $\chi^2_7$ и $\chi^2_9$. На уровне 0,95 ширина полосы $16{,}9 - 14{,}1 = 2{,}8$:

Смысл по-пацански: оценки по сырым данным «знают больше», чем частоты, и подгоняют теорию к данным чуть сильнее, чем положено по правилу «минус одна степень свободы на параметр», но не настолько сильно, как ОМП по частотам. Поэтому и закон получается «межеумочный».

Проверь себя