Esercizio - Stima Monte Carlo di pi greco e di una probabilità rara con campionamento per importanza

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 →; limiti di coda in 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 →.

Testo.

1. Si stima $\pi$ lanciando $N$ punti uniformi in $[0,1]^2$ e contando quelli con $u^2+v^2<1$. Quanti campioni servono perché $|\hat\pi-\pi|\le0{,}005$ con probabilità almeno $0{,}95$? Calcolare il numero con Hoeffding, con Chebyshev e con il teorema del limite centrale.
2. Con $N=10^6$ si ottengono $784\,674$ punti dentro. Stima di $\pi$, deviazione standard stimata e intervallo di confidenza al $95\%$.
3. Si vuole stimare $p=P(Z>4)$, $Z\sim\mathcal N(0,1)$ ($p=3{,}167\cdot10^{-5}$). Quanti campioni servono con il metodo ingenuo per un errore relativo del $10\%$?
4. Campionamento per importanza con la proposta $q=\mathcal N(\mu,1)$: scrivere il peso, la varianza relativa per campione in funzione di $\mu$ e scegliere $\mu$. Stima con $N=10^5$.

Soluzione

1. Sia $p=\tfrac\pi4=0{,}7854$ e $\hat\pi=4\hat p$. L'errore $|\hat\pi-\pi|\le0{,}005$ equivale a $|\hat p-p|\le\varepsilon=0{,}00125$ con $\delta=0{,}05$.

• Hoeffding: $N\ge\dfrac{\ln(2/\delta)}{2\varepsilon^2}=\dfrac{\ln40}{2\cdot1{,}5625\cdot10^{-6}}=\dfrac{3{,}6889}{3{,}125\cdot10^{-6}}=\boxed{1\,180\,441}$.
• Chebyshev: $N\ge\dfrac{p(1-p)}{\delta\varepsilon^2}=\dfrac{0{,}1686}{0{,}05\cdot1{,}5625\cdot10^{-6}}=\boxed{2\,157\,413}$.
• TLC: $\text{sd}(\hat\pi)=4\sqrt{p(1-p)/N}\le0{,}005/1{,}96$, quindi $N\ge\dfrac{(1{,}96\cdot4)^2\,p(1-p)}{0{,}005^2}=\dfrac{61{,}47\cdot0{,}1686}{2{,}5\cdot10^{-5}}=\boxed{414\,396}$.

Il TLC dà il numero più piccolo ma è un'approssimazione, non una garanzia: in $20\,000$ ripetizioni simulate con $N=414\,396$ la frazione di stime entro $0{,}005$ è stata $0{,}948$, appena sotto il $0{,}95$ richiesto; con $N=1\,180\,441$ è stata $0{,}999$ (Hoeffding è prudente).

2. $\hat p=0{,}784674$, quindi $\hat\pi=4\hat p=3{,}138696$. La deviazione standard stimata è $4\sqrt{\hat p(1-\hat p)/N}=4\sqrt{0{,}784674\cdot0{,}215326/10^6}=4\cdot4{,}110\cdot10^{-4}=0{,}001644$. Intervallo al $95\%$: $3{,}138696\pm1{,}96\cdot0{,}001644=[3{,}1355;\ 3{,}1419]$. Contiene $\pi=3{,}14159$ (appena, vicino al bordo: succede nel $5\%$ delle ripetizioni che l'intervallo non contenga $\pi$).

3. Il metodo ingenuo ha successo con probabilità $p=3{,}167\cdot10^{-5}$ e l'errore relativo della media è $\sqrt{(1-p)/(Np)}$. Per $0{,}10$: $N=\dfrac{1-p}{p\cdot0{,}01}\approx\boxed{3\,157\,000}$ campioni. Con $N=10^5$ si vedono in media solo $3{,}2$ successi: la stima è inservibile (nella prova del codice sotto se ne sono visti $6$, stima $6\cdot10^{-5}$, quasi il doppio del vero).

4. Il campionamento per importanza stima $p=E_\varphi[\mathbf 1[Z>4]]$ campionando $Y\sim\mathcal N(\mu,1)$ e pesando con $w(y)=\dfrac{\varphi(y)}{\varphi(y-\mu)}=e^{-y^2/2+(y-\mu)^2/2}=e^{-\mu y+\mu^2/2}$: $$\hat p=\frac1N\sum_j\mathbf 1[Y_j>4]\,e^{-\mu Y_j+\mu^2/2}.$$ È non distorto perché $E_q[\mathbf 1[Y>4]w(Y)]=\int_4^\infty\varphi=p$. Il secondo momento per campione è $$M_2(\mu)=\int_4^\infty\frac{\varphi(y)^2}{\varphi(y-\mu)}dy=\int_4^\infty\frac{e^{-y^2/2-\mu y+\mu^2/2}}{\sqrt{2\pi}}dy,\qquad\text{varianza relativa}=\frac{M_2}{p^2}-1.$$ Calcolando numericamente:

Conviene $\mu\approx4$, cioè centrare la proposta sulla soglia (con $\mu$ troppo grande si campiona oltre la soglia con pesi troppo variabili). Con $\mu=4$ il peso è $e^{-4y+8}$. L'errore relativo con $N$ campioni è $\sqrt{4{,}51/N}$: per $N=10^5$, $0{,}67\%$ (contro $3\,157\,000$ campioni del metodo ingenuo per il $10\%$: guadagno di circa $7000$ volte in varianza). La prova del codice sotto, con $N=10^5$, ha dato $\hat p=3{,}136\cdot10^{-5}$ con intervallo al $95\%$ $[3{,}095;\ 3{,}177]\cdot10^{-5}$, che contiene il valore vero $3{,}167\cdot10^{-5}$ (varianza relativa stimata $4{,}43$, contro $4{,}51$ esatta).

Codice di verifica

import numpy as np
from math import erfc, log, pi, sqrt
rng = np.random.default_rng(11)
p = pi / 4

# 1. numero di campioni per |pi_hat - pi| <= 0.005 al 95%
eps, delta = 0.005 / 4, 0.05
print("Hoeffding:", log(2 / delta) / (2 * eps**2), " Chebyshev:", p * (1 - p) / (delta * eps**2),
      " TLC:", (1.96 * 4)**2 * p * (1 - p) / 0.005**2)

# 2. intervallo di confidenza
N, k = 10**6, 784674
ph = k / N; sd = 4 * sqrt(ph * (1 - ph) / N)
print(4 * ph, sd, 4 * ph - 1.96 * sd, 4 * ph + 1.96 * sd)

# 3-4. evento raro P(Z > 4)
vero = 0.5 * erfc(4 / sqrt(2)); N = 10**5
y = rng.standard_normal(N) + 4                     # proposta N(4, 1)
w = np.exp(-4 * y + 8) * (y > 4)
print("vero", vero, " IS:", w.mean(), "+-", 1.96 * w.std() / sqrt(N), " varianza relativa:", w.var() / vero**2)
print("ingenuo: successi su 1e5 =", (rng.standard_normal(N) > 4).sum())