Logistisk regresjon og klassifikasjon

En kommune vil prioritere vedlikehold av gatelys og forutsi hvilke lyktestolper som svikter i løpet av vinteren. Hvert eksempel har label 1 (svikter) eller 0 (holder), så dette er binær klassifikasjon. Lineær regresjon passer dårlig: den kan predikere $-0,4$ eller $1,7$, og det er ikke sannsynligheter.

Logistisk regresjon starter likevel lineært. ƒLineær klassifikasjonsscore gir én score $z=w^{\top }x+b$, der positive vekter trekker scoren opp. Scoren kan være et hvilket som helst reelt tall. Deretter kommer sigmoiden, ƒsigmoidfunksjonen:

Sigmoiden gir $\sigma (0)=0,5$, nærmer seg 1 for store positive scorer og 0 for store negative, og er symmetrisk: $\sigma (-z)=1-\sigma (z)$. Den er brattest rundt $z=0$, der den deriverte $\sigma (z)(1-\sigma (z))$ er $0,25$. Ute i halene endrer selv en stor scoreendring sannsynligheten nesten ikke.

Hvorfor ikke en trappefunksjon som gir 1 når $z\ge 0$ og ellers 0? Fordi trappen har derivert 0 overalt unntatt i selve hoppet, og da har gradient descent ingenting å gå etter. Sigmoiden er glatt, så hver liten vektendring gir en liten, målbar endring i tapet.

Tallet $p$ er modellens estimerte sannsynlighet for klasse 1, ikke en garantert sann sannsynlighet. Det kan bare tas bokstavelig hvis modellen er godt kalibrert. Sannsynlighetskalibrering betyr at blant alle stolper som får $p\approx 0,3$, svikter omtrent 30 prosent. Det sjekkes på valideringsdata, gruppe for gruppe. Klassevekter, oversampling (kapittel 8), sterk regularisering og endrede forhold etter treningen kan alle forskyve sannsynlighetene, selv når rangeringen av stolpene er god.

Sannsynligheter kan også uttrykkes som odds, ƒodds: forholdet mellom sjansen for at noe skjer og sjansen for at det ikke skjer. $p=0,8$ gir odds $4$ (fire mot én), $p=0,5$ gir odds $1$, og $p=0,2$ gir odds $0,25$. Odds går fra 0 til uendelig.

Tar vi logaritmen, får vi en logit, eller log-odds, som går fra $-\infty$ til $\infty$, akkurat som scoren. Det er nøkkelen til modellen: logistisk regresjon antar at log-oddsen er lineær i featurene, ƒlog-odds:

Logit og sigmoid er hverandres inverse. Det er derfor modellen heter logistisk regresjon: den modellerer log-oddsen som en lineær funksjon av featurene.

Dette gir en presis tolkning av vektene. Øker feature $x_j$ med én enhet mens resten holdes fast, øker log-oddsen med $w_j$, og oddsen ganges med $e^{w_j}$. For gatelysene er vekten for tidligere feil $0,9$, så én feil til ganger oddsen med $e^{0,9}\approx 2,46$. Det betyr ikke at sannsynligheten øker med $0,9$ eller med 90 prosent. Hvor mye $p$ endres, avhenger av hvor på sigmoidkurven stolpen ligger: mest nær $p=0,5$, minst ute i halene.

Symmetrien $\sigma (-z)=1-\sigma (z)$ gir også en numerisk stabil sigmoid: for negative $z$ kan du bruke $\sigma (z)=e^z/(1+e^z)$, slik at eksponenten aldri blir et stort positivt tall. Det er kodesteget på slutten av dette steget.

import numpy as np

w_feil = 0.9  # @param -1.5..1.5 0.1

for z in [-3.0, -0.7, 1.0, 3.0]:            # score før én feil til
    p0 = 1 / (1 + np.exp(-z))
    p1 = 1 / (1 + np.exp(-(z + w_feil)))    # én feil til: scoren øker med w
    odds0, odds1 = p0 / (1 - p0), p1 / (1 - p1)
    print(f"z = {z:4.1f}:  p {p0:.3f} -> {p1:.3f}  (endring {p1 - p0:+.3f}),"
          f"  oddsforhold {odds1 / odds0:.3f}")
print(f"e^w = {np.exp(w_feil):.3f}")

Modellen gir en sannsynlighet, men kommunen må ta en beslutning: skal stolpen inspiseres eller ikke? En beslutningsterskel $\tau$ gjør sannsynligheten om til klasse med ƒbeslutningsregel: predikér 1 når $p\ge \tau$. Terskelen er ikke en del av treningen. Endrer du $\tau$, er score og sannsynlighet de samme; bare beslutningen flytter seg.

$\tau =0,5$ er en konvensjon, ikke en optimal regel. Senker du terskelen, blir flere stolper positive: du fanger flere av dem som faktisk svikter, men inspiserer også flere friske. Hever du den, skjer det motsatte. Riktig $\tau$ avhenger av hva en oversett feil koster mot en unødvendig inspeksjon, og det kommer vi tilbake til i kapittel 6.

En beslutningsgrense er alle punkter der $p=\tau$ akkurat. Fordi sigmoiden er voksende, er $p\ge \tau$ det samme som at scoren er minst logiten til terskelen, ƒbeslutningsgrense for en terskel:

Med $\tau =0,5$ er høyresiden 0. Med to features er grensen en rett linje, med tre et plan og generelt et hyperplan. Logistisk regresjon gir altså en lineær grense selv om sannsynlighetskurven er S-formet. Å endre $\tau$ forskyver grensen parallelt, mens vektene bestemmer retningen. Trenger klassene en buet grense, må du lage nye features, for eksempel $x_1^2$ eller $x_1{x}_2$, eller velge en annen modell.

For å lære vektene trenger vi et tap som vurderer sannsynlighetene. Binary cross-entropy (BCE) er gjennomsnittet av $-\ln$ til sannsynligheten modellen ga den riktige klassen, ƒbinary cross-entropy:

For $y=1$ er bidraget $-\ln p$: $0,105$ for $p=0,9$, $0,693$ for $p=0,5$ og $2,303$ for $p=0,1$. Når $p$ går mot 0, går tapet mot uendelig. For $y=0$ er bidraget $-\ln (1-p)$. BCE er den negative log-likelihooden til en Bernoulli-modell, så å minimere den er maksimum likelihood, akkurat som MSE svarer til normalfordelt støy i kapittel 2.

Hvorfor ikke bare bruke MSE på sannsynlighetene? Tre grunner:

• •Gradienten forsvinner. Med MSE er den deriverte med hensyn på scoren $2(p-y)p(1-p)$. Er modellen sikker og feil, for eksempel $p=0,01$ når $y=1$, er $p(1-p)$ nesten 0, og punktet drar nesten ikke i vektene. Med BCE er den deriverte $p-y\approx -0,99$, altså full styrke.
• •Konveksitet. BCE er konveks i vektene for logistisk regresjon. MSE gjennom en sigmoid er det ikke, og kan ha flate platåer.
• •Straffen. MSE kan aldri gi mer enn 1 per punkt. BCE straffer selvsikre feil langt hardere.

I NumPy skrives BCE slik, med sannsynlighetene klippet bort fra 0 og 1 så logaritmen ikke blir uendelig:

Python
p = np.clip(p, 1e-7, 1 - 1e-7)
bce = -np.mean(y * np.log(p) + (1 - y) * np.log(1 - p))

Vanlige feil er å glemme minustegnet (tapet blir negativt, og minimering gjør modellen verre), å bytte om $p$ og $1-p$, og å droppe $(1-y)$-leddet slik at de negative eksemplene ikke teller. Biblioteker regner gjerne BCE direkte fra scoren $z$, som er numerisk enda mer stabilt.

import numpy as np

def sigmoid(z):
    return 1 / (1 + np.exp(-z))

feilsikkerhet = 3.0  # @param 0..8 0.5
y = np.array([1, 0, 1, 0, 1])
z = np.array([2.0, -1.5, 0.3, 0.8, -feilsikkerhet])   # siste punkt: y = 1, men negativ score
p = np.clip(sigmoid(z), 1e-7, 1 - 1e-7)              # klipp bort fra 0 og 1 før log
tap = -(y * np.log(p) + (1 - y) * np.log(1 - p))

print(f"siste punkt: z = {z[-1]:.1f}, p = {p[-1]:.4f}, tap = {tap[-1]:.2f}")
print(f"BCE = {tap.mean():.3f}   MSE = {np.mean((p - y) ** 2):.3f}")
# Gradienten med hensyn på scoren z for det siste punktet
print(f"BCE-gradient p - y:            {p[-1] - 1:+.4f}")
print(f"MSE-gradient 2(p - y)p(1 - p): {2 * (p[-1] - 1) * p[-1] * (1 - p[-1]):+.4f}")

Vektene læres med gradient descent fra kapittel 3. Med kjerneregelen blir den deriverte av BCE for ett punkt med hensyn på scoren overraskende enkel:

Siden $z=w^{\top }x+b$, er $\partial z/\partial w_j=x_j$, og vi får ƒgradienten til BCE:

Det er samme mønster som for lineær regresjon: prediksjon minus fasit, ganger input. Et punkt med $y=1$ og lav $p$ drar vektene mot høyere score, mens et punkt som allerede er riktig og sikkert, nesten ikke drar.

Logistisk regresjon har ingen lukket formel for de beste vektene, så den trenes iterativt. Fordi BCE er konveks her, finner gradient descent det globale minimumet med passende $\eta$. Ett unntak: er klassene perfekt lineært separerbare, kan tapet alltid presses lavere ved å skalere vektene opp, og vektene vokser uten grense. Regularisering, i scikit-learn styrt av parameteren C, stopper det.

Læringsraten virker som i kapittel 3, med én forskjell: tapsflaten er ikke en ren parabel, så for stor $\eta$ gir et tap som hopper opp og ned fra epoke til epoke i stedet for å synke jevnt. Uskalerte features gjør det verre, fordi én bratt retning i flaten krever liten $\eta$ for alle vektene. Standardiser derfor featurene før gradientbasert trening.

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(7)
x = np.concatenate([rng.normal(2, 1, 30), rng.normal(5, 1, 30)])  # uskalert feature
y = np.concatenate([np.zeros(30), np.ones(30)])
eta = 0.5  # @param 0.1..3 0.1
standardiser = 0  # @param 0..1
if standardiser:
    x = (x - x.mean()) / x.std()

w = b = 0.0
bce = []
for epoke in range(30):
    p = 1 / (1 + np.exp(-np.clip(w * x + b, -30, 30)))
    bce.append(-np.mean(y * np.log(p) + (1 - y) * np.log(1 - p)))
    w -= eta * np.mean((p - y) * x)    # gradienten: (p - y) ganger x
    b -= eta * np.mean(p - y)
print(f"BCE etter {len(bce)} epoker: {bce[-1]:.3f}   (w = {w:.2f}, b = {b:.2f})")

plt.figure()
plt.plot(bce)
plt.xlabel("epoke")
plt.ylabel("BCE")
plt.tight_layout()
plt.show()

Gatelysene kan svikte på ulike måter: pæra, kabelen eller styringen. Med $K$ klasser gir modellen én score per klasse, $z_k=w_k^{\top }x+b_k$, og softmax gjør scorene om til sannsynligheter, ƒsoftmax:

Eksponentialen gjør alle tallene positive, og divisjonen får dem til å summere til 1. Den største scoren får den største sannsynligheten, så argmax er det samme før og etter softmax. Legger du samme konstant til alle scorene, endres ingenting, og det utnyttes i kode: trekk fra den største scoren før du tar eksponentialen, så unngår du overflyt. Med to klasser er softmax det samme som sigmoiden av differansen $z_1-z_0$, så logistisk regresjon er et spesialtilfelle.

Tapet er kategorisk kryssentropi, ƒkategorisk kryssentropi: $-\ln$ av sannsynligheten for den sanne klassen. Med fasiten som en one-hot-vektor $y$ er gradienten med hensyn på scorene igjen $p-y$.

Et alternativ er én-mot-resten: $K$ separate binære modeller, én per klasse. Da summerer ikke sannsynlighetene til 1 av seg selv. Scikit-learns LogisticRegression bruker softmax-varianten for flere klasser med standardinnstillingene.

import numpy as np

def softmax_naiv(z):
    e = np.exp(z)
    return e / e.sum()

def softmax(z):
    e = np.exp(z - z.max())      # trekk fra maks: samme svar, ingen overflyt
    return e / e.sum()

forskyvning = 0  # @param 0..1000 50
z = np.array([2.0, 1.0, -1.0]) + forskyvning   # score for pære, kabel og styring
print(f"naiv:   {np.round(softmax_naiv(z), 4)}")
print(f"stabil: {np.round(softmax(z), 4)}")

y = 1                            # indeksen til den sanne klassen (kabel)
print(f"kryssentropi = {-np.log(softmax(z)[y]):.3f}")
print(f"gradient p - y: {np.round(softmax(z) - np.eye(3)[y], 4)}")