Розділ 9

Чому саме квадрати: статистика

У розділі 2 ми обрали квадрати за зручність: гладкість, формула, єдиний мінімум. Теорія ймовірностей дає глибшу відповідь. Квадрати — найкращий вибір, коли похибки гауссові, а гауссові вони дуже часто. Вона ж пояснює, звідки беруться «±» біля результатів.

9.1. Два питання без відповіді

Проблема

У розділі 2 кожен критерій давав свою «найкращу» оцінку: середнє, медіану, середину розмаху. Квадрати ми обрали за зручність. Але чи є причина вважати, що вони дають найточнішу відповідь?
У розділі 8 ми писали $H_C = 102{,}0119 \pm 1{,}3$ мм. Що саме означає це «±», і чому воно обчислюється формулою $\sigma_0\sqrt{Q_{jj}}$?

Щоб відповісти, треба сказати, що ми знаємо про похибки. Будемо вважати похибку кожного виміру випадковою величиною: до повторення експерименту невідомо, якою вона буде, але відомо, як вона розподілена.

9.2. Максимальна правдоподібність

Нехай похибки незалежні й мають нормальний (гауссів) розподіл з нульовим середнім і стандартним відхиленням $\sigma$. Імовірність (точніше, щільність) отримати саме такі нев'язки $r_1, \dots, r_n$ за параметрів $\boldsymbol\theta$:

$$\mathcal L(\boldsymbol\theta) = \prod_{i=1}^n \frac{1}{\sigma\sqrt{2\pi}}\,e^{-r_i^2/(2\sigma^2)} = \text{const}\cdot\exp\Bigl(-\frac{1}{2\sigma^2}\sum_{i=1}^n r_i^2\Bigr).$$

Принцип максимальної правдоподібності: беремо ті параметри, за яких спостережені дані найімовірніші. Експонента найбільша, коли показник найменший, тобто коли мінімальна сума квадратів. МНК — це метод максимальної правдоподібності для гауссового шуму. Так і обґрунтував його Гаусс у 1809 році.

А якщо шум інший? Той самий принцип дає інші критерії, і тут знову з'являються знайомі з розділу 2 кандидати:

Розподіл похибки	Щільність ∝	Максимальна правдоподібність	Оцінка для одного невідомого
Гаусса	$e^{-r^2/2\sigma^2}$	$\min\sum r_i^2$ (L2)	середнє
Лапласа	$e^{-\|r\|/b}$	$\min\sum \|r_i\|$ (L1)	медіана
рівномірний на $[-a, a]$	стала всередині, 0 зовні	$\min\max\|r_i\|$ (L∞)	середина розмаху

Отже, критерій у розділі 2 обирали не за смаком: кожен з них оптимальний для «свого» шуму. Перевіримо це експериментом. Повторимо той самий вимір 2000 разів і подивимося, наскільки розкидаються три оцінки.

Гаусса Лапласа рівномірний

За гауссового шуму середнє справді перемагає. Медіана розкидається приблизно на 25% більше (для великих $n$ — у $\sqrt{\pi/2} \approx 1{,}25$ раза), тобто щоб медіаною досягти тієї ж точності, треба в $\pi/2 \approx 1{,}57$ раза більше вимірів. Але варто перемкнути шум, і переможець зміниться.

9.3. Звідки береться гауссів шум

Чому ж так часто вважають, що шум саме гауссів? Похибка реального виміру — це сума багатьох дрібних незалежних причин: рефракція, центрування приладу, відлік по рейці, коливання температури, дрижання рук. Центральна гранична теорема каже: сума великої кількості незалежних дрібних доданків має приблизно гауссів розподіл, які б розподіли не мали окремі доданки. У геодезії цю думку називають гіпотезою елементарних похибок (Гаген, 1837).

Уже за $k \approx 5$ гістограму важко відрізнити від гауссової кривої. Тому квадрати — розумний вибір «за замовчуванням». Але саме «за замовчуванням»: промахи (розділ 1) — не сума дрібних причин, а одна велика. Вони роблять хвости розподілу товстішими, як у Лапласа, і тоді сума квадратів перестає бути найкращою. До цього повернемося в розділі 10.

9.4. Теорема Гаусса — Маркова

Є й друге обґрунтування, яке взагалі не вимагає гауссовості. Припустімо лише, що похибки мають нульове середнє, однакову дисперсію $\sigma^2$ і некорельовані. Тоді серед усіх оцінок, які

є лінійними функціями вимірів (як середнє, але не як медіана), і
є незміщеними (у середньому по багатьох експериментах дають правильну відповідь),

оцінка МНК має найменшу дисперсію для кожного параметра. Англійською її називають BLUE: Best Linear Unbiased Estimator. Якщо дисперсії різні, найкращим стає зважений МНК з вагами $p_i \propto 1/\sigma_i^2$. Ось і відповідь, чому ваги в розділі 8 були саме такими.

Зверніть увагу на слово «лінійними». За шуму Лапласа медіана обіграла середнє, і суперечності тут немає: медіана — нелінійна оцінка, тож теорема про неї нічого не каже.

9.5. Точність оцінок і еліпс похибок

Оцінка МНК — лінійна функція вимірів: $\hat{\boldsymbol\theta} = (A^\mathsf{T}A)^{-1}A^\mathsf{T}\mathbf y$. Тому шум у $\mathbf y$ передається в параметри за простим правилом (виведення — у подробицях):

$$\operatorname{Cov}(\hat{\boldsymbol\theta}) = \sigma^2\,(A^\mathsf{T}A)^{-1}, \qquad\text{а для зваженого МНК}\quad \sigma_0^2\,(A^\mathsf{T}PA)^{-1}.$$

На діагоналі стоять дисперсії параметрів. Звідси «±» у розділі 8: $\sigma_{H_j} = \sigma_0\sqrt{Q_{jj}}$. Невідоме $\sigma$ замінюють оцінкою $\hat\sigma^2 = S/(n - k)$ з розділу 5. Недіагональні елементи — коваріації: вони кажуть, що похибки параметрів пов'язані. Для двох параметрів матриця коваріацій описує еліпс похибок: область, куди з певною ймовірністю потрапляє оцінка.

Формула не бреше. Хмара розв'язків лягає точно в передбачений еліпс. Суцільний еліпс (1σ) містить близько 39% точок, пунктирний — 95%. Для двох вимірів «один сигма» означає 39%, а не 68%, як на прямій.
Точність залежить від геометрії, а не лише від приладу. Поставте станції з одного боку: та сама похибка відстаней у 1 м дає еліпс у кілька метрів завдовжки. Відношення розміру похибки точки до похибки виміру — геометричний фактор. У GNSS його називають DOP (dilution of precision), і приймач показує його в кожен момент часу.
Точність можна спланувати заздалегідь. Матриця $(J^\mathsf{T}J)^{-1}$ залежить лише від розташування станцій, а не від результатів вимірів. Геодезисти так і проєктують мережі: спершу розраховують очікувані еліпси, потім ідуть у поле.

9.6. Статистичний погляд на машинне навчання

Зв'язок з ML

Функція втрат = мінус логарифм правдоподібності. MSE відповідає гауссовому шуму, MAE — шуму Лапласа. Для класифікації та сама ідея дає крос-ентропію (cross-entropy): це мінус логарифм правдоподібності для ймовірностей класів. Вибір функції втрат — це неявне припущення про те, як влаштовано шум.
Регуляризація — це апріорне знання. Якщо до правдоподібності додати гауссове апріорне припущення «параметри швидше малі, ніж великі», максимум апостеріорної ймовірності (MAP) дає МНК зі штрафом $\lambda|\boldsymbol\theta|^2$. Це ridge-регресія з розділу 10.
Незміщеність — не все. Теорема Гаусса — Маркова обирає найкращу серед незміщених оцінок. Трохи зміщена оцінка (ridge, обрізаний SVD з розділу 6) може мати значно меншу дисперсію й загалом меншу помилку. Це той самий компроміс зміщення й дисперсії з розділу 7.
Невизначеність прогнозів. $\operatorname{Cov}(\hat{\boldsymbol\theta})$ дає не лише «±» для параметрів, а й довірчі інтервали для передбачень моделі. statsmodels друкує їх у таблиці summary(): стандартні похибки, t-статистики, p-значення. Уся ця таблиця — формули цього розділу.

🔍 Хочете більше подробиць? Виведення коваріації, n − k, доведення Гаусса — Маркова, довірчі інтервали, χ²-тест

Звідки $\sigma^2(A^\mathsf{T}A)^{-1}$

Позначимо $B = (A^\mathsf{T}A)^{-1}A^\mathsf{T}$, тоді $\hat{\boldsymbol\theta} = B\mathbf y$. Для будь-якої лінійної функції випадкового вектора $\operatorname{Cov}(B\mathbf y) = B\operatorname{Cov}(\mathbf y)B^\mathsf{T}$. Якщо $\operatorname{Cov}(\mathbf y) = \sigma^2 I$:

$$\operatorname{Cov}(\hat{\boldsymbol\theta}) = \sigma^2 (A^\mathsf{T}A)^{-1}A^\mathsf{T}A(A^\mathsf{T}A)^{-1} = \sigma^2 (A^\mathsf{T}A)^{-1}.$$

Звідси, зокрема, випливає, що дисперсія середнього з $n$ вимірів дорівнює $\sigma^2/n$ (там $A^\mathsf{T}A = n$), тобто похибка спадає як $1/\sqrt n$. Саме це показував інтерактив «Скільки разів міряти?» у розділі 1.

Чому ділимо на $n - k$

Нев'язки $\mathbf r = (I - P)\mathbf y$ — проєкція шуму на $(n-k)$-вимірний підпростір (розділ 5). Звідси $\mathbb E\,|\mathbf r|^2 = \sigma^2 \operatorname{tr}(I - P) = \sigma^2(n - k)$, тож $\hat\sigma^2 = S/(n - k)$ — незміщена оцінка. Ділення на $n$ систематично занижувало б похибку: модель «з'їла» $k$ напрямків шуму, підлаштувавши під них параметри.

Ідея доведення теореми Гаусса — Маркова

Будь-яку іншу лінійну незміщену оцінку можна записати як $\tilde{\boldsymbol\theta} = (B + D)\mathbf y$ з $DA = 0$ (інакше вона була б зміщеною). Тоді $\operatorname{Cov}(\tilde{\boldsymbol\theta}) = \sigma^2\bigl(BB^\mathsf{T} + DD^\mathsf{T}\bigr)$, бо перехресні доданки $BD^\mathsf{T} = (A^\mathsf{T}A)^{-1}(DA)^\mathsf{T} = 0$. Матриця $DD^\mathsf{T}$ невід'ємно визначена, тож дисперсія будь-якої іншої оцінки не менша. Знову Піфагор: «зайва» частина $D$ перпендикулярна до корисної.

Довірчі інтервали й перевірка гіпотез

За гауссового шуму величина $(\hat\theta_j - \theta_j)/\hat\sigma_{\theta_j}$ має розподіл Стьюдента з $n - k$ ступенями свободи. Звідси довірчий інтервал $\hat\theta_j \pm t_{n-k}\,\hat\sigma_{\theta_j}$ і тест значущості: якщо нуль лежить в інтервалі, параметр (наприклад, коефіцієнт при $t^3$ з розділу 7) можна прибрати.

Для контролю мережі перевіряють саму величину $\hat\sigma_0$. Статистика $(n-k)\,\hat\sigma_0^2/\sigma_0^2$ має розподіл $\chi^2$ з $n - k$ ступенями свободи. У мережі з розділу 8 з промахом ($\hat\sigma_0 = 11{,}2$ мм/√км при очікуваних близько 2) вона дорівнює приблизно $4\cdot(11{,}2/2)^2 \approx 125$. Для $\chi^2$ з чотирма ступенями свободи це неймовірно велике значення, і тест однозначно сигналізує про промах. Далі шукають, який вимір винен (тест Баарди, розділ 10).

Еліпси й імовірності

Для двовимірного гауссового розподілу квадрат «відстані Махаланобіса» $\boldsymbol\delta^\mathsf{T}\operatorname{Cov}^{-1}\boldsymbol\delta$ має розподіл $\chi^2$ з двома ступенями свободи. Імовірність потрапити в еліпс «$k$ сигм» дорівнює $1 - e^{-k^2/2}$: 39,3% для $k = 1$, 86,5% для $k = 2$, 95% для $k \approx 2{,}45$. У тривимірному просторі (GNSS) відповідні числа ще менші, тож «1σ» там означає лише близько 20%.

9.7. Код

Python: змагання трьох оцінок за трьох видів шуму, перевірка формули $\sigma^2(A^\mathsf{T}A)^{-1}$ повторними експериментами та порівняння ділення на $n$ і на $n - k$. JavaScript: центральна гранична теорема «гістограмою з символів» і максимальна правдоподібність перебором для вимірів з промахом.

9.8. Перевірте себе

Виміри: 1, 2, 9. Похибки вважаємо розподіленими за Лапласом. Яка оцінка максимальної правдоподібності?

Вимір має похибку $\sigma = 5$ мм. Яка похибка (у мм) середнього з 25 таких вимірів?

Пряму підібрали за 10 точками, сума квадратів нев'язок $S = 12$. Чому дорівнює незміщена оцінка $\hat\sigma^2$?

Який відсоток розв'язків (приблизно) потрапляє у двовимірний еліпс похибок «1σ»?

Далі: МНК і машинне навчання

Тепер у нас є все, щоб подивитися на МНК очима машинного навчання. У розділі 10 кулька з розділу 3 нарешті покотиться всерйоз (градієнтний спуск), до суми квадратів додасться штраф (ridge-регресія), а промахи, які МНК «розмазує», ми навчимося ігнорувати (функція Губера, RANSAC).