Esercizio - Generazione di variabili aleatorie con trasformazione inversa e accettazione-rifiuto

Teoria: Metodi Monte Carlo - stima, generazione di variabili aleatorie e MCMCMonte Carlo stima una grandezza come media di campioni casuali: errore tipico σ/√N, numero di campioni N ≥ ln(2/δ)/(2ε²) per Hoeffding (probabilità), N ≥ 3 ln(2/δ)/(pε²) per Chernoff (errore relativo). Stima di π con p=π/4; integrali come media di g(U); eventi rari con campionamento per importanza (pesi φ/q); conteggio DNF con Karp-Luby in tempo polinomiale perché P(Y=1) ≥ 1/m. Generazione: trasformazione inversa X=F^{-1}(U) e accettazione-rifiuto (f ≤ M g, efficienza 1/M). MCMC: Metropolis-Hastings costruisce una catena con stazionaria π (nota a meno di costante) imponendo il bilancio dettagliato con probabilità di accettazione min(1, π(y)Q(y,x)/(π(x)Q(x,y))).Metodi Monte Carlo - stima, generazione di variabili aleatorie e MCMC →; 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 →, 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 →.

Testo. Si dispone di $U\sim\text{U}(0,1)$.

1. Densità di Rayleigh $f(x)=xe^{-x^2/2}$ per $x\ge0$: generarla con la trasformazione inversa; calcolare $P(X>2)$, media e varianza e confrontarle con una simulazione.
2. Variabile discreta con $P(X=k)=k/10$ per $k=1,2,3,4$: algoritmo di generazione e controllo.
3. $\text{Beta}(3,2)$, $f(x)=12x^2(1-x)$ su $[0,1]$: accettazione-rifiuto con $g$ uniforme. Quale $M$ minimo? Probabilità di accettazione e numero medio di tentativi. Cosa cambia con $M=2$?
4. Dimostrare che l'algoritmo del punto 3 restituisce la densità $f$.

Soluzione

1. La funzione di distribuzione è $F(x)=\int_0^xte^{-t^2/2}dt=1-e^{-x^2/2}$. Invertendo $u=1-e^{-x^2/2}$: $x=\sqrt{-2\ln(1-u)}$; poiché $1-U\sim U$, si può usare $$\boxed{X=\sqrt{-2\ln U}}.$$ Coda: $P(X>2)=e^{-2}=0{,}1353$. Media: $E[X]=\int_0^\infty x^2e^{-x^2/2}dx=\sqrt{\pi/2}=1{,}2533$ (metà dell'integrale gaussiano $\int_{-\infty}^\infty x^2e^{-x^2/2}dx=\sqrt{2\pi}$). Secondo momento: $E[X^2]=\int_0^\infty x^3e^{-x^2/2}dx=2$ (con $s=x^2/2$: $\int_0^\infty2s\,e^{-s}ds=2$), quindi $\text{Var}=2-\tfrac\pi2=0{,}4292$. Simulazione con $10^6$ campioni (codice sotto): media $1{,}2534$, varianza $0{,}4293$, $P(X>2)=0{,}1354$.

2. Le somme cumulative sono $F(1)=0{,}1$, $F(2)=0{,}3$, $F(3)=0{,}6$, $F(4)=1$. Si restituisce il più piccolo $k$ con $U\le F(k)$: $U\in[0;0{,}1]\to1$, $(0{,}1;0{,}3]\to2$, $(0{,}3;0{,}6]\to3$, $(0{,}6;1]\to4$; ogni intervallo ha lunghezza uguale alla probabilità cercata. Frequenze su $10^6$ prove: $0{,}1003;\ 0{,}1996;\ 0{,}3000;\ 0{,}4001$.

3. Si cerca $M$ minimo con $12x^2(1-x)\le M\cdot1$ su $[0,1]$. La derivata è $12(2x-3x^2)=12x(2-3x)$, nulla in $x=\tfrac23$ dove $f$ è massima: $f(\tfrac23)=12\cdot\tfrac49\cdot\tfrac13=\tfrac{16}9=1{,}778$. Quindi $\boxed{M=\tfrac{16}9}$. L'algoritmo: estrai $Y,U$ uniformi; accetta se $U\le\dfrac{12Y^2(1-Y)}{16/9}=\tfrac{27}4Y^2(1-Y)$ (il rapporto vale $1$ in $\tfrac23$). La probabilità di accettare a ogni tentativo è $\tfrac1M=\tfrac9{16}=0{,}5625$, i tentativi sono geometrici di media $M=1{,}778$ (simulazione: $1{,}774$). Con $M=2$ l'algoritmo resta corretto (basta $f\le Mg$) ma accetta con probabilità $\tfrac12$ e servono in media $2$ tentativi (simulazione: $1{,}999$). I campioni di $\text{Beta}(3,2)$ hanno media $\tfrac35=0{,}6$ e varianza $\tfrac{3\cdot2}{25\cdot6}=0{,}04$ (simulazione: $0{,}5995$ e $0{,}0401$ con $M=\tfrac{16}9$, $0{,}6003$ e $0{,}0400$ con $M=2$).

4. Sia $Y\sim\text{U}(0,1)$ e $U\sim\text{U}(0,1)$ indipendenti. $P(Y\in dy,\ \text{accetta})=1\cdot\dfrac{f(y)}{M}dy$. Integrando su $[0,1]$: $P(\text{accetta})=\tfrac1M\int_0^1f=\tfrac1M$ (con $M=\tfrac{16}9$: $0{,}5625$, come sopra). Dato "accetta", la densità di $Y$ è $\dfrac{f(y)/M}{1/M}=f(y)$. I tentativi sono indipendenti, quindi i valori restituiti hanno esattamente densità $f$. $\blacksquare$

Codice di verifica

import numpy as np
rng = np.random.default_rng(1)

# 1. Rayleigh con la trasformazione inversa
X = np.sqrt(-2 * np.log(rng.random(10**6)))
print("1.", X.mean(), X.var(), (X > 2).mean(), " attesi", np.sqrt(np.pi / 2), 2 - np.pi / 2, np.exp(-2))

# 2. discreta con P(k) = k/10
k = np.searchsorted(np.cumsum([.1, .2, .3, .4]), rng.random(10**6), side="right") + 1
print("2.", np.bincount(k)[1:] / 10**6)

# 3. Beta(3,2) per accettazione-rifiuto
def beta32(N, M):
    out, tent = [], 0
    while len(out) < N:
        y = rng.random(); tent += 1
        if rng.random() <= 12 * y * y * (1 - y) / M: out.append(y)
    return np.array(out), tent / N
for M in (16 / 9, 2):
    b, t = beta32(200_000, M)
    print("3. M =", round(M, 3), "media", b.mean(), "var", b.var(), "tentativi per campione", t)