Esercizio - Dimostrazioni su esponenziale, Erlang e processo di Poisson (domande da orale)

Teoria: Distribuzione esponenziale e processo di Poisson per l'informaticaL'esponenziale Exp(λ) è l'unica legge continua senza memoria: P(X>s+t | X>s)=P(X>t). Il minimo di esponenziali indipendenti è Exp(Σλ_i) e l'indice del minimo è i con probabilità λ_i/Σλ (indipendente dal valore del minimo); il massimo di n esponenziali Exp(λ) ha media H_n/λ. Il processo di Poisson di intensità λ ha interarrivi Exp(λ) i.i.d., N(t) ~ Poisson(λt) e incrementi indipendenti; la scomposizione (ogni evento tenuto con probabilità p) dà due processi di Poisson indipendenti, la sovrapposizione dà un processo di intensità somma. Applicazioni: carichi Poisson delle urne (poissonizzazione, con costo e√m), server replicati, arrivi durante un servizio.Distribuzione esponenziale e processo di Poisson per l'informatica →.

Testo. Dimostrare e verificare numericamente:

1. Il massimo di tre variabili $\text{Exp}(1)$ indipendenti ha media $\tfrac{11}6$ e varianza $\tfrac{49}{36}$.
2. Se $X\sim\text{Exp}(\lambda)$ allora $\lceil X\rceil$ è geometrica di parametro $1-e^{-\lambda}$. Per $\lambda=\ln2$, scrivere la legge.
3. Dualità Erlang-Poisson: per $S_k\sim\text{Erlang}(k,\lambda)$ vale $P(S_k>t)=P(N(t)\le k-1)$, con $N(t)\sim\text{Poisson}(\lambda t)$. Verificarlo per $k=2$, $\lambda=2$, $t=1$.
4. In un processo di Poisson di intensità $\lambda=3$, dato $N(2)=1$, l'istante dell'evento è uniforme su $[0,2]$: calcolare $P(S_1\le0{,}5\mid N(2)=1)$.

Soluzione

1. Il massimo è la somma di tre tempi indipendenti: il primo completamento $T_1\sim\text{Exp}(3)$, poi (assenza di memoria) $T_2\sim\text{Exp}(2)$, poi $T_3\sim\text{Exp}(1)$, tutti indipendenti. Quindi $$E[\max]=\tfrac13+\tfrac12+1=\tfrac{11}6=1{,}8333,\qquad\text{Var}[\max]=\tfrac19+\tfrac14+1=\tfrac{4+9+36}{36}=\tfrac{49}{36}=1{,}3611.$$ (La varianza della somma di variabili indipendenti è la somma delle varianze, $\tfrac1{\mu^2}$ per ciascuna.) Simulazione con $10^6$ terne (codice sotto): media $1{,}834$, varianza $1{,}361$.

2. Per $k\ge1$, $P(\lceil X\rceil=k)=P(k-1<X\le k)=e^{-\lambda(k-1)}-e^{-\lambda k}=\bigl(e^{-\lambda}\bigr)^{k-1}\bigl(1-e^{-\lambda}\bigr)$: è la geometrica (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 →) con probabilità di successo $p=1-e^{-\lambda}$. Per $\lambda=\ln2$: $e^{-\lambda}=\tfrac12$, $p=\tfrac12$ e $P(\lceil X\rceil=k)=2^{-k}$: $\tfrac12,\tfrac14,\tfrac18,\tfrac1{16},\dots$ La simulazione dà $0{,}4992;\ 0{,}2505;\ 0{,}1251;\ 0{,}0626;\ 0{,}0314$. Questo è l'ingrediente di un campionatore: un'esponenziale discretizzata è geometrica (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 →).

3. Dimostrazione. $S_k>t$ significa che il $k$-esimo evento avviene dopo $t$, cioè che in $[0,t]$ ci sono al più $k-1$ eventi: $\{S_k>t\}=\{N(t)\le k-1\}$ (stesso evento, per definizione di $N$ a partire da $S_k$). Quindi $P(S_k>t)=\sum_{j=0}^{k-1}e^{-\lambda t}\frac{(\lambda t)^j}{j!}$. Numeri: $k=2$, $\lambda=2$, $t=1$: $P(S_2>1)=e^{-2}(1+2)=3e^{-2}=0{,}4060$. Direttamente, con la densità di Erlang: $P(S_2>1)=\int_1^\infty4se^{-2s}ds=\bigl[-(2s+1)e^{-2s}\bigr]_1^\infty=3e^{-2}=0{,}4060$ ✓ (simulazione: $0{,}4064$).

4. Con $N(2)=1$, la densità congiunta di $S_1=s$ e $\{N(2)=1\}$ è $\lambda e^{-\lambda s}\cdot e^{-\lambda(2-s)}=\lambda e^{-2\lambda}$, costante in $s\in[0,2]$; dividendo per $P(N(2)=1)=2\lambda e^{-2\lambda}$ si ottiene la densità condizionata $\tfrac12$: uniforme su $[0,2]$. Quindi $$P(S_1\le0{,}5\mid N(2)=1)=\frac{0{,}5}{2}=\boxed{0{,}25}\quad\text{per ogni }\lambda.$$ La simulazione con $\lambda=3$ ($8856$ casi su $600\,000$ con esattamente un evento, contro i $600\,000\cdot6e^{-6}\approx8924$ attesi) dà $0{,}2511$. Il risultato non dipende da $\lambda$: l'intensità decide quanti eventi cadono in $[0,t]$, non dove.

Codice di verifica

import numpy as np
from math import exp, log
rng = np.random.default_rng(3)

M = rng.exponential(1, size=(1_000_000, 3)).max(axis=1)
print("1.", M.mean(), M.var(), "attesi", 11 / 6, 49 / 36)

C = np.ceil(rng.exponential(1 / log(2), 1_000_000))
print("2.", [round((C == k).mean(), 4) for k in range(1, 6)])

S2 = rng.exponential(1 / 2, size=(1_000_000, 2)).sum(axis=1)
print("3.", (S2 > 1).mean(), 3 * exp(-2))

lam, t, tot, hit = 3, 2, 0, 0
for _ in range(600_000):
    a = rng.exponential(1 / lam)
    if a <= t and a + rng.exponential(1 / lam) > t:      # esattamente un evento in [0, t]
        tot += 1; hit += a <= 0.5
print("4.", tot, hit / tot)