Massima verosimiglianza

Gli stimatori visti finora, come la media campionaria o i minimi quadrati, non fanno ipotesi sulla distribuzione dei dati. La stima a massima verosimiglianza (maximum likelihood, ML) fa il contrario: si sceglie un modello probabilistico per i dati e se ne stimano i parametri in modo che siano il più possibile coerenti con ciò che si è osservato.

Esempio: moneta truccata

Lanciamo una moneta 8 volte e otteniamo:

Cinque teste e tre croci. Il parametro da stimare è $p$, la probabilità di testa. Se i lanci sono indipendenti, la probabilità di osservare esattamente questa sequenza è il prodotto delle probabilità dei singoli lanci:

$$p^5\,(1-p)^3$$

Se $p$ fosse noto, questa formula darebbe la probabilità dei dati. Ma noi conosciamo i dati, e $p$ no. Letta come funzione di $p$, con i dati fissati, la stessa formula si chiama verosimiglianza (likelihood), $\mathcal{L}(p \mid \text{dati})$.

Per capire quale $p$ spiega meglio i dati, immaginiamo diverse monete, truccate o no, e chiediamoci, per ciascuna, quanto sarebbe stato probabile ottenere proprio la nostra sequenza. Per ogni candidato si calcolano due fattori: $p^5$, la probabilità delle cinque teste, e $(1-p)^3$, quella delle tre croci. Il loro prodotto è la verosimiglianza:

Per esempio, con una moneta onesta ($p = 0{,}5$): $0{,}5^5 \cdot 0{,}5^3 = 1/256 \approx 0{,}0039$. Circa 4 volte su mille una moneta onesta produrrebbe esattamente questa sequenza.

Una moneta con $p = 0{,}9$ renderebbe facili le teste, ma tre croci su otto molto improbabili ($0{,}1^3 = 0{,}001$); una con $p = 0{,}3$ farebbe il contrario con le cinque teste. Il valore che rende i dati osservati più verosimili è $\hat p = 5/8 = 0{,}625$.

Questo esempio si basa su pochi dati (8 lanci), quindi il risultato osservato potrebbe essere dovuto semplicemente al caso. Tuttavia, questo riguarda l'affidabilità della stima, cioè quanto il valore stimato sia vicino a quello reale, e non il procedimento di stima, che consiste nel trovare il valore del parametro che rende massima la verosimiglianza dei dati osservati.

Stima a massima verosimiglianza

La stima a massima verosimiglianza è il valore dei parametri che massimizza la verosimiglianza:

$$\hat\theta_{\text{ML}} = \arg\max_{\theta}\ \mathcal{L}(\theta \mid \text{dati})$$

Con dati indipendenti e identicamente distribuiti (i.i.d.), $\mathcal{L}$ è il prodotto delle probabilità (o, con dati continui, delle densità) dei singoli dati. In pratica, però, non si massimizza $\mathcal{L}$ direttamente ma il suo logaritmo, $\ln \mathcal{L}$; o, equivalentemente, si minimizza $-\ln \mathcal{L}$, la meno-log-verosimiglianza (negative log-likelihood), come si fa con una funzione di costo.

Il massimo non si sposta

Il logaritmo è una funzione crescente: se $a > b$, allora $\ln a > \ln b$. Cambia i valori, ma non l'ordine. Riprendendo la tabella della moneta:

Il $p$ che massimizza $\ln \mathcal{L}$ è lo stesso che massimizza $\mathcal{L}$.

I prodotti diventano somme

Il logaritmo ha due proprietà che fanno al caso nostro: $\ln(a \cdot b) = \ln a + \ln b$ e $\ln(a^k) = k \ln a$. Applicate alla moneta:

$$\ln \mathcal{L} = \ln\bigl(p^5 (1-p)^3\bigr) = 5 \ln p + 3 \ln(1-p)$$

Un prodotto di potenze è diventato una somma di termini semplici. Per trovare il massimo si deriva e si pone la derivata uguale a zero; derivare una somma significa derivare un termine alla volta:

$$\frac{d \ln \mathcal{L}}{dp} = \frac{5}{p} - \frac{3}{1-p} = 0 \quad\Longrightarrow\quad 5(1-p) = 3p \quad\Longrightarrow\quad \hat p = \frac{5}{8}$$

Derivare direttamente $p^5(1-p)^3$ richiede la regola del prodotto e porta allo stesso risultato con più fatica.

Numeri troppo piccoli per il computer

C'è anche un motivo pratico. Con 1100 lanci di una moneta onesta, la verosimiglianza di una sequenza è $0{,}5^{1100} \approx 10^{-331}$: troppo piccola per i numeri in virgola mobile di un computer, che non scendono sotto circa $5 \cdot 10^{-324}$ e la arrotondano a zero. A quel punto tutti i candidati hanno verosimiglianza 0 e non c'è più niente da confrontare. Il logaritmo invece vale $1100 \cdot \ln 0{,}5 \approx -762$, un numero del tutto ordinario.

Esempio: media di una gaussiana

L'idea è la stessa della moneta: si sceglie il valore del parametro che rende i dati osservati più compatibili con il modello. Cambiano solo i dati e il parametro:

• con la moneta, i dati erano teste e croci e il parametro era la probabilità di testa $p$;
• con la gaussiana, i dati sono misure numeriche e il parametro è la media $\mu$, cioè il centro della campana.

Vediamo perché, in questo caso, la stima migliore è proprio la media dei dati.

Il problema: dove sta il centro della campana?

Abbiamo misurato due valori, $y(1) = 4$ e $y(2) = 6$. Sappiamo che vengono da una gaussiana con varianza 1, ma non ne conosciamo la media $\mu$, che fissa la posizione della campana sull'asse:

• con $\mu = 4$ il centro coincide con il primo dato;
• con $\mu = 6$ coincide con il secondo;
• con $\mu = 5$ sta esattamente a metà.

Quale scegliere? Quella che rende i dati osservati più verosimili.

Quanto è compatibile un dato con la media?

Per la moneta ogni lancio aveva probabilità $p$ o $1-p$. Con dati continui quel ruolo lo gioca la densità: per una gaussiana di media $\mu$ e varianza 1,

$$f(y \mid \mu) = \frac{1}{\sqrt{2\pi}}\, e^{-\frac{1}{2}(y - \mu)^2}$$

Spiegazione della formula

La formula ha tre pezzi:

• $(y - \mu)^2$ è la distanza del dato dal centro, al quadrato. Il quadrato toglie il segno: 4 e 6 sono ugualmente lontani da 5.
• $e^{-\frac{1}{2}(\dots)}$ trasforma la distanza in un numero tra 0 e 1: vale 1 al centro e cala man mano che ci si allontana. Il $\tfrac{1}{2}$ viene dalla varianza: in generale è $\tfrac{1}{2\sigma^2}$, e qui $\sigma^2 = 1$.
• $\frac{1}{\sqrt{2\pi}}$ serve solo a far valere 1 l'area sotto la campana. Non dipende da $\mu$, quindi non cambia quale media vince.

Distanza dal centro	0	1	2	3
Densità $f$	0,40	0,24	0,05	0,004

In breve: più il dato è vicino a $\mu$, più la densità è alta.

La densità è alta se il dato è vicino al centro $\mu$ e bassa se è lontano. Il grafico mette a confronto due campane candidate:

Con $\mu = 5$ entrambi i dati sono a distanza 1 dal centro e ricevono la stessa densità, 0,24. Con $\mu = 4$ il primo dato è proprio sul picco (0,40), ma il secondo finisce sulla coda (0,05).

Per trovare la media migliore si procede come per la moneta. I dati sono indipendenti, quindi si moltiplicano le densità, e si passa al logaritmo, che trasforma il prodotto in una somma ed elimina gli esponenziali:

$$\ln \mathcal{L}(\mu) = \text{cost} - \tfrac{1}{2}(4 - \mu)^2 - \tfrac{1}{2}(6 - \mu)^2$$

La costante viene dai fattori $\frac{1}{\sqrt{2\pi}}$ e non dipende da $\mu$, quindi si può ignorare. Restano due termini, uno per dato: sono penalità, che abbassano la log-verosimiglianza tanto più quanto $\mu$ è lontana dal dato.

Proviamo diversi valori della media

Calcoliamo la log-verosimiglianza, senza la costante, per alcune medie candidate:

Si cerca il valore più alto, cioè il meno negativo: è −1, con $\mu = 5$. Con $\mu = 4$ il primo dato è al centro, ma il secondo è a distanza 2 e paga 2; con $\mu = 6$ succede il contrario. Con $\mu = 5$ entrambi distano solo 1 e pagano 0,5 ciascuno: è il compromesso migliore.

Il risultato esatto con la derivata

La tabella prova solo alcuni valori. Per trovare il massimo esatto si deriva la funzione, senza la costante:

$$g(\mu) = -\tfrac{1}{2}(4 - \mu)^2 - \tfrac{1}{2}(6 - \mu)^2 \qquad\Longrightarrow\qquad g'(\mu) = (4 - \mu) + (6 - \mu)$$

e si pone la derivata uguale a zero:

$$(4 - \mu) + (6 - \mu) = 0 \quad\Longrightarrow\quad 10 - 2\mu = 0 \quad\Longrightarrow\quad \hat\mu_{\text{ML}} = 5$$

È lo stesso valore della tabella, ma ora sappiamo che è il massimo esatto e non solo il migliore tra quelli provati. Che sia un massimo, e non un minimo, lo conferma la derivata seconda: $g''(\mu) = -2 < 0$, la funzione è una parabola rivolta verso il basso. L'equazione si legge anche così: $(4 - \mu)$ e $(6 - \mu)$ sono gli scarti dei dati dal centro, e il massimo è dove si compensano, uno a $-1$ e l'altro a $+1$.

Il risultato, $\hat\mu_{\text{ML}} = \frac{4 + 6}{2}$, è la media campionaria dei dati. Con $N$ dati il ragionamento è identico e porta a $\hat\mu_{\text{ML}} = \frac{1}{N}\sum_{i=1}^{N} y(i)$: la media, scelta di solito "a intuito", è la stima a massima verosimiglianza quando i dati sono gaussiani.

Risultato: i minimi quadrati

Riprendiamo il modello lineare $y(i) = \varphi(i)^\top \theta + \varepsilon(i)$ e aggiungiamo un'ipotesi: il rumore è gaussiano, con media zero e varianza $\lambda^2$, indipendente da un'osservazione all'altra. Allora ogni $y(i)$ è una gaussiana centrata su $\varphi(i)^\top\theta$, e la meno-log-verosimiglianza diventa:

$$-\ln \mathcal{L}(\theta) = \text{cost} + \frac{1}{2\lambda^2} \sum_{i=1}^{N} \bigl( y(i) - \varphi(i)^\top \theta \bigr)^2$$

La costante e il fattore $1/(2\lambda^2)$ non dipendono da $\theta$, quindi non spostano il minimo. Resta la somma dei quadrati degli errori:

Con rumore gaussiano, la stima a massima verosimiglianza di un modello lineare coincide con quella a minimi quadrati.

È una giustificazione probabilistica di una scelta che finora sembrava arbitraria: perché proprio il quadrato degli errori? Perché è ciò che esce da un rumore gaussiano. Con altre ipotesi sul rumore escono altre funzioni di costo. E con dati che valgono 0 o 1, modellati come lanci di una moneta con $p$ che dipende dagli ingressi, esce la funzione di costo della regressione logistica.

I passaggi per il modello lineare

Chiamiamo $\varepsilon(i) = y(i) - \varphi(i)^\top\theta$ l'errore della singola osservazione. La densità di ciascuna osservazione è

$$f\bigl(y(i)\bigr) = \frac{1}{\sqrt{2\pi\lambda^2}}\, e^{-\varepsilon(i)^2 / (2\lambda^2)}$$

Per l'indipendenza, $\mathcal{L}$ è il prodotto delle $N$ densità, e il logaritmo trasforma il prodotto in una somma:

$$\begin{aligned} \ln \mathcal{L}(\theta) &= \sum_{i=1}^{N} \ln \frac{1}{\sqrt{2\pi\lambda^2}} \\ &\quad - \sum_{i=1}^{N} \frac{\varepsilon(i)^2}{2\lambda^2} \\ &= -\frac{N}{2}\ln(2\pi\lambda^2) \\ &\quad - \frac{1}{2\lambda^2} \sum_{i=1}^{N} \varepsilon(i)^2 \end{aligned}$$

Il primo termine non dipende da $\theta$. Massimizzare $\ln \mathcal{L}$ equivale quindi a minimizzare la somma dei quadrati.

Le proprietà dello stimatore ML

Sotto ipotesi di regolarità abbastanza generali, lo stimatore a massima verosimiglianza è:

1. asintoticamente corretto: con pochi dati può essere distorto, ma il bias sparisce al crescere di $N$;
2. consistente: al crescere di $N$ la stima converge al valore vero;
3. asintoticamente efficiente: con molti dati la sua varianza si avvicina al limite di Cramér-Rao, la minima possibile per uno stimatore corretto;
4. asintoticamente normale: con molti dati la stima ha una distribuzione gaussiana centrata sul valore vero.

Riepilogo

• La verosimiglianza è la probabilità (o la densità) dei dati osservati, vista come funzione dei parametri.
• La stima ML la massimizza; in pratica si minimizza $-\ln \mathcal{L}$, che trasforma i prodotti in somme.
• Per una moneta dà la frazione di teste; per una gaussiana, la media campionaria.
• Con rumore gaussiano, ML su un modello lineare equivale ai minimi quadrati.