Metoda najmniejszych kwadratów — kompletny traktat: historia, teoria i dowody

Metoda najmniejszych kwadratów nie narodziła się jako ćwiczenie z algebry. Wyrosła z najtrudniejszych problemów rachunkowych XVIII i XIX wieku — wyznaczania kształtu Ziemi i orbit ciał niebieskich z pomiarów obarczonych błędem. Ten artykuł prześledzi całą tę drogę: od prekursorów, przez spór Legendre’a z Gaussem, po dwa różne dowody Gaussa i narodziny współczesnej regresji. Z datami, nazwiskami, oryginalnymi tekstami i dowodami.

Część I. Problem: niezgodne obserwacje (XVIII wiek)

W XVIII wieku astronomia i geodezja postawiły przed matematykami kłopot nowego rodzaju. Aby wyznaczyć kilka nieznanych wielkości (np. parametry orbity albo spłaszczenie Ziemi), wykonywano wiele pomiarów — więcej, niż było niewiadomych. Każdy pomiar był nieco błędny, więc układ równań stawał się nadokreślony i sprzeczny: nie istniało rozwiązanie pasujące do wszystkich obserwacji jednocześnie.

Kształt Ziemi: spór, który napędził geodezję

Newton przewidywał, że Ziemia jest spłaszczona na biegunach (elipsoida obrotacyjna). Aby to rozstrzygnąć, Francuska Akademia Nauk wysłała dwie wyprawy mierzące długość łuku południka: do Laponii (Maupertuis, 1736–37) i do Peru (Bouguer i La Condamine, 1735–44). Z wielu pomiarów geodezyjnych trzeba było wyznaczyć dwa parametry elipsoidy — klasyczny układ nadokreślony.

Prekursorzy: pierwsze sposoby łączenia obserwacji

• Roger Cotes (1682–1716) zaproponował ważoną średnią obserwacji (opublikowane pośmiertnie w Harmonia mensurarum, 1722) — pomysł, by dokładniejszym pomiarom dawać większą wagę.
• Tobias Mayer (1723–1762) przy badaniu libracji Księżyca (ok. 1750) zastosował „metodę uśredniania": dzielił równania na grupy i sumował je, redukując nadokreślony układ do rozwiązywalnego.
• Ruđer Bošković (1711–1787), mierząc z Christopherem Maire łuk południka pod Rzymem (De litteraria expeditione…, 1755), zaproponował kryterium najbliższe duchowi MNK: minimalizować sumę wartości bezwzględnych odchyleń przy warunku, że suma odchyleń (ze znakiem) wynosi zero. To regresja oparta na normie $L_1$ — alternatywa, do której statystyka wróciła dopiero w XX wieku.
• Leonhard Euler i Pierre-Simon Laplace próbowali innych „metod kombinacji", ale żadna nie była ani ogólna, ani prosta.

Brakowało jednej, uniwersalnej zasady. Pojawiła się na początku XIX wieku — i to niemal jednocześnie u dwóch ludzi.

Część II. Legendre: pierwsza publikacja (1805)

Pierwszą publikacją metody był dodatek do dzieła Adriena-Marie Legendre’a Nouvelles méthodes pour la détermination des orbites des comètes (1805), zatytułowany Sur la méthode des moindres carrés. To Legendre nadał metodzie nazwę — méthode des moindres carrés („metoda najmniejszych kwadratów") — i sformułował ją z rozbrajającą jasnością:

Legendre podał metodę jako przepis rachunkowy — skuteczny i elegancki, ale bez uzasadnienia, dlaczego akurat kwadraty błędów, a nie inne ich potęgi. To pytanie zostawił otwarte. Odpowiedź dał Gauss.

Część III. Gauss: geneza i pierwsze uzasadnienie (1809)

Zaginiona planeta Ceres (1801)

1 stycznia 1801 roku Giuseppe Piazzi z Palermo dostrzegł słaby, ruchomy obiekt — planetę karłowatą Ceres. Obserwował ją przez 41 nocy, po czym Ceres zbliżyła się do Słońca i zniknęła w jego blasku. Astronomowie mieli garść niezgodnych pomiarów i jedno pytanie: gdzie szukać planety, gdy wyłoni się zza Słońca?

24-letni wówczas Gauss podjął się zadania. Korzystając z własnych metod — w tym minimalizacji sumy kwadratów błędów dla nadmiarowych obserwacji — obliczył orbitę i wskazał, gdzie szukać. W grudniu 1801 roku Franz von Zach i Heinrich Olbers odnaleźli Ceres blisko przewidywań Gaussa. Sukces rozsławił go w całej Europie.

Theoria Motus (1809) i spór o pierwszeństwo

W Theoria Motus Corporum Coelestium in Sectionibus Conicis Solem Ambientium (1809) Gauss przedstawił metodę wraz z głębokim uzasadnieniem i napisał, że stosował ją od 1795 roku — na długo przed publikacją Legendre’a. Wybuchł jeden z najsłynniejszych sporów o pierwszeństwo w historii nauki. W klasycznym angielskim przekładzie Davisa (1857) Gauss pisze:

Dziś ocenia się go zwykle tak:

• Legendre pierwszy opublikował metodę i nadał jej nazwę;
• Gauss pierwszy jej użył (bez publikacji) i — co ważniejsze — pierwszy wyjaśnił, dlaczego kwadraty są właściwe.

Oryginalne rozumowanie Gaussa

Gauss zadał pytanie, którego Legendre nie postawił: dlaczego suma kwadratów? Punktem wyjścia uczynił zasadę, którą astronomowie i tak stosowali:

Genialność polegała na odwróceniu pytania: zamiast „jaka jest najlepsza wartość przy danym prawie błędów?", Gauss zapytał — jakie prawo błędów $\varphi(\varepsilon)$ sprawia, że najlepszą wartością jest właśnie średnia? Odpowiedź wymusza kształt rozkładu błędów, a tym kształtem okazuje się rozkład normalny.

Ostatni krok jest natychmiastowy: skoro $\Omega \propto \exp(-h^2\sum\varepsilon_i^2)$, to maksymalizacja wiarygodności = minimalizacja sumy kwadratów $\sum\varepsilon_i^2$. Oto najmniejsze kwadraty jako estymacja największej wiarygodności.

Część IV. Laplace: fundament probabilistyczny

Argument Gaussa z 1809 r. opierał się na postulacie średniej — eleganckim, lecz arbitralnym. Niezależną podporę dał Pierre-Simon Laplace. W 1810 roku udowodnił centralne twierdzenie graniczne: suma wielu drobnych, niezależnych wpływów ma rozkład w przybliżeniu normalny.

To wyjaśniło, dlaczego błąd pomiaru — będący sumą wielu maleńkich, niezależnych zaburzeń — jest w praktyce gaussowski. Laplace połączył to z metodą najmniejszych kwadratów (Théorie analytique des probabilités, 1812), dając jej uzasadnienie wielkopróbkowe, niezależne od postulatu Gaussa. Tę zbieżność dwóch dróg nazywa się czasem syntezą Gaussa-Laplace’a — to ona ugruntowała rolę rozkładu normalnego i MNK w całej nauce. Więcej w artykule o centralnym twierdzeniu granicznym.

Część V. Gauss raz jeszcze: twierdzenie Gaussa-Markowa (1823)

Sam Gauss uznał argument z 1809 r. za niezadowalający (zależał od założenia normalności). W rozprawie Theoria combinationis observationum erroribus minimis obnoxiae (1821–1823) podał zupełnie nowe uzasadnienie, niewymagające żadnego konkretnego prawa błędów.

W pełnej wersji (model liniowy $\mathbf{y}=\mathbf{X}\boldsymbol\beta+\boldsymbol\varepsilon$) ten sam argument daje $\hat{\boldsymbol\beta}=(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}$ jako estymator o najmniejszej macierzy kowariancji wśród liniowych nieobciążonych.

Część VI. Od błędów pomiarowych do regresji (XIX–XX w.)

Przez cały XIX wiek najmniejsze kwadraty służyły astronomii i geodezji — do szacowania stałych fizycznych. Przełom w stronę statystyki zależności przyszedł pod koniec stulecia:

• Francis Galton (1822–1911) badając wzrost rodziców i dzieci zauważył „regresję ku przeciętnej" (1885) i ukuł termin regresja. To on przeniósł metodę z „błędów" na związki między zmiennymi.
• Karl Pearson (1857–1936) sformalizował współczynnik korelacji (1896) i rozwinął aparat regresji.
• George Udny Yule (1871–1951) pokazał (1897), że regresja wieloraka to po prostu najmniejsze kwadraty zastosowane do wielu zmiennych objaśniających.
• Ronald A. Fisher (1890–1962) ujednolicił teorię: metoda największej wiarygodności, analiza wariancji, własności estymatorów. W jego ujęciu MNK to estymator największej wiarygodności przy błędach normalnych — domykając koło zapoczątkowane przez Gaussa.

Tak narzędzie astronomów stało się fundamentem ekonometrii: tym samym rachunkiem, którym Gauss łapał planetę, dziś szacujemy zwrot z edukacji czy elastyczność popytu.

Część VII. Dwa filary, jeden wniosek

Uzasadnienie	Rok	Autor	Założenie	Wniosek
Probabilistyczne	1809	Gauss	błędy normalne	MNK = estymator największej wiarygodności
Graniczne (CLT)	1810	Laplace	suma wielu drobnych przyczyn	błędy są w przybliżeniu normalne
Wariancyjne (Gauss-Markow)	1823	Gauss	$\mathbb{E}[\varepsilon]=0$, stała wariancja, brak korelacji	MNK = najlepszy liniowy nieobciążony estymator

Pierwszy filar wyjaśnia, skąd wziął się rozkład normalny. Drugi — dlaczego błędy bywają normalne. Trzeci — że metoda jest optymalna nawet bez normalności. Razem tłumaczą, czemu najmniejsze kwadraty przetrwały dwa stulecia jako domyślne narzędzie nauk ilościowych.

Część VIII. Sformułowanie macierzowe modelu

Współczesna ekonometria zapisuje model liniowy w postaci macierzowej, co czyni cały rachunek zwartym i ogólnym. Mamy $n$ obserwacji i $k$ parametrów:

$$ \mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}, $$

gdzie $\mathbf{y}$ to wektor $n\times 1$ obserwacji zmiennej objaśnianej, $\mathbf{X}$ to macierz $n\times k$ regresorów (pierwsza kolumna zwykle jedynkowa — wyraz wolny), $\boldsymbol{\beta}$ to wektor $k\times 1$ nieznanych parametrów, a $\boldsymbol{\varepsilon}$ to wektor $n\times 1$ składników losowych.

Założenia A1–A4 to dokładnie warunki twierdzenia Gaussa-Markowa. Normalność A5 nie jest potrzebna do optymalności MNK — służy tylko do dokładnych testów $t$ i $F$ w małej próbie.

Część IX. Wyprowadzenie estymatora (postać macierzowa)

Minimalizujemy sumę kwadratów reszt jako iloczyn skalarny:

$$ S(\boldsymbol{\beta}) = \boldsymbol{\varepsilon}^\top\boldsymbol{\varepsilon} = (\mathbf{y}-\mathbf{X}\boldsymbol{\beta})^\top(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}) = \mathbf{y}^\top\mathbf{y} - 2\boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{y} + \boldsymbol{\beta}^\top\mathbf{X}^\top\mathbf{X}\boldsymbol{\beta}. $$

Część X. Geometria: macierz kapeluszowa i rzut prostopadły

Wartości dopasowane to $\hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol\beta} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y} = \mathbf{H}\mathbf{y}$, gdzie

$$ \mathbf{H} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top $$

to macierz kapeluszowa (ang. hat matrix) — rzut prostopadły na przestrzeń kolumn $\mathbf{X}$. Reszty to $\mathbf{e} = \mathbf{y}-\hat{\mathbf{y}} = (\mathbf{I}-\mathbf{H})\mathbf{y} = \mathbf{M}\mathbf{y}$, gdzie $\mathbf{M} = \mathbf{I}-\mathbf{H}$ to macierz anihilująca.

Część XI. Własności estymatora MNK

Podstawiając $\mathbf{y}=\mathbf{X}\boldsymbol\beta+\boldsymbol\varepsilon$ do wzoru na $\hat{\boldsymbol\beta}$, otrzymujemy kluczową tożsamość:

$$ \hat{\boldsymbol\beta} = \boldsymbol\beta + (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\boldsymbol\varepsilon. $$

• Nieobciążoność. $\mathbb{E}[\hat{\boldsymbol\beta}\mid\mathbf{X}] = \boldsymbol\beta + (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\,\mathbb{E}[\boldsymbol\varepsilon\mid\mathbf{X}] = \boldsymbol\beta$ (z A3).
• Macierz kowariancji. Przy sferyczności (A4): $$ \mathrm{Var}(\hat{\boldsymbol\beta}\mid\mathbf{X}) = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top(\sigma^2\mathbf{I})\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1} = \sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}. $$
• Estymator wariancji błędu. $s^2 = \dfrac{\mathbf{e}^\top\mathbf{e}}{n-k}$ jest nieobciążony: $\mathbb{E}[s^2]=\sigma^2$ (dzielnik $n-k$, nie $n$ — analogia poprawki Bessela; $k$ stopni swobody „zużywają" parametry).

Część XII. Twierdzenie Gaussa-Markowa — dowód macierzowy

To uogólnienie dowodu z części V na pełny model liniowy.

Część XIII. Rozkłady i wnioskowanie (przy normalności)

Przy dodatkowym założeniu A5 ($\boldsymbol\varepsilon\sim\mathcal{N}$) dostajemy dokładne rozkłady, na których opiera się wnioskowanie:

$$ \hat{\boldsymbol\beta}\sim\mathcal{N}\!\big(\boldsymbol\beta,\ \sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}\big),\qquad \frac{(n-k)s^2}{\sigma^2}\sim\chi^2_{n-k}. $$

Ponieważ $\sigma^2$ jest nieznane, w statystyce $t$ zastępujemy je przez $s^2$. Dla pojedynczego współczynnika:

$$ \frac{\hat\beta_j-\beta_j}{\operatorname{SE}(\hat\beta_j)}\sim t_{n-k},\qquad \operatorname{SE}(\hat\beta_j)=s\sqrt{[(\mathbf{X}^\top\mathbf{X})^{-1}]_{jj}}. $$

Dla $q$ łącznych ograniczeń liniowych $\mathbf{R}\boldsymbol\beta=\mathbf{r}$ używamy statystyki $F$:

$$ F=\frac{(\mathbf{R}\hat{\boldsymbol\beta}-\mathbf{r})^\top\big[\mathbf{R}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{R}^\top\big]^{-1}(\mathbf{R}\hat{\boldsymbol\beta}-\mathbf{r})}{q\,s^2}\sim F_{q,\,n-k}. $$

To formalna podstawa testów hipotez, p-wartości i przedziałów ufności. Poniżej wyprowadzamy te rozkłady od podstaw — z jednego założenia o normalności składników losowych.

Część XIV. Twierdzenie Frischa-Waugha-Lovella

Jeden z najważniejszych wyników o interpretacji regresji wielorakiej (Frisch i Waugh 1933, Lovell 1963).

To twierdzenie tłumaczy, dlaczego współczynniki regresji wielorakiej interpretujemy jako efekty ceteris paribus, i jest narzędziem dowodowym w teorii danych panelowych i estymatorach efektów stałych. Z FWL wynika też wzór na pojedynczy współczynnik regresji wielorakiej i jego odporność: dodanie nieskorelowanego regresora nie zmienia pozostałych oszacowań.

Część XV. Naruszenia założeń — mapa

Każde założenie A1–A5 może zostać złamane; każde łamanie ma inną diagnozę i lekarstwo. Szczegóły w osobnych artykułach:

Złamane założenie	Skutek	Diagnoza	Lekarstwo
A1 liniowość	obciążenie	wzór reszt	człony kwadratowe, logarytmy
A2 pełny rząd	brak/niestabilność	VIF	usuń/połącz zmienne, regularyzacja
A3 egzogeniczność	obciążenie i niezgodność	argument merytoryczny	zmienne instrumentalne, panel
A4 homoskedastyczność	błędne SE	Breusch-Pagan	odporne SE, WLS
A4 brak autokorelacji	błędne SE	Breusch-Godfrey	błędy HAC
A5 normalność	testy w małej próbie	Q-Q, Shapiro	większa próba (CTG), bootstrap

Część XVI. Rozszerzenia metody

• Ważona MNK (WLS): przy znanej heteroskedastyczności minimalizujemy $\sum w_i\varepsilon_i^2$; $\hat{\boldsymbol\beta}=(\mathbf{X}^\top\mathbf{W}\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{W}\mathbf{y}$.
• Uogólniona MNK (GLS): dla dowolnej $\mathrm{Var}(\boldsymbol\varepsilon)=\sigma^2\boldsymbol\Omega$: $\hat{\boldsymbol\beta}=(\mathbf{X}^\top\boldsymbol\Omega^{-1}\mathbf{X})^{-1}\mathbf{X}^\top\boldsymbol\Omega^{-1}\mathbf{y}$ (Aitken).
• Zmienne instrumentalne / 2SLS: przy endogeniczności, $\hat{\boldsymbol\beta}_{IV}=(\mathbf{Z}^\top\mathbf{X})^{-1}\mathbf{Z}^\top\mathbf{y}$ — patrz IV.
• Regresja grzbietowa (ridge) i lasso: dodajemy karę za wielkość współczynników; ridge: $\hat{\boldsymbol\beta}=(\mathbf{X}^\top\mathbf{X}+\lambda\mathbf{I})^{-1}\mathbf{X}^\top\mathbf{y}$ — stabilizuje przy współliniowości i w wysokim wymiarze.
• Nieliniowa MNK (NLS): minimalizacja $\sum(y_i-f(\mathbf{x}_i,\boldsymbol\beta))^2$ dla nieliniowej $f$, metodami iteracyjnymi (Gauss-Newton).

Część XVII. Metody obliczeniowe

W praktyce nie odwraca się jawnie $\mathbf{X}^\top\mathbf{X}$ — to numerycznie niestabilne, bo wskaźnik uwarunkowania podnosi się do kwadratu: $\kappa(\mathbf{X}^\top\mathbf{X})=\kappa(\mathbf{X})^2$. Lepsze są:

• Rozkład QR: $\mathbf{X}=\mathbf{Q}\mathbf{R}$ (ortogonalno-trójkątny), wtedy $\mathbf{R}\hat{\boldsymbol\beta}=\mathbf{Q}^\top\mathbf{y}$ — rozwiązanie przez podstawianie wsteczne; standard w lm() (R) i lstsq.
• Rozkład SVD: $\mathbf{X}=\mathbf{U}\boldsymbol\Sigma\mathbf{V}^\top$ — najstabilniejszy, radzi sobie nawet z (prawie) zdegenerowanym $\mathbf{X}$; daje pseudoodwrotność Moore’a-Penrose’a.
• Rozkład Choleskiego $\mathbf{X}^\top\mathbf{X}=\mathbf{L}\mathbf{L}^\top$ — szybki, gdy uwarunkowanie jest dobre.
• Eliminacja Gaussa — algorytm, który Gauss rozwinął właśnie po to, by ręcznie rozwiązywać równania normalne; do dziś podstawa rozwiązywania układów liniowych.

Część XVIII. Związki z innymi metodami estymacji

• Największa wiarygodność: przy błędach normalnych MNW pokrywa się dokładnie z MNK (część III) — to ten sam estymator.
• Bayes: z płaskim priorem rozkład a posteriori $\boldsymbol\beta$ jest skupiony wokół $\hat{\boldsymbol\beta}_{\text{MNK}}$; z normalnym priorem sprzężonym otrzymujemy regresję grzbietową (ridge).
• Metoda momentów: warunek $\mathbb{E}[\mathbf{X}^\top\boldsymbol\varepsilon]=\mathbf{0}$ w wersji próbkowej to $\mathbf{X}^\top\mathbf{e}=\mathbf{0}$ — znów równania normalne. MNK jest więc estymatorem metody momentów.

Część XIX. Przykład liczbowy: regresja krok po kroku

Cała teoria sprowadza się do kilku działań arytmetycznych. Przejdźmy je ręcznie na małym zbiorze ($n=5$, jeden regresor i wyraz wolny), aby zobaczyć, skąd biorą się wszystkie liczby raportowane przez program statystyczny.

Krok 1. Sumy i macierze

Macierz planu $\mathbf{X}$ ma w pierwszej kolumnie same jedynki (wyraz wolny), w drugiej wartości $x_i$. Stąd

$$ \mathbf{X}^\top\mathbf{X} = \begin{bmatrix} n & \sum x_i \\ \sum x_i & \sum x_i^2 \end{bmatrix} = \begin{bmatrix} 5 & 15 \\ 15 & 55 \end{bmatrix}, \qquad \mathbf{X}^\top\mathbf{y} = \begin{bmatrix} \sum y_i \\ \sum x_i y_i \end{bmatrix} = \begin{bmatrix} 20 \\ 66 \end{bmatrix}. $$

Wyznacznik $\det(\mathbf{X}^\top\mathbf{X}) = 5\cdot 55 - 15^2 = 275 - 225 = 50$, więc

$$ (\mathbf{X}^\top\mathbf{X})^{-1} = \frac{1}{50}\begin{bmatrix} 55 & -15 \\ -15 & 5 \end{bmatrix} = \begin{bmatrix} 1{,}1 & -0{,}3 \\ -0{,}3 & 0{,}1 \end{bmatrix}. $$

Krok 2. Estymator

$$ \hat{\boldsymbol\beta} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y} = \begin{bmatrix} 1{,}1 & -0{,}3 \\ -0{,}3 & 0{,}1 \end{bmatrix}\begin{bmatrix} 20 \\ 66 \end{bmatrix} = \begin{bmatrix} 22 - 19{,}8 \\ -6 + 6{,}6 \end{bmatrix} = \begin{bmatrix} 2{,}2 \\ 0{,}6 \end{bmatrix}. $$

Ten sam wynik daje wzór skalarny dla regresji prostej (część skupiona w sumach odchyleń od średniej): z $S_{xx} = \sum(x_i-\bar x)^2 = 10$ i $S_{xy} = \sum(x_i-\bar x)(y_i-\bar y) = 6$,

$$ \hat\beta_1 = \frac{S_{xy}}{S_{xx}} = \frac{6}{10} = 0{,}6, \qquad \hat\beta_0 = \bar y - \hat\beta_1\bar x = 4 - 0{,}6\cdot 3 = 2{,}2. $$

Dopasowana prosta: $\hat y = 2{,}2 + 0{,}6\,x$.

Krok 3. Wartości dopasowane i reszty

$i$	$x_i$	$y_i$	$\hat y_i = 2{,}2+0{,}6x_i$	$e_i = y_i-\hat y_i$	$e_i^2$
1	1	2	2{,}8	$-0{,}8$	0{,}64
2	2	4	3{,}4	$0{,}6$	0{,}36
3	3	5	4{,}0	$1{,}0$	1{,}00
4	4	4	4{,}6	$-0{,}6$	0{,}36
5	5	5	5{,}2	$-0{,}2$	0{,}04
			suma	$0{,}0$	2{,}40

Reszty sumują się do zera ($\sum e_i = 0$) — to bezpośrednia konsekwencja równań normalnych przy obecnym wyrazie wolnym. Suma kwadratów reszt wynosi $\text{RSS} = \mathbf{e}^\top\mathbf{e} = 2{,}40$.

Krok 4. Wariancja resztowa i błędy standardowe

Mamy $n=5$, $k=2$, więc $n-k=3$ stopni swobody:

$$ s^2 = \frac{\text{RSS}}{n-k} = \frac{2{,}40}{3} = 0{,}80. $$

Macierz kowariancji estymatora: $\widehat{\mathrm{Var}}(\hat{\boldsymbol\beta}) = s^2(\mathbf{X}^\top\mathbf{X})^{-1}$, której przekątna daje błędy standardowe:

$$ \operatorname{SE}(\hat\beta_0) = \sqrt{0{,}80\cdot 1{,}1} = \sqrt{0{,}88} \approx 0{,}938, \qquad \operatorname{SE}(\hat\beta_1) = \sqrt{0{,}80\cdot 0{,}1} = \sqrt{0{,}08} \approx 0{,}283. $$

Krok 5. Statystyki $t$ i istotność

Przy hipotezie zerowej $\beta_j = 0$:

$$ t_0 = \frac{2{,}2}{0{,}938} \approx 2{,}345, \qquad t_1 = \frac{0{,}6}{0{,}283} \approx 2{,}121. $$

Wartość krytyczna $t$ dla $3$ stopni swobody i poziomu $0{,}05$ (dwustronnie) wynosi $3{,}182$. Obie statystyki są od niej mniejsze, więc przy tak małej próbie nie odrzucamy hipotezy o zerowych współczynnikach na poziomie 5% — mimo wyraźnego dopasowania. To dobra ilustracja, jak mało informacji niesie $n=5$. Zobacz: testy hipotez, wartość p.

Krok 6. Dopasowanie modelu

Całkowita suma kwadratów $\text{TSS} = \sum(y_i-\bar y)^2 = 4+0+1+0+1 = 6$. Stąd

$$ R^2 = 1 - \frac{\text{RSS}}{\text{TSS}} = 1 - \frac{2{,}4}{6} = 0{,}60, $$

skorygowane $\bar R^2 = 1 - (1-R^2)\dfrac{n-1}{n-k} = 1 - 0{,}4\cdot\dfrac{4}{3} \approx 0{,}467$, oraz statystyka

$$ F = \frac{(\text{TSS}-\text{RSS})/(k-1)}{\text{RSS}/(n-k)} = \frac{3{,}6/1}{2{,}4/3} = 4{,}5 = t_1^2, $$

co potwierdza tożsamość $F = t^2$ dla pojedynczego regresora. Model wyjaśnia 60% zmienności $y$. Zobacz: współczynnik $R^2$.

Część XX. Przedział ufności dla E[y₀] a przedział predykcji

Po dopasowaniu modelu zadajemy w punkcie $\mathbf{x}_0$ dwa różne pytania. Pierwsze: gdzie leży średnia warunkowa $\mathbb{E}[y\mid\mathbf{x}_0]=\mathbf{x}_0^\top\boldsymbol\beta$? Drugie: jaką wartość przyjmie nowa, pojedyncza obserwacja $y_0=\mathbf{x}_0^\top\boldsymbol\beta+\varepsilon_0$? Oba mają ten sam estymator punktowy $\hat y_0=\mathbf{x}_0^\top\hat{\boldsymbol\beta}$, lecz różną niepewność — i stąd rozróżnienie na przedział ufności i przedział predykcji.

Zastępując nieznane $\sigma^2$ przez $s^2$ i korzystając z rozkładu $t_{n-k}$, otrzymujemy dwa przedziały (na poziomie ufności $1-\alpha$):

$$ \underbrace{\hat y_0\pm t_{n-k,\,\alpha/2}\;s\sqrt{\mathbf{x}_0^\top(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{x}_0}}_{\text{przedział ufności dla }\mathbb{E}[y_0]},\qquad \underbrace{\hat y_0\pm t_{n-k,\,\alpha/2}\;s\sqrt{1+\mathbf{x}_0^\top(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{x}_0}}_{\text{przedział predykcji dla nowego }y_0}. $$

Część XXI. Dekompozycja wariancji: R², skorygowane R² i tablica ANOVA

Jak dobrze model objaśnia zmienność danych? Odpowiedź daje rozłożenie całkowitej zmienności $\mathbf{y}$ na część wyjaśnioną przez regresję i część resztową — fundament analizy wariancji (ANOVA) oraz współczynnika determinacji $R^2$. Trzy sumy kwadratów (model z wyrazem wolnym; $\bar y$ to średnia $y$):

$$ \underbrace{\textstyle\sum_i (y_i-\bar y)^2}_{\text{TSS — całkowita}}, \qquad \underbrace{\textstyle\sum_i (\hat y_i-\bar y)^2}_{\text{ESS — wyjaśniona}}, \qquad \underbrace{\textstyle\sum_i (y_i-\hat y_i)^2}_{\text{RSS — resztowa}}. $$

Stąd współczynnik determinacji:

$$ R^2=\frac{\text{ESS}}{\text{TSS}}=1-\frac{\text{RSS}}{\text{TSS}}\in[0,1], $$

interpretowany jako udział zmienności $y$ wyjaśniony przez model. Równoważnie $R^2$ jest kwadratem współczynnika korelacji między $y_i$ a $\hat y_i$.

Problem w tym, że $R^2$ nigdy nie maleje po dodaniu regresora (każda nowa kolumna co najwyżej zmniejsza RSS). By karać za nadmiar parametrów, używa się skorygowanego $R^2$:

$$ \bar R^2=1-\frac{\text{RSS}/(n-k)}{\text{TSS}/(n-1)}=1-(1-R^2)\frac{n-1}{n-k}, $$

który rośnie tylko wtedy, gdy nowy regresor poprawia dopasowanie bardziej, niż wynikałoby z samego losu; $\bar R^2$ może być ujemne. Te same sumy kwadratów układają się w tablicę ANOVA, z której wprost wynika statystyka $F$ testu łącznej istotności regresji:

Źródło	Suma kwadratów	df	Średni kwadrat	$F$
Regresja	ESS	$k-1$	$\text{ESS}/(k-1)$	$\dfrac{\text{ESS}/(k-1)}{\text{RSS}/(n-k)}$
Reszty	RSS	$n-k$	$\text{RSS}/(n-k)=s^2$
Całkowita	TSS	$n-1$

Część XXII. Regresja wieloraka i obciążenie przez zmienną pominiętą — przykład liczbowy

Część XIV pokazała, jak twierdzenie FWL odzyskuje współczynnik oczyszczony z wpływu innych zmiennych. Tu policzymy drugą stronę medalu: co się dzieje, gdy istotną zmienną pominiemy. Użyjemy tych samych danych co w części XIV:

$x_1$	1	2	3	4	5
$x_2$	2	1	5	3	4
$y$	3	2	6	5	7

Pełna regresja $y$ na $\{1,x_1,x_2\}$ — rozwiązanie równań normalnych $\hat{\boldsymbol\beta}=(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}$ (dokładnie, w ułamkach):

$$ \hat\beta_0=\tfrac{23}{80}=0{,}2875,\qquad \hat\beta_1=\tfrac{19}{32}=0{,}59375,\qquad \hat\beta_2=\tfrac{27}{32}=0{,}84375. $$

Regresja skrócona, z pominięciem $x_2$ (tylko $y$ na $\{1,x_1\}$):

$$ \tilde\beta_0=\tfrac{13}{10}=1{,}3,\qquad \tilde\beta_1=\tfrac{11}{10}=1{,}1. $$

Współczynnik przy $x_1$ skoczył z $0{,}594$ do $1{,}1$ — niemal dwukrotnie. To obciążenie przez zmienną pominiętą: skrócony współczynnik przypisuje $x_1$ także część efektu $x_2$, bo zmienne są skorelowane.

Aneks A. Zestawienie najważniejszych wzorów

Wielkość	Wzór
Estymator	$\hat{\boldsymbol\beta}=(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}$
Wartości dopasowane	$\hat{\mathbf{y}}=\mathbf{H}\mathbf{y}$, $\mathbf{H}=\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top$
Reszty	$\mathbf{e}=\mathbf{M}\mathbf{y}$, $\mathbf{M}=\mathbf{I}-\mathbf{H}$
Kowariancja	$\mathrm{Var}(\hat{\boldsymbol\beta})=\sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}$
Estymator $\sigma^2$	$s^2=\mathbf{e}^\top\mathbf{e}/(n-k)$
Statystyka $t$	$(\hat\beta_j-\beta_j)/\operatorname{SE}(\hat\beta_j)\sim t_{n-k}$
$R^2$	$1-\mathbf{e}^\top\mathbf{e}/\sum(y_i-\bar y)^2$
Skorygowane $R^2$	$1-(1-R^2)\dfrac{n-1}{n-k}$
Statystyka $F$ (ANOVA)	$\dfrac{\text{ESS}/(k-1)}{\text{RSS}/(n-k)}\sim F_{k-1,\,n-k}$

Aneks B. Oznaczenia

Symbol	Znaczenie
$n,\ k$	liczba obserwacji, liczba parametrów
$\mathbf{y},\ \mathbf{X},\ \boldsymbol\beta,\ \boldsymbol\varepsilon$	wektor obserwacji, macierz regresorów, parametry, błędy
$\hat{\boldsymbol\beta},\ \mathbf{e},\ \hat{\mathbf{y}}$	estymator, reszty, wartości dopasowane
$\mathbf{H},\ \mathbf{M}$	macierz kapeluszowa, anihilująca
$\sigma^2,\ s^2$	wariancja błędu i jej estymator

Aneks C. Sylwetki

Krótkie noty o postaciach, których prace złożyły się na metodę najmniejszych kwadratów i regresję.

Adrien-Marie Legendre (1752–1833) — francuski matematyk, autor fundamentalnych prac z teorii liczb (Essai sur la théorie des nombres), funkcji eliptycznych i geometrii (jego Éléments de géométrie przez stulecie wypierały Euklidesa w szkołach). To on opublikował metodę najmniejszych kwadratów jako pierwszy (1805) i nadał jej nazwę. W ostatnich latach życia stracił rządową pensję, gdy odmówił poparcia kandydata władz w głosowaniu Instytutu; zmarł w ubóstwie. Jedyną zachowaną jego podobizną jest karykatura Boilly’ego.

Carl Friedrich Gauss (1777–1855) — niemiecki matematyk, astronom i fizyk, nazywany princeps mathematicorum („księciem matematyków"). Pochodził z ubogiej rodziny w Brunszwiku; jego talent dostrzeżono w dzieciństwie. W Disquisitiones Arithmeticae (1801) zbudował nowoczesną teorię liczb; obliczył orbitę Ceres; sformułował twierdzenie zasadnicze algebry (rozprawa doktorska, 1799); współtworzył geodezję, magnetyzm i elektrostatykę (jednostka indukcji magnetycznej nosi jego imię). Dla MNK kluczowe są Theoria Motus (1809) i Theoria combinationis (1821–1823).

Pierre-Simon de Laplace (1749–1827) — francuski matematyk i astronom, twórca mechaniki nieba (Traité de mécanique céleste) i nowoczesnej teorii prawdopodobieństwa (Théorie analytique des probabilités, 1812). Jego centralne twierdzenie graniczne wyjaśniło, dlaczego błędy pomiarowe bywają normalne, i połączyło MNK z rachunkiem prawdopodobieństwa. Sprawnie przetrwał rewolucję, cesarstwo i restaurację, zachowując wpływowe stanowiska naukowe.

Ruđer Josip Bošković (1711–1787) — uczony z Dubrownika, jezuita, astronom i fizyk działający we Włoszech. Zaproponował (1755–1760) kryterium najbliższe duchem MNK: minimalizację sumy wartości bezwzględnych odchyleń przy zerowej sumie odchyleń ze znakiem — dzisiejszą regresję $L_1$, do której statystyka wróciła dopiero w XX wieku.

Andriej Markow (1856–1922) — rosyjski matematyk, uczeń Czebyszewa, twórca teorii łańcuchów Markowa. Jego nazwisko trafiło do nazwy twierdzenia Gaussa-Markowa za sprawą późniejszych podręczników, choć sam wynik miał już Gauss w 1823 r.

Francis Galton (1822–1911) — brytyjski uczony, kuzyn Darwina, badacz dziedziczności. Wprowadził pojęcie regresji („powrotu do średniej") i korelacji, przenosząc metodę z błędów pomiarowych na związki między zmiennymi. (Był też propagatorem eugeniki — idei dziś odrzuconej.)

Karl Pearson (1857–1936) — brytyjski matematyk i statystyk, który sformalizował współczynnik korelacji, test $\chi^2$ i znaczną część aparatu statystyki matematycznej; założyciel pierwszej katedry statystyki (University College London) i czasopisma Biometrika.

Ronald Aylmer Fisher (1890–1962) — brytyjski statystyk i genetyk, jeden z twórców współczesnej statystyki. Wprowadził metodę największej wiarygodności, analizę wariancji, planowanie eksperymentu i pojęcie istotności statystycznej, ujmując MNK jako estymator największej wiarygodności przy błędach normalnych.

Aneks D. Teksty oryginalne i źródła

Zestawienie pierwotnych prac, w których metoda powstała, wraz z miejscem kluczowych fragmentów i wskazaniem przekładów. Cytaty w częściach II i III pochodzą wprost z tych wydań.

Legendre (1805). A. M. Legendre, Nouvelles méthodes pour la détermination des orbites des comètes, Paris: Firmin Didot, 1805. Pierwsza publikacja metody znajduje się w dodatku Sur la méthode des moindres carrés (s. 72–80). To stamtąd pochodzi słynne zdanie „…rendre minimum la somme des carrés des erreurs" (część II).

Gauss (1809). C. F. Gauss, Theoria Motus Corporum Coelestium in Sectionibus Conicis Solem Ambientium, Hamburg: Perthes & Besser, 1809. Łaciński oryginał; rozważania o najmniejszych kwadratach mieszczą się w Księdze II, dziale 3, artykuły 175–186. Standardowy przekład angielski: C. H. Davis, Theory of the Motion of the Heavenly Bodies Moving about the Sun in Conic Sections, Boston: Little, Brown, 1857 — z niego pochodzi cytat o pierwszeństwie (część III).

Gauss (1821–1823). C. F. Gauss, Theoria combinationis observationum erroribus minimis obnoxiae, Getynga (część pierwsza 1821, druga 1823) — tu znajduje się dojrzałe uzasadnienie wariancyjne (twierdzenie Gaussa-Markowa). Przekład angielski z komentarzem: G. W. Stewart, Theory of the Combination of Observations Least Subject to Errors, SIAM, 1995.

Laplace (1812). P. S. Laplace, Théorie analytique des probabilités, Paris: Courcier, 1812 — probabilistyczne ujęcie błędów i centralne twierdzenie graniczne.

Przekłady i opracowania. Angielskie przekłady wielu prac historycznych (m.in. Legendre’a, Laplace’a) udostępnił Richard J. Pulskamp (probabilityandfinance.com). Klasyczne ujęcie historyczne: S. M. Stigler, The History of Statistics: The Measurement of Uncertainty before 1900, Harvard University Press, 1986.

Aneks E. Bibliografia i dalsze lektury

Poniżej współczesne opracowania i klasyczne prace źródłowe, uzupełniające teksty pierwotne z aneksu D. Dla pełnej historii metody niezastąpione pozostaje dzieło Stiglera (1986), cytowane już wyżej.

Współczesne podręczniki ekonometrii i statystyki.

• W. H. Greene, Econometric Analysis, Pearson, 8. wyd., 2018 — encyklopedyczny standard akademicki; pełny wykład macierzowego modelu liniowego, twierdzenia Gaussa-Markowa, FWL i rozszerzeń.
• J. M. Wooldridge, Introductory Econometrics: A Modern Approach, Cengage, 7. wyd., 2019 — najpopularniejszy wstęp; oraz tegoż Econometric Analysis of Cross Section and Panel Data, MIT Press, 2. wyd., 2010 (poziom zaawansowany).
• F. Hayashi, Econometrics, Princeton University Press, 2000 — rygorystyczne, oparte na metodzie momentów ujęcie MNK i GMM.
• R. Davidson, J. G. MacKinnon, Econometric Theory and Methods, Oxford University Press, 2004 — silny nacisk na geometrię rzutów i twierdzenie FWL.
• C. R. Rao, Linear Statistical Inference and Its Applications, Wiley, 2. wyd., 1973 — klasyka teorii estymacji liniowej.
• H. Theil, Principles of Econometrics, Wiley, 1971 — historyczny punkt odniesienia współczesnej ekonometrii.

Klasyczne prace o metodzie i jej historii.

• R. Frisch, F. V. Waugh, „Partial Time Regressions as Compared with Individual Trends", Econometrica 1(4), 1933, s. 387–401 — źródło twierdzenia FWL.
• M. C. Lovell, „Seasonal Adjustment of Economic Time Series and Multiple Regression Analysis", Journal of the American Statistical Association 58(304), 1963, s. 993–1010 — uogólnienie i nazwa twierdzenia.
• A. C. Aitken, „On Least Squares and Linear Combination of Observations", Proceedings of the Royal Society of Edinburgh 55, 1935, s. 42–48 — uogólniona metoda najmniejszych kwadratów (GLS).
• R. L. Plackett, „The Discovery of the Method of Least Squares", Biometrika 59(2), 1972, s. 239–251 — szczegółowa rekonstrukcja sporu o pierwszeństwo.
• S. M. Stigler, „Gauss and the Invention of Least Squares", The Annals of Statistics 9(3), 1981, s. 465–474 — analiza twierdzenia Gaussa o pierwszeństwie z 1809 r.

Kalendarium

Rok	Wydarzenie
1722	Cotes — ważona średnia obserwacji (pośmiertnie)
1736–44	Wyprawy do Laponii i Peru — pomiar kształtu Ziemi
~1750	Tobias Mayer — metoda uśredniania (libracja Księżyca)
1755–60	Bošković — minimalizacja sumy wartości bezwzględnych ($L_1$)
1795	Gauss zaczyna (wg własnych słów) używać MNK
1801	Piazzi odkrywa Ceres; Gauss wyznacza orbitę; planeta odnaleziona
1805	Legendre — pierwsza publikacja metody najmniejszych kwadratów
1808	Robert Adrain — niezależne wyprowadzenie prawa błędów
1809	Gauss, Theoria Motus — wyprowadzenie z postulatu średniej
1810	Laplace — centralne twierdzenie graniczne
1823	Gauss, Theoria combinationis — twierdzenie Gaussa-Markowa
1885	Galton — „regresja", przeniesienie metody na związki zmiennych
1896	Pearson — współczynnik korelacji
1922	Fisher — metoda największej wiarygodności

Dziedzictwo

Z jednego astronomicznego problemu wyrosło narzędzie leżące u podstaw regresji liniowej i ekonometrii, teorii estymacji, uczenia maszynowego (minimalizacja błędu kwadratowego) oraz numerycznej algebry liniowej (eliminacja Gaussa, rozkłady QR i SVD). A rozkład wyprowadzony przez Gaussa „przy okazji" nosi dziś jego imię i jest najważniejszym rozkładem statystyki — między innymi dzięki centralnemu twierdzeniu granicznemu.

---

Powiązane: Metoda najmniejszych kwadratów — wyprowadzenie krok po kroku · Założenia KMNK i twierdzenie Gaussa-Markowa · Rozkład normalny · Centralne twierdzenie graniczne