20 - Statistica: stimare, confrontare, decidere

Corso Machine Learning Engineer · Modulo 02 · Matematica e statistica · Lezione 5

Obiettivo. Passare dalla probabilità alla statistica: dai dati osservati risalire ai parametri sconosciuti che li hanno generati, quantificando quanto ci si può fidare. Alla fine saprai riassumere un campione senza lasciarti ingannare, spiegare (con la dimostrazione) perché la varianza campionaria divide per n − 1, stimare un parametro con la massima verosimiglianza, costruire intervalli di confidenza (anche con il bootstrap), leggere correttamente un test di ipotesi, e derivare da zero la retta di regressione ai minimi quadrati.

Falso mito da sfatare. «La media del mio campione è la media vera.» È solo la migliore stima che hai, e cambia da un campione all'altro. La statistica non elimina l'incertezza: la misura.

Introduzione

La probabilità va dal modello ai dati: «se la moneta è onesta, come saranno 100 lanci?». La statistica fa il percorso inverso: «ho visto 62 teste su 100 lanci, che cosa posso dire della moneta?». I termini di base, spiegati qui al primo uso:

• Popolazione: l'insieme di tutto ciò di cui vuoi parlare (tutti i clienti, tutte le misure possibili).
• Campione: la parte che hai osservato davvero, x₁, …, xₙ. Deve essere scelto a caso (e indipendente) perché i risultati valgano.
• Parametro: un numero che descrive la popolazione (la media vera μ); è fisso e sconosciuto.
• Statistica: un numero calcolato dal campione (la media campionaria x̄); è casuale, perché dipende dal campione estratto.
• Stimatore: la ricetta (la formula) con cui da un campione si ricava una stima di un parametro.

Le cinque voci si dispongono su due piani: parametro e popolazione sono fissi e ignoti, statistica e campione sono casuali perché cambiano a ogni estrazione. Tutta la statistica inferenziale vive nel passaggio da un piano all'altro, e ogni formula di questa lezione è un modo di dire quanto la statistica può discostarsi dal parametro.

L'analogia: la minestra. Per sapere se la minestra è salata non la bevi tutta: ne assaggi un cucchiaio. Il cucchiaio è il campione, la pentola la popolazione. Funziona se hai mescolato bene (campione casuale) e se il cucchiaio è abbastanza grande. E sai che due cucchiai diversi daranno risposte leggermente diverse: quanto, lo dice la statistica.

Il percorso parte dalla statistica descrittiva, passa per stimatori, distorsione e correzione n − 1, massima verosimiglianza, intervalli di confidenza e bootstrap, test di ipotesi e regressione ai minimi quadrati, e finisce con la differenza fra correlazione e causalità. Si appoggia sul teorema del limite centrale e sulla probabilità condizionata.

Statistica descrittiva

Centro: la media x̄ = (1/n) Σ xᵢ (il centro di massa), la mediana (il valore in mezzo dopo l'ordinamento), la moda (il più frequente). Dispersione: varianza s² e deviazione standard s; lo scarto interquartile (IQR, la distanza tra il primo e il terzo quartile, cioè tra il valore sotto cui sta il 25% dei dati e quello sotto cui sta il 75%) che ignora le code.

La differenza fra media e mediana è il test più rapido per capire la forma: se la media è molto più grande della mediana, c'è una coda a destra. Mostriamolo con un esempio sintetico di redditi mensili (distribuzione log-normale, a coda lunga), a cui aggiungiamo un singolo valore estremo.

import numpy as np
import pandas as pd

rng = np.random.default_rng(10)
redditi = rng.lognormal(mean=7.4, sigma=0.5, size=1000)          # valori tipici di qualche migliaio
con_estremo = np.append(redditi, 5_000_000.0)                       # un solo valore anomalo

def riassunto(v):
    q1, q3 = np.percentile(v, [25, 75])
    return {"media": v.mean(), "mediana": np.median(v), "dev.std": v.std(ddof=1), "IQR": q3 - q1}

tab = pd.DataFrame({"senza estremo": riassunto(redditi), "con un estremo": riassunto(con_estremo)}).round(0)
print(tab)
print("asimmetria (skewness) senza estremo:", round(float(pd.Series(redditi).skew()), 2))

Output ottenuto:

         senza estremo  con un estremo
media           1786.0          6779.0
mediana         1581.0          1582.0
dev.std          923.0        157981.0
IQR             1082.0          1082.0
asimmetria (skewness) senza estremo: 1.48

Cosa fa il codice. Genera 1000 redditi log-normali (l'esponenziale di una variabile normale, quindi sempre positivi e con coda a destra), poi ne aggiunge uno da cinque milioni. La funzione riassunto calcola media, mediana, deviazione standard (con ddof=1, la correzione che spieghiamo tra poco) e IQR. Il confronto delle due colonne mostra chi regge: la media e la deviazione standard esplodono per colpa di una sola osservazione, mediana e IQR si muovono appena. Sono indicatori robusti. Quale usare dipende dalla domanda: la spesa totale si conosce dalla media, il caso tipico dalla mediana.

Stimatori e la correzione n − 1

Distorsione, varianza, errore quadratico

Uno stimatore θ̂ di un parametro θ si valuta con tre criteri:

• Distorsione (bias): E[θ̂] − θ. Uno stimatore non distorto sbaglia «a destra e a sinistra» in egual misura, in media.
• Varianza: quanto θ̂ oscilla da un campione all'altro.
• Errore quadratico medio: MSE = E[(θ̂ − θ)²] = distorsione² + varianza (si dimostra aggiungendo e togliendo E[θ̂] dentro il quadrato: il doppio prodotto si annulla).

La media campionaria è non distorta (E[x̄] = μ, per la linearità del valore atteso) e ha varianza σ²/n. La sua deviazione standard, σ/√n, si chiama errore standard: dice di quanto la media campionaria oscilla.

Perché la varianza campionaria divide per n − 1

Il candidato naturale per stimare σ² è (1/n) Σ (xᵢ − x̄)². Ma è distorto verso il basso. L'intuizione: x̄ è il centro «su misura» di questo campione, è il punto che minimizza gli scarti quadratici; i dati sono dunque sempre più vicini a x̄ che alla vera μ. Quindi gli scarti attorno a x̄ sottostimano quelli attorno a μ.

Dimostrazione. Aggiungiamo e togliamo μ: xᵢ − x̄ = (xᵢ − μ) − (x̄ − μ). Elevando al quadrato e sommando su i (il doppio prodotto vale −2(x̄−μ)·Σ(xᵢ−μ) = −2n(x̄−μ)², perché Σ(xᵢ−μ) = n(x̄−μ)):

Σ (xᵢ − x̄)² = Σ (xᵢ − μ)² − n (x̄ − μ)².

Prendiamo il valore atteso: E[Σ(xᵢ−μ)²] = nσ² (ogni termine ha media σ²) e E[n(x̄−μ)²] = n · Var(x̄) = n · σ²/n = σ². Quindi

E[Σ (xᵢ − x̄)²] = nσ² − σ² = (n − 1)σ².

Dividendo per n − 1 (e non per n) si ottiene uno stimatore non distorto: s² = (1/(n−1)) Σ (xᵢ − x̄)². La differenza 1 è il «grado di libertà» consumato per stimare x̄ dagli stessi dati.

rng = np.random.default_rng(1)
sigma2_vera = 4.0                       # varianza vera della popolazione: N(0, 4)
n, ripetizioni = 5, 200_000
campioni = rng.normal(0, np.sqrt(sigma2_vera), size=(ripetizioni, n))

stima_n = campioni.var(axis=1, ddof=0)              # divide per n
stima_n1 = campioni.var(axis=1, ddof=1)             # divide per n - 1
print("varianza vera:", sigma2_vera)
print("media delle stime, divisore n    :", round(float(stima_n.mean()), 3), "| teoria (n-1)/n * σ² =", (n - 1) / n * sigma2_vera)
print("media delle stime, divisore n - 1:", round(float(stima_n1.mean()), 3))

print("MSE divisore n    :", round(float(((stima_n - sigma2_vera) ** 2).mean()), 3))
print("MSE divisore n - 1:", round(float(((stima_n1 - sigma2_vera) ** 2).mean()), 3))

Output ottenuto:

varianza vera: 4.0
media delle stime, divisore n    : 3.19 | teoria (n-1)/n * σ² = 3.2
media delle stime, divisore n - 1: 3.988
MSE divisore n    : 5.764
MSE divisore n - 1: 7.983

Cosa fa il codice. Ripete 200 000 volte l'esperimento «estrai 5 valori da una normale con varianza vera 4 e stima la varianza»: ogni riga di campioni è un campione, var(axis=1, ...) calcola la varianza di ciascuna riga con i due divisori (in NumPy ddof = delta degrees of freedom, il numero sottratto a n). La media delle stime con divisore n cade vicino a 3.2 = (4/5)·4, confermando la distorsione calcolata; con n − 1 cade vicino a 4. Le ultime due righe riportano un dettaglio da non dimenticare: lo stimatore non distorto non è quello con MSE minore. Esiste spesso un compromesso tra distorsione e varianza: ritornerà con forza in tutto il machine learning (è il compromesso tra sotto- e sovra-adattamento).

Massima verosimiglianza

L'analogia: il detective

Vedi sul pavimento certe impronte e devi indovinare il numero di scarpe del sospettato. Per ogni numero di scarpe candidato puoi chiederti: «se il colpevole portasse quel numero, quanto sarebbe probabile lasciare impronte così?» Scegli il numero che rende le impronte osservate più plausibili. Questo è il principio della massima verosimiglianza (Fisher, 1922).

Formalmente. Dato il campione x₁, …, xₙ indipendente con densità (o massa) f(x; θ), la verosimiglianza è, come funzione di θ con i dati fissati,

L(θ) = Π f(xᵢ; θ),

il prodotto perché le osservazioni sono indipendenti. Si massimizza L, o meglio il suo logaritmo ℓ(θ) = Σ ln f(xᵢ; θ) (il logaritmo trasforma il prodotto in somma, non cambia dove sta il massimo e evita che un prodotto di molti numeri piccoli diventi 0 in virgola mobile: lezione 4 del modulo 01). Il massimo di una funzione liscia si trova dove la sua pendenza (la derivata, che definiremo bene nella lezione di analisi) è zero.

Esempio 1: la moneta (Bernoulli)

Osservi k teste in n lanci. La log-verosimiglianza è ℓ(p) = k ln p + (n − k) ln(1 − p). Pendenza zero:

k/p − (n − k)/(1 − p) = 0 ⇒ k(1 − p) = (n − k)p ⇒ p̂ = k/n.

La stima è la frequenza osservata: il metodo giustifica l'intuizione.

Esempio 2: la normale

Per N(μ, σ²) la log-verosimiglianza è ℓ = −(n/2) ln(2πσ²) − Σ(xᵢ−μ)²/(2σ²). Rispetto a μ, annullare la pendenza dà Σ(xᵢ − μ) = 0, quindi μ̂ = x̄. Rispetto a σ², dà σ̂² = (1/n) Σ(xᵢ − x̄)²: lo stimatore con divisore n, quello distorto di prima. La massima verosimiglianza non garantisce uno stimatore non distorto, ma garantisce buone proprietà per n grande (consistenza).

from scipy.optimize import minimize_scalar

k, n = 62, 100
griglia = np.linspace(0.01, 0.99, 981)
loglik = k * np.log(griglia) + (n - k) * np.log(1 - griglia)
print("massimo sulla griglia in p =", round(float(griglia[np.argmax(loglik)]), 3), "| formula k/n =", k / n)

res = minimize_scalar(lambda p: -(k * np.log(p) + (n - k) * np.log(1 - p)), bounds=(0.001, 0.999), method="bounded")
print("ottimizzatore numerico:", round(float(res.x), 4))

x = np.random.default_rng(2).normal(5.0, 2.0, size=300)
print("μ̂ (media)           :", round(float(x.mean()), 4))
print("σ̂² con divisore n   :", round(float(x.var(ddof=0)), 4), "| con n-1:", round(float(x.var(ddof=1)), 4))

Output ottenuto:

massimo sulla griglia in p = 0.62 | formula k/n = 0.62
ottimizzatore numerico: 0.62
μ̂ (media)           : 4.893
σ̂² con divisore n   : 4.1309 | con n-1: 4.1447

Cosa fa il codice. Calcola la log-verosimiglianza della moneta su una griglia di valori di p e ne cerca il massimo: cade su 0.62 = 62/100, come la formula. Poi lo stesso problema è risolto da un ottimizzatore numerico (si minimizza −ℓ, perché gli ottimizzatori cercano minimi): è come si lavora quando una formula chiusa non esiste, ad esempio per la regressione logistica. L'ultima parte mostra le stime di un campione normale: la media coincide con μ̂, e le due varianze differiscono di un fattore n/(n−1) (quasi nulla a n = 300: il problema è serio per campioni piccoli).

Intervalli di confidenza

Una stima puntuale (x̄ = 4.87) non dice quanto vale. Un intervallo di confidenza al 95% è un intervallo [a, b] costruito con una procedura tale che, ripetendola su molti campioni, il 95% degli intervalli ottenuti contiene il parametro vero.

Attenzione alla lettura. Il 95% descrive la procedura, non il singolo intervallo: una volta calcolato, [a, b] contiene μ oppure no (nel disegno, 19 su 20 sì). Dire «μ ha il 95% di probabilità di stare in [a, b]» è, nell'impostazione frequentista di questo corso, una scorrettezza di linguaggio.

Costruzione per la media

Per il TLC (lezione 4), x̄ è circa normale con deviazione standard σ/√n: quindi P(|x̄ − μ| ≤ 1.96 σ/√n) ≈ 0.95, e l'intervallo è x̄ ± 1.96 · σ/√n. Ma σ è quasi sempre sconosciuta: si sostituisce con s. Questa sostituzione aggiunge incertezza (specie con n piccolo), e William Gosset, firmando «Student» (1908), ha mostrato che la quantità (x̄ − μ)/(s/√n) segue una distribuzione con code più pesanti della normale, la t di Student con n − 1 gradi di libertà. L'intervallo corretto è x̄ ± t · s/√n, con t il quantile della t.

from scipy import stats

rng = np.random.default_rng(3)
n, prove, mu = 10, 40_000, 0.0
campioni = rng.normal(mu, 1.0, size=(prove, n))
medie = campioni.mean(axis=1)
s = campioni.std(axis=1, ddof=1)
se = s / np.sqrt(n)

z = 1.96
t = stats.t.ppf(0.975, df=n - 1)
copertura_z = (np.abs(medie - mu) <= z * se).mean()
copertura_t = (np.abs(medie - mu) <= t * se).mean()
print(f"quantile t con {n - 1} gradi di libertà: {t:.3f} (contro 1.960 della normale)")
print(f"copertura con 1.96 (sbagliato per n piccolo): {copertura_z:.3f}")
print(f"copertura con t di Student                 : {copertura_t:.3f}")

Output ottenuto:

quantile t con 9 gradi di libertà: 2.262 (contro 1.960 della normale)
copertura con 1.96 (sbagliato per n piccolo): 0.919
copertura con t di Student                 : 0.951

Cosa fa il codice. Ripete 40 000 volte l'esperimento «campione di 10 osservazioni da una normale con μ = 0», calcola per ciascuna la media e l'errore standard s/√n e verifica in quanti casi l'intervallo contiene μ: la copertura. Con il moltiplicatore normale 1.96 si copre meno del 95% promesso (circa il 92%): l'intervallo è troppo stretto perché ignora l'incertezza di s. Con il quantile della t (2.262 per 9 gradi di libertà) la copertura torna al 95%. Per campioni grandi i due quantili coincidono, per quelli piccoli no.

Il bootstrap: quando non c'è la formula

Come costruisci un intervallo per la mediana, o per un rapporto, o per una metrica di un modello? Per molte statistiche non esiste una formula semplice. Il bootstrap (Efron, 1979) offre una ricetta universale:

1. Hai un campione di n osservazioni. Trattalo come se fosse la popolazione.
2. Estrai con ripetizione n osservazioni da esso: un campione bootstrap.
3. Calcola su di esso la statistica che ti interessa.
4. Ripeti migliaia di volte: ottieni una distribuzione della statistica.
5. L'intervallo percentile al 95% è fra il 2.5° e il 97.5° percentile di quella distribuzione.

L'idea è che il campione è la miglior descrizione della popolazione che hai: ricampionarlo simula «cosa sarebbe successo con altri campioni». L'ordine dei cinque passi è obbligato, perché ciascuno usa il risultato del precedente, e il passo 2 è quello che fa la differenza: estraendo con ripetizione ogni campione bootstrap ha la stessa dimensione dell'originale ma una composizione diversa, ed è la variabilità fra queste composizioni a stimare l'incertezza.

rng = np.random.default_rng(4)
campione = rng.lognormal(7.4, 0.5, size=60)                  # 60 redditi sintetici

mediane = np.array([np.median(rng.choice(campione, size=campione.size, replace=True)) for _ in range(10_000)])
lo, hi = np.percentile(mediane, [2.5, 97.5])
print(f"mediana del campione: {np.median(campione):.0f}")
print(f"intervallo bootstrap 95 % per la mediana: [{lo:.0f}, {hi:.0f}]")

# frazione di osservazioni distinte in un campione bootstrap
n = campione.size
distinti = np.mean([len(np.unique(rng.integers(0, n, n))) / n for _ in range(5_000)])
print(f"frazione media di osservazioni distinte: {distinti:.4f} | 1 - 1/e = {1 - np.exp(-1):.4f}")

Output ottenuto:

mediana del campione: 1515
intervallo bootstrap 95 % per la mediana: [1361, 1789]
frazione media di osservazioni distinte: 0.6359 | 1 - 1/e = 0.6321

Cosa fa il codice. rng.choice(campione, size=n, replace=True) è il campione bootstrap (estrazione con ripetizione); una list comprehension ne produce 10 000 e ne calcola la mediana. I percentili 2.5 e 97.5 delimitano l'intervallo. L'ultima parte verifica un fatto elegante: in un campione bootstrap la frazione di osservazioni distinte è circa il 63.2%. Perché: la probabilità che una data osservazione non venga mai estratta in n estrazioni è (1 − 1/n)ⁿ, che per n grande tende a e⁻¹ ≈ 0.368 (un limite notevole che vedremo nella lezione di analisi); quindi sono presenti circa 1 − 0.368 = 63.2% di osservazioni distinte, le altre sono ripetute. Questa quota è quella delle osservazioni «fuori sacco» (out-of-bag) usate dalle foreste casuali per valutarsi.

Trade-off. Il bootstrap è potente ma non magico: con campioni molto piccoli, o per statistiche sensibili agli estremi (il massimo), può funzionare male; e non corregge un campione non rappresentativo.

Test di ipotesi e p-value

Hai già incontrato l'idea nel modulo 00 (lezione 6, A/B test). La ripetiamo con una domanda diversa: due gruppi hanno la stessa media?

• Ipotesi nulla H₀: «nessuna differenza» (le medie vere sono uguali).
• Statistica del test: un numero che misura quanto i dati si discostano da H₀. Per due medie: t = (x̄₁ − x̄₂) / √(s₁²/n₁ + s₂²/n₂) (test di Welch).
• p-value: la probabilità, se H₀ fosse vera, di osservare una statistica estrema almeno quanto quella ottenuta. Non è la probabilità che H₀ sia vera.

Le tre voci sono ordinate dalla più semplice alla più fraintesa. L'ipotesi nulla è una dichiarazione che si può rifiutare ma non dimostrare, la statistica del test è un numero che si calcola dai dati, e il p-value è la probabilità di quel numero (o di uno più estremo) calcolata supponendo vera H₀: non risponde alla domanda «quanto è probabile che io abbia ragione».

rng = np.random.default_rng(5)
a = rng.normal(100, 15, size=50)                     # gruppo di controllo
b = rng.normal(106, 15, size=50)                     # gruppo trattato: differenza vera di 6 punti

t_stat, p = stats.ttest_ind(b, a, equal_var=False)   # test di Welch
diff = b.mean() - a.mean()
se_diff = np.sqrt(a.var(ddof=1) / 50 + b.var(ddof=1) / 50)
d_cohen = diff / np.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
print(f"differenza osservata: {diff:.2f} | errore standard: {se_diff:.2f}")
print(f"statistica t = {t_stat:.3f} | p-value = {p:.4f}")
print(f"effetto standardizzato (d di Cohen) = {d_cohen:.2f}")
print(f"intervallo di confidenza 95 % della differenza: [{diff - 1.98 * se_diff:.2f}, {diff + 1.98 * se_diff:.2f}]")

Output ottenuto:

differenza osservata: 9.02 | errore standard: 2.66
statistica t = 3.393 | p-value = 0.0010
effetto standardizzato (d di Cohen) = 0.68
intervallo di confidenza 95 % della differenza: [3.76, 14.28]

Cosa fa il codice. Costruisce due gruppi con una differenza vera di 6 punti (nota perché la abbiamo inserita noi), calcola il test di Welch (equal_var=False non assume varianze uguali) e riporta, oltre al p-value, la differenza con il suo errore standard, l'effetto standardizzato d (la differenza in unità di deviazione standard, indipendente dalla scala) e un intervallo di confidenza (con il moltiplicatore 1.98 approssimato per circa 98 gradi di libertà). Nota che la differenza osservata (circa 9) non è esattamente quella vera (6): è la variabilità del campione, e l'intervallo di confidenza, infatti, contiene 6. Questo è il modo corretto di riportare un risultato: il p-value dice se la differenza è distinguibile dal rumore, d e l'intervallo dicono quanto è grande e quanto è incerta.

Che cosa fa davvero un p-value, per simulazione

rng = np.random.default_rng(6)
p_values = np.array([stats.ttest_ind(rng.normal(0, 1, 30), rng.normal(0, 1, 30)).pvalue for _ in range(5_000)])
print("H0 vera, 5000 test: frazione con p < 0.05 =", round(float((p_values < 0.05).mean()), 4))
print("decili dei p-value (uniformi se H0 è vera):", np.round(np.percentile(p_values, np.arange(10, 100, 10)), 2))

Output ottenuto:

H0 vera, 5000 test: frazione con p < 0.05 = 0.0498
decili dei p-value (uniformi se H0 è vera): [0.11 0.2  0.3  0.4  0.5  0.59 0.7  0.8  0.9 ]

Cosa fa il codice. Ripete 5000 volte un test su due gruppi estratti dalla stessa popolazione (H₀ vera per costruzione). Il risultato mostra la proprietà più importante dei p-value: se H₀ è vera, sono uniformemente distribuiti tra 0 e 1, quindi il 5% sta sotto 0.05 per puro caso. Questi sono i falsi positivi. Se fai 20 test su dati senza alcun effetto, ne aspetti circa uno «significativo» (1 − 0.95²⁰ ≈ 64% di averne almeno uno). La dichiarazione dell'American Statistical Association sui p-value (2016) ricorda che il p-value non misura la probabilità che l'ipotesi sia vera né l'importanza dell'effetto, e che le conclusioni non vanno basate su una soglia sola.

Regressione lineare ai minimi quadrati

L'analogia: l'elastico

Hai una nuvola di punti (xᵢ, yᵢ) e vuoi la retta ŷ = a + b·x che meglio li descrive. Immagina di collegare ciascun punto alla retta con una molla verticale: la retta si assesta dove la tensione complessiva è minima. La tensione di una molla cresce col quadrato dell'allungamento, quindi la retta di equilibrio minimizza la somma dei quadrati dei residui rᵢ = yᵢ − (a + b·xᵢ). È il metodo dei minimi quadrati (Legendre 1805, Gauss 1809).

Derivazione

Vogliamo minimizzare S(a, b) = Σ (yᵢ − a − b·xᵢ)². Una funzione di due variabili ha il minimo dove le derivate parziali (la pendenza lungo a tenendo ferma b, e viceversa; le definiamo nella lezione di analisi) sono entrambe zero.

Rispetto ad a: −2 Σ (yᵢ − a − b·xᵢ) = 0, cioè n·ȳ − n·a − b·n·x̄ = 0, quindi a = ȳ − b·x̄: la retta passa sempre per il punto delle medie (x̄, ȳ).

Rispetto a b: −2 Σ xᵢ (yᵢ − a − b·xᵢ) = 0. Sostituendo a = ȳ − b x̄ e riordinando: Σ xᵢ (yᵢ − ȳ) = b · Σ xᵢ (xᵢ − x̄). Poiché Σ x̄ (yᵢ − ȳ) = 0 e Σ x̄ (xᵢ − x̄) = 0, si può sostituire xᵢ con (xᵢ − x̄), e si arriva a:

b = Σ (xᵢ − x̄)(yᵢ − ȳ) / Σ (xᵢ − x̄)² = Cov(x, y) / Var(x).

La pendenza è il rapporto tra covarianza e varianza di x (quella campionaria, con lo stesso divisore al numeratore e al denominatore). Si definisce anche il coefficiente di determinazione R² = 1 − Σ rᵢ² / Σ (yᵢ − ȳ)²: la quota della variabilità di y che la retta spiega; per la regressione semplice coincide con il quadrato della correlazione.

rng = np.random.default_rng(8)
x = rng.uniform(0, 10, size=80)
y = 2.0 + 0.5 * x + rng.normal(0, 1.0, size=80)         # retta vera: a = 2, b = 0.5, rumore σ = 1

b = ((x - x.mean()) * (y - y.mean())).sum() / ((x - x.mean()) ** 2).sum()
a = y.mean() - b * x.mean()
residui = y - (a + b * x)
r2 = 1 - (residui ** 2).sum() / ((y - y.mean()) ** 2).sum()
print(f"a = {a:.3f}, b = {b:.3f}, R² = {r2:.3f}")
print("b = cov/var:", round(float(np.cov(x, y)[0, 1] / np.var(x, ddof=1)), 3), "| R² = corr²:", round(float(np.corrcoef(x, y)[0, 1] ** 2), 3))
print("controllo np.polyfit:", np.round(np.polyfit(x, y, 1), 3), "| scipy linregress:", round(float(stats.linregress(x, y).slope), 3))
print("somma dei residui (≈ 0):", round(float(residui.sum()), 10))

# incertezza sulla pendenza con il bootstrap sulle coppie (x, y)
pendenze = []
for _ in range(5_000):
    idx = rng.integers(0, x.size, x.size)
    xb, yb = x[idx], y[idx]
    pendenze.append(((xb - xb.mean()) * (yb - yb.mean())).sum() / ((xb - xb.mean()) ** 2).sum())
print("intervallo bootstrap 95 % della pendenza:", np.round(np.percentile(pendenze, [2.5, 97.5]), 3))

Output ottenuto:

a = 1.849, b = 0.524, R² = 0.664
b = cov/var: 0.524 | R² = corr²: 0.664
controllo np.polyfit: [0.524 1.849] | scipy linregress: 0.524
somma dei residui (≈ 0): 0.0
intervallo bootstrap 95 % della pendenza: [0.45  0.604]

Cosa fa il codice. Genera 80 punti da una retta nota (a = 2, b = 0.5) più rumore, e applica le formule derivate. Le stime sono vicine ai valori veri ma non identiche (è l'incertezza del campione). Poi verifica la stessa pendenza con np.cov/np.var (covarianza su varianza), R² come quadrato della correlazione, np.polyfit e scipy.stats.linregress. La somma dei residui è zero (conseguenza della derivata rispetto ad a). L'ultimo blocco ricampiona le coppie (x, y) con ripetizione e ricalcola la pendenza, ottenendo un intervallo per la pendenza: se la procedura è corretta, lo contiene in circa il 95% degli esperimenti (in questo contiene la pendenza vera, 0.5; in un esperimento su venti non succederebbe).

Trade-off. I minimi quadrati sono molto sensibili agli outlier (il quadrato amplifica gli scarti grandi); assumono una relazione lineare, residui con dispersione costante e indipendenti. Sempre: disegna i residui (lezione sulla visualizzazione: il quartetto di Anscombe ha la stessa retta e R² per quattro situazioni molto diverse).

Correlazione e causalità

Due grandezze possono essere correlate perché una causa l'altra, perché una causa l'altra al contrario, per caso, o perché una terza le influenza entrambe (un confondente). Simuliamo il terzo caso: la temperatura fa aumentare sia i gelati venduti sia le persone che fanno il bagno (e quindi, purtroppo, certi incidenti).

rng = np.random.default_rng(9)
temp = rng.normal(25, 5, size=5_000)
gelati = 10 + 2.0 * temp + rng.normal(0, 4, 5_000)
bagnanti = -5 + 1.5 * temp + rng.normal(0, 4, 5_000)

print("corr(gelati, bagnanti) =", round(float(np.corrcoef(gelati, bagnanti)[0, 1]), 3))

def residuo(y, x):
    pend, interc = np.polyfit(x, y, 1)
    return y - (interc + pend * x)

print("corr dopo aver tolto l'effetto della temperatura =", round(float(np.corrcoef(residuo(gelati, temp), residuo(bagnanti, temp))[0, 1]), 3))

Output ottenuto:

corr(gelati, bagnanti) = 0.822
corr dopo aver tolto l'effetto della temperatura = -0.007

Cosa fa il codice. Costruisce due variabili che dipendono soltanto dalla temperatura (più rumore indipendente): non c'è alcun legame diretto tra loro, ma la correlazione è alta. Poi calcola i residui di ciascuna rispetto alla temperatura (ciò che resta dopo aver tolto l'effetto lineare della temperatura) e correla i residui: la correlazione quasi sparisce. Questa è la logica del «controllo per una variabile» e il motivo per cui gli esperimenti randomizzati (lezione 6 del modulo 00) sono il metodo più solido per parlare di cause.

Errori tipici (sintomo, causa, rimedio)

• Riassumere con la sola media. Rimedio: mediana, IQR e un istogramma.
• Divisore n per la varianza di un campione. Rimedio: ddof=1.
• Leggere l'intervallo di confidenza come probabilità sul parametro. Rimedio: «procedura che copre il 95% delle volte».
• Usare 1.96 con campioni piccoli. Rimedio: quantili della t.
• «p < 0.05 quindi l'effetto è grande/vero». Rimedio: riportare effetto e intervallo; sapere che con H₀ vera il 5% dei test è «significativo».
• Test multipli senza correzione. Rimedio: una metrica primaria dichiarata, oppure correzione di Bonferroni (α/m).
• Bootstrap su dati dipendenti (serie temporali, righe dello stesso soggetto): ricampionare le singole righe spezza la struttura. Rimedio: bootstrap a blocchi o per gruppo.
• Estrapolare la retta fuori dall'intervallo dei dati. Rimedio: dichiarare il dominio di validità.
• Confondere correlazione e causalità. Rimedio: confondenti, esperimenti.

Gli errori si raggruppano in tre famiglie: quelli di descrizione (la sola media, il divisore sbagliato), quelli di lettura (il 95%, il p-value, 1.96 su pochi dati) e quelli di metodo (test multipli, bootstrap su dati dipendenti, estrapolazione, causalità). Per le ultime due famiglie non c'è una riga di codice che protegga: serve dichiarare prima del calcolo che cosa si vuole concludere e con quali ipotesi.

Come lavorare nella pratica

Descrivi prima, con una tabella di riassunti robusti e un grafico. Riporta sempre stima, incertezza e dimensione del campione, e preferisci gli intervalli ai soli p-value, dichiarando in anticipo il test e la soglia. Simula per verificare coperture e falsi positivi, e documenta le ipotesi (indipendenza, normalità approssimata, varianze) controllandole sui dati.

Il lessico si riassume in poche righe. Popolazione e campione sono il tutto e la parte osservata; parametro, statistica e stimatore sono il numero vero ignoto, il numero calcolato dai dati e la formula che lo produce. La distorsione è l'errore medio sistematico di uno stimatore, l'errore standard la sua deviazione standard (per la media, σ/√n) e i gradi di libertà sono i valori liberi di variare dopo i vincoli imposti dalle stime. La verosimiglianza è la probabilità dei dati vista come funzione del parametro; l'intervallo di confidenza nasce da una procedura che copre il parametro con frequenza dichiarata; il residuo è lo scarto fra osservato e previsto; il confondente è una terza variabile che influenza entrambe le grandezze.

Mini-progetto

1. Simula 10 000 campioni da una Poisson con λ = 4 e n = 15. Per ciascuno calcola un intervallo di confidenza al 95% per λ con x̄ ± 1.96·√(x̄/n) e misura la copertura. Cambia n in 5, 15, 100: come cambia?
2. Con il bootstrap, stima un intervallo al 95% per la differenza tra le mediane dei due gruppi del test di Welch e confrontalo con quello per le medie.
3. Costruisci un esempio in cui R² della retta è alto ma i residui mostrano una curvatura, e disegnali.

Il primo punto è la verifica del criterio «copertura»: l'intervallo di Poisson usa un'approssimazione normale e aspettarsi che funzioni peggio per n piccolo è la parte da misurare, non da assumere. Il secondo riusa il codice del bootstrap con una statistica diversa, il terzo mostra perché un R² alto non basta a dichiarare lineare una relazione.

Challenge

• Easy. Calcola con la formula e con np.polyfit la retta per i punti (1,2), (2,3), (3,5), (4,4) e verifica che la somma dei residui sia zero.
• Medium. Dimostra che MSE = distorsione² + varianza aggiungendo e togliendo E[θ̂].
• Hard. Dimostra che per la regressione semplice R² = r² (il quadrato del coefficiente di correlazione).
• Expert. Ricava lo stimatore di massima verosimiglianza di λ per una Poisson (ℓ(λ) = Σ(xᵢ ln λ − λ − ln xᵢ!)) e dimostra che è x̄. È non distorto? Qual è il suo errore standard?
• Enterprise. Prendi una metrica di una tua dashboard (conversione, tempo di risposta). Scrivi stima, intervallo bootstrap al 95%, e una frase su cosa non puoi concludere con quel campione.

Le prime quattro sfide sono dimostrazioni e calcoli che puoi controllare con il codice della lezione; la quinta è la più formativa perché ti obbliga a scrivere la frase conclusiva con i suoi limiti. Un risultato statistico non è completo finché non dice anche che cosa il campione non permette di affermare.

Che cosa sai fare adesso

Sai riassumere un campione con indicatori robusti, spiegare e verificare per simulazione perché la varianza campionaria divide per n − 1, stimare un parametro con la massima verosimiglianza, costruire un intervallo di confidenza con la t di Student o con il bootstrap e misurarne la copertura. Sai leggere un p-value per quello che è, derivare la retta ai minimi quadrati e riconoscere un confondente. Il passo successivo è l'analisi matematica: la derivata che qui abbiamo usato «a occhio» per trovare i massimi della verosimiglianza e i minimi dei minimi quadrati; la regressione con molte variabili ritorna in Regressione lineare.

Nella stessa serie trovi anche: 15 - scikit-learn: stimatori, pipeline e validazione, 17 - Calcolo combinatorio: contare senza elencare, 24 - Algebra lineare 2: autovalori, SVD e PCA, 25 - Calcolo matriciale e differenziazione automatica: dal Jacobiano a un autograd scritto da zero, 28 - Statistica bayesiana: prior, posterior, coniugati e perché la regolarizzazione è una credenza, 29 - Campionamento e metodi Monte Carlo: stimare quello che non si riesce a calcolare, 74 - Reinforcement learning: l'agente che impara per tentativi (e il dilemma esplorazione-sfruttamento) (che richiamano questa lezione).