06 - Inferenza statistica: stimatori e proprietà

Hai misurato venti volte la stessa grandezza e vuoi dire di quanto le misure si disperdono. La formula della varianza campionaria ha un denominatore, e il libro dice 19, mentre la calcolatrice ha due tasti, uno che divide per 20 e uno per 19. La risposta breve è che nessuno dei due è «quello giusto»: con 19 lo stimatore è non distorto, con 20 coincide con quello di massima verosimiglianza, e con 21 ha l'errore quadratico medio più piccolo fra i tre, come calcolo qui sotto. Scegliere vuol dire decidere quale proprietà ti interessa.

Questo articolo, il sesto della serie Probabilità e Statistica, ti dà il vocabolario per farlo: che cos'è uno stimatore, come se ne misura la qualità, quanto bene si può fare al massimo e come si costruisce uno stimatore buono. Servono variabile aleatoria, valore atteso e varianza, e conviene aver letto l'articolo su teorema centrale del limite e legge dei grandi numeri.

Parametro, stimatore, stima

Immagina di pesare venti volte lo stesso oggetto con una bilancia imprecisa. Ammetti che ogni pesata X_i sia una variabile aleatoria Normale con media μ (il peso vero) e varianza σ² (l'imprecisione della bilancia), e che le pesate siano indipendenti. I numeri μ e σ² sono i parametri: descrivono la distribuzione e tu non li conosci. Una statistica è una qualunque funzione delle venti pesate, come la media o la più grande. Uno stimatore è una statistica scelta per avvicinarsi a un parametro; la stima è il numero che lo stimatore assume sul campione che hai davvero in mano.

La distinzione conta perché lo stimatore, prima di raccogliere i dati, è una variabile aleatoria: ripetendo le venti pesate otterresti ogni volta una media diversa. La sua distribuzione si chiama distribuzione campionaria, ed è lei, non la singola stima, che permette di giudicare lo stimatore.

Se chiamo T uno stimatore di un parametro θ, la sua qualità si descrive con tre numeri: il bias (distorsione) b = E[T] − θ, cioè di quanto sbaglia in media; la varianza Var(T), cioè quanto oscilla da campione a campione; e l'errore quadratico medio (MSE, mean squared error) E[(T − θ)²], che misura la distanza media al quadrato dal bersaglio.

I tre numeri sono legati da una identità: MSE = Var(T) + b². La dimostrazione è di una riga: scrivi T − θ come (T − E[T]) + (E[T] − θ), eleva al quadrato e prendi il valore atteso; il doppio prodotto vale zero perché E[T − E[T]] = 0. L'identità dice che uno stimatore può compensare un po' di distorsione con una varianza più piccola, ed è il motivo per cui «non distorto» non coincide con «migliore». L'identità, con il nome di «rischio quadratico = bias² + varianza», è nelle lezioni di inferenza parametrica del corso MIT 18.650.

Dividere per n−1, per n o per n+1

Per una Normale con n osservazioni, sia Q = Σ(X_i − X̄)² la somma dei quadrati degli scarti dalla media campionaria. Un risultato classico per i campioni Normali dice che Q/σ² segue una distribuzione chi-quadro con n − 1 gradi di libertà, quindi ha media n − 1 e varianza 2(n − 1). Da qui segue tutto, per lo stimatore Q/d con divisore d: la media è σ²(n − 1)/d e la varianza è 2σ⁴(n − 1)/d².

Con d = n − 1 il bias è nullo; con d = n la media è σ²(n − 1)/n, quindi lo stimatore sottostima sempre di σ²/n. La ragione intuitiva è che gli scarti sono misurati dalla media del campione, che è il punto più vicino possibile ai dati, e non dalla media vera: Q è per costruzione più piccola di quanto sarebbe misurando da μ.

Il programma verifica queste formule con una simulazione: genera 200.000 campioni di n = 20 pesate Normali con σ² = 4, calcola Q per ciascuno e confronta, per i tre divisori, media e MSE ottenuti dalla simulazione con quelli esatti. Il seme fisso rende l'esecuzione ripetibile, e l'errore standard della media dei quadrati degli errori è stampato accanto, perché una simulazione ha senso solo se dichiari di quanto può sbagliare.

# Tre stimatori della varianza: divisore n-1, n, n+1 (campione Normale, n = 20)
import numpy as np

n, sigma2 = 20, 4.0
R = 200_000                      # campioni simulati
rng = np.random.default_rng(42)
X = rng.normal(0.0, np.sqrt(sigma2), size=(R, n))
somma_quadrati = ((X - X.mean(axis=1, keepdims=True)) ** 2).sum(axis=1)

print(f"n = {n}, sigma^2 = {sigma2}, campioni simulati = {R}")
print("divisore  media esatta  media sim.  MSE esatto  MSE sim. (+/- errore std)")
for d in (n - 1, n, n + 1):
    stima = somma_quadrati / d
    # Esatto: somma_quadrati / sigma^2 ~ chi-quadro con n-1 gradi di liberta',
    # quindi media = sigma^2 (n-1)/d e varianza = 2 sigma^4 (n-1)/d^2
    media = sigma2 * (n - 1) / d
    var = 2 * sigma2**2 * (n - 1) / d**2
    mse = var + (media - sigma2) ** 2
    err2 = (stima - sigma2) ** 2
    print(f"{d:8d}  {media:12.4f}  {stima.mean():10.4f}  {mse:10.4f}  "
          f"{err2.mean():8.4f} +/- {err2.std() / np.sqrt(R):.4f}")

# Output:
# n = 20, sigma^2 = 4.0, campioni simulati = 200000
# divisore  media esatta  media sim.  MSE esatto  MSE sim. (+/- errore std)
#       19        4.0000      3.9999      1.6842    1.6896 +/- 0.0061
#       20        3.8000      3.7999      1.5600    1.5649 +/- 0.0052
#       21        3.6190      3.6189      1.5238    1.5283 +/- 0.0046

Leggi l'uscita riga per riga. Con divisore 19 la media è 4,0000, cioè non c'è distorsione, ma l'MSE vale 1,6842. Con divisore 20 la media scende a 3,8, il bias è −0,2 e il suo quadrato è 0,04, ma la varianza passa da 1,6842 a 1,52 e l'MSE complessivo cala a 1,56. Con 21 il bias è −0,381, la varianza è 1,3787 e l'MSE arriva a 1,5238, il minimo dei tre.

Le colonne simulate cadono entro circa un errore standard dai valori esatti. Il guadagno è piccolo (circa il 10% di MSE fra 19 e 21) e si riduce al crescere di n; per questo nella pratica si usa il divisore 19 perché è non distorto, e le varianze di gruppi diversi si possono mediare senza accumulare un errore sistematico, non perché sia il più preciso.

Che cosa succede quando n cresce

Un buon stimatore migliora con più dati. Si dice consistente se, per n che tende all'infinito, converge in probabilità al parametro: la probabilità che sbagli di più di una tolleranza fissata tende a zero. La media campionaria lo è per la legge dei grandi numeri. Una contro-prova utile è lo stimatore «prendo solo la prima pesata, X₁»: è non distorto, ma la sua varianza resta σ² qualunque sia n, quindi non converge a nulla. La consistenza e l'assenza di distorsione sono proprietà indipendenti: lo stimatore con divisore n è distorto e consistente, X₁ è non distorto e non consistente.

La consistenza dice dove va lo stimatore, non quanto si avvicina in fretta. Per questo si studia anche la normalità asintotica: per molti stimatori √n(T_n − θ) si comporta, per n grande, come una Normale di media zero e varianza V(θ).

Per lo stimatore di massima verosimiglianza, sotto condizioni di regolarità, la varianza limite è l'inverso dell'informazione di Fisher, di cui parlo fra poco: il teorema con le condizioni (modello identificabile, supporto che non dipende da θ, parametro non sul bordo dello spazio dei parametri) è enunciato nella lezione sulla massima verosimiglianza del corso MIT 18.650. La normalità asintotica è ciò che permette di costruire intervalli di confidenza, come vedrai nell'articolo sui test di ipotesi.

Le proprietà che si chiedono a uno stimatore sono quattro, e ciascuna risponde a una domanda diversa:

• Non distorto: in media, su molti campioni, centra il parametro?
• Consistente: con più dati si avvicina al parametro?
• Efficiente: ha la varianza più piccola possibile tra quelli della sua classe?
• Sufficiente: conserva tutta l'informazione sul parametro contenuta nel campione?

Nessuna delle quattro implica le altre. Le prime due riguardano il comportamento medio e quello al limite, la terza confronta stimatori fra loro, la quarta riguarda l'informazione e non l'errore. Il resto dell'articolo le mostra una alla volta su un esempio che puoi rifare, cominciando dal limite inferiore alla varianza, che è ciò che dà un senso preciso alla parola «efficiente».

Informazione di Fisher e limite di Cramér-Rao

Quanta informazione su θ porta una osservazione? Si risponde guardando la log-verosimiglianza, cioè il logaritmo della densità (o della probabilità) vista come funzione di θ. La sua derivata rispetto a θ si chiama score, u(x, θ). L'informazione di Fisher è I(θ) = Var(u) = −E[u′]: se la log-verosimiglianza è a punta (la derivata cambia in fretta), i dati distinguono bene un valore di θ dai vicini e I è grande. La definizione, con la sua equivalenza fra le due forme, è nella stessa lezione MIT sulla massima verosimiglianza.

Per una Bernoulli di parametro p, la log-probabilità di x ∈ {0, 1} è x log p + (1 − x) log(1 − p), lo score è x/p − (1 − x)/(1 − p), e −u′ = x/p² + (1 − x)/(1 − p)². Il valore atteso è 1/p + 1/(1 − p) = 1/(p(1 − p)). L'informazione è minima per p = 1/2, dove la moneta è più imprevedibile, e cresce senza limite verso 0 e 1.

Il limite di Cramér-Rao collega informazione e varianza: per ogni stimatore non distorto di θ basato su n osservazioni indipendenti, e sotto condizioni di regolarità, Var(T) ≥ 1/(n I(θ)). La condizione più semplice da controllare è che il supporto della distribuzione non dipenda da θ: la Uniforme su [0, θ] non la soddisfa, e per lei il limite non vale. Per la Bernoulli il limite è p(1 − p)/n, che è proprio la varianza della frequenza campionaria: quello stimatore è efficiente.

Per la distribuzione esponenziale di parametro λ (il tasso, con media 1/λ) l'informazione è I(λ) = 1/λ², e il limite è λ²/n. Qui il quadro è meno ordinato di quello dei libri più sbrigativi. Chiamo T la somma delle n osservazioni; T segue una distribuzione Gamma e 1/T una Gamma inversa, di cui scipy.stats.invgamma restituisce i momenti esatti (la documentazione di SciPy ne dà la definizione).

Gli stimatori della forma c/T comprendono il massimo verosimile (c = n, cioè 1/X̄) e quello non distorto (c = n − 1). Il programma calcola per entrambi, e per c = n − 2, media, varianza e MSE esatti, e li confronta con una simulazione.

# Cramer-Rao per Exp(lambda): il limite non e' raggiunto a n finito, nemmeno dallo stimatore non distorto
import numpy as np
from scipy import stats

lam, n, R = 2.5, 20, 200_000
rng = np.random.default_rng(100)
X = rng.exponential(scale=1 / lam, size=(R, n))
T = X.sum(axis=1)                      # statistica sufficiente: T ~ Gamma(n, scala 1/lam)

# Informazione di Fisher di una osservazione: varianza dello score u(x) = 1/lam - x
I_sim = (1 / lam - X[:, 0]).var()
print(f"I(lambda): teorica {1 / lam**2:.4f}, dallo score simulato {I_sim:.4f}")
bound = 1 / (n * (1 / lam**2))
print(f"limite di Cramer-Rao per non distorti: {bound:.4f}")

# 1/T segue una Gamma inversa(n, scala lam): momenti esatti con scipy
inv = stats.invgamma(a=n, scale=lam)
m1, m2 = inv.mean(), inv.var() + inv.mean() ** 2
print("stimatore c/T  media esatta  var esatta  MSE esatto  var sim.   media sim.")
for c, nome in ((n, "MLE c=n"), (n - 1, "non distorto c=n-1"), (n - 2, "MSE minimo c=n-2")):
    media = c * m1
    var = c**2 * (m2 - m1**2)
    mse = var + (media - lam) ** 2
    s = c / T
    print(f"{nome:20s} {media:9.4f} {var:11.4f} {mse:11.4f} {s.var():9.4f} {s.mean():10.4f}")

# Output:
# I(lambda): teorica 0.1600, dallo score simulato 0.1586
# limite di Cramer-Rao per non distorti: 0.3125
# stimatore c/T  media esatta  var esatta  MSE esatto  var sim.   media sim.
# MLE c=n                 2.6316      0.3847      0.4020    0.3862     2.6304
# non distorto c=n-1      2.5000      0.3472      0.3472    0.3485     2.4988
# MSE minimo c=n-2        2.3684      0.3116      0.3289    0.3128     2.3673

Il limite per n = 20 e λ = 2,5 è 0,3125. Lo stimatore di massima verosimiglianza ha media 2,6316, quindi un bias di +0,1316, e varianza 0,3847: il limite non lo riguarda, perché non è non distorto. Quello non distorto (c = 19) ha varianza 0,3472, più alta del limite: è il migliore fra i non distorti, ma a n finito non raggiunge il limite; la varianza vale λ²/(n − 2), e il rapporto col limite, n/(n − 2), tende a 1 solo per n grande.

Con c = 18 la varianza (0,3116) scende persino sotto il limite, e non c'è contraddizione, perché lo stimatore è distorto, e il suo MSE (0,3289) è il più basso dei tre. La riga «I(lambda)» mostra anche la simulazione dello score: 0,1586 contro 0,16 teorico, con errore standard di circa 0,001.

Dall'informazione all'intervallo di confidenza

Se lo stimatore di massima verosimiglianza è approssimativamente Normale con varianza 1/(n I(θ)), allora un intervallo che va da T − z/√(n I) a T + z/√(n I), con z = 1,96, dovrebbe contenere θ nel 95% dei campioni. Con I stimata sui dati è l'intervallo di Wald.

Per una proporzione p diventa f ± 1,96 √(f(1 − f)/n), con f la frequenza campionaria. «Dovrebbe» è la parola chiave: il programma calcola la copertura vera, cioè la quota di campioni per cui l'intervallo contiene davvero p = 0,10, in due modi. Con la formula esatta, sommando le probabilità binomiali dei valori di k per cui l'intervallo contiene p, e con 200.000 campioni simulati per ciascun n.

# Copertura reale dell'intervallo di Wald al 95% per una proporzione p (Bernoulli)
import numpy as np
from scipy import stats

p, R, z = 0.10, 200_000, stats.norm.ppf(0.975)
rng = np.random.default_rng(2024)
print(f"p = {p}, z = {z:.4f}, campioni per riga = {R}")
print("   n  copertura esatta  copertura simulata (+/- err. std)")
for n in (20, 50, 100, 500, 2000):
    k = rng.binomial(n, p, size=R)          # numero di successi in ciascun campione
    ph = k / n
    se = np.sqrt(ph * (1 - ph) / n)         # 1/sqrt(n I(p)) con I stimata in ph
    copre = ((ph - z * se <= p) & (p <= ph + z * se)).mean()
    # Esatta: somma delle probabilita' binomiali dei k per cui l'intervallo contiene p
    kk = np.arange(n + 1)
    pk = kk / n
    sek = np.sqrt(pk * (1 - pk) / n)
    contiene = (pk - z * sek <= p) & (p <= pk + z * sek)
    esatta = stats.binom.pmf(kk, n, p)[contiene].sum()
    err = np.sqrt(copre * (1 - copre) / R)
    print(f"{n:4d}  {esatta:16.4f}  {copre:12.4f} +/- {err:.4f}")

# Output:
# p = 0.1, z = 1.9600, campioni per riga = 200000
#    n  copertura esatta  copertura simulata (+/- err. std)
#   20            0.8760        0.8766 +/- 0.0007
#   50            0.8789        0.8783 +/- 0.0007
#  100            0.9324        0.9325 +/- 0.0006
#  500            0.9427        0.9433 +/- 0.0005
# 2000            0.9462        0.9463 +/- 0.0005

Con venti osservazioni l'intervallo «al 95%» copre il parametro nell'87,6% dei campioni, e con cinquanta nell'87,9%; ci vogliono circa 2000 osservazioni per arrivare vicino al 95% (94,6%). La simulazione conferma la formula esatta entro un errore standard di 0,0005-0,0007. Una proprietà asintotica, cioè vera al limite, non garantisce nulla per un n finito: quando p è vicino a 0 la Normale è un cattivo modello della distribuzione di f, e quando k = 0 l'intervallo collassa in un punto che non contiene mai p.

La lettura corretta del «95%» è questa: è la frequenza con cui la procedura contiene il parametro nel lungo periodo, e la tabella mostra che a volte è molto minore del dichiarato. Una volta calcolato, l'intervallo contiene p oppure no: non ha senso dire che lo contiene «con probabilità 95%».

Sufficienza e Rao-Blackwell

Una statistica T è sufficiente per θ se, una volta noto il valore di T, la distribuzione del campione non dipende più da θ: il resto dei dati è rumore che non dice niente sul parametro. Per un campione di Bernoulli, la somma dei successi T = ΣX_i è sufficiente, perché data la somma tutte le disposizioni dei successi sono equiprobabili qualunque sia p. Si riconosce con il teorema di fattorizzazione: T è sufficiente quando la densità congiunta si scrive come prodotto di un termine che dipende dai dati solo attraverso T (e da θ) e di uno che non dipende da θ.

Il teorema di Rao-Blackwell dà un uso pratico alla sufficienza: se W è un qualunque stimatore non distorto e T è sufficiente, lo stimatore E[W | T] è ancora non distorto e ha varianza minore o uguale. Condizionare a T equivale a fare la media su tutti i campioni che hanno lo stesso valore di T, e la media non può aumentare la varianza.

Nell'esempio Bernoulli, parti da W = X₁ (la prima moneta sola) e condizioni alla somma: per simmetria E[X₁ | T = t] = t/n, quindi lo stimatore migliorato è la frequenza campionaria. Il programma lo verifica senza fidarsi della simmetria: fra i campioni con somma 3, 6 e 10 calcola la media di X₁ e la confronta con t/n.

# Rao-Blackwell su Bernoulli(p): da X_1 alla media campionaria, condizionando alla somma
import numpy as np

p, n, R = 0.3, 20, 1_000_000
rng = np.random.default_rng(3)
X = (rng.random(size=(R, n)) < p)
W = X[:, 0].astype(float)        # stimatore non distorto ma rumoroso: usa una sola osservazione
T = X.sum(axis=1)                # statistica sufficiente: numero di successi
RB = T / n                       # E[W | T] = T/n

print(f"W = X_1     : media {W.mean():.4f} (+/- {W.std() / np.sqrt(R):.4f}), varianza {W.var():.4f} (esatta {p * (1 - p):.4f})")
print(f"E[W|T] = T/n: media {RB.mean():.4f} (+/- {RB.std() / np.sqrt(R):.4f}), varianza {RB.var():.5f} (esatta {p * (1 - p) / n:.5f})")
# Controllo diretto di E[X_1 | T = t] = t/n, sui campioni con somma t
for t in (3, 6, 10):
    sel = T == t
    err = W[sel].std() / np.sqrt(sel.sum())
    print(f"T = {t:2d}: campioni {sel.sum():6d}, media di X_1 {W[sel].mean():.4f} (+/- {err:.4f}), t/n = {t / n:.4f}")

# Output:
# W = X_1     : media 0.3000 (+/- 0.0005), varianza 0.2100 (esatta 0.2100)
# E[W|T] = T/n: media 0.2999 (+/- 0.0001), varianza 0.01051 (esatta 0.01050)
# T =  3: campioni  71963, media di X_1 0.1490 (+/- 0.0013), t/n = 0.1500
# T =  6: campioni 191630, media di X_1 0.2996 (+/- 0.0010), t/n = 0.3000
# T = 10: campioni  30833, media di X_1 0.5016 (+/- 0.0028), t/n = 0.5000

La varianza di X₁ è 0,21 e quella di T/n è 0,0105: esattamente il fattore n = 20, come previsto dalle formule p(1 − p) e p(1 − p)/n. Entrambe hanno media 0,3 entro l'errore stampato. Le tre righe finali confermano la identità condizionale: 0,1490 contro 0,15, 0,2996 contro 0,30, 0,5016 contro 0,50, con errori standard di 0,001-0,003. Dei venti dati, la sola somma basta per ottenere lo stimatore migliore; la posizione dei successi non aggiunge nulla.

Esercizio svolto: la probabilità di zero eventi in un processo di Poisson

Il numero di guasti al giorno di una macchina è, per ipotesi, Poisson di parametro λ, e vuoi stimare la probabilità di una giornata senza guasti, che vale e^−λ. Hai osservato n giorni. Un primo stimatore non distorto è W = 1 se il primo giorno non ci sono guasti e 0 altrimenti: ha valore atteso e^−λ, ma usa un solo giorno.

La somma T dei guasti nei n giorni è sufficiente e segue una Poisson di parametro nλ. Dato T = t, i t guasti si distribuiscono sui n giorni come t lanci indipendenti, ciascuno nel giorno 1 con probabilità 1/n; la probabilità che il giorno 1 resti vuoto è (1 − 1/n)^t. Rao-Blackwell dà quindi lo stimatore (1 − 1/n)^T.

Poiché la famiglia di Poisson è una famiglia esponenziale, T è anche completa (nessuna funzione non banale di T ha valore atteso nullo per ogni λ), e il teorema di Lehmann-Scheffé garantisce che questo stimatore è l'unico non distorto di varianza minima. La sua varianza si ottiene dalla funzione generatrice di una Poisson: con a = 1 − 1/n, E[a^2T] = exp(nλ(a² − 1)), da cui Var = exp(nλ(a² − 1)) − e^−2λ. Il programma confronta questo stimatore, la sua varianza esatta e lo stimatore «naturale» e^−X̄ (si sostituisce λ con la media campionaria), per λ = 2 e n = 10.

# Esercizio svolto: UMVUE di exp(-lambda) per un campione Poisson(lambda), con Rao-Blackwell
import numpy as np

lam, n, R = 2.0, 10, 200_000
rng = np.random.default_rng(11)
X = rng.poisson(lam, size=(R, n))
T = X.sum(axis=1)                       # sufficiente e completa: T ~ Poisson(n lam)

W = (X[:, 0] == 0).astype(float)        # stimatore non distorto di P(X=0) = exp(-lam)
a = 1 - 1 / n
RB = a ** T                             # E[W | T = t] = (1 - 1/n)^t
plug = np.exp(-X.mean(axis=1))          # stimatore plug-in exp(-media), per confronto

target = np.exp(-lam)
var_w = target * (1 - target)
var_rb = np.exp(n * lam * (a**2 - 1)) - target**2     # dalla funzione generatrice di Poisson(n lam)
print(f"valore da stimare exp(-lambda) = {target:.5f}")
print(f"W = 1{{X_1=0}}   : media {W.mean():.5f}, varianza {W.var():.5f} (esatta {var_w:.5f})")
print(f"(1-1/n)^T      : media {RB.mean():.5f}, varianza {RB.var():.6f} (esatta {var_rb:.6f})")
print(f"exp(-media)    : media {plug.mean():.5f}, varianza {plug.var():.6f}, MSE {((plug - target) ** 2).mean():.6f}")

# Output:
# valore da stimare exp(-lambda) = 0.13534
# W = 1{X_1=0}   : media 0.13445, varianza 0.11638 (esatta 0.11702)
# (1-1/n)^T      : media 0.13521, varianza 0.004050 (esatta 0.004055)
# exp(-media)    : media 0.14895, varianza 0.004407, MSE 0.004593

Il bersaglio è e⁻² = 0,13534. Lo stimatore grezzo W ha varianza 0,117 (esatta 0,11702), enorme rispetto al valore da stimare; lo stimatore di Rao-Blackwell ha varianza 0,004055 (simulata 0,004050), circa 29 volte più piccola, e media 0,13521, coerente col non distorto.

Lo stimatore plug-in e^−X̄, il primo a cui verrebbe da pensare, è distorto: la sua media simulata è 0,14895, una sovrastima di circa 0,0136 dovuta al fatto che e^−x è una funzione convessa (disuguaglianza di Jensen), e ha MSE 0,004593, peggiore di quello dello stimatore non distorto (0,004055). Qui, a differenza del caso della varianza, l'unico stimatore non distorto di varianza minima vince su entrambi i fronti.

Momenti contro massima verosimiglianza

Il metodo dei momenti costruisce uno stimatore uguagliando i momenti teorici, che dipendono da θ, ai momenti calcolati sul campione, e risolvendo per θ; la giustificazione di questa procedura è nel capitolo del corso MIT sul metodo dei momenti. Per una Gamma di forma α e scala s, media e varianza sono αs e αs²; risolvendo, α = media²/varianza, e basta sostituire i valori del campione.

La massima verosimiglianza sceglie invece il parametro che rende più probabili i dati osservati, cioè massimizza la somma dei logaritmi delle densità. Per la Gamma non c'è una formula chiusa e serve un ottimizzatore numerico: scipy.stats.gamma.fit(x, floc=0) lo fa, e floc=0 fissa a zero il parametro di posizione, come previsto dalla documentazione di fit; senza quell'argomento la stima avrebbe un parametro in più.

Il programma genera 2000 campioni di n = 100 da una Gamma con α = 3 e scala 2, stima α in entrambi i modi, e confronta le varianze con il limite asintotico di Cramér-Rao, che si ottiene invertendo la matrice di informazione di Fisher dei due parametri, la cui voce (α, α) richiede la derivata seconda del logaritmo della funzione gamma, la trigamma scipy.special.polygamma(1, α). Il limite è un riferimento asintotico, non un limite stretto a n = 100.

# Gamma(alfa, scala): metodo dei momenti contro massima verosimiglianza per la forma alfa
import numpy as np
from scipy import stats
from scipy.special import polygamma

alfa, scala, n, R = 3.0, 2.0, 100, 2000
rng = np.random.default_rng(50)
mm = np.empty(R)
ml = np.empty(R)
for r in range(R):
    x = rng.gamma(alfa, scala, size=n)
    mm[r] = x.mean() ** 2 / x.var(ddof=1)            # momenti: media = alfa*scala, varianza = alfa*scala^2
    ml[r] = stats.gamma.fit(x, floc=0)[0]            # massima verosimiglianza, posizione fissata a 0

# Limite asintotico: (I^-1)_{alfa,alfa}/n con I = [[psi'(alfa), 1/scala], [1/scala, alfa/scala^2]]
I = np.array([[polygamma(1, alfa), 1 / scala], [1 / scala, alfa / scala**2]])
cr = np.linalg.inv(I)[0, 0] / n
print(f"vero alfa = {alfa}, n = {n}, ripetizioni = {R}")
print(f"limite asintotico di Cramer-Rao per alfa: varianza {cr:.4f}")
print(f"momenti : media {mm.mean():.4f}, varianza {mm.var():.4f}, rapporto col limite {mm.var() / cr:.3f}")
print(f"MLE     : media {ml.mean():.4f}, varianza {ml.var():.4f}, rapporto col limite {ml.var() / cr:.3f}")

# Output:
# vero alfa = 3.0, n = 100, ripetizioni = 2000
# limite asintotico di Cramer-Rao per alfa: varianza 0.1623
# momenti : media 3.0932, varianza 0.2603, rapporto col limite 1.604
# MLE     : media 3.0841, varianza 0.1905, rapporto col limite 1.173

Il limite asintotico per la varianza di α è 0,1623. Il metodo dei momenti ha varianza 0,2603, circa 1,6 volte il limite; la massima verosimiglianza 0,1905, circa 1,17 volte. Con 2000 ripetizioni una varianza è nota con un errore relativo di circa il 3%, quindi la differenza fra i due è ben fuori dal rumore.

Entrambi hanno una piccola sovrastima media (3,09 e 3,08 su 3). Il prezzo del guadagno è un ottimizzatore al posto di una formula, con le sue possibili non convergenze. In pratica i due si usano insieme: i momenti danno un punto di partenza e un controllo di ragionevolezza, la verosimiglianza la stima finale.

Tre errori che si commettono davvero

• «Non distorto vuol dire migliore»: è falso già per la varianza di una Normale, dove il divisore 21 ha l'MSE più basso.
• «Nessuno stimatore può avere varianza sotto il limite di Cramér-Rao»: il limite vale per gli stimatori non distorti di un modello regolare, e nell'esempio esponenziale uno stimatore distorto lo supera.
• «Un intervallo al 95% contiene il parametro con probabilità 95%»: il 95% descrive la procedura su molti campioni, non il singolo intervallo già calcolato, e a n finito può essere molto meno di 95%.

I tre errori hanno la stessa origine: usare una proprietà asintotica o media come se fosse una garanzia sul caso singolo. Ciascuno si smonta con un programma di poche righe, e la simulazione, quando c'è una formula esatta con cui confrontarla, è un modo economico per accorgersi dell'errore prima di scriverlo in una relazione.

Adesso sai distinguere stimatore e stima, calcolare bias, varianza e MSE di uno stimatore (e scegliere fra due guardando l'MSE), calcolare l'informazione di Fisher e il limite di Cramér-Rao di un modello semplice, applicare Rao-Blackwell a uno stimatore grezzo, e controllare con una simulazione se un intervallo di confidenza copre quanto dichiara. Nel prossimo articolo userai questi strumenti per decidere con un test di ipotesi, e nell'articolo sull'ANOVA li estenderai a più gruppi.

Nella stessa serie trovi anche Variabili aleatorie discrete, Distribuzioni multivariate e covarianza, Regressione lineare semplice e multipla, Stima bayesiana: prior, posterior, MAP, Bootstrap e resampling e Analisi dati esplorativa (EDA) con Python.