05 - Teorema centrale del limite e legge dei grandi numeri

Un istituto di ricerca vuole stimare la quota di elettori favorevoli a una proposta, con un margine di errore di tre punti percentuali e una confidenza del 95%. Quante persone deve intervistare? La risposta è 1.068, e arrivarci richiede due risultati che sono il cuore della teoria: la legge dei grandi numeri, che garantisce che la frequenza osservata si avvicini alla probabilità vera, e il teorema centrale del limite, che dice quanto si avvicina e con quale forma di errore.

Qui li enunci, li metti alla prova su esempi in cui l'errore si calcola esattamente (monete, dadi, somme di variabili esponenziali), vedi dove falliscono (la distribuzione di Cauchy) e chiudi con il sondaggio, controllando con un calcolo esatto se quei 1.068 intervistati bastano davvero. Ti servono le variabili continue, in particolare la normale, e il concetto di indipendenza dell'articolo 01.

Che cosa vuol dire «la media converge»

Siano X₁, X₂, … variabili indipendenti e con la stessa distribuzione (si dice i.i.d.), con media μ e varianza σ², e sia X̄_n la media delle prime n. Questa media è essa stessa una variabile aleatoria, con media μ e varianza σ²/n (la varianza di una somma di indipendenti è la somma delle varianze, poi si divide per n²). «Converge a μ» può significare cose diverse.

La convergenza in probabilità richiede che, per ogni margine ε > 0, la probabilità che X̄_n disti da μ più di ε tenda a zero. La convergenza quasi certa è più forte: riguarda la singola successione di osservazioni e dice che con probabilità 1 le medie correnti finiscono per restare vicine a μ. La convergenza in distribuzione, infine, non dice nulla sui valori ma solo sulla loro legge: la funzione di ripartizione di X_n tende a quella del limite nei punti in cui questa è continua.

La legge debole dei grandi numeri afferma la convergenza in probabilità di X̄_n a μ, e con varianza finita si dimostra in una riga con la disuguaglianza di Chebyshev, P(|Y − E[Y]| ≥ ε) ≤ Var(Y)/ε² (teorema 8.1 del capitolo 8 di Grinstead e Snell): applicata a X̄_n dà P(|X̄_n − μ| ≥ ε) ≤ σ²/(nε²), che tende a zero.

La legge forte, dovuta a Kolmogorov, afferma la convergenza quasi certa ed è vera con la sola ipotesi che la media esista; ne esiste una dimostrazione elementare di Etemadi (1981), che qui non riporto. La differenza fra le due leggi è concreta: la debole garantisce che a ogni n grande sei probabilmente vicino a μ, la forte che la tua singola serie di dati non si allontanerà più, da un certo punto in poi.

Il primo programma, che usa il conteggio dei successi della distribuzione binomiale, confronta tre numeri per la moneta equa (μ = 1/2, σ² = 1/4): la probabilità esatta che la frequenza di testa disti dal valore vero di almeno 0,05, calcolata sommando la distribuzione binomiale, il maggiorante di Chebyshev e l'approssimazione normale che il teorema centrale del limite fornisce. Poi segue una singola traiettoria, la media corrente di un milione di lanci di un dado (μ = 3,5, σ² = 35/12), confrontata con l'errore tipico σ/√n.

import numpy as np
from scipy.stats import binom, norm

# Media di n lanci di moneta equa: quanto e' probabile che la frequenza di testa si discosti da 1/2 di almeno eps?
p, eps = 0.5, 0.05
print(f"{'n':>5} {'esatta':>9} {'Chebyshev':>10} {'CLT':>9}")
for n in (100, 400, 1600):
    k = np.arange(n + 1)
    esatta = binom.pmf(k, n, p)[np.abs(k / n - p) >= eps - 1e-12].sum()   # somma esatta sulla binomiale
    cheb = p * (1 - p) / (n * eps**2)                                    # sigma^2 / (n eps^2)
    clt = 2 * (1 - norm.cdf(eps * n**0.5 / (p * (1 - p)) ** 0.5))        # approssimazione normale
    print(f"{n:5d} {esatta:9.5f} {min(cheb, 1):10.5f} {clt:9.5f}")

# Una traiettoria: media corrente dei lanci di un dado, con errore tipico sigma/sqrt(n), sigma^2 = 35/12
rng = np.random.default_rng(101)
lanci = rng.integers(1, 7, size=10**6)
media = np.cumsum(lanci) / np.arange(1, 10**6 + 1)
sigma = (35 / 12) ** 0.5
print(f"{'n':>8} {'media corrente':>15} {'|media-3,5|':>12} {'sigma/sqrt(n)':>14}")
for n in (10, 100, 10**3, 10**4, 10**5, 10**6):
    print(f"{n:8d} {media[n - 1]:15.4f} {abs(media[n - 1] - 3.5):12.4f} {sigma / n**0.5:14.4f}")

# Output:
#     n    esatta  Chebyshev       CLT
#   100   0.36820    1.00000   0.31731
#   400   0.05104    0.25000   0.04550
#  1600   0.00007    0.06250   0.00006
#        n  media corrente  |media-3,5|  sigma/sqrt(n)
#       10          3.5000       0.0000         0.5401
#      100          3.4000       0.1000         0.1708
#     1000          3.4360       0.0640         0.0540
#    10000          3.4993       0.0007         0.0171
#   100000          3.5047       0.0047         0.0054
#  1000000          3.5023       0.0023         0.0017

Con 100 lanci la probabilità esatta di uno scarto di almeno cinque punti percentuali è 0,368, e il maggiorante di Chebyshev vale 1: è vero ma non dice nulla. Con 400 lanci l'esatta è 0,0510 contro un Chebyshev di 0,25, con 1600 è 0,00007 contro 0,0625. Il maggiorante è una garanzia valida per qualsiasi distribuzione con varianza finita, e proprio per questo è lasco (qui, fino a quasi mille volte); l'approssimazione normale (0,317, 0,0455, 0,00006) segue l'esatta molto da vicino, e il suo errore è l'argomento della prossima sezione.

La traiettoria del dado mostra la legge forte all'opera: lo scarto dalla media 3,5 oscilla (0,100 a n = 100, 0,064 a 1.000, 0,0007 a 10.000) ma resta dell'ordine di σ/√n, che dimezza ogni volta che n quadruplica, cioè scende di un fattore 10 solo moltiplicando le osservazioni per cento.

Quando la legge dei grandi numeri non vale: la Cauchy

Le due leggi richiedono che la media esista. La distribuzione di Cauchy standard, con densità 1/(π(1 + x²)), ha code così pesanti che l'integrale di |x| diverge: non ha valore atteso. In questo caso la media campionaria di n osservazioni è ancora una Cauchy standard, qualunque sia n: mediare non serve a niente.

I quartili della Cauchy standard sono −1 e +1, per cui il programma stima i quartili e la mediana delle medie di n osservazioni, su 2.000 ripetizioni, per n uguale a 1, 100 e 10.000, e li confronta con quelli di una distribuzione che invece ha media, l'esponenziale di media 1.

import numpy as np

# Cauchy standard: nessuna media. La media di n osservazioni e' ancora Cauchy standard, per ogni n.
# I quartili della Cauchy standard sono -1 e +1 (esatti).
rng = np.random.default_rng(103)
R = 2000                                          # replicati per ogni n
print(f"{'n':>7} {'q25 delle medie':>16} {'mediana':>9} {'q75':>8}   (attesi -1, 0, +1; errore tipico circa 0,06)")
for n in (1, 100, 10_000):
    medie = rng.standard_cauchy((R, n)).mean(axis=1)
    q25, q50, q75 = np.quantile(medie, [0.25, 0.5, 0.75])
    print(f"{n:7d} {q25:16.3f} {q50:9.3f} {q75:8.3f}")

# Confronto con una distribuzione che ha media: Exp(1), media 1. I quartili delle medie si stringono attorno a 1.
for n in (1, 100, 10_000):
    medie = rng.exponential(1.0, (R, n)).mean(axis=1)
    q25, q50, q75 = np.quantile(medie, [0.25, 0.5, 0.75])
    print(f"Exp  n={n:6d}: q25 {q25:.4f}  mediana {q50:.4f}  q75 {q75:.4f}")

# Output:
#       n  q25 delle medie   mediana      q75   (attesi -1, 0, +1; errore tipico circa 0,06)
#       1           -0.996    -0.046    0.907
#     100           -1.058    -0.034    1.040
#   10000           -1.070    -0.036    1.011
# Exp  n=     1: q25 0.2980  mediana 0.7037  q75 1.3764
# Exp  n=   100: q25 0.9278  mediana 0.9951  q75 1.0655
# Exp  n= 10000: q25 0.9942  mediana 1.0006  q75 1.0071

Per la Cauchy i quartili delle medie restano a circa −1 e +1 e la mediana a circa zero con n = 1, 100 e 10.000 (per esempio −1,070 e 1,011 a n = 10.000), senza alcun restringimento: l'errore di stima a 10.000 osservazioni è grande quanto quello di una sola.

Per l'esponenziale, invece, l'intervallo fra i quartili si riduce da (0,298; 1,376) a (0,928; 1,066) e infine a (0,994; 1,007), cioè la larghezza scende come 1/√n. Il messaggio è operativo: prima di mediare grandezze con code estreme (pagamenti di assicurazioni, tempi di risposta, rendimenti nei crolli) controlla che la media esista; se non esiste, o ha varianza infinita, il teorema che sta per seguire non si applica.

Il teorema centrale del limite e il suo errore

La legge dei grandi numeri dice che X̄_n si avvicina a μ; il teorema centrale del limite (di Lindeberg-Lévy) dice come: se le X_i sono i.i.d. con varianza finita e positiva, la somma standardizzata (S_n − nμ)/(σ√n), dove S_n = X₁ + … + X_n, ha una funzione di ripartizione che tende a Φ, quella della normale standard, in ogni punto.

La forma di partenza della distribuzione non conta, purché la varianza sia finita. Il teorema è un limite e non dice quanto è grande l'errore a n fissato. Lo dice la disuguaglianza di Berry-Esseen (Berry 1941): la massima distanza fra la funzione di ripartizione della somma standardizzata e Φ è al più Cρ/(σ³√n), dove ρ = E[|X − μ|³] e C è una costante universale, per cui Shevtsova (2011) ha dimostrato il valore 0,4748 nel caso i.i.d.

Per misurare l'errore senza simulare serve un caso in cui la legge della somma sia nota: la somma di n esponenziali di media 1 ha esattamente distribuzione Gamma di forma n, quindi la distanza di Kolmogorov si calcola con la funzione di ripartizione esatta, su una griglia fitta. Il programma la calcola per n da 1 a 1.000, calcola ρ per integrazione numerica (la formula chiusa è 12/e − 2), la confronta con il maggiorante di Berry-Esseen, e infine simula 200.000 somme con n = 30 per vedere la stessa distanza da un campione.

import numpy as np
from scipy.stats import gamma, norm, kstest
from scipy.integrate import quad

# X ~ Exp(1): media 1, varianza 1, asimmetria 2. S_n = X_1+...+X_n ha legge Gamma(n, 1) ESATTA.
rho, _ = quad(lambda x: abs(x - 1) ** 3 * np.exp(-x), 0, np.inf, limit=200)   # E|X - mu|^3
print(f"rho = E|X-1|^3 = {rho:.6f}  (formula 12/e - 2 = {12 / np.e - 2:.6f})")
C = 0.4748                                       # costante di Berry-Esseen per somme i.i.d. (Shevtsova 2011)
z = np.linspace(-6, 6, 120001)
print(f"{'n':>5} {'distanza esatta':>16} {'maggiorante B-E':>16} {'dist*sqrt(n)':>13}")
for n in (1, 5, 30, 200, 1000):
    # F_n(z) = P((S_n - n)/sqrt(n) <= z) = Gamma(n,1).cdf(n + z sqrt(n))
    Fn = gamma.cdf(n + z * n**0.5, a=n)
    dist = np.max(np.abs(Fn - norm.cdf(z)))
    print(f"{n:5d} {dist:16.5f} {C * rho / n**0.5:16.5f} {dist * n**0.5:13.4f}")

# Conferma per simulazione a n = 30: 200000 somme, distanza di Kolmogorov-Smirnov dalla normale standard
rng = np.random.default_rng(107)
n, M = 30, 200_000
exact30 = np.max(np.abs(gamma.cdf(n + z * n**0.5, a=n) - norm.cdf(z)))
S = rng.exponential(1.0, size=(M, n)).sum(axis=1)
Zs = (S - n) / n**0.5
D = kstest(Zs, norm.cdf).statistic
print(f"n=30, {M} somme simulate: D = {D:.5f}, distanza esatta {exact30:.5f}, differenza {D - exact30:+.5f} (rumore tipico 1/sqrt(M) = {1 / M**0.5:.5f})")
terzo = Zs**3
print(f"terzo momento delle somme standardizzate: {terzo.mean():.3f} +/- {terzo.std() / M**0.5:.3f} (esatto 2/sqrt(30) = {2 / 30**0.5:.3f})")

# Output:
# rho = E|X-1|^3 = 2.414553  (formula 12/e - 2 = 2.414553)
#     n  distanza esatta  maggiorante B-E  dist*sqrt(n)
#     1          0.15866          1.14643        0.1587
#     5          0.05963          0.51270        0.1333
#    30          0.02429          0.20931        0.1330
#   200          0.00940          0.08106        0.1330
#  1000          0.00421          0.03625        0.1330
# n=30, 200000 somme simulate: D = 0.02781, distanza esatta 0.02429, differenza +0.00352 (rumore tipico 1/sqrt(M) = 0.00224)
# terzo momento delle somme standardizzate: 0.341 +/- 0.010 (esatto 2/sqrt(30) = 0.365)

La distanza esatta scende come 1/√n, da 0,159 a n = 1 a 0,0243 a n = 30 e a 0,0042 a n = 1.000, e il maggiorante di Berry-Esseen la copre sempre ma con un margine di quasi nove volte (per n = 1 vale 1,146, quindi non dice nulla). Il prodotto della distanza per √n si stabilizza a 0,1330: è il primo termine della correzione di Edgeworth, γ/(6√(2π)) con γ = 2 l'asimmetria dell'esponenziale, un risultato che qui non derivo, e che spiega perché la lentezza dell'approssimazione viene dall'asimmetria.

La simulazione a n = 30 dà D = 0,0278 contro una distanza esatta di 0,0243: la differenza (0,0035) è circa una volta e mezza il rumore tipico, perché la distanza campionaria è il massimo di scarti rumorosi e in media non può stare sotto quella vera.

Il terzo momento delle somme simulate è 0,341 ± 0,010 contro l'esatto 2/√30 = 0,365: è due volte e mezza l'errore standard, una fluttuazione di questo seme (altri cinque semi, provati a parte, danno valori fra 0,353 e 0,391 attorno allo stesso valore esatto). La regola «con n ≥ 30 vale la normale» è dunque un'abitudine, non un teorema: per una distribuzione asimmetrica come questa il 2,4% di errore nella funzione di ripartizione a n = 30 può essere troppo.

Esercizio svolto: quanti intervistati servono

Torniamo al sondaggio. Ogni intervistato risponde sì con probabilità p sconosciuta; la frequenza p̂ dei sì su n persone ha media p e varianza p(1 − p)/n, e per il teorema centrale del limite è approssimativamente normale. L'intervallo p̂ ± m contiene p con probabilità circa 95% quando m = 1,96√(p(1 − p)/n), dove 1,96 è il quantile 97,5% della normale standard, cioè n = 1,96²p(1 − p)/m².

Poiché p non si conosce, si usa il caso peggiore p = 1/2, che massimizza p(1 − p), con m = 0,03: n = 1.067,07, da arrotondare per eccesso a 1.068.

Il programma verifica questa risposta con il calcolo esatto della probabilità di copertura, sommando la distribuzione binomiale, per quattro valori veri di p e per due tipi di intervallo: quello a margine fisso di ±0,03 e quello di Wald, la cui ampiezza è stimata dalla frequenza osservata, p̂ ± 1,96√(p̂(1 − p̂)/n). Infine simula 200.000 sondaggi con p = 1/2 per confermare la copertura.

import numpy as np
from math import ceil
from scipy.stats import binom

# Margine di errore +/-3 punti, confidenza 95%: n >= z^2 p(1-p) / m^2, nel caso peggiore p = 1/2.
z, m = 1.959964, 0.03                             # z = quantile 97,5% della normale standard
n = ceil(z**2 * 0.25 / m**2)
print(f"n minimo (caso peggiore p=1/2): {z**2 * 0.25 / m**2:.2f} -> {n}")

k = np.arange(n + 1)
phat = k / n
print(f"{'p vero':>7} {'copertura +/-3%':>16} {'copertura Wald':>15}")
for p in (0.5, 0.3, 0.1, 0.02):
    pesi = binom.pmf(k, n, p)
    fisso = pesi[np.abs(phat - p) <= m].sum()                              # intervallo phat +/- 0,03
    semi = z * np.sqrt(phat * (1 - phat) / n)                              # Wald: ampiezza stimata da phat
    wald = pesi[(phat - semi <= p) & (p <= phat + semi)].sum()
    print(f"{p:7.2f} {fisso:16.4f} {wald:15.4f}")

rng = np.random.default_rng(109)
R = 200_000
x = rng.binomial(n, 0.5, size=R) / n
cop = np.mean(np.abs(x - 0.5) <= m)
esatta = binom.pmf(k, n, 0.5)[np.abs(phat - 0.5) <= m].sum()
tol = 3 * (esatta * (1 - esatta) / R) ** 0.5
print(f"p=0.5, {R} sondaggi simulati: copertura {cop:.4f} (esatta {esatta:.4f} +/- {tol:.4f})")
assert abs(cop - esatta) < tol

# Output:
# n minimo (caso peggiore p=1/2): 1067.07 -> 1068
#  p vero  copertura +/-3%  copertura Wald
#    0.50           0.9533          0.9462
#    0.30           0.9674          0.9511
#    0.10           0.9988          0.9445
#    0.02           1.0000          0.9290
# p=0.5, 200000 sondaggi simulati: copertura 0.9532 (esatta 0.9533 +/- 0.0014)

Con 1.068 intervistati e p = 1/2 la copertura esatta del margine fisso è 0,9533, in linea con il 95% richiesto (l'eccedenza è dovuta all'arrotondamento di n), e la simulazione dà 0,9532, dentro la tolleranza di ±0,0014. Il margine fisso è prudente: per p = 0,3 la copertura sale a 0,9674 e per p = 0,1 a 0,9988, perché a quel p la varianza è minore.

L'intervallo di Wald, invece, mostra dove l'approssimazione normale si indebolisce: 0,9462 a p = 0,5, 0,9511 a 0,3, 0,9445 a 0,1 e solo 0,9290 a p = 0,02, perché quando il successo è raro la distribuzione di p̂ è asimmetrica e la stima della varianza dalla frequenza osservata è inaffidabile. Resta un limite del modello, non del calcolo: tutto ciò vale per un campione casuale semplice con risposte sincere; un sondaggio reale ha altre fonti di errore, come la mancata risposta e la selezione del campione, che questo calcolo non misura.

Che cosa resta fuori

Questo articolo tratta il caso di variabili i.i.d. con varianza finita. Quando le variabili non sono identicamente distribuite serve la condizione di Lindeberg (1922), che chiede che nessun termine domini la somma; il teorema di Lindeberg-Feller ne fa una condizione necessaria e sufficiente.

Non tratto neppure le correzioni di Edgeworth oltre il primo termine, le versioni per vettori (che portano a una normale multivariata, come nell'articolo 04), né le leggi dei grandi numeri per variabili dipendenti. La dimostrazione del teorema, per funzioni caratteristiche, richiede strumenti di analisi che qui non servono. Il passo successivo è l'inferenza statistica, dove il teorema centrale del limite giustifica gli intervalli di confidenza.

Domande frequenti

La legge dei grandi numeri dice che dopo molti testa deve uscire croce?

No. I lanci sono indipendenti e la moneta non ha memoria: la frequenza si avvicina a 1/2 perché i lanci futuri diluiscono lo scarto accumulato, non perché lo compensino. Lo scarto assoluto può continuare a crescere come √n pur con frequenza che converge.

A partire da quante osservazioni vale la normale?

Dipende dalla forma della distribuzione. Per l'esponenziale, la distanza dalla normale a 30 osservazioni è 0,024; per distribuzioni simmetriche è molto minore. Il maggiorante di Berry-Esseen dà una garanzia, ma prudente.

Perché la media di dati Cauchy non converge?

Perché la distribuzione non ha valore atteso: i valori estremi sono abbastanza frequenti da dominare la somma. La media di n osservazioni ha esattamente la stessa distribuzione di una sola.

Che cosa sai fare adesso

Sai distinguere le convergenze in probabilità, quasi certa e in distribuzione, dimostrare la legge debole con Chebyshev e riconoscere un caso, la Cauchy, in cui nessuna delle due leggi vale. Sai enunciare il teorema centrale del limite, misurarne l'errore con un calcolo esatto e confrontarlo con il maggiorante di Berry-Esseen, e dimensionare un campione verificando poi la copertura invece di fidarti della sola formula. Il passo successivo è sfruttare questi risultati per stimare un parametro da dati: stimatori, loro proprietà e intervalli di confidenza.

Nella stessa serie trovi anche Test di ipotesi: t-test, z-test, p-value e Analisi dati esplorativa (EDA) con Python. L'elenco completo è nella pagina della serie, Probabilità e Statistica.