Eddig trendszámítás/regresszió témakörében két blogbejegyzésem született, amelyeket alább olvashattok. Mindkét bejegyzés a matematikai alapok mély megértésére koncentrál.
A többváltozós regressziónál nem egyetlen független (magyarázó) változó van, hanem több. Ez a gyakorlatban annyit tesz, hogy több faktor írja le a független és függő változók közötti kapcsolatot. A többváltozós lineáris regresszió egyenlete a következőképpen írható fel:
Ahol
a tengelymetszet, a többi
pedig az n-edik dimenzióhoz tartozó meredekséget írja le.
Fontos, hogy a többváltozós regresszió során többdimenziós függvényeket állítunk elő, amelyek nem keverendők össze a polinomiális trenddel/regresszióval, ami alapvetően kétdimenziós!
Az
a hibatag, tehát a modellünk
-al fog eltérni a megfigyelt értékektől.
Lényeges még az is, hogy a trendszámítás és a regressziószámítás matematikailag ekvivalensek, ezért ami az egyikre megáll, a legtöbb esetben a másikra is igaz.
Paraméterek kiszámítása
Ahhoz, hogy megkapjuk a paramétereket, tehát az -ekhez tartozó meredekségeket, a következő mátrixműveleteket szükséges elvégezni:
ahol:
-
: az együtthatóvektor (a becsült paraméterek vektora), tehát a számítás eredménye:
: a magyarázó (független) változók mátrixa, az első oszlopban a tengelymetszethez szükséges egyesekkel
: az eredményváltozó (függő változó) megfigyelt értékeinek vektora
X nyilván egy mátrix lesz, ha több mint kettő dimenzióban gondolkodunk, Y pedig egy vektor.
Szeretnék visszautalni az első trendszámítással kapcsolatos cikkemre, amelyben kifejtem, hogy az úgynevezett normálegyenletek az SSE (Sum of Squared Errors) deriváltjának nullára rendezésével kaphatók vissza. Ez igaz a fenti képletre is, csak ez a többdimenziós eset.
Az sajnos hátrányban van, aki nem ismeri a lineáris algebra alapjait. Tervezek írni ezzel kapcsolatban egy cikket és linkelni ehhez a bejegyzéshez.
Fontos, hogy ki kell egészítenünk az X-mátrixot egy csupa egyesekből álló oszloppal, hogy megkapjuk a tengelymetszetet. A végső mátrix valahogy így fog kinézni:
Centrálás a nagyobb numerikus pontosság érdekében
Ha az adatsorban nagyon nagy értékek is szerepelnek, akkor a transzponált és az eredeti mátrix szorzása miatt létrejövő négyzetszámok akár túl is léphetik a 64 bites lebegőpontos aritmetika korlátait. Tehát elképzelhető, hogy túl nagy számot kapunk és úgynevezett túlcsordulás következik be.
A centrálás semmi más nem jelent, csak azt, hogy a függő és független változók összes értékéből levonjuk az adott változó átlagát. Így minden egyes érték a nulla körül fog szóródni, és jóval kisebb számokkal dolgozunk, viszont pontosan ugyanazt az eredményt kapjuk vissza. Ebben az esetben az egyesek eltűnnek az első sorból, és a tengelymetszetet csak trükkösen kaphatjuk vissza.
A csupa egyesekből álló oszlop elrejtését és a centrálást követően az X-mátrix így fog kinézni:
Az Y mátrix pedig így:
Mivel ilyenkor hiányozni fog a tengelymetszet változó, ezért ezt külön szükséges kiszámolni, melyre a képlet a következő:
A Python kód eddig
import pandas as pd
import numpy as np
def linear_regression(x_df:pd.DataFrame, y_df:pd.DataFrame):
"""
Átalakítjuk az x és az y értékek DataFrame-jeit NumPy tömbökké.
Ügyelünk a típushelyességre. (dtype=np.float64)
"""
x = x_df.to_numpy(dtype=np.float64)
y = y_df.to_numpy(dtype=np.float64)
"""
Az axis=0 miatt ezek az oszlopok átlagai,
tehát annyi átlagot kapunk, ahány oszlop
van és nem annyit, hány sor.
"""
x_means = x.mean(axis=0)
"""
Kiszámoljuk az y átlagát is.
"""
y_mean = y.mean()
"""
A NumPy array okos, így mindig az oszlopoknak
megfelelő értéket vonja ki az adott oszlop
n-edik sorából. Teszi mindezt automatikusan
anélkül, hogy az indexekkel kellene törődnünk,
vagy ciklusban egyenként számításokat végeznünk.
(A nulladik x_means értéket az összes nulladik
oszlopbeli x_0_i-ből levonja.)
"""
x_centered = x - x_means
"""
Az y egy vektor (tehát egyetlen oszlop).
Minden egyes értékből le lesz vonva az y átlag.
"""
y_centered = y - y_mean
"""
Itt már csak követjük a képletet. A slopes
a centrálással, és a csupa egyes oszlop nélnül
ténylegesen csak az x_n szorzókat (meredekségeket)
adja vissza, a tengelymetszetet nem!
Ne feledjük, mátrix szorzáshoz a @ jelet,
Hadamard szorzáshoz pedig a * jelet alkalmazzuk!
A dot product (skaláris szorzat) kiszámításához
az np.dot() metódust alkalmazzuk.
"""
xt = x_centered.T #transzponálás
xt_x_inv = np.linalg.inv(xt @ x_centered) #inverz
xty = xt @ y_centered #transzponált és y szorzata
slopes = xt_x_inv @ xty #béták
"""
Az intercept kiszámítása a képlet alapján.
"""
intercept = y_mean - np.sum(slopes * x_means)
return slopes, intercept
A model helyességének mérése
Nézzünk meg három mutatót, amelyek fontosak a modell helyessége és magyarázóereje mérésének szempontjából.
SSE
Az SSE (Sum of Squared Errors) azt mondja meg, hogy mennyi a kumulált négyzetes eltérés a valós értékek valamint a modell által számított értékek között. A képlete a következő:
SST
Az SST (Total Sum of Squares) kiszámítja, hogy mekkora a négyzetes eltérés a tényleges eredményváltozó értékek és az értékek átlaga között. Ez lényegében az Y-változó varianciája (oké, eltérés-négyzetösszege). Képlet:
SSR
Az SSR (Sum of Squares due to Regression) megadja, hogy a modell által számított értékek mennyire térnek el az Y-átlagától. Képlet:
Mindhárom mutató négyzetes eltérések összegét jelenti! Nézzük, hogy mire tudjuk ezeket használni!
R-négyzet
Az (determinációs tényező) annak az indikátora, hogy a független (magyarázó) változók a függő (eredmény) változó varianciájának (szóródásának) mekkora hányadát (hány százalékát) magyarázzák meg. A képlete a következő:
A mutató megértése rendkívül egyszerű, ha egy kicsit belegondolsz az előzőek alapján, hogy mi az SSR és az SST. Az SSR a modell által számított értékek eltérése az átlagtól, míg az SST a tényleges értékek eltérése az átlagtól. Hogyha a modell tökéletes, akkor a kettő érték egyforma, tehát ugyanúgy a függő változó varianciáját adják vissza, így a végeredmény egy lesz. Ha az értékek nagy mértékben eltérnek, akkor feltehetőleg egy kisebb számot fogunk kapni.
Korrigált R-négyzet (Adjusted R-squared)
Az remek mutató, de van egy kisebb szépséghibája: ha újabb és újabb magyarázó változókat adsz hozzá a modellhez, az értéke szinte soha nem fog csökkenni, sőt az esetek túlnyomó többségében inkább nőni fog, még akkor is, ha az új változónak égvilágon semmi köze a függő változóhoz (pl. a részvényárakat próbálod megjósolni a budapesti esős napok számával).
Miért van ez? Mert a legkisebb négyzetek módszere (OLS) arra törekszik, hogy a lehető legkisebbre csökkentse a hibaösszeget (SSE). Ha adsz a modellnek egy újabb magyarázó változót, az OLS találni fog benne valamilyen minimális, véletlenszerű egybeesést a függő változóval, még akkor is, ha az a változó valójában teljesen lényegtelen. A modell ezt a véletlen ingadozást (a zajt) is felhasználja arra, hogy egy picit közelebb tornázza magát a pontokhoz. Így az nőni fog, de valójában nem a valódi összefüggést értettük meg jobban, hanem a mintánk egyedi torzulásait tanultuk meg. Ezt hívjuk túlillesztésnek (overfitting).
Ezért van szükségünk a korrigált R-négyzetre ( ), amely bünteti a modellt, ha feleslegesen sok magyarázó változót pakolunk bele.
A képlet
Vagy ha a négyzetes eltérés-összegekkel írjuk fel:
Mi mit jelent a képletben?
- : A determinációs tényező (sima R-négyzet).
- : A megfigyelések (mintaelemek) száma.
- : A magyarázó (független) változók száma (a tengelymetszet nélkül).
- : A teljes eltérés-négyzetösszeg (SST) szabadságfoka.
- : A maradvány eltérés-négyzetösszeg (SSE) szabadságfoka.
Miért működik a büntetés?
A képletben lévő szorzó értéke mindig nagyobb, mint 1. Ha újabb magyarázó változót ( ) adunk a modellhez, a nevező csökken, így a tört értéke nő. Ez a büntetőtag megnöveli a értékét, amit a végén kivonunk 1-ből.
A korrigált csak akkor fog növekedni, ha az új változó által hozott magyarázóerő-javulás nagyobb, mint a plusz változó bevonása miatt járó büntetés. Ha a változó nem hoz érdemi javulást, a korrigált értéke csökkenni fog.
F-mutató
Az értékek szignifikanciáját is tudjuk tesztelni hipotézisvizsgálat segítségével. Az F-érték kiszámítása itt a következőképpen alakul:
ahol:
- : a függő változó értékeinek száma (mintaelemszám)
- : a független (magyarázó) változók száma (az mátrix oszlopainak száma a tengelymetszet nélkül)
- : a regressziós átlagos négyzetösszeg ( )
- : a hiba átlagos négyzetösszege ( )
A regressziós modell tesztelése mindig jobb oldali próba, tehát a kritikus tartomány jobb oldali. A hipotézisvizsgálatban a
nullhipotézis az, hogy minden magyarázóváltozó inszignifikáns, tehát nincs magyarázóereje, a
ellenponthipotézis pedig az, hogy legalább egynek van szignifikáns magyarázóereje.
Nade, hogy jön ide az MSR/MSE mutató? Létezik egy ilyen összefüggés:
Ha MSR értéke nagy, MSE értéke pedig kicsi, akkor a hányadosuk egy nagy szám lesz. Ez szignifikáns kapcsolatra utalhat. Ha MSR értéke kicsi és MSE értéke is viszonylag alacsony, akkor lehetséges, hogy a kisebb reziduális szórás ellenére sem szignifikáns a modell, mert kevés varianciát magyaráz. Fontos megérteni, hogy a kettő aránya a lényeges, és nem az abszolútértékük!
A kritikus F-eloszlás értékének kikeresése már nem egy nagy szám, hiszen k és (n - k - 1) szabadságfoknál lévő F-értéket kell kikeresni a megfelelő táblázatból, vagy használhatjuk akár a Python-vonatkozó könyvtárát is. Mutatom is hogyan:
from scipy import stats
n = 50 # Megfigyelések száma
k = 4 # Magyarázó változók száma
alpha = 0.05 # Szignifikanciaszint
"""
Az F-eloszlás inverz eloszlásfüggvénye (PPF) alapján
meghatározzuk a kritikus F-értéket.
A stats.f.ppf() paraméterei:
1. paraméter: Kvantilis valószínűségi szintje (1 - alpha).
Ez nem a szignifikanciaszint, hanem a
kumulatív valószínűség, amelyhez a kritikus
értéket keressük.
2. paraméter: Számláló szabadságfoka (k).
3. paraméter: Nevező szabadságfoka (n - k - 1).
Az alpha a szignifikanciaszint, amely a jobboldali
kritikus tartomány valószínűségét határozza meg.
Az 1 - alpha a kumulatív valószínűségi szint, amelyhez
a kritikus F-értéket meghatározzuk.
A ppf a Percent Point Function rövidítése.
"""
f_critical = stats.f.ppf(
1 - alpha,
k,
n - k - 1
)
Fontos, hogyha ennél az értéknél kisebbet ad a saját F-statisztika értékünk, akkor elfogadtuk a H0 nullpont hipotézist, miszerint a modell magyarázóereje NEM szignifikáns.
A saját Python OLS (Ordinary Least Squares) kódom
A dataclass ne zavarjon meg benneteket! Ez csak azt határozza meg, hogy egy regressziószámításnak pontosan milyen kimeneti paraméterekkel kell rendelkeznie. Azért előnyös, mert a Python kódod ennek hatására automatikus javaslatokkal fog élni, tehát ha beírod, hogy result.s, akkor felajánlja az összes s-kezdőbetűs mutatót. Az __str__ arra jó, hogyha printeljük az eredményt, akkor formázottan, táblázatban jelennek meg az értékek. Más programozási nyelveknél ez a tostring metódus felülírása lenne.
from dataclasses import dataclass
import numpy as np
import pandas as pd
from scipy import stats
@dataclass
class RegressionResult:
slopes: np.ndarray
intercept: float
sst: float
ssr: float
sse: float
r_squared: float
r_adjusted: float
f: float
h0_accepted: bool
significant: bool
def __str__(self) -> str:
lines = [
f" slopes │ {np.round(self.slopes, 4)} ",
f" intercept │ {self.intercept:.4f} ",
f" SST │ {self.sst:.4f} ",
f" SSR │ {self.ssr:.4f} ",
f" SSE │ {self.sse:.4f} ",
f" R^2 │ {self.r_squared:.4f} ",
f" Adj. R^2 │ {self.r_adjusted:.4f} ",
f" F-stat │ {self.f:.4f} ",
f" H0 accepted │ {str(self.h0_accepted)} ",
f" Significant │ {str(self.significant)} ",
]
return "\n".join(lines)
def linear_regression(
x_df: pd.DataFrame, y_df: pd.DataFrame, alpha: float = 0.05
) -> RegressionResult:
x = x_df.to_numpy(dtype=np.float64)
#A ravel() egydimenziós tömbbé (vektorrá) lapítja ki a DataFrame-ből kinyert 2D-s NumPy tömböt.
y = y_df.to_numpy(dtype=np.float64).ravel()
n = x.shape[0]
df_1 = x.shape[1]
df_2 = n - df_1 - 1
y_mean = y.mean()
y_centered = y - y_mean
x_means = x.mean(axis=0)
x_centered = x - x_means
xt = x_centered.T
xtx_inv = np.linalg.inv(xt @ x_centered)
xty = xt @ y_centered
slopes = xtx_inv @ xty
intercept = y_mean - np.sum(slopes * x_means)
y_hat = intercept + x @ slopes
sst = np.sum((y - y_mean) ** 2)
ssr = np.sum((y_hat - y_mean) ** 2)
sse = np.sum((y - y_hat) ** 2)
r_squared = ssr / sst
r_adjusted = 1 - ((1 - r_squared) * (n - 1) / df_2)
f = (ssr / df_1) / (sse / df_2)
f_critical = stats.f.ppf(1 - alpha, df_1, df_2)
h0_accepted = bool(f < f_critical)
significant = not h0_accepted
return RegressionResult(
slopes=slopes,
intercept=intercept,
ssr=ssr,
sse=sse,
sst=sst,
f=f,
r_squared=r_squared,
r_adjusted=r_adjusted,
h0_accepted=h0_accepted,
significant=significant,
)
