Lineární regrese: nejmenší čtverce geometricky i maticově

Nechť p ∈ ℕ, p ≥ 1, je počet příznaků a w = (w₀, w₁, …, w_p) ∈ ℝ^p+1 je vektor parametrů. Lineární model s parametry w je funkce f_w : ℝ^p → ℝ,

f_w(x₁, …, x_p) = w₀ + ∑_j=1^p w_jx_j.

Číslo w₀ je úsek, čísla w_j pro j = 1, …, p jsou koeficienty. Hodnota ŷ = f_w(x) je předpověď pro vektor příznaků x ∈ ℝ^p.

Nechť jsou dále dána data — n ≥ 1 pozorování (x_i, y_i), i = 1, …, n, kde x_i = (x_i1, …, x_ip) ∈ ℝ^p jsou příznaky a y_i ∈ ℝ naměřená hodnota. Reziduum i-tého pozorování je číslo

e_i = y_i − ŷ_i,   kde ŷ_i = f_w(x_i).

Lineární regrese je úloha najít z dat parametry w, pro které jsou rezidua v daném smyslu co nejmenší (nejčastěji metodou nejmenších čtverců, viz níže). Model tedy předpovídá úsek plus vážený součet příznaků a reziduum je, o kolik se skutečnost od předpovědi liší.

Nechť (x_i, y_i), i = 1, …, n, jsou data lineárního modelu s p příznaky (x_i ∈ ℝ^p, y_i ∈ ℝ). Součet čtverců reziduí je funkce SSE : ℝ^p+1 → ℝ,

SSE(w) = ∑_i=1ⁿ (y_i − f_w(x_i))².

Metoda nejmenších čtverců volí parametry ŵ ∈ ℝ^p+1 jako bod minima této funkce:

SSE(ŵ) ≤ SSE(w)  pro všechna w ∈ ℝ^p+1.

Pro jeden příznak (p = 1, x_i ∈ ℝ) označme průměry x̄ = 1/n ∑_i=1ⁿ x_i a ȳ = 1/n ∑_i=1ⁿ y_i. Nejsou-li všechna x_i stejná, má úloha právě jedno řešení (důkaz ve větě níže; dál ho píšeme bez stříšky, w₀, w₁):

w₁ = ∑_i=1ⁿ (x_i − x̄)(y_i − ȳ)/∑_i=1ⁿ (x_i − x̄)²,   w₀ = ȳ − w₁x̄.

Sklon je tedy kovariance x a y dělená rozptylem x a přímka prochází těžištěm (x̄, ȳ).

Nechť n ≥ 2, (x_i, y_i) ∈ ℝ² pro i = 1, …, n, a nechť čísla x₁, …, x_n nejsou všechna stejná. Nechť x̄, ȳ jsou jejich průměry. Pak funkce

SSE(w₀, w₁) = ∑_i=1ⁿ (y_i − w₀ − w₁x_i)²

má na ℝ² právě jeden bod minima, a to w₁ = ∑ (x_i − x̄)(y_i − ȳ)/∑ (x_i − x̄)², w₀ = ȳ − w₁x̄.

Důkaz. SSE je součet čtverců afinních funkcí parametrů, tedy konvexní kvadratická funkce. Její Hessova matice je 2 · [n, ∑ x_i; ∑ x_i, ∑ x_i²] s determinantem 4(n ∑ x_i² − (∑ x_i)²) = 4n ∑ (x_i − x̄)² > 0 (protože x_i nejsou všechna stejná) a kladnou diagonálou, takže je pozitivně definitní a funkce je ostře konvexní. Bod minima je proto jediný a je to právě bod, kde jsou obě parciální derivace nulové:

První rovnice je přesně podmínka „přímka prochází těžištěm“. ∎

Průměry: x̄ = 525/8 = 65,625 m², ȳ = 53,4/8 = 6,675 mil. Kč. Odchylky a jejich součiny (barvy podle rozboru vzorce pro sklon výše):

Všech osm součinů je kladných: každý byt je buď větší a dražší, nebo menší a levnější než průměr — proto tak jasný kladný sklon.

Nechť n, p ∈ ℕ, n, p ≥ 1, a nechť (x_i, y_i), i = 1, …, n, jsou data, x_i = (x_i1, …, x_ip) ∈ ℝ^p, y_i ∈ ℝ. Označme y = (y₁, …, y_n)^T ∈ ℝⁿ. Matice plánu je matice X ∈ ℝ^n×(p+1) se sloupci číslovanými 0, …, p a prvky

X_i0 = 1,   X_ij = x_ij   (i = 1, …, n;  j = 1, …, p),

tj. řádek na každé pozorování, první sloupec samé jedničky (pro úsek) a další sloupce příznaky. Pro w = (w₀, …, w_p)^T ∈ ℝ^p+1 je vektor předpovědí všech pozorování ŷ = Xw ∈ ℝⁿ a SSE(w) = ‖y − Xw‖².

Normální rovnice jsou soustava p + 1 lineárních rovnic pro neznámý vektor w ∈ ℝ^p+1:

X^TXw = X^Ty.

Podle věty níže je w bodem minima SSE, právě když řeší normální rovnice; je-li matice X^TX regulární, je řešení jediné, w = (X^TX)⁻¹X^Ty.

Nechť X ∈ ℝ^n×(p+1), y ∈ ℝⁿ a S : ℝ^p+1 → ℝ, S(w) = ‖y − Xw‖². Pak

1. w ∈ ℝ^p+1 je bodem minima S, právě když X^TXw = X^Ty;
2. jsou-li sloupce X lineárně nezávislé, je X^TX regulární a bod minima je jediný, w = (X^TX)⁻¹X^Ty.

Důkaz. Rozepíšeme čtverec normy: ‖y − Xw‖² = y^Ty − 2w^TX^Ty + w^TX^TXw. Gradient podle w je −2X^Ty + 2X^TXw; je nulový, právě když X^TXw = X^Ty. Hessova matice 2X^TX je pozitivně semidefinitní (v^TX^TXv = ‖Xv‖² ≥ 0 pro každé v), takže S je konvexní a diferencovatelná — a pro takovou funkci je bod minimem, právě když je v něm gradient nulový. To je bod 1. Jsou-li sloupce X lineárně nezávislé, je Xv ≠ 0 pro každé v ≠ 0, tedy ‖Xv‖² > 0; matice X^TX je pozitivně definitní, a proto regulární, a soustava má jediné řešení. ∎

Tři měření x = (1, 2, 3), y = (1, 3, 2), model ŷ = w₀ + w₁x (budou se hodit i v příští kapitole):

Barvy odpovídají členům vzorce w = (X^TX)⁻¹X^Ty výše.

V matici X^TX je 3 počet bodů, 6 součet x a 14 součet x²; v X^Ty součet y a součet x·y. Data se opravdu „sesypala“ do pěti čísel.

Byt o 70 m², 5 km od centra, v domě starém 30 let. S koeficienty zaokrouhlenými jako v tabulce:

Barvy podle rozboru lineárního modelu v první kapitole: úsek, koeficienty, hodnoty příznaků a předpověď. Záporné koeficienty jsou zapsané odečtením.

S nezaokrouhlenými koeficienty vyjde 6,00 mil. Kč — u plochy se zaokrouhlení násobí sedmdesáti, proto ten rozdíl.

Nechť X ∈ ℝ^n×(p+1) je matice plánu a y ∈ ℝⁿ. Prostor sloupců matice X je množina všech možných vektorů předpovědí

Col(X) = {Xw : w ∈ ℝ^p+1} ⊆ ℝⁿ;

je to lineární podprostor ℝⁿ rozepnutý sloupci X. Kolmý průmět vektoru y do Col(X) je vektor ŷ ∈ ℝⁿ splňující

1. ŷ ∈ Col(X),
2. v^T(y − ŷ) = 0 pro každé v ∈ Col(X).

Takový vektor existuje a je jediný a je to bod Col(X) nejbližší k y: ‖y − ŷ‖ ≤ ‖y − v‖ pro všechna v ∈ Col(X). Podmínku 2 stačí ověřit pro sloupce X, takže pro vektor reziduí e = y − ŷ zní X^Te = 0; po dosazení ŷ = Xw to jsou přesně normální rovnice X^TXw = X^Ty. Předpověď nejmenších čtverců je tedy kolmý průmět y do prostoru sloupců a reziduum je kolmé na všechno, co model umí vyrobit.

Nechť n ≥ 2, y = (y₁, …, y_n) ∈ ℝⁿ jsou naměřené hodnoty, které nejsou všechny stejné, ȳ = 1/n ∑_i=1ⁿ y_i jejich průměr a ŷ = (ŷ₁, …, ŷ_n) ∈ ℝⁿ předpovědi nějakého modelu. Definujeme

1. celkový součet čtverců SST = ∑_i=1ⁿ (y_i − ȳ)² (je kladný, protože y_i nejsou všechna stejná),
2. vysvětlený součet čtverců SSR = ∑_i=1ⁿ (ŷ_i − ȳ)²,
3. reziduální součet čtverců SSE = ∑_i=1ⁿ (y_i − ŷ_i)².

Koeficient determinace předpovědí ŷ je číslo

R² = 1 − SSE/SST.

Je-li ŷ předpověď lineárního modelu s úsekem odhadnutého nejmenšími čtverci z těchto dat, platí SST = SSR + SSE (věta níže; celkový = vysvětlený + nevysvětlený), a tedy R² = SSR/SST ∈ [0, 1]. Pro jiné předpovědi (třeba na nových datech) rozklad platit nemusí a R² může být i záporné.

Nechť X ∈ ℝ^n×(p+1) je matice plánu, jejíž první sloupec je vektor jedniček 1 = (1, …, 1)^T ∈ ℝⁿ, y ∈ ℝⁿ, ŷ je kolmý průmět y do Col(X), e = y − ŷ a ȳ průměr složek y. Pak SST = SSR + SSE.

Důkaz. Označme ȳ1 vektor se všemi složkami rovnými ȳ. Pak y − ȳ1 = (ŷ − ȳ1) + e. Vektor ŷ − ȳ1 leží v Col(X) (ŷ v něm leží a 1 je sloupec X) a e je na Col(X) kolmý, takže (ŷ − ȳ1)^Te = 0 a

Bez úseku v modelu sloupec jedniček chybí, kolmost na něj neplatí a rozklad (i R² v téhle podobě) se rozpadne.

Tři body z minulé kapitoly: y = (1, 3, 2), ȳ = 2, ŷ = (1,5; 2; 2,5), e = (−0,5; 1; −0,5).

Osm bytů z kapitoly o nejmenších čtvercích: SSE = 3,807 a SST = ∑ (y_i − 6,675)² = 31,755, tedy R² = 1 − 3,807/31,755 ≈ 0,880. Přímka vysvětlí 88 % rozptylu cen osmi bytů. (Barvy podle rozboru vzorce pro R² výše.)

Nechť je lineární model s úsekem a p ≥ 1 příznaky odhadnut nejmenšími čtverci z n pozorování, kde n > p + 1, a nechť R² je jeho koeficient determinace na těchto (trénovacích) datech. Upravený koeficient determinace je číslo

R²_upr = 1 − (1 − R²) · n − 1/n − p − 1.

Protože n − 1/n − p − 1 > 1 a s rostoucím p roste, penalizuje počet příznaků p vzhledem k počtu pozorování n; vždy platí R²_upr ≤ R².

Šedesát bytů, n = 60 (barvy podle vzorce výše):

Oprava u 43 příznaků nafoukne nevysvětlený podíl skoro čtyřikrát, ale trénovací R² se přiblížilo k jedničce tolik, že to nestačí.

Nechť n ≥ 2 a u = (u₁, …, u_n), v = (v₁, …, v_n) ∈ ℝⁿ jsou hodnoty dvou příznaků na n pozorováních, z nichž ani jeden není konstantní (nejsou všechna u_i stejná ani všechna v_i stejná). Nechť ū = 1/n ∑_i=1ⁿ u_i a v̄ = 1/n ∑_i=1ⁿ v_i. Pearsonův korelační koeficient příznaků u, v je číslo

r = ∑_i=1ⁿ (u_i − ū)(v_i − v̄)/√∑_i=1ⁿ (u_i − ū)² · ∑_i=1ⁿ (v_i − v̄)².

Z Cauchyho–Schwarzovy nerovnosti plyne r ∈ [−1, 1]; |r| = 1, právě když všechny body (u_i, v_i) leží na jedné přímce, a r = 0, právě když je sklon přímky nejmenších čtverců nulový — mezi příznaky pak není žádný lineární vztah (jiný být může).

Čitatel je stejný jako u vzorce pro sklon; jmenovatel místo rozptýlenosti jednoho příznaku bere geometrický průměr obou, takže r nemá jednotky a je symetrický.

Nechť X ∈ ℝ^n×(p+1) je matice plánu se sloupci X_·0 = 1, X_·1, …, X_·p. Příznaky jsou (přesně) kolineární, právě když jsou sloupce X lineárně závislé, tj.

∃ c ∈ ℝ^p+1, c ≠ 0 :   Xc = c₀1 + c₁X_·1 + … + c_pX_·p = 0,

neboli některý sloupec je lineární kombinací ostatních (včetně sloupce jedniček). Pak je X^TX singulární a normální rovnice nemají jediné řešení. Pro j = 1, …, p označme R_j² koeficient determinace regrese příznaku X_·j (s úsekem) na ostatní příznaky. Příznaky jsou přibližně kolineární (multikolinearita), když je pro některé j číslo R_j² blízké 1 — j-tý příznak jde ostatními vyjádřit skoro přesně. Ostrá hranice neexistuje; v praxi se jako varovný signál bere inflační faktor 1/1 − R_j² nad 5 či 10. Matice X^TX je pak skoro singulární a odhady koeficientů mají velký rozptyl, i když předpovědi zůstávají dobré.

Nechť D = (z₁, …, z_n) jsou data o n ≥ 1 pozorováních. Převýběr (bootstrapový výběr) z D je posloupnost

D^* = (z_J1, …, z_Jn),

kde indexy J₁, …, J_n se losují nezávisle a každý rovnoměrně z {1, …, n} (každé pozorování má v každém tahu pravděpodobnost 1/n). Jde tedy o n tahů s vracením — některý byt se vytáhne dvakrát, jiný vůbec.

Nechť t je odhad (funkce dat, třeba koeficient w_plocha z nejmenších čtverců), B ∈ ℕ, B ≥ 2, a D^*₁, …, D^*_B nezávislé převýběry. Označme t_b hodnotu odhadu na D^*_b a m = 1/B ∑_b=1^B t_b. Bootstrapový odhad směrodatné odchylky odhadu t je

s_boot = √1/B ∑_b=1^B (t_b − m)².

Rozptyl odhadu přes mnoho převýběrů tedy napodobuje, jak by se odhad lišil, kdybychom měli jiná data ze stejného zdroje. Čísla výše vznikla takhle (B = 200):

for s in 1 … 200:
    D = 60 bytů vytažených náhodně s vracením z 60 původních
    w = nejmenší čtverce na D (plocha, vzdálenost, stáří[, pokoje])
    zapiš w_plocha (a w_pokoje)
směrodatná odchylka = √( průměr (w_plocha − průměr w_plocha)² )

Barvy podle vzorce pro rozptyl odhadu výše.

Vzorec tedy čeká zhruba 2,3× větší kolísání, převýběry naměřily 2,8× — 200 převýběrů ze 60 bytů je samo jen odhad a vzorec navíc předpokládá pevná data. Řád ale sedí: pokoje „ukradly“ 83 % rozptylu plochy.

Nechť X ∈ ℝ^n×(p+1) je matice plánu s lineárně nezávislými sloupci (prvním sloupcem jedniček) a y ∈ ℝⁿ. Projekční matice (anglicky hat matrix) je

H = X(X^TX)⁻¹X^T ∈ ℝ^n×n.

Předpovědi nejmenších čtverců jsou ŷ = Xw = Hy, tedy ŷ_i = ∑_k=1ⁿ h_iky_k — lineární kombinace naměřených hodnot. Páka i-tého pozorování (i = 1, …, n) je diagonální prvek h_ii: váha, s jakou vlastní y_i vstupuje do vlastní předpovědi ŷ_i. Platí

1. 1/n ≤ h_ii ≤ 1 pro každé i (dolní mez díky sloupci jedniček),
2. ∑_i=1ⁿ h_ii = p + 1, tj. součet pák je počet parametrů (u přímky 2).

Pro přímku (p = 1, x_i ∈ ℝ nejsou všechna stejná) vychází

h_ii = 1/n + (x_i − x̄)²/∑_k=1ⁿ (x_k − x̄)².

Devět bytů včetně překlepu: x̄ = 645/9 ≈ 71,67 m², ∑ (x_k − x̄)² = 5918 (barvy podle rozboru výše).

Polovinu předpovědi pro 120 m² určuje sám překlep. Běžné orientační pravidlo (Hoaglin a Welsch, 1978) radí zkontrolovat body s pákou nad dvojnásobkem průměru — tady nad 0,44.

Nechť δ ∈ ℝ, δ > 0, je práh. Huberova ztráta s prahem δ je funkce L_δ : ℝ → ℝ, která reziduu e ∈ ℝ přiřadí

L_δ(e) = e²/2,   je-li |e| ≤ δ,

L_δ(e) = δ · (|e| − δ/2),   je-li |e| > δ.

Odhad s Huberovou ztrátou volí parametry minimalizující ∑_i L_δ(e_i) místo SSE. Malá rezidua trestá jako čtverce (hladce, s minimem jako u nejmenších čtverců), velká jen lineárně, takže jeden překlep o 5,5 milionu nepřetáhne přímku k sobě. V bodě |e| = δ obě větve navazují hodnotou i derivací.