Lineær regresjon og tapsfunksjoner

En lineær modell beregner en vektet sum av featurene og legger til et konstantledd. Det er ƒlineær regresjonsmodell:

Vektene $w_j$ og konstantleddet $b$ er hver en modellparameter: tall som treningen tilpasser fra data. Et biasledd $b$ er prediksjonen når alle featurene er null, og det lar linjen flyttes opp og ned. Uten $b$ må linjen gå gjennom origo.

Med én feature er modellen $\hat{y}=wx+b$, og $w$ er stigningen: øker $x$ med én enhet, øker prediksjonen med $w$. Med flere features er $w_j$ endringen i $\hat{y}$ når $x_j$ øker med én og alle de andre featurene holdes faste. Det er en egenskap ved modellen og ikke nødvendigvis en årsakssammenheng. Korrelerte features og variabler som mangler i datasettet, kan gjøre en vekt vanskelig å tolke.

For hvert eksempel er en residual avviket $r_i=y_i-{\hat{y}}_i$. En positiv residual betyr at modellen predikerte for lavt, og en negativ at den predikerte for høyt. Residualene er de loddrette avstandene mellom punktene og linjen i figuren.

Et residualplott, med residual mot prediksjon, avslører ofte det modellen ikke fanger: en bue tyder på en ikke-lineær sammenheng, og en vifte som åpner seg tyder på at støyen vokser med nivået.

import numpy as np

# Features: areal (10 m²), antall soverom, avstand til sentrum (km)
X = np.array([[6.0, 2.0, 3.0],
              [8.5, 3.0, 7.5],
              [4.0, 1.0, 1.0]])
y = np.array([4.1, 4.6, 3.3])            # salgspris i millioner kroner
w = np.array([0.35, 0.20, -0.12])
b = 1.0  # @param 0..3 0.05

y_hat = X @ w + b
print(f"prediksjoner:   {np.round(y_hat, 3)}")
print(f"residualer y-ŷ: {np.round(y - y_hat, 3)}")

x_ny = X[0].copy()
x_ny[0] += 1                              # én enhet (10 m²) mer areal, alt annet likt
print(f"endring i ŷ for 10 m² mer: {(x_ny @ w + b) - y_hat[0]:.3f}")

En tapsfunksjon samler alle residualene i ett tall som treningen kan gjøre minst mulig. De to vanligste er ƒmean squared error og ƒmean absolute error:

Forskjellen ligger i hvordan én stor feil teller. I MSE teller en feil på 10 hundre ganger så mye som en feil på 1, i MAE bare ti ganger. MSE straffer derfor store avvik hardt og lar noen få ekstreme punkter dra modellen mot seg, mens MAE er mer robust mot uteliggere. Roten av MSE, RMSE, har samme enhet som målvariabelen og er lettere å tolke.

En nyttig måte å se forskjellen på er å spørre hvilket konstant tall som er best å predikere for alle punktene. For MSE er svaret gjennomsnittet av $y$, og for MAE er det medianen. Gjennomsnittet flytter seg mye når én verdi blir ekstrem, medianen nesten ikke.

MSE har også en praktisk fordel: den er glatt og har en enkel derivert overalt, mens MAE har en knekk der residualen er 0. Det gjør MSE enkel å optimere. Tap og evalueringsmetrikk trenger likevel ikke være like. Man kan trene med MSE og rapportere både MAE og RMSE, fordi de belyser ulike sider av feilene.

import numpy as np

X = np.array([[1.0], [2.0], [3.0], [4.0], [5.0]])
y = np.array([2.8, 5.1, 7.2, 8.9, 11.1])
w = np.array([2.0])
b = 1.0
malefeil = 0.0  # @param 0..10 0.5

def tap(y):
    r = y - (X @ w + b)
    return np.mean(r ** 2), np.mean(np.abs(r))

y_feil = y.copy()
y_feil[-1] += malefeil                    # den siste målingen blir feil
mse0, mae0 = tap(y)
mse, mae = tap(y_feil)
print(f"uten feil: MSE = {mse0:.3f}, MAE = {mae0:.3f}")
print(f"med feil:  MSE = {mse:.3f}, MAE = {mae:.3f}, RMSE = {np.sqrt(mse):.3f}")
print(f"MSE ble {mse / mse0:.1f} ganger større, MAE {mae / mae0:.1f} ganger")

MSE har enheten til målvariabelen i andre, og tallet alene sier lite: er en MSE på 4 bra eller dårlig? Det avhenger av hvor mye $y$ varierer. ƒForklart varians løser dette ved å sammenligne modellen med den enkleste tenkelige baselinen, en modell som alltid predikerer snittet $\bar{y}$:

$R^2=1$ betyr perfekte prediksjoner, og $R^2=0$ betyr at modellen ikke er bedre enn snittet. $R^2$ kan også bli negativ på nye data, og det betyr at modellen er dårligere enn å gjette snittet. Det skjer typisk når en overtilpasset modell møter testdata, eller når prediksjonene har en systematisk forskyvning. Til tross for navnet er $R^2$ altså ikke nødvendigvis et kvadrat, og den har ingen nedre grense.

Tolk tallet med forsiktighet. En høy $R^2$ på treningsdata sier lite om generalisering, så den skal måles på validerings- eller testdata. Den sier heller ikke at modellantakelsene holder, eller at feilene er små nok for beslutningen. Og den avhenger av hvor spredt $y$ er i datasettet: samme modell kan få ulik $R^2$ på to testsett med like store feil.

I scikit-learn returnerer model.score(X, y) for regresjonsmodeller nettopp $R^2$, og r2_score(y_true, y_pred) regner det ut fra prediksjoner. Rekkefølgen på argumentene betyr noe, fordi $S{S}_{\text{tot}}$ regnes ut fra de sanne verdiene.

import numpy as np

y = np.array([12.0, 15.0, 11.0, 18.0, 14.0])          # sanne verdier
forskyvning = 0.0  # @param 0..10 0.5
modell = np.array([12.5, 14.0, 11.5, 17.0, 15.0]) + forskyvning   # god modell, ev. forskjøvet
baseline = np.full_like(y, y.mean())                   # predikerer alltid snittet
daarlig = np.array([17.0, 11.0, 16.0, 12.0, 13.0])     # verre enn snittet

def r2(y, y_hat):
    return 1 - np.sum((y - y_hat) ** 2) / np.sum((y - y.mean()) ** 2)

for navn, p in [("modell", modell), ("baseline", baseline), ("dårlig", daarlig)]:
    print(f"{navn:9s} R² = {r2(y, p):7.3f}")

Hvorfor er det nettopp kvadratet av residualene vi minimerer? Svaret kommer fra en antakelse om støyen. Anta at dataene er generert som

der støyen ${\epsilon }_i$ er uavhengig for hvert punkt. Det kalles Gaussisk støy. Da er hver $y_i$ normalfordelt rundt prediksjonen ${\hat{y}}_i$, og tettheten for alle de observerte dataene er produktet av $n$ slike tettheter. Maksimum likelihood betyr å velge parameterne som gjør de observerte dataene mest sannsynlige. Det er lettere å maksimere logaritmen, fordi produktet da blir en sum. Det gir ƒGaussisk log-likelihood:

Det første leddet avhenger ikke av $w$ og $b$. Det andre er en negativ konstant ganger summen av kvadrerte residualer. Å maksimere $\log L$ er derfor nøyaktig det samme som å minimere ${\sum }_i(y_i-{\hat{y}}_i{)}^2$, og dermed MSE. Støyvariansen ${\sigma }^2$ skalerer bare leddet, så den endrer ikke hvilke parametere som er best.

Bytter vi antakelsen, bytter vi tap. Laplace-fordelt støy med skala $s$ har spissere topp og tyngre haler, og da blir ƒLaplace-log-likelihood en konstant minus $\frac{1}{s}{\sum }_i|y_i-{\hat{y}}_i|$. Maksimum likelihood blir det samme som å minimere MAE. Tunge haler betyr at store avvik ikke er «umulige», så noen få uteliggere får mindre makt. Å velge tapsfunksjon er altså å velge en antakelse om hvordan feilene ser ut.

En lineær modell er additiv: hver feature bidrar med sitt eget ledd, uavhengig av de andre. Matematisk betyr det at

uansett verdien av $x_2$. Effekten av $x_1$ er den samme overalt. Tenk på plantevekst med vann ($x_1$) og gjødsel ($x_2$). En lineær modell sier at ett gram gjødsel gir like mye vekst om planten blir vannet eller står i tørr jord. I virkeligheten hjelper gjødsel nesten ikke uten vann. Når effekten av én feature avhenger av en annen, kaller vi det en vekselvirkning eller et samspill, og det kan den rene lineære modellen ikke uttrykke.

Løsningen er å lage et interaksjonsledd, en ny feature $x_1{x}_2$. Det gir ƒmodell med samspillsledd:

Nå avhenger effekten av $x_1$ av $x_2$. Modellen er fortsatt lineær i vektene, så den tilpasses med samme metode og samme MSE. Vi har bare lagt til en kolonne i $X$. Er $x_2$ binær, som helg eller hverdag, gir leddet hver gruppe sin egen stigning.

Uten samspill gir modellen parallelle linjer: gjødsel flytter linjen opp, men endrer ikke stigningen. Med samspill vifter linjene ut. Det samme prinsippet forklarer hvorfor en lineær modell aldri kan lære XOR-mønsteret uten et produktledd eller et skjult lag (kapittel 11).

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(5)
n = 200
vann = rng.uniform(0, 10, n)                # liter per uke
gjodsel = rng.uniform(0, 10, n)             # gram per uke
samspill = 0.08  # @param 0..0.15 0.01
vekst = 2 + 0.3 * vann + 0.02 * gjodsel + samspill * vann * gjodsel + rng.normal(0, 0.5, n)

X_lin = np.column_stack([vann, gjodsel])
X_int = np.column_stack([vann, gjodsel, vann * gjodsel])   # ny kolonne x1·x2
for navn, X in [("uten x1·x2", X_lin), ("med x1·x2", X_int)]:
    m = LinearRegression().fit(X, vekst)
    print(f"{navn:10s} vekter = {np.round(m.coef_, 3)}, R² = {m.score(X, vekst):.3f}")

# m er nå modellen med x1·x2: effekten av +1 g gjødsel er w2 + w12 · vann
for v in [0, 8]:
    print(f"effekt av +1 g gjødsel når vann = {v}: {m.coef_[1] + m.coef_[2] * v:.2f}")

Samme triks som for samspill gir kurver. Vil vi modellere en sammenheng som bøyer seg, legger vi til polynomiske features som $x^2$ og $x^3$:

Kurven i $x$ er ikke rett, men modellen er fortsatt lineær i vektene $w$, og det er det som teller for metoden. Vi tilpasser vanlig lineær regresjon på kolonnene $x$, $x^2$ og $x^3$. Med flere features gir PolynomialFeatures(degree=2) både kvadratene og alle parvise produkter, så samspillsleddene fra forrige steg er et spesialtilfelle. Antall kolonner vokser raskt med graden, og høye potenser får svært ulike skalaer, så skalering er lurt. For høy grad gir en kurve som slynger seg gjennom støyen. Det er overtilpasning, og det tar kapittel 7 for seg.

Hvordan finner vi vektene? For MSE finnes en lukket løsning. Samler vi en kolonne med enere (for $b$) og featurene i $X$, er minimumet gitt av ƒnormallikningen:

Løsningen krever at $X^{\top }X$ er inverterbar, altså at ingen kolonne er en lineær kombinasjon av de andre, og den blir dyr når antall features er stort. Gradient descent er alternativet som skalerer til store datasett og modeller uten lukket løsning, og det har sitt eget kapittel (kapittel 3).