29 - Campionamento e metodi Monte Carlo: stimare quello che non si riesce a calcolare

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

Obiettivo. Trasformare la simulazione da strumento di verifica a strumento di calcolo: stimare un'aspettativa che non si sa integrare, e sapere quanto vale la stima. Alla fine sai dimostrare che la media campionaria non è distorta e che il suo errore scende come 1/√n, sai generare campioni da una distribuzione qualsiasi con la trasformata inversa o per rifiuto, sai spostare i campioni dove servono con il campionamento per importanza, e sai spiegare perché il trucco della riparametrizzazione rende derivabile un'aspettativa che altrimenti non lo è.

Falso mito da sfatare. «Monte Carlo è un metodo approssimato: se vuoi il numero giusto, usa una formula.» No: Monte Carlo è un metodo esatto in media con un errore quantificabile, e in dimensione alta è l'unico che funziona. La formula chiusa, quando esiste, vince sempre; il punto è che quasi non esiste mai.

Introduzione

Nella lezione sulle variabili aleatorie la simulazione è servita come verifica: si generavano numeri per vedere «succedere» la legge dei grandi numeri e il teorema del limite centrale, cioè per convincersi di due teoremi già dimostrati. Qui la simulazione cambia ruolo e diventa uno stimatore: lo strumento con cui si calcola un integrale che nessuno sa risolvere a mano.

Non è un argomento di contorno. Tre lezioni successive useranno questi strumenti dandoli per noti: il variational autoencoder si addestra solo grazie alla riparametrizzazione, le instabilità dell'ELBO sono instabilità di uno stimatore Monte Carlo, e tutto il modulo sull'apprendimento per rinforzo, a partire dall'equazione di Bellman, stima ritorni attesi per simulazione. Chi arriva lì senza questa lezione copia formule.

Risultato finale di questa lezione. Scrivere uno stimatore per simulazione, dichiararne l'errore standard, scegliere il numero di campioni a partire dall'errore voluto, e riconoscere i modi in cui uno stimatore corretto produce comunque numeri inservibili.

Roadmap della lezione.

1. Lo stimatore Monte Carlo, la dimostrazione che non è distorto e il suo errore σ/√n.
2. Il confronto con una quadratura su griglia, in dimensione 1, 3 e 10.
3. Come si generano numeri casuali, e perché il seme fa parte dell'esperimento.
4. La trasformata inversa: derivazione, esempio continuo e discreto, limiti.
5. Il campionamento per rifiuto: dimostrazione di correttezza e crollo in alta dimensione.
6. Il campionamento per importanza, la sua varianza, e il caso in cui esplode.
7. Variabili antitetiche e variate di controllo, con la riduzione misurata.
8. Il trucco della riparametrizzazione, confrontato numericamente con REINFORCE.
9. Un esperimento completo su una coda a cinque stati, con intervalli di confidenza.

La scaletta va dal metodo alle sue fragilità: prima lo stimatore e il suo errore, poi i tre modi di produrre campioni (trasformata inversa, rifiuto, importanza), poi come ridurre la varianza, infine la riparametrizzazione e un esperimento completo. Ogni tecnica ha una sezione in cui si rompe: leggi quelle prima di fidarti di una stima.

L'idea di Monte Carlo, e il prezzo della sua precisione

L'analogia: contare i pesci nel lago

Per sapere quanti pesci ci sono in un lago non si svuota il lago. Si pesca un centinaio di pesci, li si marca e li si rimette in acqua; qualche giorno dopo si ripesca e si guarda che frazione dei pescati porta il marchio. Se su cento ripescati dodici sono marcati, i cento marcati sono circa il 12% della popolazione: ottocento pesci.

Di quella stima vanno notate due cose, che sono esattamente le due cose di Monte Carlo. Il metodo non guarda tutto il lago, guarda un campione, e funziona proprio perché il campione è preso a caso. E il risultato non è un numero, è un numero con un'incertezza: ripescando si troverebbero undici marcati, o quattordici. Un biologo che riportasse «ottocento pesci» senza dire «più o meno quanti» non avrebbe fatto una misura. In termini di ingegneria: quando una quantità si scrive come valore atteso, si può stimarla estraendo campioni e facendone la media.

Lo stimatore, e la dimostrazione che non è distorto

Sia X una variabile aleatoria con densità p, e f una funzione di cui interessa il valore medio. La quantità che si vuole è

I = E_p[f(X)] = ∫ f(x) p(x) dx.

Lo stimatore Monte Carlo prende n campioni indipendenti x₁, …, xₙ estratti da p e ne fa la media:

Î_n = (1/n) · Σᵢ f(xᵢ).

Che questo stimatore sia non distorto (unbiased) si dimostra in due righe, e la dimostrazione usa solo la linearità del valore atteso: E[Î_n] = E[(1/n) Σᵢ f(xᵢ)] = (1/n) Σᵢ E[f(xᵢ)] = (1/n) · n · I = I. In parole: la media delle medie è la media. Non c'è nessuna approssimazione in questo passaggio, e non c'è nessuna ipotesi su f oltre al fatto che I esista. ∎

«Non distorto» significa che in media lo stimatore dà la risposta giusta, non che una esecuzione la dia. Sono due affermazioni diverse e confonderle è l'errore numero uno di questa materia.

L'errore: σ/√n, e perché un decimale costa cento volte

Quanto si discosta Î_n da I? Gli f(xᵢ) sono indipendenti e hanno tutti varianza σ² = Var_p[f(X)]; la varianza di una somma è la somma delle varianze, e dividendo per n il divisore si eleva al quadrato:

Var[Î_n] = (1/n²) · Σᵢ Var[f(xᵢ)] = (1/n²) · n · σ² = σ²/n,

quindi la deviazione standard dello stimatore, che si chiama errore standard, vale σ/√n. Leggila come una frase: l'incertezza scende come la radice quadrata del numero di campioni. La conseguenza è brutale e va detta esplicitamente: per guadagnare una cifra decimale servono cento volte i campioni. Da n a 4n l'errore si dimezza, da n a 100n si divide per dieci. Chi promette «basta simulare un po' di più» non ha fatto questo conto.

C'è però una proprietà che rende Monte Carlo insostituibile, e sta in ciò che σ/√n non contiene: la dimensione dello spazio. Integrare in una dimensione o in mille non cambia il tasso di convergenza, cambia solo σ. Nessun metodo deterministico ha questa proprietà, ed è il motivo per cui statistica bayesiana, fisica computazionale e apprendimento per rinforzo girano tutti su simulazioni.

Verifichiamo entrambe le affermazioni su un caso con risposta nota: X ∼ N(0, 1) e f(x) = x², per cui I = 1 e σ² = Var[X²] = 2, quindi l'errore standard teorico è √2/√n.

import numpy as np

def stima_monte_carlo(n, rng):
    """Stima E[f(X)] con f(x) = x² e X ~ N(0, 1): il valore esatto è 1."""
    x = rng.standard_normal(n)
    f = x ** 2
    return f.mean(), f.std(ddof=1) / np.sqrt(n)      # σ̂/√n, NON σ̂

rng = np.random.default_rng(42)
print("      n      stima   errore stimato   errore vero   σ/√n teorico")
for n in (10, 100, 1000, 10000, 100000):
    stima, errore_stimato = stima_monte_carlo(n, rng)
    print(f"{n:>7} {stima:>10.5f} {errore_stimato:>16.5f}"
          f" {abs(stima - 1.0):>13.5f} {np.sqrt(2.0 / n):>13.5f}")

# l'errore tipico si misura ripetendo l'esperimento, non eseguendolo una volta
rng = np.random.default_rng(7)
print("\n      n   errore quadratico medio   σ/√n teorico   rapporto")
for n in (10, 100, 1000, 10000, 100000):
    errori = np.array([(rng.standard_normal(n) ** 2).mean() - 1.0 for _ in range(1000)])
    rmse = float(np.sqrt((errori ** 2).mean()))
    print(f"{n:>7} {rmse:>25.6f} {np.sqrt(2.0 / n):>14.6f} {rmse / np.sqrt(2.0 / n):>10.4f}")

Output ottenuto (Python 3.14.4, NumPy 2.3.5, semi 42 e 7):

      n      stima   errore stimato   errore vero   σ/√n teorico
     10    0.89688          0.36726       0.10312       0.44721
    100    0.59126          0.07694       0.40874       0.14142
   1000    1.00197          0.04563       0.00197       0.04472
  10000    1.01770          0.01452       0.01770       0.01414
 100000    1.00663          0.00452       0.00663       0.00447

      n   errore quadratico medio   σ/√n teorico   rapporto
     10                  0.431696       0.447214     0.9653
    100                  0.140293       0.141421     0.9920
   1000                  0.044682       0.044721     0.9991
  10000                  0.014107       0.014142     0.9975
 100000                  0.004492       0.004472     1.0043

Cosa fa il codice. La prima tabella esegue una simulazione per ogni n: stima, errore standard dai dati (σ̂/√n, dove il /√n è la parte che quasi tutti dimenticano), errore vero — che conosciamo solo perché la risposta è nota — ed errore teorico. La seconda ripete l'esperimento mille volte per ogni n e riporta la radice dell'errore quadratico medio: è così che si misura l'errore tipico di uno stimatore.

Che cosa ci dice. La seconda tabella conferma la legge su quattro ordini di grandezza: il rapporto fra errore misurato e σ/√n teorico sta fra 0,965 e 1,004. La prima dice qualcosa di più interessante alla riga n = 100, dove l'errore vero (0,409) è più di cinque volte quello stimato (0,077). Non è un bug: è un campione sfortunato in cui la stima è bassa (0,591) e la varianza campionaria è sottostimata insieme a lei, perché x² è asimmetrica e chi perde la coda destra perde insieme media e dispersione. Un errore standard è esso stesso una stima, su n piccolo è rumoroso quanto la quantità che sorveglia, e per di più è correlato con l'errore che dovrebbe misurare.

Monte Carlo contro una griglia: quando conviene davvero

L'alternativa deterministica a Monte Carlo è la quadratura su griglia: si divide il dominio in celle, si valuta la funzione al centro di ciascuna e si somma. In una dimensione la regola del punto medio ha errore O(h²) con h = 1/m il passo, cioè O(m⁻²): raddoppiando i nodi l'errore si divide per quattro. Monte Carlo, a confronto, con il suo O(n⁻¹/²) fa una figura pessima.

Il problema è che in d dimensioni una griglia con m nodi per asse ne richiede m^d. A parità di budget di valutazioni N, la griglia può permettersi solo m ≈ N^{1/d} nodi per asse, quindi il suo errore è O(N^{−2/d}). Confrontalo con O(N^{−1/2}): la griglia vince finché 2/d > 1/2, cioè finché d < 4, e perde dopo. È la maledizione della dimensionalità vista in una disuguaglianza di una riga.

Misuriamolo su ∫ exp(−‖x‖²) dx sul cubo [0,1]^d, che è separabile e ha quindi valore esatto ((√π/2)·erf(1))^d, con un budget fisso di 10 000 valutazioni.

import math
import numpy as np

def integrale_esatto(d):
    """∫ exp(−‖x‖²) su [0,1]^d: l'integrando è separabile, quindi è una potenza."""
    return (math.sqrt(math.pi) / 2.0 * math.erf(1.0)) ** d

def quadratura_punto_medio(d, bilancio):
    """Regola del punto medio: m nodi per asse, quindi m^d valutazioni."""
    m = int(bilancio ** (1.0 / d))
    nodi = (np.arange(m) + 0.5) / m
    una_dimensione = float(np.exp(-nodi ** 2).sum() / m)
    return una_dimensione ** d, m, m ** d

def quadratura_monte_carlo(d, n, rng):
    u = rng.random((n, d))
    f = np.exp(-(u ** 2).sum(axis=1))
    return f.mean(), f.std(ddof=1) / np.sqrt(n)

rng = np.random.default_rng(42)
bilancio = 10_000
print("  d | m per asse  valutazioni  errore griglia | errore MC  errore standard MC")
for d in (1, 3, 10):
    esatto = integrale_esatto(d)
    griglia, m, usate = quadratura_punto_medio(d, bilancio)
    mc, errore_standard = quadratura_monte_carlo(d, bilancio, rng)
    print(f" {d:>2} | {m:>10} {usate:>12} {abs(griglia - esatto):>15.3e} |"
          f" {abs(mc - esatto):>9.3e} {errore_standard:>19.3e}")

Output ottenuto (seme 42):

  d | m per asse  valutazioni  errore griglia | errore MC  errore standard MC
  1 |      10000        10000       3.066e-10 | 2.105e-03           2.005e-03
  3 |         21         9261       1.163e-04 | 3.780e-03           1.998e-03
 10 |          2         1024       5.889e-03 | 3.253e-05           5.457e-04

Cosa fa il codice. quadratura_punto_medio sfrutta la separabilità: l'integrale su d dimensioni è la potenza d-esima di quello su una. int(bilancio ** (1.0 / d)) è il numero di nodi per asse che il budget consente: 10 000 in una dimensione, 21 in tre (21³ = 9261), e soltanto 2 in dieci, perché 3¹⁰ = 59 049 sarebbe già sei volte fuori budget.

Che cosa ci dice, e va letto con onestà. In dimensione 1 la griglia stravince: 3,07·10⁻¹⁰ contro 2,11·10⁻³, sette ordini di grandezza con le stesse valutazioni. Chi integra una funzione liscia di una variabile usando Monte Carlo spreca un fattore dieci milioni. In dimensione 3 la griglia è ancora avanti di trenta volte. In dimensione 10 il rapporto si rovescia: strozzata a due nodi per asse, la griglia sbaglia di 5,9·10⁻³ contro un errore standard Monte Carlo di 5,5·10⁻⁴.

Un avvertimento sulla riga d = 10. L'errore Monte Carlo lì riportato è 3,25·10⁻⁵, cioè più piccolo del suo stesso errore standard 5,46·10⁻⁴: è un tiro fortunato, lo stesso fenomeno della riga n = 100 di sopra con il segno rovesciato. Il numero da confrontare con la griglia è 5,46·10⁻⁴; riportare l'altro significherebbe dichiarare una precisione che il metodo non ha. La regola da portare via: Monte Carlo non è «meglio» della quadratura, è indifferente alla dimensione. Sotto le quattro dimensioni, con funzioni lisce, usa una griglia; sopra, la griglia non esiste come opzione.

Generare numeri casuali: il seme è parte dell'esperimento

Tutto quanto precede presuppone di saper estrarre campioni, e alla base c'è una sola primitiva: un generatore di numeri uniformi in [0, 1). Tutto il resto — normali, esponenziali, distribuzioni discrete, permutazioni — si costruisce da lì, ed è ciò che fanno le due sezioni successive.

Quei numeri non sono casuali: sono prodotti da una funzione deterministica dello stato interno del generatore, il cui valore iniziale è il seme. Si chiamano pseudocasuali perché superano i test statistici di casualità pur essendo riproducibili, e questa non è una debolezza da tollerare: è la proprietà che rende un esperimento Monte Carlo un esperimento invece di un aneddoto. Il seme fa parte del protocollo di misura e va dichiarato come si dichiara la versione di una libreria.

In NumPy si usa np.random.default_rng(seme), che restituisce un generatore con il proprio stato, basato sull'algoritmo PCG64 di O'Neill. La vecchia interfaccia np.random.seed() più np.random.rand() è dichiarata congelata dalla documentazione ufficiale, che la riserva ai casi in cui servono gli stessi valori delle versioni vecchie, e va evitata nel codice nuovo perché agisce su uno stato globale condiviso con ogni libreria del processo: basta che una dipendenza estragga un numero per spostare la tua sequenza, e il risultato «riproducibile» smette di esserlo senza che niente lo segnali.

import numpy as np

rng_a = np.random.default_rng(42)
rng_b = np.random.default_rng(42)
print("stesso seme, due generatori:", rng_a.random(3), "|", rng_b.random(3))

rng = np.random.default_rng(42)
print("primo blocco: ", rng.random(3))
print("secondo blocco:", rng.random(3))       # lo stato è avanzato: numeri diversi
print("algoritmo:", type(rng.bit_generator).__name__)

# due flussi indipendenti e riproducibili, senza toccare nessuno stato globale
radice = np.random.SeedSequence(42)
flussi = [np.random.default_rng(s) for s in radice.spawn(2)]
print("flusso 0:", flussi[0].random(2))
print("flusso 1:", flussi[1].random(2))

Output ottenuto (seme 42):

stesso seme, due generatori: [0.77395605 0.43887844 0.85859792] | [0.77395605 0.43887844 0.85859792]
primo blocco:  [0.77395605 0.43887844 0.85859792]
secondo blocco: [0.69736803 0.09417735 0.97562235]
algoritmo: PCG64
flusso 0: [0.91674416 0.91098667]
flusso 1: [0.46749078 0.0464489 ]

Cosa fa il codice. Due generatori con lo stesso seme producono la stessa sequenza cifra per cifra, ma lo stato avanza: la seconda chiamata allo stesso generatore dà numeri diversi. L'ultimo blocco usa SeedSequence.spawn, il modo corretto di ottenere flussi indipendenti da un unico seme quando si simula in parallelo — mettere default_rng(42) in ogni processo ripeterebbe lo stesso identico lavoro, con riduzione di varianza nulla e nessun errore. La regola da portare via: un generatore si passa come argomento. Una funzione che chiama np.random.rand() dentro di sé non è testabile, e due funzioni che lo fanno non sono componibili.

La trasformata inversa

La derivazione

Data una variabile aleatoria con funzione di ripartizione (CDF) F(x) = P(X ≤ x) continua e strettamente crescente, si vuole generare campioni di X partendo da U ∼ U(0,1). L'affermazione è che basta applicare F⁻¹:

se U ∼ U(0,1) allora X = F⁻¹(U) ha funzione di ripartizione F.

La dimostrazione è tre passaggi e usa solo la definizione. Si calcola la CDF di F⁻¹(U): P(F⁻¹(U) ≤ x). Siccome F è crescente, applicarla a entrambi i membri non cambia il verso della disuguaglianza, quindi l'evento {F⁻¹(U) ≤ x} coincide con {U ≤ F(x)}. Infine, per una uniforme su [0,1] vale P(U ≤ t) = t per ogni t ∈ [0,1], e F(x) sta proprio lì dentro. Mettendo insieme: P(F⁻¹(U) ≤ x) = P(U ≤ F(x)) = F(x), che è la tesi. ∎

L'intuizione geometrica: dove la densità è alta la CDF sale ripida, quindi un intervallo di u di ampiezza fissa viene mappato in un intervallo di x stretto, e lì cadono molti campioni. La trasformata inversa redistribuisce i punti uniformi secondo la pendenza di F.

Due esempi svolti fino in fondo

Esponenziale. Per X ∼ Exp(λ) si ha F(x) = 1 − e^{−λx} per x ≥ 0. Si pone u = 1 − e^{−λx} e si risolve in x: e^{−λx} = 1 − u, quindi −λx = ln(1 − u), quindi x = −ln(1 − u)/λ. Una riga di codice, nessun rifiuto, nessun ciclo.

Discreta. Con esiti a₁, …, a_k e probabilità p₁, …, p_k la CDF è una scala: si calcolano le cumulate c_j = p₁ + … + p_j e si restituisce il primo a_j con u ≤ c_j. L'inversa generalizzata è F⁻¹(u) = min{ a_j : F(a_j) ≥ u }, in NumPy una sola chiamata a searchsorted.

import numpy as np

def esponenziale_per_trasformata_inversa(n, lam, rng):
    """F⁻¹(u) = −ln(1 − u)/λ. `rng.random()` sta in [0, 1): 1 − u non è mai 0."""
    u = rng.random(n)
    return -np.log(1.0 - u) / lam

rng = np.random.default_rng(42)
lam = 0.4
campioni = esponenziale_per_trasformata_inversa(200_000, lam, rng)
print(f"media    {campioni.mean():.5f}   attesa 1/λ    = {1 / lam:.5f}")
print(f"mediana  {np.median(campioni):.5f}   attesa ln2/λ  = {np.log(2) / lam:.5f}")
print(f"dev.std. {campioni.std(ddof=1):.5f}   attesa 1/λ    = {1 / lam:.5f}")

# caso discreto: l'esito di una richiesta di servizio, tre valori possibili
esiti = ("servita", "in coda", "rifiutata")
probabilita = np.array([0.6, 0.3, 0.1])
cumulate = np.cumsum(probabilita)
indici = np.searchsorted(cumulate, rng.random(200_000))
for k, nome in enumerate(esiti):
    print(f"{nome:>10}: frequenza {np.mean(indici == k):.5f}   attesa {probabilita[k]:.5f}")

Output ottenuto (seme 42):

media    2.50114   attesa 1/λ    = 2.50000
mediana  1.73343   attesa ln2/λ  = 1.73287
dev.std. 2.50291   attesa 1/λ    = 2.50000
   servita: frequenza 0.60023   attesa 0.60000
   in coda: frequenza 0.30019   attesa 0.30000
 rifiutata: frequenza 0.09959   attesa 0.10000

Cosa fa il codice. La funzione applica letteralmente la formula derivata sopra. Il dettaglio di 1.0 - u invece di u non è cosmetico: rng.random() dà valori in [0, 1), quindi lo zero è possibile e l'uno no, e 1 − u non si annulla mai; scrivendo −ln(u)/λ si otterrebbe un inf con probabilità piccolissima ma non nulla, cioè un guasto che compare in produzione dopo settimane. Nel caso discreto searchsorted restituisce l'indice della prima cumulata che supera u: è l'inversa generalizzata in una riga, in O(log k).

Che cosa ci dice. Tutte e tre le statistiche dell'esponenziale e le tre frequenze discrete cadono sulla terza cifra. Che media e deviazione standard coincidano entrambe con 1/λ è una proprietà dell'esponenziale, non un errore di stampa.

Il limite del metodo, ed è grosso. Serve F⁻¹ in forma chiusa. Per esponenziale, uniforme, Weibull, Cauchy e discrete finite c'è; per la normale no (la sua CDF è una funzione speciale, invertibile solo numericamente), e per quasi tutte le distribuzioni interessanti in statistica bayesiana nemmeno.

Il campionamento per rifiuto

L'idea geometrica

Campionare da una densità p equivale a estrarre un punto uniforme nell'area sotto il suo grafico e tenerne l'ascissa. Il campionamento per rifiuto (rejection sampling) dice: se non so estrarre punti sotto p, ne estraggo sotto una curva più comoda che le sta sopra, e butto via quelli che finiscono nello spazio fra le due. Servono una proposta q da cui so campionare e una costante M tali che

p(x) ≤ M · q(x) per ogni x.

L'algoritmo è: estrai x ∼ q; estrai u ∼ U(0,1); accetta x se u ≤ p(x) / (M·q(x)), altrimenti ricomincia.

La dimostrazione che i campioni accettati hanno distribuzione p

Calcoliamo la probabilità congiunta di estrarre un valore ≤ x e di accettarlo. Condizionando sul valore estratto: P(X ≤ x, accetto) = ∫_{−∞}^{x} q(t) · [p(t)/(M q(t))] dt, perché la probabilità di accettare dato t è proprio il rapporto. I due q(t) si semplificano e resta (1/M) ∫_{−∞}^{x} p(t) dt = F_p(x)/M.

Ponendo x → +∞ nella stessa formula si ottiene la probabilità di accettare, qualunque sia il valore: P(accetto) = 1/M. Dividendo le due si ha la distribuzione condizionata dei soli campioni accettati: P(X ≤ x | accetto) = [F_p(x)/M] / (1/M) = F_p(x). ∎

Due conseguenze operative. I campioni accettati sono esatti, qualunque M valido si usi; e il tasso di accettazione è 1/M, quindi servono in media M tentativi per campione. Un M inutilmente grande non introduce distorsione, costa solo tempo — ed è un guaio sottile, perché un metodo lento ma corretto non lancia eccezioni.

Verifichiamolo su una Beta(2, 5), la cui densità 30·x·(1−x)⁴ non ha inversa della CDF in forma chiusa, con proposta uniforme. Il massimo sta in x = 1/5 e vale M = 2,4576, quindi il tasso teorico è 1/M = 0,406901.

import numpy as np
from scipy.stats import beta

def densita_bersaglio(x):
    """Densità Beta(2, 5) scritta a mano: 30·x·(1 − x)⁴ su [0, 1]."""
    return 30.0 * x * (1.0 - x) ** 4

M = densita_bersaglio(0.2)        # massimo della densità, con proposta uniforme q(x) = 1

def campiona_per_rifiuto(tentativi, rng):
    x = rng.random(tentativi)                     # proposta q = U(0, 1)
    altezza = rng.random(tentativi) * M           # punto uniforme sotto M·q(x)
    accettati = altezza <= densita_bersaglio(x)
    return x[accettati], accettati.mean()

rng = np.random.default_rng(42)
campioni, tasso = campiona_per_rifiuto(500_000, rng)
print(f"M = {M:.6f}   tasso di accettazione teorico 1/M = {1 / M:.6f}")
print(f"tasso misurato {tasso:.6f}   campioni accettati {campioni.size}")
print(f"media    {campioni.mean():.6f}   attesa 2/7 = {2 / 7:.6f}")
print(f"varianza {campioni.var(ddof=1):.6f}   attesa     = {beta.var(2, 5):.6f}")
for q in (0.1, 0.5, 0.9):
    print(f"quantile {q:.1f}: misurato {np.quantile(campioni, q):.6f}"
          f"   esatto {beta.ppf(q, 2, 5):.6f}")

Output ottenuto (seme 42, SciPy 1.16.3):

M = 2.457600   tasso di accettazione teorico 1/M = 0.406901
tasso misurato 0.407112   campioni accettati 203556
media    0.285713   attesa 2/7 = 0.285714
varianza 0.025444   attesa     = 0.025510
quantile 0.1: misurato 0.093102   esatto 0.092595
quantile 0.5: misurato 0.264270   esatto 0.264450
quantile 0.9: misurato 0.509763   esatto 0.510316

Cosa fa il codice. Invece di estrarre u ∼ U(0,1) e confrontarlo con p(x)/(M·q(x)), estrae un'altezza uniforme in [0, M] e la confronta con p(x): è la stessa disuguaglianza moltiplicata per M, ma si legge per quello che è — un punto uniforme nel rettangolo, di cui si tiene l'ascissa se cade sotto la curva. La versione vettoriale genera tutti i tentativi in un colpo e maschera, invece di usare un while: si paga in memoria ciò che si guadagna in velocità, e il numero di campioni ottenuti non è noto in anticipo.

Che cosa ci dice. Tasso misurato 0,407112 contro teorico 0,406901, cinque cifre di accordo, e i 203 556 campioni accettati su 500 000 sono esattamente quelli che il tasso prevede. Media, varianza e tre quantili coincidono con i valori esatti di SciPy entro la terza o quarta cifra.

Perché in alta dimensione il rifiuto smette di funzionare

Il tasso di accettazione 1/M è un rapporto di volumi, e i rapporti di volumi in dimensione alta si comportano in modo controintuitivo. L'esperimento più nitido è il più semplice: estrarre punti uniformi nel cubo [−1, 1]^d e accettare quelli dentro la palla unitaria. Il tasso è il rapporto fra il volume della palla, π^{d/2}/Γ(d/2 + 1), e quello del cubo, 2^d.

import numpy as np
from scipy.special import gammaln

def tasso_teorico(d):
    """Volume della palla unitaria diviso il volume del cubo [−1, 1]^d."""
    log_palla = (d / 2) * np.log(np.pi) - gammaln(d / 2 + 1)
    return float(np.exp(log_palla - d * np.log(2.0)))

rng = np.random.default_rng(42)
n = 200_000
print("  d | tasso misurato   tasso teorico   tentativi per 1 campione")
for d in (1, 2, 5, 10, 20):
    punti = rng.uniform(-1.0, 1.0, size=(n, d))
    dentro = (punti ** 2).sum(axis=1) <= 1.0
    print(f" {d:>2} | {dentro.mean():>14.6f} {tasso_teorico(d):>15.8f}"
          f" {1 / tasso_teorico(d):>24.1f}")

Output ottenuto (seme 42):

  d | tasso misurato   tasso teorico   tentativi per 1 campione
  1 |       1.000000      1.00000000                      1.0
  2 |       0.786720      0.78539816                      1.3
  5 |       0.165875      0.16449341                      6.1
 10 |       0.002700      0.00249039                    401.5
 20 |       0.000000      0.00000002               40631627.6

Cosa fa il codice. gammaln calcola il logaritmo della funzione Gamma: si lavora in scala logaritmica perché π^{d/2} e Γ(d/2+1) diventano entrambi enormi e il loro rapporto piccolissimo, e in doppia precisione il numeratore straripa molto prima che il risultato perda senso. Il confronto è sul quadrato della norma, così si evita una radice inutile.

Che cosa ci dice. Si accetta il 79% dei tentativi in due dimensioni, il 16% in cinque, lo 0,25% in dieci — già 400 tentativi per campione. In venti il tasso teorico è 2,5·10⁻⁸ e la simulazione, su 200 000 tentativi, ne accetta zero: ne servirebbero in media 40 milioni per un solo campione. Nessuna ottimizzazione lo risolve: è il volume del cubo che si concentra negli angoli. La regola da portare via: il rifiuto è esatto, semplice ed eccellente in una o due dimensioni; oltre le cinque è morto, ed è per questo che esistono le catene di Markov Monte Carlo, che non estraggono campioni indipendenti ma costruiscono un percorso che visita le regioni ad alta densità.

Il campionamento per importanza

La riscrittura, e il peso

Nel lago dell'analogia iniziale, se i pesci che interessano stanno tutti in una insenatura, pescare a caso in tutto il lago è uno spreco: conviene pescare dove sono e poi correggere il conteggio per quanto li si è sovracampionati. È il campionamento per importanza, e la riscrittura che lo giustifica sta in una riga: si moltiplica e si divide per una densità di proposta q, scelta da noi, purché positiva ovunque f·p non sia nullo.

E_p[f(X)] = ∫ f(x) p(x) dx = ∫ f(x) · [p(x)/q(x)] · q(x) dx = E_q[ f(X) · w(X) ], con w(x) = p(x)/q(x).

In parole: si può stimare una media rispetto a p estraendo da una qualsiasi altra distribuzione q, purché si corregga ogni contributo con il peso di importanza w = p/q, grande dove q sottocampiona e piccolo dove sovracampiona. La condizione di positività non è un cavillo: se q è nulla dove f·p non lo è, quella regione non viene mai visitata e lo stimatore diventa distorto senza nessun sintomo.

La varianza è Var_q[f·w] = ∫ f(x)² p(x)²/q(x) dx − I², e quell'integrale è la ricetta per scegliere q: va messa massa dove |f|·p è grande. La proposta ottima q*(x) ∝ |f(x)|·p(x) azzera la varianza per f ≥ 0, ed è inutilizzabile perché la sua costante di normalizzazione è l'integrale cercato; serve però da bussola.

Il guadagno è spettacolare sugli eventi rari. Stimiamo P(X > 4) per X ∼ N(0,1), che vale 3,167·10⁻⁵: con Monte Carlo ingenuo, su centomila campioni, se ne attendono tre sopra la soglia. Con proposta q = N(4, 1) metà dei campioni la supera, e il peso vale p(x)/q(x) = exp(8 − 4x).

import numpy as np
from scipy.stats import norm

soglia = 4.0
esatto = float(norm.sf(soglia))          # P(X > 4) con X ~ N(0, 1)

def monte_carlo_ingenuo(n, rng):
    x = rng.standard_normal(n)
    indicatore = (x > soglia).astype(float)
    return indicatore.mean(), indicatore.std(ddof=1) / np.sqrt(n)

def monte_carlo_per_importanza(n, rng):
    """Proposta q = N(4, 1): il peso p/q vale exp(8 − 4x)."""
    x = rng.standard_normal(n) + soglia
    contributo = np.exp(8.0 - 4.0 * x) * (x > soglia)
    return contributo.mean(), contributo.std(ddof=1) / np.sqrt(n)

rng = np.random.default_rng(42)
n = 100_000
print(f"valore esatto Φ(−4) = {esatto:.6e}")
for nome, stimatore in (("ingenuo", monte_carlo_ingenuo),
                        ("per importanza", monte_carlo_per_importanza)):
    stima, errore_standard = stimatore(n, rng)
    relativo = errore_standard / stima if stima > 0 else float("inf")
    print(f"{nome:>15}: stima {stima:.6e}   errore standard {errore_standard:.3e}"
          f"   errore relativo {relativo:>8.2%}")

Output ottenuto (seme 42):

valore esatto Φ(−4) = 3.167124e-05
        ingenuo: stima 9.000000e-05   errore standard 3.000e-05   errore relativo   33.33%
 per importanza: stima 3.173136e-05   errore standard 2.127e-07   errore relativo    0.67%

Cosa fa il codice. Lo stimatore ingenuo conta quanti campioni superano la soglia: è la frequenza di un evento raro, cioè una binomiale con pochissimi successi. Quello per importanza sposta la proposta sulla soglia (rng.standard_normal(n) + soglia è un campione da N(4,1)) e pesa ogni contributo con exp(8 − 4x), ottenuto semplificando il rapporto delle due densità gaussiane.

Che cosa ci dice. Errore relativo dello 0,67% contro il 33%, cioè una varianza cinquanta volte più piccola a parità di campioni. E c'è di più: la stima ingenua vale 9,0·10⁻⁵, quasi tre volte il valore vero, perché ha pescato nove campioni oltre la soglia invece dei tre attesi; il suo intervallo al 95% [3,0·10⁻⁵, 1,5·10⁻⁴] contiene il valore vero per un soffio. È uno stimatore corretto che, su un evento raro, produce un numero su cui nessuno dovrebbe decidere.

Quando la varianza esplode: code troppo leggere

Il modo in cui il campionamento per importanza fallisce è particolarmente sgradevole: lo stimatore resta non distorto ma acquista varianza infinita, e una varianza infinita non si vede guardando la stima. La condizione sta in ∫ f²p²/q: se q ha code più leggere di f·p, il rapporto p/q diverge nelle code più in fretta di quanto q le smorzi, e l'integrale non converge.

In numeri: bersaglio p = N(0,1), f(x) = x², valore esatto 1, proposta q = N(0, σ_q²) con peso σ_q · exp(x²(1/(2σ_q²) − 1/2)). Per σ_q = 1,5 l'esponente è negativo e tutto va bene; per σ_q = 0,5 vale +1,5x² e il secondo momento contiene ∫ x⁴ e^{x²} dx, che diverge.

import numpy as np

def contributi_per_importanza(n, sigma_proposta, rng):
    """Bersaglio p = N(0, 1) e f(x) = x², quindi E_p[f] = 1 esattamente."""
    x = rng.standard_normal(n) * sigma_proposta
    log_peso = -x ** 2 / 2.0 + x ** 2 / (2.0 * sigma_proposta ** 2) + np.log(sigma_proposta)
    return x ** 2 * np.exp(log_peso)

rng = np.random.default_rng(42)
for sigma in (1.5, 0.5):
    contributi = contributi_per_importanza(200_000, sigma, rng)
    corse = np.cumsum(contributi) / np.arange(1, contributi.size + 1)
    print(f"proposta N(0, {sigma}²)   code {'più pesanti' if sigma > 1 else 'più leggere'} del bersaglio")
    for n in (1_000, 10_000, 50_000, 100_000, 200_000):
        print(f"   n = {n:>7}: stima {corse[n - 1]:.6f}")
    print(f"   contributo massimo osservato {contributi.max():.3e}"
          f"   varianza campionaria {contributi.var(ddof=1):.3e}")

Output ottenuto (seme 42):

proposta N(0, 1.5²)   code più pesanti del bersaglio
   n =    1000: stima 0.967729
   n =   10000: stima 0.990041
   n =   50000: stima 0.996901
   n =  100000: stima 0.997347
   n =  200000: stima 0.998902
   contributo massimo osservato 1.987e+00   varianza campionaria 4.913e-01
proposta N(0, 0.5²)   code più leggere del bersaglio
   n =    1000: stima 0.904171
   n =   10000: stima 0.921567
   n =   50000: stima 0.812826
   n =  100000: stima 0.903204
   n =  200000: stima 0.847854
   contributo massimo osservato 9.829e+03   varianza campionaria 6.363e+02

Cosa fa il codice. Il peso si calcola in scala logaritmica e si esponenzia una volta sola: con σ_q = 0,5 l'esponente arriva a 1,5·x², che per x ≈ 2,3 vale già e⁸. np.cumsum diviso per gli indici produce la media corrente: il grafico della convergenza, letto come tabella.

Che cosa ci dice, ed è il punto della sezione. Con la proposta larga la media corrente sale ordinatamente verso 1 (0,9677 → 0,9900 → 0,9969 → 0,9973 → 0,9989), contributo massimo 1,99 e varianza 0,49. Con la proposta stretta non converge: 0,904 → 0,922 → 0,813 → 0,903 → 0,848. Guarda il passaggio da 10 000 a 50 000: la stima peggiora quintuplicando i campioni. Il colpevole è il contributo massimo, 9,8·10³, cinquemila volte quello della proposta buona: un singolo campione vale più di tutti gli altri insieme, e la media salta ogni volta che ne arriva uno.

Il sintomo e il rimedio. Una stima che sembra convergere e poi salta, o una varianza campionaria che cresce con n, significa quasi sempre code della proposta più leggere del bersaglio. Stampa sempre il peso massimo e la dimensione campionaria efficace ESS = (Σwᵢ)²/Σwᵢ²: se su duecentomila campioni l'ESS è nell'ordine delle decine, non hai duecentomila campioni, ne hai qualche decina. In caso di dubbio si sceglie una proposta con code più pesanti del bersaglio: un peso che tende a zero nelle code costa efficienza, un peso che esplode costa la validità della stima.

Riduzione della varianza: antitetiche e variate di controllo

Siccome l'errore scende come 1/√n, aumentare n è il modo più caro di diventare precisi. L'alternativa è abbassare σ: costruire uno stimatore della stessa quantità con varianza minore. Due tecniche costano quasi nulla.

Variabili antitetiche. Invece di n campioni indipendenti si usano n/2 coppie (U, 1−U). La varianza di una media di coppia è Var[(f(U) + f(1−U))/2] = [Var f + Cov]/2, perché i due termini hanno la stessa distribuzione. Se f è monotona, f(U) e f(1−U) sono anticorrelate, Cov < 0, e la varianza scende — gratis, perché la coppia costa un numero casuale in meno, non in più.

Variate di controllo. Si sfrutta una g di cui si conosce esattamente la media, stimando f(x) − c·(g(x) − E[g]). Il termine aggiunto ha media zero, quindi lo stimatore resta non distorto per qualunque c, che agisce solo sulla varianza Var f − 2c·Cov(f,g) + c²·Var g. È una parabola in c, minima in c* = Cov(f,g)/Var(g); sostituendo si ottiene Var f · (1 − ρ²). Più g somiglia a f, più si guadagna.

Misuriamole su f(u) = e^u con u uniforme in [0,1] — il costo di un'operazione che cresce esponenzialmente con il carico, dati sintetici inventati da noi — il cui valore esatto è e − 1. Come variata di controllo g(u) = u, con E[g] = 1/2. Il confronto è a parità di valutazioni di f: altrimenti le antitetiche barerebbero.

import math
import numpy as np

def costo(u):
    """Costo di un'operazione che cresce esponenzialmente con il carico u ∈ [0, 1]."""
    return np.exp(u)

esatto = math.e - 1.0
rng = np.random.default_rng(42)
n = 200_000                           # valutazioni di f: le stesse per i tre stimatori

semplice = costo(rng.random(n))

v = rng.random(n // 2)
antitetico = (costo(v) + costo(1.0 - v)) / 2.0

u = rng.random(n)
c = float(np.cov(costo(u), u, ddof=1)[0, 1] / u.var(ddof=1))
controllo = costo(u) - c * (u - 0.5)

print(f"valore esatto e − 1 = {esatto:.6f}")
print("stimatore                   stima   varianza per valutazione   errore standard")
teoriche = {"semplice": 0.2420356, "antitetico": 0.0078250, "variata di controllo": 0.0039401}
for nome, campione, valutazioni in (("semplice", semplice, 1),
                                    ("antitetico", antitetico, 2),
                                    ("variata di controllo", controllo, 1)):
    varianza = campione.var(ddof=1) * valutazioni
    print(f"{nome:>20} {campione.mean():>11.6f} {varianza:>26.6f}"
          f" {campione.std(ddof=1) / np.sqrt(campione.size):>17.2e}"
          f"   (teorica {teoriche[nome]:.6f})")
print(f"coefficiente stimato c = {c:.6f}   ottimo teorico 12·Cov(e^U, U) = {12 * (1 - 0.5 * esatto):.6f}")

Output ottenuto (seme 42):

valore esatto e − 1 = 1.718282
stimatore                   stima   varianza per valutazione   errore standard
            semplice    1.718252                   0.242112          1.10e-03   (teorica 0.242036)
          antitetico    1.718174                   0.007763          1.97e-04   (teorica 0.007825)
variata di controllo    1.718315                   0.003945          1.40e-04   (teorica 0.003940)
coefficiente stimato c = 1.689440   ottimo teorico 12·Cov(e^U, U) = 1.690309

Cosa fa il codice. Le tre varianze sono riportate per valutazione di f: lo stimatore antitetico ne usa due per ogni media di coppia, quindi la sua varianza va moltiplicata per due prima del confronto, altrimenti il confronto è truccato a suo favore di un fattore esattamente 2. Il coefficiente c è stimato dai dati con np.cov, come si fa quando Cov(f,g) non si sa calcolare a mano.

Che cosa ci dice. A parità di 200 000 valutazioni: 0,2421 per lo stimatore semplice, 0,00776 per l'antitetico (31 volte meno), 0,00395 per la variata di controllo (61 volte meno). Tradotto: per eguagliare con la forza bruta la precisione che la variata di controllo ottiene con 200 000 valutazioni ne servirebbero dodici milioni. Le tre varianze misurate centrano quelle teoriche alla terza cifra, e il coefficiente stimato (1,68944) centra l'ottimo (1,690309).

Trade-off. Le antitetiche richiedono f monotona, o almeno prevalentemente tale: su una f simmetrica rispetto a 1/2 la correlazione diventa positiva e la varianza aumenta. Le variate di controllo richiedono una g correlata di media nota, che non sempre esiste; e stimare c sugli stessi campioni usati per la media introduce una distorsione trascurabile su n grande ma non su n piccolo, dove conviene un campione pilota separato.

Il trucco della riparametrizzazione

Il problema: derivare rispetto a chi genera i campioni

Finora la distribuzione da cui si campiona era data. In quasi tutto l'apprendimento automatico moderno non lo è: la si vuole ottimizzare. Serve cioè il gradiente di

L(θ) = E_{x ∼ p_θ}[ f(x) ] rispetto a θ.

E qui lo scambio che si fa sempre — derivata dentro l'aspettativa — non è lecito. ∇_θ E_{x∼p_θ}[f(x)] non è E_{x∼p_θ}[∇_θ f(x)], e il motivo si vede scrivendo ∇_θ ∫ f(x) p_θ(x) dx: la dipendenza da θ non sta in f, che spesso non contiene θ affatto, sta nella misura. Derivando f si deriva la cosa sbagliata, e si ottiene zero.

Prima soluzione: lo stimatore score function (REINFORCE)

Si deriva sotto il segno di integrale rispetto alla densità: ∇_θ ∫ f p_θ dx = ∫ f ∇_θ p_θ dx. Ma ∇_θ p_θ non è un'aspettativa, perché non è qualcosa · p_θ. Si rimedia con l'identità ∇_θ p_θ = p_θ · ∇_θ log p_θ, che è la derivata del logaritmo riletta al contrario:

∇_θ L(θ) = E_{x ∼ p_θ}[ f(x) · ∇_θ log p_θ(x) ].

È lo stimatore score function, noto nell'apprendimento per rinforzo come REINFORCE (Williams, 1992). Funziona sempre, anche con f non derivabile o discreta, perché f non viene mai derivata. Il prezzo è pesante: f(x) moltiplica un termine a media nulla e rumoroso, quindi la varianza è alta e cresce col valore assoluto di f. Si mitiga sottraendo a f una costante b (il baseline) — che è esattamente una variata di controllo — ma resta alta.

Seconda soluzione: la riparametrizzazione

L'idea è spostare la casualità fuori dal percorso del gradiente. Se x ∼ N(μ, σ²) si può scrivere x = μ + σ·ε con ε ∼ N(0,1): la distribuzione di ε non dipende da θ = (μ, σ), e tutta la dipendenza dal parametro è in una funzione deterministica. Allora

∇_θ E_{x∼p_θ}[f(x)] = ∇_θ E_{ε∼N(0,1)}[ f(μ + σε) ] = E_{ε}[ ∇_θ f(μ + σε) ],

e adesso lo scambio è lecito, perché la misura su cui si integra non dipende più da θ: per μ si ottiene E_ε[f′(μ + σε)], per σ si ottiene E_ε[ε·f′(μ + σε)]. È il trucco della riparametrizzazione (Kingma e Welling, 2014), ed è ciò che rende addestrabile un variational autoencoder: senza, il campionamento nello spazio latente interromperebbe la retropropagazione. Il prezzo è che serve f derivabile e una riparametrizzazione esplicita — condizioni che la gaussiana soddisfa e una distribuzione discreta no.

Confrontiamo le varianze dei due stimatori su ∇_μ E_{x∼N(μ,1)}[x²] con μ = 1,5, il cui valore vero è 2μ = 3. Le varianze teoriche si calcolano a mano: μ⁴ + 14μ² + 15 = 51,5625 per lo score function, 28 per lo score con baseline, 4σ² = 4 per la riparametrizzazione.

import numpy as np

mu, sigma = 1.5, 1.0
gradiente_vero = 2.0 * mu                # d/dμ E[x²] = d/dμ (μ² + σ²) = 2μ

def stimatore_score_function(n, rng):
    """REINFORCE: f(x) · ∂log p(x)/∂μ, con ∂log p/∂μ = (x − μ)/σ²."""
    x = mu + sigma * rng.standard_normal(n)
    return x ** 2 * (x - mu) / sigma ** 2

def stimatore_score_con_baseline(n, rng):
    """Lo stesso, con E[f] sottratta a f: è una variata di controllo."""
    x = mu + sigma * rng.standard_normal(n)
    return (x ** 2 - (mu ** 2 + sigma ** 2)) * (x - mu) / sigma ** 2

def stimatore_riparametrizzato(n, rng):
    """x = μ + σε con ε ∼ N(0,1): ∂f/∂μ = f′(x) · ∂x/∂μ = 2(μ + σε)."""
    eps = rng.standard_normal(n)
    return 2.0 * (mu + sigma * eps)

rng = np.random.default_rng(42)
teoriche = {"score function": mu ** 4 + 14 * mu ** 2 + 15,
            "score + baseline": 28.0,
            "riparametrizzazione": 4.0 * sigma ** 2}
print(f"gradiente vero = {gradiente_vero:.6f}")
print("stimatore                  media su 2·10⁶   varianza misurata   varianza teorica")
for nome, f in (("score function", stimatore_score_function),
                ("score + baseline", stimatore_score_con_baseline),
                ("riparametrizzazione", stimatore_riparametrizzato)):
    campione = f(2_000_000, rng)
    print(f"{nome:>20} {campione.mean():>17.5f} {campione.var(ddof=1):>19.4f}"
          f" {teoriche[nome]:>18.4f}")

print("\nerrore standard del gradiente, a parità di n")
print("      n   score function   riparametrizzazione")
for n in (10, 100, 1000, 10000):
    a = stimatore_score_function(n, rng)
    b = stimatore_riparametrizzato(n, rng)
    print(f"{n:>7} {a.std(ddof=1) / np.sqrt(n):>16.5f}"
          f" {b.std(ddof=1) / np.sqrt(n):>21.5f}")

Output ottenuto (seme 42):

gradiente vero = 3.000000
stimatore                  media su 2·10⁶   varianza misurata   varianza teorica
      score function           3.00276             51.7719            51.5625
    score + baseline           2.99788             27.9164            28.0000
 riparametrizzazione           3.00128              3.9978             4.0000

errore standard del gradiente, a parità di n
      n   score function   riparametrizzazione
     10          1.57648               0.51661
    100          1.02247               0.20037
   1000          0.25905               0.06371
  10000          0.07301               0.01981

Cosa fa il codice. Nessuna libreria di differenziazione automatica: i due gradienti sono scritti a mano, ed è il modo migliore per vedere che sono due formule diverse per lo stesso numero. Lo score function moltiplica f(x) = x² per lo score (x − μ)/σ²; la riparametrizzazione deriva f e la valuta in μ + σε. Entrambe stimano 3, e lo fanno: 3,00276 e 3,00128.

Che cosa ci dice. Le varianze misurate centrano quelle teoriche alla terza cifra (51,77 contro 51,5625; 3,998 contro 4,000), il che conferma insieme la derivazione e il codice. Il rapporto è 12,9 a favore della riparametrizzazione, e il baseline recupera meno della metà del divario (51,6 → 28,0): utile, non sostitutivo. La seconda tabella traduce la varianza in quello che si sente addestrando: con dieci campioni lo score function ha un errore standard di 1,58 su un valore vero di 3, il 53% — il segno è giusto, la grandezza è un'opinione; la riparametrizzazione sta allo 0,52, cioè al 17%.

Trade-off, ed è netto. La riparametrizzazione ha varianza molto minore ma richiede f derivabile e una distribuzione riparametrizzabile. Lo score function non richiede nulla di ciò, ed è per questo l'unico utilizzabile su azioni discrete. Non è che una sia migliore: coprono casi diversi, e la seconda si paga con molti più campioni.

L'esperimento: una coda a 5 stati

Mettiamo insieme tutto su un sistema piccolo ma non banale, interamente sintetico: uno sportello con una sala d'attesa da quattro posti, quindi cinque stati possibili (0, 1, 2, 3, 4 utenti nel sistema). Il tempo è discreto; in ogni passo, prima il servizio e poi l'arrivo: se c'è almeno un utente, con probabilità s = 0,5 viene servito e se ne va; dopodiché, con probabilità a = 0,1 arriva un nuovo utente, che entra se c'è posto e viene respinto se il sistema è pieno. Le due probabilità sono inventate da noi e dichiarate tali: servono a far comparire i fenomeni, non a descrivere uno sportello reale.

Il vantaggio di un sistema così piccolo è che ha una risposta esatta con cui confrontarsi, e un sistema vero non ce l'ha. La catena è di nascita e morte: da uno stato n ≥ 1 si sale con probabilità (1−s)·a = 0,05 e si scende con s·(1−a) = 0,45, mentre da 0 si sale con a = 0,1. Il bilancio dettagliato π_n · (su) = π_{n+1} · (giù) dà i rapporti π₁/π₀ = 0,1/0,45 e π_{n+1}/π_n = 0,05/0,45 = 1/9 per n ≥ 1; normalizzando si ottiene la distribuzione stazionaria, e da lì il numero medio nel sistema L = Σ n·π_n. La legge di Little L = λ·W restituisce infine l'attesa media W, dove λ è il tasso di arrivi accettati, cioè a meno quelli respinti.

import numpy as np
from collections import deque

ARRIVO, SERVIZIO, CAPIENZA = 0.1, 0.5, 4

def riferimento_esatto():
    """Distribuzione stazionaria della catena a 5 stati, e legge di Little."""
    su = (1 - SERVIZIO) * ARRIVO                   # da n a n+1, per n ≥ 1
    giu = SERVIZIO * (1 - ARRIVO)                  # da n a n−1, per n ≥ 1
    rapporti = [1.0, ARRIVO / giu]                 # il primo salto parte da 0
    for _ in range(3):
        rapporti.append(rapporti[-1] * su / giu)
    pi = np.array(rapporti) / sum(rapporti)
    numero_medio = float(np.arange(5) @ pi)
    arrivi_accettati = ARRIVO - pi[4] * su         # gli arrivi respinti non entrano
    return pi, numero_medio, numero_medio / arrivi_accettati

pi, numero_medio, attesa_esatta = riferimento_esatto()
print("distribuzione stazionaria:", np.round(pi, 6))
print(f"numero medio nel sistema L = {numero_medio:.6f}"
      f"   attesa media esatta W = L/λ = {attesa_esatta:.6f}")

Output ottenuto:

distribuzione stazionaria: [8.00024e-01 1.77783e-01 1.97540e-02 2.19500e-03 2.44000e-04]
numero medio nel sistema L = 0.224851   attesa media esatta W = L/λ = 2.248780

Cosa fa il codice. Costruisce i rapporti fra probabilità stazionarie consecutive invece di risolvere il sistema lineare: per una catena di nascita e morte il bilancio dettagliato li dà uno alla volta, e la normalizzazione finale si fa dividendo per la somma. np.arange(5) @ pi è il prodotto scalare fra i valori di stato e le loro probabilità, cioè il valore atteso.

Che cosa ci dice. Lo sportello è vuoto l'80% del tempo e pieno lo 0,024%. L'attesa media esatta è 2,2488 passi, da confrontare con i 2 passi che costerebbe il solo servizio se non ci fosse mai coda (la durata del servizio è geometrica di media 1/s = 2): la coda aggiunge appena un quarto di passo. Questo è il numero che la simulazione deve ritrovare, e non ritrovarlo significa che la simulazione è sbagliata.

Ora la stima per simulazione, con l'intervallo di confidenza al 95% dello stimatore — cioè stima ± 1,96 · errore standard, dove l'errore standard è calcolato sulla dispersione delle medie di ogni ripetizione, non su quella delle singole attese. Due configurazioni: corse brevi da 400 passi, e corse lunghe da 5000 passi di cui i primi 200 scartati come riscaldamento.

def simula_coda(passi, rng, scarto=0):
    """Coda a 5 stati (0..4): dentro ogni passo prima il servizio, poi l'arrivo.
    L'attesa di un cliente è il numero di fine-passo in cui risulta nel sistema.
    Si contano solo i clienti arrivati dopo `scarto` e usciti entro la fine:
    quelli ancora dentro all'ultimo passo sono troncati, e vanno scartati."""
    in_sistema, attese = deque(), []
    for t in range(passi):
        if in_sistema and rng.random() < SERVIZIO:
            arrivato = in_sistema.popleft()
            if arrivato >= scarto:
                attese.append(t - arrivato)
        if rng.random() < ARRIVO and len(in_sistema) < CAPIENZA:
            in_sistema.append(t)
    return np.array(attese, dtype=float)

rng = np.random.default_rng(42)
for etichetta, passi, scarto, lista in (
        ("corse brevi, nessun riscaldamento (400 passi)", 400, 0, (25, 100, 400, 1600)),
        ("corse lunghe, 200 passi scartati (5000 passi)", 5000, 200, (25, 50, 100, 400))):
    print(f"\n{etichetta}")
    print("ripetizioni   attesa media   errore standard   IC 95%            copre W?")
    for ripetizioni in lista:
        medie = np.array([simula_coda(passi, rng, scarto).mean() for _ in range(ripetizioni)])
        stima = medie.mean()
        errore = medie.std(ddof=1) / np.sqrt(ripetizioni)
        basso, alto = stima - 1.96 * errore, stima + 1.96 * errore
        print(f"{ripetizioni:>11} {stima:>14.5f} {errore:>17.5f}"
              f"   [{basso:.5f}, {alto:.5f}]  {'sì' if basso <= attesa_esatta <= alto else 'NO'}")

Output ottenuto (seme 42):

corse brevi, nessun riscaldamento (400 passi)
ripetizioni   attesa media   errore standard   IC 95%            copre W?
         25        2.24170           0.06869   [2.10706, 2.37634]  sì
        100        2.25600           0.04361   [2.17052, 2.34147]  sì
        400        2.22397           0.01757   [2.18952, 2.25841]  sì
       1600        2.22711           0.00893   [2.20961, 2.24461]  NO

corse lunghe, 200 passi scartati (5000 passi)
ripetizioni   attesa media   errore standard   IC 95%            copre W?
         25        2.24896           0.01883   [2.21206, 2.28586]  sì
         50        2.24007           0.01480   [2.21106, 2.26908]  sì
        100        2.23429           0.01117   [2.21240, 2.25617]  sì
        400        2.24940           0.00542   [2.23878, 2.26001]  sì

Cosa fa il codice. La deque contiene gli istanti di arrivo dei clienti presenti, in ordine: il servizio è a disciplina FIFO, quindi esce sempre il primo della fila con popleft. L'attesa registrata è t − arrivato, cioè il numero di conteggi di fine passo in cui quel cliente risulta nel sistema — ed è esattamente la grandezza che compare nella legge di Little. Il parametro scarto esclude dal conteggio i clienti arrivati durante il riscaldamento. I clienti ancora in coda all'ultimo passo non vengono mai registrati: è un troncamento inevitabile, e la sua conseguenza si vede nella prima tabella.

Che cosa ci dice, ed è il momento centrale della lezione. Guarda la colonna «copre W?» della prima tabella. A 25, 100 e 400 ripetizioni l'intervallo di confidenza contiene il valore vero 2,24878. A 1600 ripetizioni smette di contenerlo: l'intervallo è [2,20961, 2,24461] e il valore vero cade appena fuori. Non è sfortuna e non è un errore di calcolo: le corse da 400 passi producono una stima distorta verso il basso di circa 0,02, per due motivi che agiscono nella stessa direzione — si parte dal sistema vuoto e i primi clienti non aspettano (transitorio), e i clienti ancora dentro alla fine non vengono contati, benché siano proprio quelli con l'attesa più lunga (troncamento).

Leggila bene, perché è controintuitiva. Finché le ripetizioni erano poche, l'intervallo era largo e la distorsione ci stava comodamente dentro: il protocollo sembrava corretto. Aumentando le ripetizioni l'intervallo si è stretto attorno al numero sbagliato, e la distorsione, che c'era sempre stata, è diventata visibile solo quando la precisione è cresciuta abbastanza da smascherarla. Più dati non correggono una distorsione: la rendono solo più evidente. La seconda tabella mostra il rimedio — corse dodici volte più lunghe con i primi 200 passi buttati via — e lì l'intervallo copre il valore vero a tutte e quattro le numerosità, compresa quella più stretta (±0,0054).

La saturazione: un evento raro, e due modi di misurarlo

Seconda domanda sullo stesso sistema, e stavolta su un evento raro: partendo dal sistema vuoto, qual è la probabilità che la sala d'attesa si riempia (stato 4) entro otto passi? Servono almeno quattro arrivi in otto passi, ciascuno con probabilità 0,1, e quasi nessun servizio nel frattempo. Anche qui esiste la risposta esatta, per programmazione dinamica sulla catena con lo stato 4 reso assorbente.

import math
import numpy as np

ORIZZONTE = 8

def nucleo(peso_arrivo, peso_non_arrivo):
    """Matrice 4×4 sugli stati 0..3 e peso della fuga verso lo stato 4.
    Con (ARRIVO, 1−ARRIVO) dà la catena vera; con (a²/a′, (1−a)²/(1−a′))
    dà la somma dei p²/q, che è ciò che serve per il secondo momento."""
    Q = np.zeros((4, 4))
    Q[0, 0], Q[0, 1] = peso_non_arrivo, peso_arrivo
    for n in (1, 2, 3):
        Q[n, n - 1] = SERVIZIO * peso_non_arrivo
        Q[n, n] = SERVIZIO * peso_arrivo + (1 - SERVIZIO) * peso_non_arrivo
        if n + 1 <= 3:
            Q[n, n + 1] = (1 - SERVIZIO) * peso_arrivo
    return Q, (1 - SERVIZIO) * peso_arrivo

def somma_sui_percorsi(peso_arrivo, peso_non_arrivo):
    """Somma sui percorsi che toccano lo stato 4 entro ORIZZONTE passi.
    `coda` è il fattore dei passi che restano DOPO l'assorbimento: il peso
    continua ad accumularsi anche quando l'evento è già deciso."""
    Q, fuga = nucleo(peso_arrivo, peso_non_arrivo)
    coda = peso_arrivo + peso_non_arrivo
    u, totale = np.array([1.0, 0.0, 0.0, 0.0]), 0.0
    for k in range(ORIZZONTE):
        totale += fuga * u[3] * coda ** (ORIZZONTE - k - 1)
        u = u @ Q
    return totale

def varianza_teorica(arrivo_proposto):
    """Var_q[w·1_A] = E_q[w²·1_A] − P(A)², con E_q[w²·1_A] = Σ p²/q sui percorsi."""
    secondo = somma_sui_percorsi(ARRIVO ** 2 / arrivo_proposto,
                                 (1 - ARRIVO) ** 2 / (1 - arrivo_proposto))
    return secondo - esatto ** 2

def episodi(n, arrivo_proposto, rng):
    """n episodi di ORIZZONTE passi da sistema vuoto, in parallelo.
    Restituisce i contributi 1{ha toccato lo stato 4} · peso del percorso."""
    stato = np.zeros(n, dtype=np.int64)
    saturato, log_peso = np.zeros(n, dtype=bool), np.zeros(n)
    sopra = math.log(ARRIVO / arrivo_proposto)
    sotto = math.log((1 - ARRIVO) / (1 - arrivo_proposto))
    for _ in range(ORIZZONTE):
        stato -= (stato > 0) & (rng.random(n) < SERVIZIO)
        arrivo = rng.random(n) < arrivo_proposto
        log_peso += np.where(arrivo, sopra, sotto)
        entra = arrivo & (stato < CAPIENZA)
        stato += entra
        saturato |= entra & (stato == CAPIENZA)
    return saturato * np.exp(log_peso)

esatto = somma_sui_percorsi(ARRIVO, 1 - ARRIVO)      # pesi veri: è la probabilità
rng = np.random.default_rng(42)
print(f"probabilità esatta di saturazione in {ORIZZONTE} passi = {esatto:.6e}")
print("stimatore                   n        stima   errore standard"
      "   errore relativo   episodi utili   varianza teorica")
misurate = {}
for nome, n, arrivo_proposto in (("ingenuo", 20_000, ARRIVO),
                                 ("ingenuo", 2_000_000, ARRIVO),
                                 ("per importanza", 20_000, 0.5)):
    contributi = episodi(n, arrivo_proposto, rng)
    stima, varianza = contributi.mean(), contributi.var(ddof=1)
    misurate[nome] = varianza
    print(f"{nome:>15} {n:>11} {stima:>12.6e} {math.sqrt(varianza / n):>17.3e}"
          f" {math.sqrt(varianza / n) / stima if stima > 0 else float('nan'):>17.1%}"
          f" {int((contributi > 0).sum()):>15} {varianza_teorica(arrivo_proposto):>18.6e}")
print(f"rapporto fra le varianze, misurato  = "
      f"{misurate['ingenuo'] / misurate['per importanza']:.1f}")
print(f"rapporto fra le varianze, teorico   = "
      f"{varianza_teorica(ARRIVO) / varianza_teorica(0.5):.1f}")

Output ottenuto (seme 42):

probabilità esatta di saturazione in 8 passi = 2.217811e-04
stimatore                   n        stima   errore standard   errore relativo   episodi utili   varianza teorica
        ingenuo       20000 2.000000e-04         9.999e-05             50.0%               4       2.217319e-04
        ingenuo     2000000 2.185000e-04         1.045e-05              4.8%             437       2.217319e-04
 per importanza       20000 2.125971e-04         1.128e-05              5.3%            1899       2.670085e-06
rapporto fra le varianze, misurato  = 85.8
rapporto fra le varianze, teorico   = 83.0

Cosa fa il codice. somma_sui_percorsi è una sola funzione usata per due scopi diversi, ed è la parte da leggere con attenzione. Con i pesi veri (a, 1−a) somma le probabilità dei percorsi che toccano lo stato 4, cioè restituisce P(A) esatta. Con i pesi (a²/a′, (1−a)²/(1−a′)) somma invece p²/q, che è per definizione E_q[w²·1_A], cioè il secondo momento dello stimatore per importanza. Il fattore coda conta i passi che restano dopo l'assorbimento: il simulatore non si ferma quando l'evento si decide, e il peso continua a moltiplicarsi fino alla fine dell'orizzonte. Con i pesi veri coda = a + (1−a) = 1 e il fattore sparisce; con i pesi quadratici vale 0,02 + 1,62 = 1,64 e non sparisce. episodi simula tutti gli episodi in parallelo con operazioni su array, accumulando il logaritmo del peso per non perdere precisione.

Che cosa ci dice. Con 20 000 episodi lo stimatore ingenuo trova quattro saturazioni in tutto. La sua stima (2,0·10⁻⁴) non è nemmeno lontanissima dal vero, ma il suo errore relativo dichiarato è del 50% — e quel 50% è a sua volta calcolato su quattro eventi, quindi è una misura di cui non ci si può fidare. È la condizione peggiore in cui uno stimatore possa trovarsi: è così impreciso da non riuscire a misurare quanto è impreciso. Lo stimatore per importanza, con lo stesso numero di episodi, ne ha 1899 di utili e arriva al 5,3%. Per eguagliarlo, quello ingenuo ha bisogno di due milioni di episodi: cento volte il lavoro, per lo stesso risultato.

Il rapporto fra le varianze, e perché lo riportiamo due volte. Misurato 85,8, teorico 83,0: uno scarto del 3,4%, dentro l'incertezza con cui si può stimare una varianza su 437 eventi utili. Il valore teorico non è un'asserzione presa da un libro, lo calcola la funzione qui sopra. Vale la pena dire che una prima versione di questa derivazione dava 137,7 — sbagliato, perché sommava i pesi solo fino all'istante dell'assorbimento invece che su tutto l'orizzonte. L'errore era invisibile nella formula e si è manifestato solo come disaccordo con la misura: il tipo di controllo che si perde se si riporta un solo numero.

Quando Monte Carlo è la scelta sbagliata

Nessuna tecnica di questa lezione è universalmente migliore, e la sezione che i corsi omettono è questa.

La regola da portare via. Monte Carlo è il metodo che si usa quando gli altri non si applicano, e quasi sempre non si applicano. Ma «quasi sempre» non è «sempre», e il costo di non accorgersene è un fattore dieci milioni.

Errori tipici (sintomo, causa, rimedio)

1. Seme non fissato. Sintomo: due esecuzioni danno numeri diversi e nessuno sa quale riportare. Rimedio: rng = np.random.default_rng(42), generatore passato come argomento, seme scritto accanto al risultato.
2. Usare lo stato globale di NumPy. Sintomo: il risultato cambia aggiornando una dipendenza o spostando una riga. Causa: np.random.seed() e np.random.rand() condividono uno stato con tutto il processo. Rimedio: default_rng, e SeedSequence.spawn per i flussi paralleli.
3. Riportare la deviazione standard dei dati invece di quella dello stimatore. Sintomo: una barra d'errore che non si restringe aumentando i campioni. Causa: si è stampato σ̂ invece di σ̂/√n. Rimedio: se il numero non scende con n, non è un errore standard.
4. Fidarsi di un intervallo di confidenza su uno stimatore distorto. Sintomo: l'intervallo copre il valore vero finché le ripetizioni sono poche e smette di coprirlo quando aumentano — la riga a 1600 ripetizioni dell'esperimento sulla coda. Causa: l'intervallo misura l'errore attorno alla media dello stimatore, non la distanza dalla verità. Rimedio: scartare un riscaldamento, non troncare le osservazioni lunghe, validare su un caso con risposta nota.
5. Proposta con code più leggere del bersaglio. Sintomo: la stima sembra convergere e poi salta; la varianza campionaria cresce con n. Causa: ∫ f²p²/q diverge, lo stimatore ha varianza infinita pur restando non distorto. Rimedio: stampare peso massimo ed ESS; nel dubbio, code più pesanti.
6. Sommare i pesi di importanza solo fino al tempo d'arresto. Sintomo: la teoria predice una riduzione di varianza sistematicamente maggiore di quella misurata, e lo scarto non si chiude aumentando i campioni. Causa: il rapporto di verosimiglianza continua ad accumularsi dopo che l'evento è stato deciso, e su una mossa la somma dei p²/q non vale 1. Rimedio: sommare i pesi su tutto l'orizzonte effettivamente simulato.
7. Pesi non normalizzati, o normalizzati quando non si doveva. Sintomo: la stima è sbagliata di un fattore costante, uguale a ogni esecuzione. Causa: p o q note solo a meno della costante. Rimedio: lo stimatore grezzo (1/n)Σwᵢf(xᵢ) è non distorto ma richiede densità esatte; quello autonormalizzato Σwᵢf(xᵢ)/Σwᵢ tollera le costanti ignote ma è distorto di ordine 1/n.
8. Usare il rifiuto in alta dimensione. Sintomo: il programma gira per ore e produce pochissimi campioni, senza errori. Causa: 1/M crolla esponenzialmente — 0,25% a dieci dimensioni, 2,5·10⁻⁸ a venti. Rimedio: sotto l'1% di accettazione si cambia metodo, non macchina.
9. Scegliere n «perché sembrava abbastanza». Sintomo: nessuno sa dire che precisione ha il risultato. Rimedio: decidere prima la precisione e, stimare σ su una corsa pilota, porre n ≥ (σ/e)². È la stessa formula letta al contrario, e trasforma una simulazione in una misura.

I nove errori si dividono in due famiglie. Cinque riguardano la dichiarazione dell'incertezza (seme, stato globale, errore standard confuso con la deviazione standard, intervallo su uno stimatore distorto, n scelto a occhio) e si evitano con una verifica su un caso a risposta nota. Quattro riguardano proposte e pesi (code troppo leggere, tempo d'arresto, normalizzazione, rifiuto in alta dimensione) e si riconoscono dagli indicatori che la lezione fa stampare sempre: peso massimo, ESS, tasso di accettazione.

Best practice

• Dichiara il seme e le versioni. Senza seme un risultato non è riproducibile; senza versione non lo è a distanza di un anno.
• Riporta l'errore standard accanto alla stima. Un numero da solo non è una misura; se l'errore standard non si sa calcolare, si usa il bootstrap.
• Dimensiona n a partire dall'errore voluto, con n ≥ (σ/e)² da una corsa pilota, invece di scegliere una potenza di dieci.
• Valida su un caso con risposta nota prima di usarne uno senza. Senza riferimento non si distingue un bug da un risultato.
• Sorveglia i pesi: peso massimo, ESS e tasso di accettazione vanno stampati sempre, non solo quando qualcosa sembra storto.
• Prova la riduzione di varianza prima di aumentare n. Qui una variata di controllo da tre righe ha reso 61 volte: la stessa precisione con la forza bruta costerebbe sessantuno volte i campioni.

Le sei regole sono ordinate dalla più economica alla più impegnativa: seme, errore standard e dimensionamento di n costano una riga, la validazione su un caso noto costa una mezz'ora, la riduzione della varianza richiede di pensare allo stimatore. Se parti da zero, applica nell'ordine in cui sono scritte.

Glossario minimo

• Stimatore non distorto (unbiased estimator): il suo valore atteso è esattamente la quantità cercata. Non garantisce nulla su una singola esecuzione.
• Errore standard (standard error): deviazione standard dello stimatore, σ/√n. Non è σ, che è la dispersione dei dati.
• Trasformata inversa (inverse transform sampling): X = F⁻¹(U) con U uniforme. Esatta e senza scarti, richiede F⁻¹ in forma chiusa.
• Campionamento per rifiuto (rejection sampling): proporre da q, accettare con probabilità p/(Mq). Esatto, al costo di M tentativi per campione.
• Proposta (proposal): la distribuzione q da cui si sa estrarre, usata al posto del bersaglio p.
• Peso di importanza (importance weight): w = p/q, corregge l'aver campionato dalla distribuzione sbagliata.
• Dimensione campionaria efficace (effective sample size): ESS = (Σwᵢ)²/Σwᵢ², quanti campioni di peso uguale valgono davvero i tuoi.
• Variabili antitetiche (antithetic variates): coppie (U, 1−U), correlazione negativa e varianza minore su funzioni monotone.
• Variata di controllo (control variate): funzione correlata di media nota, sottratta con c* = Cov(f,g)/Var(g); riduce la varianza del fattore 1 − ρ².
• Riparametrizzazione (reparametrisation trick): riscrivere x ∼ p_θ come x = g_θ(ε) con ε da una distribuzione fissa, per derivare dentro l'aspettativa.

Il glossario ha un ordine nascosto: stimatore, errore, tre modi di campionare, riduzione della varianza, riparametrizzazione. Le voci che si confondono più spesso sono le prime due: uno stimatore non distorto non è uno stimatore preciso, e l'errore standard non è la deviazione standard dei dati.

Mini-progetto

Un unico script, un'unica domanda: quanto vale il tempo di risposta al 99° percentile di un servizio composto da tre chiamate in serie?

1. Simula tre latenze indipendenti, ciascuna generata per trasformata inversa da un'esponenziale con λ diverso (0,2, 0,5 e 1,0 ms⁻¹), e somma. Stima media e 99° percentile, dichiarando l'errore standard di entrambi (per il percentile serve il bootstrap: non c'è formula chiusa);
2. decidi n partendo dalla precisione: vuoi il 99° percentile a meno di 0,1 ms. Stima σ con una corsa pilota da 10 000, calcola l'n necessario, verifica a posteriori di averlo centrato;
3. aggiungi un timeout T oltre il quale la chiamata fallisce, e stima la probabilità di fallimento per T grande, prima con Monte Carlo ingenuo e poi per importanza (suggerimento: proponi le esponenziali con un λ più piccolo e pesa col rapporto delle densità). Riporta i due errori relativi e gli episodi utili di ciascuno.

Le tre richieste costruiscono un piccolo progetto di misura completo: stima con errore, dimensionamento dei campioni a partire dalla precisione, evento raro con campionamento per importanza. Se non riesci a dichiarare l'errore standard del 99° percentile, hai trovato il punto in cui serve il bootstrap.

Consegna: lo script deve stampare, per ogni punto, stima, errore standard e numero di campioni, più una riga che dice se la precisione richiesta al punto 2 è stata raggiunta.

Challenge

• Easy. Carta e penna: hai uno stimatore con errore standard 0,08 su n = 5000 campioni. Quanti ne servono per portarlo a 0,04, e perché? Soluzione: l'errore scende come 1/√n, quindi dimezzarlo richiede di quadruplicare n: servono 20 000 campioni. In generale n_nuovo = n · (e_vecchio/e_nuovo)², e l'esponente 2 è l'unica cosa da ricordare di tutta la lezione.
• Medium. Carta e penna: deriva la trasformata inversa per l'esponenziale di parametro λ partendo da F(x) = 1 − e^{−λx}, e per una discreta a tre esiti con probabilità 0,5, 0,3 e 0,2. Per la discreta, scrivi esplicitamente i tre intervalli di u e di' che cosa restituisce la formula per u = 0,5 esatto, motivando la scelta con la definizione di inversa generalizzata.
• Hard. Codice: scegli una densità con massimo noto (una Beta diversa, o una mistura di gaussiane troncata) e implementa il rifiuto con proposta uniforme. Misura il tasso su 10⁵ tentativi contro 1/M; poi rifai tutto con M quadruplicato e verifica che i campioni abbiano la stessa distribuzione (tre quantili) ma siano costati quattro volte tanto.
• Expert. Codice: studia come il divario fra score function e riparametrizzazione dipende da μ. Le varianze teoriche sono μ⁴ + 14μ² + 15 e 4σ² costante: traccia il rapporto per μ da 0 a 5, verificalo numericamente, poi aggiungi il baseline ottimo e misura quanta parte del divario recupera al crescere di μ.
• Enterprise. Prendi una quantità del tuo lavoro oggi riportata senza incertezza — un tempo di risposta medio, un tasso di conversione, il costo atteso di un batch — e stimala per simulazione. Dichiara l'errore standard, giustifica il numero di campioni con la precisione richiesta invece che col tempo disponibile, e scrivi una riga su quale distorsione potrebbe esserci nel modo in cui i dati sono stati raccolti.

I cinque livelli salgono dal conto a mano alla stima su dati tuoi: i primi due consolidano l'errore standard e la trasformata inversa, il terzo e il quarto richiedono codice e confronti con il valore teorico, l'ultimo porta la dichiarazione dell'incertezza su una quantità che oggi riporti senza errore. Fai prima quello in cui ti blocchi.

Gli esercizi di codice di questa lezione sono in un solo file, esercizi.py (circa 15 KB). Tre esercizi di codice su domini nuovi rispetto all'articolo: la trasformata inversa per una Pareto a code pesanti, dove l'errore standard esiste come numero ma non significa niente; il collezionista di figurine, con il valore esatto n·H_n da ritrovare; una variata di controllo con il coefficiente ottimo calcolato a mano. Scheletro con le parti da completare e, in fondo, le soluzioni commentate. Richiede Python 3.12+ e NumPy 2.3. Eseguilo con python3 esercizi.py: gira anche con le parti da completare lasciate vuote e ti dice quali mancano.

Recap

Mappa mentale. Î_n = (1/n)Σ f(xᵢ) stima E_p[f], è non distorta per linearità, e ha errore standard σ/√n — indipendente dalla dimensione, che è tutto il punto. Trasformata inversa: F⁻¹(U), esatta ma serve l'inversa. Rifiuto: p ≤ Mq, accettazione 1/M, muore oltre le cinque dimensioni. Importanza: E_p[f] = E_q[f·p/q], varianza ∫f²p²/q − I², esplode se q ha code leggere. Antitetiche e variate di controllo abbassano σ invece di alzare n. Riparametrizzazione: x = μ + σε sposta la casualità fuori dal percorso del gradiente.

Cheat sheet.

3 concetti: lo stimatore non distorto con errore σ/√n indipendente dalla dimensione; il campionamento per importanza come riscrittura esatta che sposta i campioni dove c'è informazione; la riparametrizzazione come spostamento della casualità fuori dal percorso del gradiente. 3 errori: confondere la dispersione dei dati con quella dello stimatore; fidarsi di un intervallo di confidenza su uno stimatore distorto; usare una proposta con code più leggere del bersaglio. 3 best practice: dichiarare seme e versioni; dimensionare n a partire dall'errore voluto; validare su un caso con risposta esatta prima di usarne uno senza. 1 esercizio: il mini-progetto di questa lezione.

Prossima lezione: le catene di Markov, cioè i processi che dimenticano il passato — e la distribuzione che alla fine non cambia più.

Domande per andare avanti

• Se il rifiuto muore in alta dimensione e l'importanza richiede una buona proposta, come si campiona da una distribuzione complicata in cento dimensioni? (Catene di Markov Monte Carlo: lezione 15 di questo modulo per la parte di teoria, e Metropolis-Hastings come applicazione.)
• La distribuzione stazionaria usata qui come riferimento esatto da dove viene, esattamente? (È l'autovettore sinistro con autovalore 1 della matrice di transizione: la lezione su autovalori e SVD dà lo strumento, la prossima lo applica.)
• Se la media stimata ha un errore standard, quanto grande deve essere la differenza fra due stime perché si possa dire che sono diverse? (È il problema del test di ipotesi: la lezione sulla statistica inquadra la domanda, qui ne abbiamo visto solo il lato della precisione.)

Le tre domande collegano la lezione alle altre: la prima porta alle catene di Markov Monte Carlo come rimedio al rifiuto in alta dimensione, la seconda alla stazionaria come autovettore sinistro della matrice di transizione, la terza al test di ipotesi, cioè al confronto fra due stime che hanno entrambe un errore standard. Non vanno risolte subito.

Nella stessa serie trovi anche: 28 - Statistica bayesiana: prior, posterior, coniugati e perché la regolarizzazione è una credenza (la lezione precedente), 30 - Catene di Markov: la proprietà che dimentica il passato, e la distribuzione che non cambia più (la lezione successiva).