Esercizio - probabilità di una somma di valori assoluti di gaussiane indipendenti

Questa pagina non ha ancora la versione ripasso: qui sotto c'è il testo completo.

Testo (Esercitazione 5, problema 1, corso Teoria dei Segnali, UniPD). Le variabili aleatorie $x$ e $y$ sono indipendenti e gaussiane $N(0,1)$. Calcolare $P[\,|x|+|y|\le 1\,]$.

Teoria usata: Vettori gaussianiX = (X₁, ..., Xₙ) è un vettore gaussiano N(m, Σ) se ogni combinazione lineare a·X è gaussiana (equivalentemente X = m + AZ con Z gaussiane standard indipendenti); se Σ è invertibile ha densità exp(−½(x−m)ᵀΣ⁻¹(x−m)) / √((2π)ⁿ det Σ). Proprietà chiave: AX + b ~ N(Am + b, AΣAᵀ), le marginali sono gaussiane e componenti non correlate sono indipendenti.Vettori gaussiani →, 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) →, Vettori aleatori assolutamente continuiUn vettore (X, Y) è assolutamente continuo se P((X, Y) ∈ A) = ∬_A f(x, y) dx dy per una densità congiunta f ≥ 0 con integrale 1; le marginali si ottengono integrando sull'altra variabile (f_X(x) = ∫ f(x, y) dy), X e Y sono indipendenti se f(x, y) = f_X(x) f_Y(y), e E[g(X, Y)] = ∬ g f. Il punto delicato degli esercizi è descrivere bene la regione dove f > 0.Vettori aleatori assolutamente continui →, Integrali doppi e teorema di FubiniL'integrale doppio di f su un rettangolo si definisce con somme inferiori e superiori su partizioni in rettangolini (per f ≥ 0 è il volume sotto il grafico); su un dominio limitato D si estende f con 0 fuori da D. Area(D) = ∬D 1. Per f continua su un dominio y-semplice {a ≤ x ≤ b, g1(x) ≤ y ≤ g2(x)} vale Fubini: ∬D f = ∫ab (∫g1(x)g2(x) f dy) dx, e simmetricamente per i domini x-semplici; scambiare l'ordine può rendere calcolabile l'integrale.Integrali doppi e teorema di Fubini →, 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 →.

(1) Densità congiunta ed evento

Poiché $x$ e $y$ sono indipendenti (Vettori aleatori assolutamente continuiUn vettore (X, Y) è assolutamente continuo se P((X, Y) ∈ A) = ∬_A f(x, y) dx dy per una densità congiunta f ≥ 0 con integrale 1; le marginali si ottengono integrando sull'altra variabile (f_X(x) = ∫ f(x, y) dy), X e Y sono indipendenti se f(x, y) = f_X(x) f_Y(y), e E[g(X, Y)] = ∬ g f. Il punto delicato degli esercizi è descrivere bene la regione dove f > 0.Vettori aleatori assolutamente continui →), la densità congiunta è il prodotto delle densità marginali $\varphi(a)=\frac1{\sqrt{2\pi}}e^{-a^2/2}$ (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) →): $$f_{xy}(a,b)=\varphi(a)\varphi(b)=\frac1{2\pi}\,e^{-\frac{a^2+b^2}2}.$$ L'evento $\{|x|+|y|\le1\}$ è la regione $D=\{(a,b):|a|+|b|\le1\}$: un rombo (un quadrato ruotato di $45^\circ$) con vertici $(\pm1,0)$ e $(0,\pm1)$, di area $2$. La probabilità è l'integrale della densità su $D$: $$P=\iint_D \varphi(a)\varphi(b)\,da\,db.$$

(2) Perché le coordinate polari non sono la strada comoda

La densità dipende solo da $a^2+b^2$, quindi per un cerchio le coordinate polari (Cambio di variabili negli integrali doppi e coordinate polariSe Φ(u, v) = (x, y) è C¹, iniettiva (salvo insiemi di area nulla) con det JΦ ≠ 0, allora ∬D f dx dy = ∬D' f(Φ(u, v)) |det JΦ(u, v)| du dv: il fattore |det J| misura come Φ dilata le aree. Trasformazioni lineari: |det A| costante (ellisse: area πab). Coordinate polari x = r cos θ, y = r sin θ: dx dy = r dr dθ, adatte a dischi, corone e settori. Applicazioni: area, massa, baricentro, momento d'inerzia.Cambio di variabili negli integrali doppi e coordinate polari →) danno subito il risultato: $P[x^2+y^2\le r^2]=\int_0^{r}\rho\,e^{-\rho^2/2}d\rho=1-e^{-r^2/2}$. Ma il bordo del rombo, in polari, è $\rho(\theta)=\dfrac1{|\cos\theta|+|\sin\theta|}$: si ottiene $$P=\frac1{2\pi}\int_0^{2\pi}\left(1-e^{-\frac1{2(|\cos\theta|+|\sin\theta|)^2}}\right)d\theta,$$ un integrale unidimensionale valido ma senza primitiva elementare (lo si può solo calcolare numericamente). Conviene usare invece la geometria del rombo, che riconduce tutto alla funzione di distribuzione $\Phi$ della gaussiana standard.

(3) Si riduce a un quarto del rombo

La densità è pari in $a$ e pari in $b$ separatamente, e anche $D$ è simmetrico rispetto ai due assi. I quattro triangoli in cui gli assi dividono il rombo hanno quindi la stessa probabilità: $$P=4\,P[\,x\ge0,\ y\ge0,\ x+y\le1\,]=4\iint_{T}\varphi(a)\varphi(b)\,da\,db,\qquad T=\{0\le a\le1,\ 0\le b\le1-a\}.$$

(4) Integrale iterato con Fubini

$T$ è un dominio semplice rispetto ad $a$: fissato $a\in[0,1]$, $b$ corre da $0$ a $1-a$ (Integrali doppi e teorema di FubiniL'integrale doppio di f su un rettangolo si definisce con somme inferiori e superiori su partizioni in rettangolini (per f ≥ 0 è il volume sotto il grafico); su un dominio limitato D si estende f con 0 fuori da D. Area(D) = ∬D 1. Per f continua su un dominio y-semplice {a ≤ x ≤ b, g1(x) ≤ y ≤ g2(x)} vale Fubini: ∬D f = ∫ab (∫g1(x)g2(x) f dy) dx, e simmetricamente per i domini x-semplici; scambiare l'ordine può rendere calcolabile l'integrale.Integrali doppi e teorema di Fubini →). L'integrale interno è un'area sotto la densità gaussiana, cioè una differenza di $\Phi$: $$\int_0^{1-a}\varphi(b)\,db=\Phi(1-a)-\Phi(0)=\Phi(1-a)-\frac12.$$ Quindi $$\boxed{P=4\int_0^1\varphi(a)\left[\Phi(1-a)-\frac12\right]da}$$ Per scrivere lo stesso risultato su tutto l'intervallo $[-1,1]$ si usa la parità: fissato $a$, $b$ varia in $[-(1-|a|),\,1-|a|]$ e l'integrale interno vale $2\Phi(1-|a|)-1$, da cui $P=\int_{-1}^{1}\varphi(a)\,[2\Phi(1-|a|)-1]\,da$.

L'integrale non ha una forma chiusa elementare (è la probabilità di una regione triangolare per una gaussiana bidimensionale), ma è un integrale di una funzione liscia su $[0,1]$. Valori dell'integrando $g(a)=\varphi(a)[\Phi(1-a)-\tfrac12]$, calcolati con Python:

Con la regola di Simpson (passo $h=0{,}25$): $\int_0^1 g\approx\frac{0{,}25}3\left[0{,}1362+4(0{,}1057)+2(0{,}0674)+4(0{,}0297)+0\right]=0{,}0677$, e quindi $P\approx4\cdot0{,}0677=0{,}271$. Il calcolo con scipy.integrate.quad dà $\int_0^1 g=0{,}067730$ e

$$P[\,|x|+|y|\le1\,]=0{,}27092.$$

(Equivalente: $P=4\int_0^1\varphi(a)\Phi(1-a)\,da-(2\Phi(1)-1)=0{,}95361-0{,}68269$.)

(5) Controllo di ragionevolezza con cerchi e quadrato

Il rombo contiene il cerchio inscritto di raggio $\frac1{\sqrt2}$ ed è contenuto nel cerchio circoscritto di raggio $1$; per i cerchi vale $1-e^{-r^2/2}$ (punto (2)): $$1-e^{-1/4}=0{,}2212\ \le\ P\ \le\ 1-e^{-1/2}=0{,}3935.$$ Inoltre il rombo sta nel quadrato $[-1,1]^2$, dove $x$ e $y$ indipendenti danno $P\le(2\Phi(1)-1)^2=0{,}4661$. Il valore $0{,}2709$ rispetta tutti i limiti.

(6) Verifica numerica e Monte Carlo

Il codice calcola l'integrale unidimensionale, l'integrale doppio sul rombo, quello in polari del punto (2) e una stima Monte Carlo (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 →): si estraggono $n$ coppie indipendenti $N(0,1)$ e si calcola la frequenza relativa delle coppie con $|x|+|y|\le1$, che per la legge dei grandi numeri tende a $P$ con errore standard $\sqrt{P(1-P)/n}$.

import numpy as np
from scipy import stats, integrate

phi, Phi = stats.norm.pdf, stats.norm.cdf          # densita' e funzione di distribuzione di N(0,1)

# integrale unidimensionale: P = 4 * int_0^1 phi(x) [Phi(1-x) - 1/2] dx
I, _ = integrate.quad(lambda x: phi(x) * (Phi(1 - x) - 0.5), 0, 1)
print("4*I =", 4 * I)

# controllo 1: integrale doppio numerico sul rombo
P2, _ = integrate.dblquad(lambda y, x: phi(x) * phi(y), -1, 1, lambda x: -(1 - abs(x)), lambda x: 1 - abs(x))
print("dblquad =", P2)

# controllo 2: coordinate polari, r_max(theta) = 1/(|cos|+|sin|); P = (1/2pi) int (1 - exp(-r_max^2/2)) dtheta
P3, _ = integrate.quad(lambda th: (1 - np.exp(-1 / (2 * (abs(np.cos(th)) + abs(np.sin(th)))**2))) / (2 * np.pi),
                       0, 2 * np.pi, points=[np.pi / 2, np.pi, 3 * np.pi / 2])
print("polari =", P3)

# controllo 3: Monte Carlo
rng = np.random.default_rng(1)
for n in (10**3, 10**5, 10**7):
    x, y = rng.standard_normal(n), rng.standard_normal(n)
    p = np.mean(np.abs(x) + np.abs(y) <= 1)
    print(f"n = {n:>8d}: stima {p:.5f}, errore standard {np.sqrt(p * (1 - p) / n):.5f}")

# limiti: quadrato circoscritto, cerchi inscritto e circoscritto
print("quadrato", (2 * Phi(1) - 1)**2, "disco r=1/sqrt2", 1 - np.exp(-1 / 4), "disco r=1", 1 - np.exp(-1 / 2))

Risultati ottenuti:

• i tre calcoli numerici deterministici (formula in $\Phi$, dblquad, polari) coincidono: $0{,}270920$;
• Monte Carlo: $n=10^3$ dà $0{,}259$ (errore standard $0{,}014$), $n=10^5$ dà $0{,}2725$ ($0{,}0014$), $n=10^7$ dà $0{,}27100$ ($0{,}00014$). Tutte le stime stanno entro due errori standard dal valore $0{,}27092$, e l'errore scala come $1/\sqrt n$: cento volte più campioni, dieci volte meno errore.

Risposta

$$P[\,|x|+|y|\le1\,]=4\int_0^1\varphi(a)\left[\Phi(1-a)-\tfrac12\right]da\approx0{,}2709.$$

Errori tipici.

• Usare per il rombo la formula del cerchio $1-e^{-r^2/2}$, oppure calcolare $P[|x|\le1]\cdot P[|y|\le1]=0{,}4661$ (che è il quadrato, un'area più grande).
• Dimenticare il fattore $4$ della simmetria, o il termine $-\tfrac12$ (cioè $-\Phi(0)$) nell'integrale interno.
• Sommare i valori assoluti come se fossero indipendenti dalla somma: $|x|+|y|$ non è gaussiana, quindi non si può standardizzare e leggere $\Phi$ in un solo punto.

Vedi anche: Esercizio - variabili aleatorie, istogrammi e momenti in Python (laboratorio 6) (generazione e controllo empirico di variabili aleatorie), 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 →.