""" Corso Machine Learning Engineer - Modulo 02, Lezione 14 "Campionamento e metodi Monte Carlo: stimare quello che non si riesce a calcolare" OBIETTIVO DI QUESTO FILE ------------------------ Tre esercizi su domini NUOVI rispetto all'articolo, scelti perché ognuno ha una risposta esatta con cui confrontare la stima: 1. trasformata inversa per una Pareto, cioè una distribuzione a CODE PESANTI. La formula si deriva in due righe; la sorpresa arriva dopo, quando si prova a mettere una barra d'errore su una media che, per alpha <= 2, ha varianza INFINITA: l'errore standard sigma/sqrt(n) esiste come numero stampato e non significa niente; 2. il collezionista di figurine: quante bustine servono per completare un album di 20 figurine? Il valore atteso è esattamente n*H_n (somma di variabili geometriche indipendenti), e la simulazione lo ritrova. Qui si esercita la parte che i corsi saltano: dichiarare l'errore standard DELLO STIMATORE e dedurne quante ripetizioni servono; 3. una variata di controllo con il coefficiente ottimo calcolato a mano: si stima E[1/(1+U)] = ln 2 e si misura di quanto scende la varianza. Il fattore di riduzione teorico è 1/(1 - rho^2) ~ 31.8. COME SI USA ----------- python3 esercizi.py Il file gira anche così com'è: stampa quali esercizi mancano ancora. AMBIENTE DICHIARATO ------------------- Python 3.12+ e NumPy 2.3. Nessun altro pacchetto. python3 -m venv .venv && source .venv/bin/activate pip install "numpy>=2.0" Il seme è fissato: con lo stesso seme i numeri sono identici a ogni esecuzione. Se cambi il seme i numeri cambiano, e DEVONO cambiare entro le barre d'errore che il file stesso stampa: è il controllo più utile dei tre. """ import numpy as np SEME = 42 # ---------------------------------------------------------------------- # ESERCIZIO 1 - Trasformata inversa per una Pareto, e la barra d'errore # che non vale niente # # Densita' di Pareto con soglia x_m > 0 ed esponente alpha > 0: # # p(x) = alpha * x_m^alpha / x^(alpha+1) per x >= x_m # # Funzione di ripartizione: # # F(x) = 1 - (x_m / x)^alpha per x >= x_m # # Invertila. Posto u = F(x), si ricava x in funzione di u: è la formula che # devi scrivere nel TODO. Tre controlli gratuiti sul risultato: # # F^-1(0) = x_m (il minimo del supporto) # F^-1(0.5) = x_m * 2^(1/alpha) (la mediana) # media = alpha * x_m / (alpha - 1) SOLO se alpha > 1 # varianza = x_m^2 * alpha / ((alpha-1)^2 * (alpha-2)) SOLO se alpha > 2 # # Il file prova DUE esponenti: alpha = 2.5 (varianza finita) e alpha = 1.5 # (media finita, varianza infinita). Guarda cosa fa l'errore standard nel # secondo caso, e nota che nessuna riga di codice si lamenta. # ---------------------------------------------------------------------- SOGLIA = 1.0 def pareto_per_trasformata_inversa(n, alpha, soglia, rng): """Campioni Pareto(alpha, soglia) ottenuti invertendo la F. Ritorna ------- np.ndarray (n,) con tutti i valori >= soglia, oppure None se non svolto. Usa `rng.random(n)`, che restituisce valori in [0, 1): lo zero può uscire, l'uno no. Scrivi la formula in modo che il valore u = 0 dia esattamente la soglia e non una divisione per zero: se ti viene `1 - u` al denominatore stai bene, se ti viene `u` hai invertito il verso e otterrai inf con probabilità circa 1e-16 per campione. Non è mai zero, ed è il tipo di bug che si manifesta in produzione dopo un mese. """ # TODO (esercizio 1): sostituisci `return None` con la formula invertita. return None # ---------------------------------------------------------------------- # ESERCIZIO 2 - Il collezionista di figurine # # Un album ha n figurine distinte. Ogni bustina contiene una figurina # scelta uniformemente a caso, con reinserimento. Quante bustine servono, # in media, per completare l'album? # # Il tempo totale T si spezza nella somma dei tempi per passare da k # figurine distinte possedute a k+1: # # T = T_0 + T_1 + ... + T_(n-1) # # Quando se ne possiedono k, la probabilità che la prossima bustina porti # una figurina nuova è (n - k)/n, quindi T_k è geometrica di media # n/(n - k). Per linearità del valore atteso: # # E[T] = somma_(k=0)^(n-1) n/(n-k) = n * (1 + 1/2 + ... + 1/n) = n * H_n # # Per n = 20: E[T] = 20 * H_20 = 71.95479... (H_20 = 3.5977396571...) # # Due cose da fare: # (a) scrivi il simulatore di UN album (quante bustine fino al completamento); # (b) ragiona sull'errore: la deviazione standard di T è circa 1.2 * n # (per n = 20 vale circa 23), quindi l'errore standard della media su # R ripetizioni è circa 23/sqrt(R). Per avere la stima a meno di # mezza bustina servono R tale che 23/sqrt(R) < 0.5, cioè R > 2100. # Il file lo verifica al posto tuo. # ---------------------------------------------------------------------- FIGURINE = 20 def bustine_per_completare(n_figurine, rng): """Numero di bustine per vedere tutte le `n_figurine`, una simulazione. Ritorna ------- int, oppure None se l'esercizio non è ancora svolto. Suggerimento di implementazione: tieni un array booleano di figurine possedute e un contatore. La versione pigra `while len(set(estratte)) < n` ricostruisce l'insieme a ogni giro ed è quadratica: su 20 figurine non te ne accorgi, su 20.000 si. Il costo di un simulatore conta, perché il numero di ripetizioni lo decide l'errore che vuoi, non la pazienza che hai. """ # TODO (esercizio 2): sostituisci `return None` con la simulazione. return None # ---------------------------------------------------------------------- # ESERCIZIO 3 - Variata di controllo con coefficiente ottimo # # Si vuole stimare # # I = E[ f(U) ] con f(u) = 1/(1+u) e U ~ Uniforme(0, 1) # # che vale esattamente ln 2 = 0.6931472 (la primitiva di 1/(1+u) è # log(1+u)). Serve come banco di prova proprio perché la risposta è nota. # # La variata di controllo è g(u) = u, con media nota E[g] = 1/2. Lo # stimatore corretto è # # f(u_i) - c * ( g(u_i) - 1/2 ) # # che è NON DISTORTO per QUALUNQUE c, perché il termine aggiunto ha media # zero. Il valore di c cambia solo la varianza, che è # # Var(f) - 2c*Cov(f,g) + c^2*Var(g) # # minima in c* = Cov(f, g) / Var(g). Sostituendo, la varianza residua è # Var(f) * (1 - rho^2) con rho la correlazione fra f e g. # # Conti esatti per questo caso (fatti a mano, verificabili integrando): # Var(f) = 1/2 - (ln 2)^2 = 0.0195470 # Var(g) = 1/12 = 0.0833333 # Cov(f,g) = (1 - ln 2) - ln2 / 2 = -0.0397208 # c* = 12 * Cov(f,g) = -0.4766496 # rho^2 = 0.9685727 # riduzione della varianza 1/(1 - rho^2) = 31.82 volte # # Il segno di c* è negativo, e ha senso: f è DECRESCENTE in u mentre g è # crescente, quindi sono anticorrelate. # # Nota pratica: c* si può anche STIMARE dai dati (`np.cov`), ma stimarlo # sugli stessi campioni che usi per la media introduce una piccola # distorsione. Su n grande è trascurabile; su n piccolo si stima c su un # campione pilota separato. # ---------------------------------------------------------------------- C_OTTIMO = -0.4766496 def stimatore_con_variata_di_controllo(n, c, rng): """Campione dello stimatore f(U) - c*(U - 1/2), lungo n. Ritorna ------- np.ndarray (n,), oppure None se l'esercizio non è ancora svolto. Il valore atteso NON dipende da c; la varianza si. Il file stampa la varianza per tre valori di c (0, c*/2, c*) così si vede la parabola. """ # TODO (esercizio 3): sostituisci `return None` con lo stimatore. return None # ---------------------------------------------------------------------- # Runner # ---------------------------------------------------------------------- def _sezione(titolo): print() print(titolo) print("-" * len(titolo)) def _da_fare(numero): print(f" [ ] Esercizio {numero} non ancora svolto: completa il TODO qui sopra.") def esegui_esercizio_1(): _sezione("Esercizio 1 - Pareto per trasformata inversa") rng = np.random.default_rng(SEME) for alpha in (2.5, 1.5): campioni = pareto_per_trasformata_inversa(200_000, alpha, SOGLIA, rng) if campioni is None: _da_fare(1) return False media_attesa = alpha * SOGLIA / (alpha - 1.0) errore_standard = campioni.std(ddof=1) / np.sqrt(campioni.size) print(f" alpha = {alpha}") print(f" minimo {campioni.min():.6f} (deve essere >= {SOGLIA})") print(f" mediana {np.median(campioni):.6f}" f" esatta {SOGLIA * 2.0 ** (1.0 / alpha):.6f}") print(f" media {campioni.mean():.6f} esatta {media_attesa:.6f}" f" errore standard dichiarato {errore_standard:.6f}") if alpha > 2.0: varianza = SOGLIA ** 2 * alpha / ((alpha - 1.0) ** 2 * (alpha - 2.0)) print(f" varianza campionaria {campioni.var(ddof=1):.6f}" f" esatta {varianza:.6f}") else: print(f" varianza campionaria {campioni.var(ddof=1):.6f}" f" esatta: NON ESISTE (è infinita)") print(" -> l'errore standard stampato sopra è un numero senza") print(" significato: non converge, cresce con n a scatti.") return True def esegui_esercizio_2(): _sezione("Esercizio 2 - il collezionista di figurine") rng = np.random.default_rng(SEME) atteso = FIGURINE * sum(1.0 / k for k in range(1, FIGURINE + 1)) print(f" valore esatto n*H_n = {atteso:.6f} bustine") print(" ripetizioni media stimata errore standard IC 95%") for ripetizioni in (50, 200, 2200, 20000): prove = np.empty(ripetizioni, dtype=float) for r in range(ripetizioni): esito = bustine_per_completare(FIGURINE, rng) if esito is None: _da_fare(2) return False prove[r] = esito stima = prove.mean() errore_standard = prove.std(ddof=1) / np.sqrt(ripetizioni) print(f" {ripetizioni:>11} {stima:>15.4f} {errore_standard:>17.4f}" f" [{stima - 1.96 * errore_standard:.4f}," f" {stima + 1.96 * errore_standard:.4f}]") print(" L'intervallo si restringe come 1/sqrt(R): quadruplicare le") print(" ripetizioni lo dimezza. Per guadagnare una cifra decimale") print(" servono cento volte le ripetizioni, non dieci.") return True def esegui_esercizio_3(): _sezione("Esercizio 3 - variata di controllo") rng = np.random.default_rng(SEME) esatto = float(np.log(2.0)) print(f" valore esatto ln 2 = {esatto:.7f}") print(" c stima varianza riduzione errore standard") base = None for c in (0.0, C_OTTIMO / 2.0, C_OTTIMO): campione = stimatore_con_variata_di_controllo(500_000, c, rng) if campione is None: _da_fare(3) return False varianza = campione.var(ddof=1) base = varianza if base is None else base print(f" {c:>7.4f} {campione.mean():>10.6f} {varianza:>10.6f}" f" {base / varianza:>11.2f}x" f" {campione.std(ddof=1) / np.sqrt(campione.size):>17.2e}") print(" Atteso dalla teoria: varianza 0.019547 a c = 0, 0.000614 a c = c*,") print(" cioè una riduzione di 31.82 volte A PARITÀ di campioni. Non si") print(" è pagato nulla: g(u) = u era già estratto.") return True def main(): print("Corso MLE - Modulo 02, Lezione 14 - esercizi di codice") print(f"NumPy {np.__version__}, seme {SEME}") fatti = sum(bool(f()) for f in (esegui_esercizio_1, esegui_esercizio_2, esegui_esercizio_3)) print() print(f"Esercizi completati: {fatti} su 3.") if fatti < 3: print("Completa i TODO e rilancia: python3 esercizi.py") if __name__ == "__main__": main() # ====================================================================== # SOLUZIONI DI RIFERIMENTO # ====================================================================== # # --- Esercizio 1 ------------------------------------------------------- # def pareto_per_trasformata_inversa(n, alpha, soglia, rng): # u = rng.random(n) # return soglia / (1.0 - u) ** (1.0 / alpha) # # Derivazione: u = 1 - (x_m/x)^alpha => (x_m/x)^alpha = 1 - u # => x_m/x = (1-u)^(1/alpha) => x = x_m * (1-u)^(-1/alpha) # # Atteso, alpha = 2.5: mediana 1.319508, media 1.666667, varianza 5.333333. # Atteso, alpha = 1.5: mediana 1.587401, media 3.000000, varianza: infinita. # Nel secondo caso la media campionaria converge ancora (vale la legge dei # grandi numeri, perché E[X] esiste), ma NON vale il teorema del limite # centrale nella forma usuale: l'intervallo sigma/sqrt(n) non ha copertura # 95%. Il sintomo in produzione: la "barra d'errore" si restringe per un po' # e poi un singolo campione la fa saltare. Rimedio: usare quantili e # bootstrap percentile, non la formula gaussiana; o stimare la mediana. # # --- Esercizio 2 ------------------------------------------------------- # def bustine_per_completare(n_figurine, rng): # possedute = np.zeros(n_figurine, dtype=bool) # distinte, bustine = 0, 0 # while distinte < n_figurine: # figurina = rng.integers(n_figurine) # bustine += 1 # if not possedute[figurina]: # possedute[figurina] = True # distinte += 1 # return bustine # # Atteso: E[T] = 20*H_20 = 71.954793. La deviazione standard esatta è # sqrt(n^2 * somma_(k=1)^(n-1) 1/k^2 * ... ) ~ 23.2 per n = 20; quindi con # R = 2200 ripetizioni l'errore standard è circa 0.49 bustine, con # R = 20000 circa 0.16. Confronta: la media di T è circa 72, cioè per # stimarla all'1% (0.72 bustine) bastano circa 1000 album; per stimarla # allo 0.1% (0.072) ne servono circa 100.000. Cento volte i campioni per una # cifra decimale: è sempre la stessa legge 1/sqrt(n). # # --- Esercizio 3 ------------------------------------------------------- # def stimatore_con_variata_di_controllo(n, c, rng): # u = rng.random(n) # return 1.0 / (1.0 + u) - c * (u - 0.5) # # Conti esatti: # E[f] = ln 2 = 0.6931472 # E[f^2] = integrale (1+u)^-2 du = 1/2 # Var(f) = 1/2 - (ln 2)^2 = 0.0195470 # E[f*g] = integrale u/(1+u) du = 1 - ln 2 = 0.3068528 # Cov(f,g) = (1 - ln2) - 0.5*ln2 = -0.0397208 # c* = Cov/Var(g) = 12*Cov = -0.4766496 # rho^2 = Cov^2/(Var(f)*Var(g)) = 0.9685727 # Var residua = Var(f)*(1-rho^2) = 0.0006143 # riduzione = 1/(1-rho^2) = 31.82 volte # # A c = c*/2 la riduzione NON è la metà: la varianza è una parabola in c, # quindi a metà strada dal vertice si recupera il 75% del guadagno. È la # ragione pratica per cui una stima grossolana di c* basta quasi sempre.