Metodi Monte Carlo - stima, generazione di variabili aleatorie e MCMC

Un metodo Monte Carlo ottiene una quantità deterministica (un'area, una probabilità, un conteggio) come media di esperimenti casuali. Qui: stima e numero di campioni, riduzione della varianza, conteggio approssimato, generazione di variabili aleatorie e campionamento da distribuzioni difficili con le catene di Markov. Prerequisiti: Legge dei grandi numeri e metodo Monte CarloSe X₁, X₂, ... sono i.i.d. con media μ, la media campionaria X̄ₙ = (X₁ + ... + Xₙ)/n converge a μ: in probabilità (legge debole, dimostrata con Chebyshev se la varianza è finita: P(|X̄ₙ − μ| > ε) ≤ σ²/(nε²)) e quasi certamente (legge forte). Metodo Monte Carlo: ∫ g = E[g(U)] si stima con la media di g(U₁), ..., g(Uₙ) per uniformi indipendenti.Legge dei grandi numeri e metodo Monte Carlo → (legge dei grandi numeri e Monte Carlo per gli integrali), Teorema del limite centrale e approssimazione normaleSe X₁, X₂, ... sono i.i.d. con media μ e varianza σ² ∈ (0, ∞), la somma standardizzata (Sₙ − nμ)/(σ√n) converge in distribuzione a N(0, 1): per n grande P(Sₙ ≤ x) ≈ Φ((x − nμ)/(σ√n)). Caso binomiale (De Moivre-Laplace): Bin(n, p) ≈ N(np, np(1 − p)), con correzione di continuità ±0,5. Si dimostra con le funzioni caratteristiche e il teorema di Lévy.Teorema del limite centrale e approssimazione normale →, Funzioni di una variabile aleatoriaPer trovare la legge di Y = g(X) si calcola la FdD: F_Y(y) = P(g(X) ≤ y) = P(X ∈ g⁻¹((−∞, y])), poi si deriva; se g è strettamente monotona e derivabile vale la formula f_Y(y) = f_X(g⁻¹(y)) · |(g⁻¹)'(y)|; se g è costante a tratti Y è discreta; se g non è iniettiva (es. X²) si spezza l'evento.Funzioni di una variabile aleatoria →, Funzione di distribuzioneLa funzione di distribuzione (FdD) F_X(x) = P(X ≤ x) è definita per ogni v.a., è crescente, continua a destra, va da 0 a 1 e determina la legge; per una v.a. discreta è a gradini, con salti in corrispondenza dei valori e di altezza pari alla densità.Funzione di distribuzione →, Disuguaglianze di Markov, Chebyshev e JensenMarkov: per X ≥ 0, P(X ≥ a) ≤ E[X]/a; Chebyshev: P(|X − μ| ≥ ε) ≤ Var(X)/ε²; Jensen: per φ convessa, φ(E[X]) ≤ E[φ(X)]. Stimano probabilità e medie conoscendo solo media e varianza.Disuguaglianze di Markov, Chebyshev e Jensen → e i limiti di coda di Limiti di Chernoff e di HoeffdingMarkov applicato a e^(tX) dà code esponenziali. Per somme di Bernoulli indipendenti con media μ: P(X ≥ (1+δ)μ) ≤ [e^δ/(1+δ)^(1+δ)]^μ ≤ e^(-δ²μ/3) per 0<δ≤1, P(X ≤ (1-δ)μ) ≤ e^(-δ²μ/2), P(X ≥ R) ≤ 2^(-R) se R ≥ 6μ. Hoeffding per X_i in [a_i,b_i]: P(|S-E[S]| ≥ t) ≤ 2exp(-2t²/Σ(b_i-a_i)²). Applicazioni: n ≥ ln(2/δ)/(2ε²) campioni per stimare una probabilità, amplificazione BPP con k = O(log 1/ε) ripetizioni, bilanciamento di insiemi con scarto √(4m ln n).Limiti di Chernoff e di Hoeffding →. Il concetto di algoritmo con errore è in Algoritmi randomizzati - Monte Carlo, Las Vegas e amplificazioneUn algoritmo randomizzato usa bit casuali indipendenti; la probabilità è sulle scelte dell'algoritmo, mai sull'input (garanzia per ogni input). Las Vegas: sempre corretto, tempo aleatorio (conta il tempo atteso). Monte Carlo: tempo fisso, risposta sbagliata con probabilità limitata; errore unilaterale (RP: accetta solo se x sta nel linguaggio, accetta con prob almeno 1/2) o bilaterale (BPP: giusto con prob almeno 2/3). ZPP = RP ∩ co-RP = Las Vegas con tempo atteso polinomiale. Amplificazione: unilaterale, k ripetizioni indipendenti, errore (1-p)^k ≤ e^(-pk), serve k ≥ ln(1/ε)/p; bilaterale, voto di maggioranza su k esecuzioni, errore binomiale esatto, esponenziale in k (Chernoff, e^(-k/18) per p=2/3) mentre Chebyshev dà solo 8/k. Las Vegas -> Monte Carlo: interrompere a c volte il tempo atteso (Markov, errore ≤ 1/c). Monte Carlo -> Las Vegas: solo se la risposta si verifica in fretta. Esempio: Miller-Rabin, errore ≤ 1/4 per esecuzione.Algoritmi randomizzati - Monte Carlo, Las Vegas e amplificazione →.

Stima Monte Carlo

Se $X_1,\dots,X_N$ sono i.i.d. con media $\mu$ e varianza $\sigma^2$, la media campionaria $\hat\mu=\frac1N\sum X_i$ ha $E[\hat\mu]=\mu$ e $\text{Var}(\hat\mu)=\sigma^2/N$: l'errore tipico è $\sigma/\sqrt N$ e per ridurlo di $10$ servono $100$ volte più campioni (legge dei grandi numeri: Legge dei grandi numeri e metodo Monte CarloSe X₁, X₂, ... sono i.i.d. con media μ, la media campionaria X̄ₙ = (X₁ + ... + Xₙ)/n converge a μ: in probabilità (legge debole, dimostrata con Chebyshev se la varianza è finita: P(|X̄ₙ − μ| > ε) ≤ σ²/(nε²)) e quasi certamente (legge forte). Metodo Monte Carlo: ∫ g = E[g(U)] si stima con la media di g(U₁), ..., g(Uₙ) per uniformi indipendenti.Legge dei grandi numeri e metodo Monte Carlo →).

Esempio (stima di $\pi$, hit-or-miss). Si lanciano punti uniformi $(U,V)$ in $[0,1]^2$; il quarto di cerchio $\{u^2+v^2<1\}$ ha area $\tfrac\pi4$. Sia $I=\mathbf 1[U^2+V^2<1]\sim\text{Bernoulli}(p)$ con $p=\tfrac\pi4=0{,}7854$. Allora $\hat\pi=4\hat p=\frac4N\sum I_j$ ha media $\pi$ e deviazione standard $4\sqrt{p(1-p)/N}=1{,}64/\sqrt N$.

Con $N=10^3,10^4,10^5,10^6$ la deviazione standard teorica è $0{,}052;\ 0{,}0164;\ 0{,}0052;\ 0{,}00164$ e le stime ottenute in cinque prove per ciascun $N$ stanno in $3{,}12$–$3{,}23$; $3{,}12$–$3{,}17$; $3{,}140$–$3{,}146$; $3{,}1405$–$3{,}1430$. Tre cifre decimali costano circa un milione di campioni: il metodo è lento, ma non dipende dalla dimensione del problema.

Quanti campioni

Problema. Quanti campioni servono perché $|\hat p-p|\le\varepsilon$ con probabilità $\ge1-\delta$?

• Hoeffding (variabili in $[0,1]$): $P(|\hat p-p|\ge\varepsilon)\le2e^{-2N\varepsilon^2}\le\delta\iff N\ge\dfrac{\ln(2/\delta)}{2\varepsilon^2}$.
• Chernoff, errore relativo ($0<\varepsilon\le1$): $P(|\hat p-p|\ge\varepsilon p)\le2e^{-Np\varepsilon^2/3}\le\delta\iff N\ge\dfrac{3\ln(2/\delta)}{p\,\varepsilon^2}$. Per $p$ piccolo il numero di campioni cresce come $1/p$: servono circa $1/p$ per vedere un singolo successo.
• Chebyshev: $P(|\hat p-p|\ge\varepsilon)\le\dfrac{p(1-p)}{N\varepsilon^2}\le\delta\iff N\ge\dfrac{p(1-p)}{\delta\varepsilon^2}$ (dipende da $1/\delta$ e non da $\ln(1/\delta)$).
• Teorema del limite centrale: $N\ge\bigl(z_{1-\delta/2}\,\sigma/\varepsilon\bigr)^2$, approssimato (nessuna garanzia per $N$ finito).

Esempio. Per ottenere $|\hat\pi-\pi|\le0{,}01$ con probabilità almeno $0{,}99$ si deve avere $|\hat p-p|\le\varepsilon=0{,}0025$ e $\delta=0{,}01$:

($p(1-p)=0{,}7854\cdot0{,}2146=0{,}1686$.) Su $20\,000$ ripetizioni simulate, con $N=423\,866$ la frazione di stime entro $0{,}01$ è $0{,}9999$ (garanzia molto prudente); con $N=178\,928$ è $0{,}990$: il TLC è il più economico e qui basta appena, ma non dà una garanzia.

Integrali

Se $U\sim\text{U}(0,1)$, $\int_0^1g(x)dx=E[g(U)]$, quindi $\hat I=\frac1N\sum g(U_j)$ (Legge dei grandi numeri e metodo Monte CarloSe X₁, X₂, ... sono i.i.d. con media μ, la media campionaria X̄ₙ = (X₁ + ... + Xₙ)/n converge a μ: in probabilità (legge debole, dimostrata con Chebyshev se la varianza è finita: P(|X̄ₙ − μ| > ε) ≤ σ²/(nε²)) e quasi certamente (legge forte). Metodo Monte Carlo: ∫ g = E[g(U)] si stima con la media di g(U₁), ..., g(Uₙ) per uniformi indipendenti.Legge dei grandi numeri e metodo Monte Carlo →). Per $g(x)=e^{-x^2}$: $I=0{,}746824$, $\sigma^2=\int_0^1e^{-2x^2}dx-I^2=0{,}04040$, $\sigma=0{,}2010$. Per un errore $\pm0{,}001$ al $95\%$: $N\ge(1{,}96\cdot0{,}2010/0{,}001)^2=155\,200$. Con $10^6$ campioni si ottiene $0{,}74705$.

Riduzione della varianza con variabili antitetiche. Si usa $\frac12\bigl(g(U)+g(1-U)\bigr)$: se $g$ è monotona, $g(U)$ e $g(1-U)$ sono negativamente correlate e la varianza della media crolla. Qui la varianza per campione scende da $0{,}0404$ a $0{,}000807$ (due valutazioni di $g$ per campione: guadagno di circa $25$ volte a parità di costo).

Eventi rari e campionamento per importanza

Se $p=P(A)$ è molto piccola, il metodo ingenuo ha errore relativo $\sqrt{(1-p)/(Np)}\approx1/\sqrt{Np}$: per un errore relativo del $10\%$ servono $N\approx100/p$. Esempio: $p=P(Z>4)=3{,}167\cdot10^{-5}$ per $Z$ normale standard: servono $3\,157\,000$ campioni per il $10\%$.

Campionamento per importanza. Per stimare $E_f[h(X)]$ si campiona da un'altra densità $q$ (con $q>0$ dove $f\,h\ne0$) e si pesa: $E_f[h(X)]=E_q\!\Bigl[h(X)\dfrac{f(X)}{q(X)}\Bigr]$, stimato da $\frac1N\sum h(Y_j)\,w(Y_j)$ con $w=f/q$ e $Y_j\sim q$. Rimane non distorto; la varianza è piccola se $q$ è concentrata dove $h\,f$ è grande.

Per $P(Z>4)$ si prende $q=\mathcal N(4,1)$ (centrata sulla zona rara): $w(y)=\dfrac{\varphi(y)}{\varphi(y-4)}=e^{-4y+8}$ e la stima è la media di $\mathbf 1[Y>4]\,e^{-4Y+8}$. Con $N=10^5$ campioni il metodo ingenuo vede $3$ successi (stima $3\cdot10^{-5}$, inaffidabile); l'importance sampling dà $3{,}17\cdot10^{-5}$ con errore relativo $0{,}7\%$ (la varianza relativa per campione è $4{,}5$ contro $3{,}2\cdot10^4$: circa $7000$ volte meno).

Conteggio approssimato: formule DNF

Una formula DNF è un OR di $m$ termini, ciascuno AND di letterali. Contare le assegnazioni che la soddisfano ($\#$DNF) è difficile in modo esatto, ma si approssima con un metodo Monte Carlo ben progettato (algoritmo di Karp e Luby). Il campionamento ingenuo (assegnazione uniforme, controllare se soddisfa) fallisce quando le soluzioni sono rare: con $n=40$ variabili e $3$ termini di $20$ letterali, la frazione di soluzioni è $\le3\cdot2^{20}/2^{40}=2{,}9\cdot10^{-6}$.

Sia $S_i$ l'insieme delle assegnazioni che soddisfano il termine $i$; se il termine ha $k_i$ letterali, $|S_i|=2^{n-k_i}$ e gli $S_i$ si possono campionare uniformemente fissando i letterali. Si vuole $|U|$ con $U=\bigcup_iS_i$.

Algoritmo (Karp-Luby). Ripeti $N$ volte: scegli $i$ con probabilità $|S_i|/\sum_j|S_j|$; scegli $x\in S_i$ uniformemente; poni $Y=1$ se $i$ è il più piccolo indice tra i termini soddisfatti da $x$, altrimenti $Y=0$. Restituisci $\bigl(\sum_j|S_j|\bigr)\cdot\bar Y$.

Dimostrazione. La coppia $(i,x)$ con $x\in S_i$ ha probabilità $\frac{|S_i|}{\sum_j|S_j|}\cdot\frac1{|S_i|}=\frac1{\sum_j|S_j|}$: la scelta è uniforme sulle coppie $(i,x)$ con $x\in S_i$. Ogni $x\in U$ ha esattamente una coppia "canonica" (con $i$ il più piccolo indice che lo soddisfa), quindi $E[Y]=\dfrac{|U|}{\sum_j|S_j|}$ e la stima è non distorta. Inoltre $|U|\ge\max_j|S_j|\ge\frac1m\sum_j|S_j|$, quindi $E[Y]\ge\frac1m$: il successo non è raro. Con Chernoff relativo, $N\ge\dfrac{3m\ln(2/\delta)}{\varepsilon^2}$ campioni danno errore relativo $\le\varepsilon$ con probabilità $\ge1-\delta$: tempo polinomiale in $m$, $n$, $1/\varepsilon$, $\ln(1/\delta)$. $\blacksquare$

Esempio. $n=4$, $\varphi=(x_1\wedge x_2)\vee(\neg x_2\wedge x_3)\vee(x_3\wedge x_4)$: $|S_1|=|S_2|=|S_3|=4$, somma $12$. L'unione è $|U|=9$ (le assegnazioni $0011$ e $1011$ stanno in $S_2\cap S_3$, e $1111$ in $S_1\cap S_3$; $12-3=9$), quindi $E[Y]=\tfrac9{12}=0{,}75$. Con $10^5$ campioni si ottiene $\bar Y=0{,}7500$ e la stima $9{,}0005$. Per $m=3$, $\varepsilon=0{,}05$, $\delta=0{,}05$: $N\ge3\cdot3\ln(40)/0{,}05^2=13\,280$. L'esercizio completo è in Esercizio - Conteggio delle soluzioni di una formula DNF con l'algoritmo di Karp e Luby.

Generazione di variabili aleatorie

Si assume di avere un generatore di $U\sim\text{U}(0,1)$ (numeri pseudocasuali).

Trasformazione inversa

Teorema. Sia $F$ una funzione di distribuzione (Funzione di distribuzioneLa funzione di distribuzione (FdD) F_X(x) = P(X ≤ x) è definita per ogni v.a., è crescente, continua a destra, va da 0 a 1 e determina la legge; per una v.a. discreta è a gradini, con salti in corrispondenza dei valori e di altezza pari alla densità.Funzione di distribuzione →) e $F^{-1}(u)=\min\{x:F(x)\ge u\}$. Se $U\sim\text{U}(0,1)$ allora $X=F^{-1}(U)$ ha funzione di distribuzione $F$.

Dimostrazione. $F^{-1}(u)\le x\iff u\le F(x)$ (definizione di minimo e monotonia di $F$). Quindi $P(X\le x)=P(U\le F(x))=F(x)$. $\blacksquare$

Esempi: esponenziale $F(x)=1-e^{-\lambda x}$, $X=-\ln(1-U)/\lambda$ (o $-\ln U/\lambda$, perché $1-U$ ha la stessa legge di $U$). Densità $f(x)=3x^2$ su $[0,1]$: $F=x^3$, $X=U^{1/3}$ (media $\tfrac34$: la simulazione dà $0{,}7506$). Pareto di esponente $3$: $F=1-x^{-3}$ per $x\ge1$, $X=(1-U)^{-1/3}$ (media $1{,}5$). Discrete: si calcolano le somme cumulative $F(1)<F(2)<\dots$ e si restituisce il primo $k$ con $F(k)\ge U$; per la geometrica $X=\lceil\ln U/\ln(1-p)\rceil$ (Distribuzione geometricaGeo(p) è il numero della prova in cui arriva il primo successo in prove indipendenti: P(X = n) = (1−p)^(n−1) p per n ≥ 1, P(X > n) = (1−p)^n (lunga attesa), media 1/p, varianza (1−p)/p², ed è senza memoria.Distribuzione geometrica →).

Accettazione-rifiuto

Quando $F^{-1}$ non è calcolabile (per esempio la normale, Distribuzione gaussiana (normale)N(μ, σ²) ha densità e^(−(x−μ)²/(2σ²)) / √(2πσ²), a campana centrata in μ con larghezza σ; media μ, varianza σ²; si standardizza con Z = (X − μ)/σ ~ N(0, 1) e si calcola P(X ≤ x) = Φ((x − μ)/σ), con Φ(−z) = 1 − Φ(z); aX + b è ancora gaussiana, N(aμ + b, a²σ²).Distribuzione gaussiana (normale) →) si usa una densità $g$ facile da campionare tale che $f(x)\le M\,g(x)$ per ogni $x$ con una costante $M\ge1$.

Algoritmo. Ripeti: estrai $Y\sim g$ e $U\sim\text{U}(0,1)$ indipendenti; se $U\le\dfrac{f(Y)}{M\,g(Y)}$ restituisci $Y$, altrimenti riprova.

Teorema. Il valore restituito ha densità $f$; la probabilità di accettare a ogni tentativo è $\tfrac1M$, quindi il numero di tentativi è geometrico di media $M$.

Dimostrazione. A ogni tentativo, per $Y=y$ la probabilità di accettazione è $\frac{f(y)}{Mg(y)}$. Quindi $P(Y\in dy,\text{accetta})=g(y)\dfrac{f(y)}{Mg(y)}dy=\dfrac{f(y)}M dy$. Integrando su $y$: $P(\text{accetta})=\dfrac1M$ ($f$ ha integrale $1$). Dato "accetta", la densità di $Y$ è $\dfrac{f(y)/M}{1/M}=f(y)$. $\blacksquare$

Esempio ($\text{Beta}(2,2)$). $f(x)=6x(1-x)$ su $[0,1]$ con $g$ uniforme: il massimo di $f$ è $f(\tfrac12)=1{,}5$, quindi $M=1{,}5$: si accetta con probabilità $6y(1-y)/1{,}5$ e servono in media $1{,}5$ tentativi (simulazione: $1{,}503$; media $0{,}5010$, varianza $0{,}0500$ contro $\tfrac12$ e $\tfrac1{20}$).

Esempio (normale). Per la parte positiva $|Z|$ ($f(x)=\sqrt{2/\pi}\,e^{-x^2/2}$ per $x\ge0$) si usa $g(x)=e^{-x}$: $\dfrac fg=\sqrt{2/\pi}\,e^{x-x^2/2}$ ha massimo in $x=1$, $M=\sqrt{2/\pi}\,e^{1/2}=1{,}3155$, accettazione $\tfrac1M=0{,}760$. Il test diventa $U\le e^{-(Y-1)^2/2}$. Si assegna poi un segno $\pm$ a caso. Su $3\cdot10^5$ campioni: $1{,}315$ tentativi per campione, $P(|Z|<1)=0{,}6834$ (atteso $0{,}6827$).

Catene di Markov Monte Carlo (MCMC)

Come campionare da una distribuzione $\pi$ su un insieme enorme, nota solo a meno di una costante, $\pi(x)=w(x)/Z$ con $Z=\sum_xw(x)$ non calcolabile (configurazioni di un modello fisico, parametri di un modello bayesiano)? Si costruisce una catena di Markov (Catene di Markov - matrice di transizione, classificazione degli stati e distribuzione stazionariaUna catena di Markov a tempo discreto e stati finiti è descritta da una matrice stocastica P: P(X_{n+1}=j | X_n=i)=P_ij, la legge al tempo n è π_0 P^n. Gli stati si classificano in classi comunicanti (irriducibilità), periodo (aperiodicità) e ricorrenza/transitorietà. Una distribuzione stazionaria soddisfa πP=π; per catene finite irriducibili esiste ed è unica, vale π_i=1/h_ii (h_ii = tempo medio di primo ritorno), e se la catena è anche aperiodica P^n converge a una matrice con tutte le righe uguali a π. Le catene reversibili (π_i P_ij = π_j P_ji) hanno stazionaria immediata.Catene di Markov - matrice di transizione, classificazione degli stati e distribuzione stazionaria →) irriducibile e aperiodica con stazionaria $\pi$: per il teorema fondamentale, dopo molti passi $X_n$ ha legge vicina a $\pi$.

Algoritmo (Metropolis-Hastings). Sia $Q(x,y)$ una proposta (una matrice stocastica che si sa campionare). Da $x$: si propone $y\sim Q(x,\cdot)$ e si accetta con probabilità $$a(x,y)=\min\Bigl(1,\ \frac{\pi(y)\,Q(y,x)}{\pi(x)\,Q(x,y)}\Bigr)=\min\Bigl(1,\ \frac{w(y)\,Q(y,x)}{w(x)\,Q(x,y)}\Bigr),$$ altrimenti si resta in $x$. La costante $Z$ si semplifica. Se $Q$ è simmetrica ($Q(x,y)=Q(y,x)$) è l'algoritmo di Metropolis: $a=\min(1,w(y)/w(x))$.

Teorema. La catena $P(x,y)=Q(x,y)\,a(x,y)$ per $y\ne x$ (e $P(x,x)=1-\sum_{y\ne x}P(x,y)$) soddisfa il bilancio dettagliato $\pi(x)P(x,y)=\pi(y)P(y,x)$, quindi $\pi$ è stazionaria.

Dimostrazione. Per $x\ne y$, $\pi(x)P(x,y)=\pi(x)Q(x,y)\min\bigl(1,\tfrac{\pi(y)Q(y,x)}{\pi(x)Q(x,y)}\bigr)=\min\bigl(\pi(x)Q(x,y),\ \pi(y)Q(y,x)\bigr)$, che è simmetrico in $(x,y)$. Quindi vale anche $\pi(y)P(y,x)$ (Catene di Markov - matrice di transizione, classificazione degli stati e distribuzione stazionariaUna catena di Markov a tempo discreto e stati finiti è descritta da una matrice stocastica P: P(X_{n+1}=j | X_n=i)=P_ij, la legge al tempo n è π_0 P^n. Gli stati si classificano in classi comunicanti (irriducibilità), periodo (aperiodicità) e ricorrenza/transitorietà. Una distribuzione stazionaria soddisfa πP=π; per catene finite irriducibili esiste ed è unica, vale π_i=1/h_ii (h_ii = tempo medio di primo ritorno), e se la catena è anche aperiodica P^n converge a una matrice con tutte le righe uguali a π. Le catene reversibili (π_i P_ij = π_j P_ji) hanno stazionaria immediata.Catene di Markov - matrice di transizione, classificazione degli stati e distribuzione stazionaria →). $\blacksquare$

Perché converga a $\pi$ occorre anche che la catena sia irriducibile e aperiodica (basta che $Q$ colleghi tra loro tutti gli stati con $\pi>0$ e che ci siano rifiuti, cioè cappi).

Esempio (cinque stati). Stati $1,\dots,5$, $w=(1,2,3,4,5)$, quindi $\pi=(\tfrac1{15},\tfrac2{15},\tfrac3{15},\tfrac4{15},\tfrac5{15})$ (qui $Z=15$ si conosce, ma non serve). Proposta simmetrica: $\pm1$ con probabilità $\tfrac12$ ciascuno, e se esce dai bordi il passo è rifiutato. Probabilità di accettazione: da $x$ a $x+1$ vale $1$ (si sale verso pesi più alti), da $x$ a $x-1$ vale $\tfrac{x-1}{x}$ ($w$ scende). La matrice è $$P=\begin{pmatrix}\tfrac12&\tfrac12&0&0&0\\\tfrac14&\tfrac14&\tfrac12&0&0\\0&\tfrac13&\tfrac16&\tfrac12&0\\0&0&\tfrac38&\tfrac18&\tfrac12\\0&0&0&\tfrac25&\tfrac35\end{pmatrix}.$$ Per esempio da $3$ a $2$: $\tfrac12\cdot\tfrac23=\tfrac13$; da $4$ a $3$: $\tfrac12\cdot\tfrac34=\tfrac38$; il resto va sul cappio ($P_{33}=1-\tfrac13-\tfrac12=\tfrac16$). Bilancio dettagliato, per esempio tra $2$ e $3$: $\pi_2P_{23}=\tfrac2{15}\cdot\tfrac12=\tfrac1{15}$ e $\pi_3P_{32}=\tfrac3{15}\cdot\tfrac13=\tfrac1{15}$ ✓. Gli autovalori sono $1;\ 0{,}801;\ 0{,}412;\ -0{,}062;\ -0{,}509$: la convergenza decade come $0{,}801^n$. Partendo da $1$, la distanza di variazione totale da $\pi$ vale $0{,}800$ ($n=1$), $0{,}322$ ($n=5$), $0{,}109$ ($n=10$), $0{,}0119$ ($n=20$), $1{,}4\cdot10^{-4}$ ($n=40$). La simulazione di $2\cdot10^5$ passi del codice sotto dà frequenze $(0{,}0678;\ 0{,}1349;\ 0{,}2023;\ 0{,}2667;\ 0{,}3283)$ contro $\pi=(0{,}0667;\ 0{,}1333;\ 0{,}2;\ 0{,}2667;\ 0{,}3333)$ e media $3{,}653$ contro $\tfrac{11}3=3{,}667$: i passi consecutivi sono correlati, quindi l'errore è più grande di quello di $2\cdot10^5$ campioni indipendenti.

Uso in pratica. Si scartano i primi passi (burn-in) e si media lungo la catena (teorema ergodico). I campioni consecutivi sono correlati: il numero di campioni equivalenti indipendenti è circa $N$ diviso per il tempo di autocorrelazione, legato al gap spettrale $1-|\lambda_2|$. La proposta incide molto: passi troppo piccoli esplorano lentamente, troppo grandi vengono rifiutati quasi sempre.

Codice

import random, math
rng = random.Random(1)

# 1) hit-or-miss per pi
N = 10**6
print(4 * sum(rng.random()**2 + rng.random()**2 < 1 for _ in range(N)) / N)

# 2) trasformazione inversa: f(x) = 3 x^2 su [0,1] -> U^(1/3)
xs = [rng.random() ** (1 / 3) for _ in range(10**5)]
print(sum(xs) / len(xs))                         # circa 0.75

# 3) accettazione-rifiuto per Beta(2,2), M = 1.5
def beta22():
    while True:
        y = rng.random()
        if rng.random() <= 6 * y * (1 - y) / 1.5: return y
print(sum(beta22() for _ in range(10**5)) / 10**5)   # circa 0.5

# 4) Metropolis su {1..5} con peso w(x) = x
w = lambda x: x
x, cont = 1, [0] * 6
for _ in range(200000):
    y = x + rng.choice((-1, 1))
    if 1 <= y <= 5 and rng.random() < min(1, w(y) / w(x)): x = y
    cont[x] += 1
print([round(c / 200000, 4) for c in cont[1:]])   # circa [0.0667, 0.1333, 0.2, 0.2667, 0.3333]

Domande d'esame

1. Quanti campioni servono per stimare con Hoeffding una probabilità con errore $0{,}01$ e confidenza $95\%$? $N\ge\dfrac{\ln(2/0{,}05)}{2\cdot0{,}01^2}=\dfrac{3{,}689}{0{,}0002}=18\,444$.
2. Dimostra che la trasformazione inversa genera la legge $F$ e che l'accettazione-rifiuto genera la densità $f$. $F^{-1}(u)\le x\iff u\le F(x)$; $P(\text{accetta})=1/M$, densità condizionata $f$.
3. Perché l'algoritmo di Karp-Luby è efficiente e quello ingenuo no? In Karp-Luby $E[Y]\ge1/m$ per ogni formula; nell'ingenuo la probabilità di successo può essere $2^{-\Omega(n)}$.
4. Dimostra che Metropolis-Hastings ha stazionaria $\pi$ e fai vedere perché $Z$ non serve. Bilancio dettagliato $\min(\pi(x)Q(x,y),\pi(y)Q(y,x))$; $a$ dipende solo dal rapporto $w(y)/w(x)$.