Esercizio - variabili aleatorie, istogrammi e momenti in Python (laboratorio 6)

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

Testo (laboratorio 6, parti 1-4, corso Teoria dei Segnali, UniPD).

1. Generazione e densità. Generare $N=1000$ campioni di tre variabili aleatorie: uniforme $U[a_1,a_2]$ con $a_1=-1$, $a_2=3$; gaussiana $N(m,\sigma^2)$ con $m=2$, $\sigma=1{,}5$; esponenziale di media $\mu=0{,}99$. Rappresentare per ciascuna l'istogramma normalizzato (stima della densità $f_x(a)$).
2. Funzione di distribuzione. Rappresentare la funzione di distribuzione empirica $F_x(a)=P[x\le a]$.
3. Momenti. Calcolare, teoricamente e dai campioni, la media $m_x=E[x]$, la potenza statistica $M_x=E[x^2]$ e la varianza $\sigma_x^2$, e verificare $M_x-m_x^2=\sigma_x^2$.
4. Variabile discreta. Alfabeto M-PAM con $M=4$ simboli equispaziati, $d=0{,}5$, $N=2^{18}$ simboli estratti in modo indipendente ed equiprobabile: istogramma delle probabilità, funzione di distribuzione, momenti e correlazione $r_x(kT)$ della sequenza.

Teoria usata: Variabili aleatorie assolutamente continueX è assolutamente continua se F_X(x) = ∫ da −∞ a x di f_X(t) dt per una densità f_X ≥ 0 con integrale 1; allora P(a < X ≤ b) = ∫ da a a b di f_X, P(X = x) = 0 per ogni x, f_X = F_X' dove F_X è derivabile, E[g(X)] = ∫ g(x) f_X(x) dx, e varianza, momenti e disuguaglianze funzionano come nel caso discreto con gli integrali al posto delle somme.Variabili aleatorie assolutamente continue →, Distribuzioni uniforme continua ed esponenzialeU(a, b) ha densità costante 1/(b − a) su [a, b], media (a + b)/2 e varianza (b − a)²/12; Exp(λ) ha densità λe^(−λx) per x ≥ 0, FdD 1 − e^(−λx), P(X > t) = e^(−λt), media 1/λ, varianza 1/λ², ed è l'unica legge continua senza memoria (versione continua della geometrica).Distribuzioni uniforme continua ed esponenziale →, 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) →, 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 →, Valore attesoIl valore atteso E[X] = Σ x p_X(x) è la media dei valori di X pesata con le loro probabilità (esiste se la serie converge assolutamente); per una funzione g vale E[g(X)] = Σ g(x) p_X(x) senza trovare la legge di g(X), ed E è lineare: E[aX + bY + c] = aE[X] + bE[Y] + c.Valore atteso →, Varianza e momentiI momenti E[X^k] e i momenti centrati E[(X − μ)^k] descrivono la forma di una legge; la varianza Var(X) = E[(X − μ)²] = E[X²] − E[X]² misura quanto X si disperde attorno alla media, vale Var(aX + b) = a² Var(X) e Var(X) = 0 solo se X è costante.Varianza e momenti →, Variabili aleatorie discrete e densità discretaUna variabile aleatoria discreta è una funzione X da Ω in R che assume un insieme finito o numerabile di valori (l'alfabeto); la sua densità discreta p_X(x) = P(X = x) basta a calcolare la probabilità di ogni evento che riguarda X.Variabili aleatorie discrete e densità discreta →, Processi aleatori - definizioni, media e autocorrelazioneUn processo aleatorio x(t), t in I (R o Z(T)), è una famiglia di variabili aleatorie sullo stesso spazio di probabilità; fissato l'esito si ottiene una realizzazione (un segnale). Si descrive con le densità di ordine N (complete), in particolare del primo e del secondo ordine, oppure solo con media m_x(t) e correlazione r_x(t,s) = E[x(t)x*(s)] (descrizione di potenza). Un processo gaussiano è determinato da media e correlazione.Processi aleatori - definizioni, media e autocorrelazione →, 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) Teoria da confrontare con i campioni

Per le tre variabili, densità, funzione di distribuzione e momenti sono ($1(a)$ è il gradino):

Con i dati del laboratorio: uniforme $m_x=1$, $M_x=\frac{1-3+9}3=\frac73=2{,}3333$, $\sigma_x^2=\frac{16}{12}=1{,}3333$; gaussiana $m_x=2$, $M_x=4+2{,}25=6{,}25$, $\sigma_x^2=2{,}25$; esponenziale $m_x=0{,}99$, $M_x=2\cdot0{,}9801=1{,}9602$, $\sigma_x^2=0{,}9801$. Il legame $M_x=\sigma_x^2+m_x^2$ vale per ogni variabile: $E[(x-m_x)^2]=E[x^2]-2m_xE[x]+m_x^2=M_x-m_x^2$ (Varianza e momentiI momenti E[X^k] e i momenti centrati E[(X − μ)^k] descrivono la forma di una legge; la varianza Var(X) = E[(X − μ)²] = E[X²] − E[X]² misura quanto X si disperde attorno alla media, vale Var(aX + b) = a² Var(X) e Var(X) = 0 solo se X è costante.Varianza e momenti →). Per l'esponenziale, per esempio, $M_x=\int_0^\infty a^2\frac1\mu e^{-a/\mu}da=2\mu^2$ (due integrazioni per parti), quindi $\sigma_x^2=2\mu^2-\mu^2=\mu^2$.

(2) Come si generano i campioni

• Uniforme. rng.random(N) dà $u\sim U[0,1]$; la trasformazione affine $x=a_1+(a_2-a_1)u$ porta l'intervallo $[0,1]$ in $[a_1,a_2]$ (Matlab rand).
• Gaussiana. rng.standard_normal(N) dà $z\sim N(0,1)$; $x=m+\sigma z$ è $N(m,\sigma^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) →: trasformazione affine di una gaussiana; Matlab randn).
• Esponenziale. rng.exponential(mu, N) ha media $\mu$ (Matlab exprnd(mu)).
• Il generatore è inizializzato con un seme fisso (default_rng(1)) per ripetere gli stessi numeri; cambiando il seme cambiano le cifre, non le conclusioni.

Istogramma normalizzato. Si divide l'asse in $n_b$ intervalli (nel laboratorio $60$, $60$, $80$) di ampiezza $\Delta$, si contano i campioni $k_j$ in ciascuno e si disegna una barra di altezza $\frac{k_j}{N\Delta}$ (density=True, Matlab 'Normalization','pdf'). L'area totale vale $1$ e l'altezza approssima $f_x(a)$ perché $P[a\in\text{intervallo}]\approx k_j/N=f_x\,\Delta$ (Variabili aleatorie assolutamente continueX è assolutamente continua se F_X(x) = ∫ da −∞ a x di f_X(t) dt per una densità f_X ≥ 0 con integrale 1; allora P(a < X ≤ b) = ∫ da a a b di f_X, P(X = x) = 0 per ogni x, f_X = F_X' dove F_X è derivabile, E[g(X)] = ∫ g(x) f_X(x) dx, e varianza, momenti e disuguaglianze funzionano come nel caso discreto con gli integrali al posto delle somme.Variabili aleatorie assolutamente continue →). Con soli $1000$ campioni e tante barre l'istogramma è molto frastagliato: poche decine di campioni per barra, errore relativo $\approx1/\sqrt{k_j}$.

Funzione di distribuzione empirica. Si ordinano i campioni, $x_{(1)}\le\dots\le x_{(N)}$, e si disegna la scalinata che passa a $\frac iN$ in $x_{(i)}$ (np.sort e np.arange(1, N+1)/N): è la frazione di campioni $\le a$, e per la legge dei grandi numeri tende a $F_x(a)$ (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 →, 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 →). Lo scarto massimo tra la scalinata e la curva teorica (statistica di Kolmogorov-Smirnov) è dell'ordine di $1/\sqrt N\approx0{,}03$.

(3) Momenti: teorici ed empirici

Ogni momento è un valore atteso, e il valore atteso si stima con la media dei campioni (Valore attesoIl valore atteso E[X] = Σ x p_X(x) è la media dei valori di X pesata con le loro probabilità (esiste se la serie converge assolutamente); per una funzione g vale E[g(X)] = Σ g(x) p_X(x) senza trovare la legge di g(X), ed E è lineare: E[aX + bY + c] = aE[X] + bE[Y] + c.Valore atteso →): $\hat m=\frac1N\sum x_i$, $\hat M=\frac1N\sum x_i^2$, $\hat\sigma^2=\frac1{N-1}\sum(x_i-\hat m)^2$. In numpy, x.var() divide per $N$ (come x.var(ddof=0)): per riprodurre il comando var di Matlab, che divide per $N-1$, si scrive x.var(ddof=1). L'errore tipico della stima della media è $\sigma_x/\sqrt N$.

(4) Alfabeto M-PAM e correlazione di una sequenza di simboli

Si costruisce l'alfabeto $A=\{(2k-(M-1))\,d\}_{k=0}^{M-1}$: con $M=4$ e $d=0{,}5$ sono $\{-1{,}5;\,-0{,}5;\,0{,}5;\,1{,}5\}$, simboli equispaziati di passo $2d=1$ e simmetrici rispetto allo zero. Ogni simbolo è estratto con probabilità $\frac1M$, indipendentemente dagli altri. I momenti teorici sono $$m_x=0,\qquad M_x=\frac{2.25+0.25+0.25+2.25}4=1{,}25=\frac{d^2(M^2-1)}3,\qquad \sigma_x^2=M_x.$$

Correlazione. La sequenza $x(kT)$ è un processo a tempo discreto a simboli indipendenti (Processi aleatori - definizioni, media e autocorrelazioneUn processo aleatorio x(t), t in I (R o Z(T)), è una famiglia di variabili aleatorie sullo stesso spazio di probabilità; fissato l'esito si ottiene una realizzazione (un segnale). Si descrive con le densità di ordine N (complete), in particolare del primo e del secondo ordine, oppure solo con media m_x(t) e correlazione r_x(t,s) = E[x(t)x*(s)] (descrizione di potenza). Un processo gaussiano è determinato da media e correlazione.Processi aleatori - definizioni, media e autocorrelazione →): per $k\ne0$ i simboli $x(nT)$ e $x(nT+kT)$ sono indipendenti, quindi $r_x(kT)=E[x(nT+kT)x(nT)]=m_x^2$, mentre per $k=0$ si ha la potenza $M_x$. Per una sequenza a media nulla: $$r_x(kT)=\begin{cases}M_x=1{,}25 & k=0\\ 0 & k\ne0\end{cases}$$ e la densità spettrale di potenza (Processi stazionari e densità spettrale di potenzaUn processo è stazionario in media, in potenza, nel primo ordine o in correlazione se la rispettiva grandezza non cambia traslando il tempo; stazionario in senso lato (sl) = media costante e r_x(t,s) = r_x(t-s); in senso stretto (ss) = tutte le densità invarianti. La densità spettrale di potenza R_x(f) = F[r_x(τ)] è non negativa e ha integrale = potenza statistica. Rumore bianco: R costante. Ciclostazionario: statistiche periodiche. Ergodico in media: la media temporale di una realizzazione converge a m_x.Processi stazionari e densità spettrale di potenza →) è costante, $R_x(f)=T\sum_kr_x(kT)e^{-i2\pi fkT}=T\,M_x$: la sequenza è un "rumore bianco" a tempo discreto. La stima dai campioni è $\hat r_x(kT)=\frac1N\sum_{n=0}^{N-1-k}x[n]x[n+k]$ (divisa per $N$, non per $N-k$: stima "biased"). Il file del corso chiama xcorr senza l'opzione 'biased': i valori che ne escono sono $N$ volte più grandi ($r(0)=N\cdot M_x\approx3{,}3\times10^5$), mentre qui si divide per $N$ per confrontare con $M_x$.

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

rng = np.random.default_rng(1)          # seme fisso: risultati riproducibili
N = 1000

# 1) generazione
a1, a2 = -1, 3
x_unif = a1 + (a2 - a1) * rng.random(N)        # x = a1 + (a2-a1) u, u ~ U[0,1]   (MATLAB: rand)
m, sigma = 2, 1.5
x_gauss = m + sigma * rng.standard_normal(N)   # x = m + sigma z, z ~ N(0,1)      (MATLAB: randn)
mu = 0.99
x_exp = rng.exponential(mu, N)                 # esponenziale di media mu         (MATLAB: exprnd)

# 2) istogrammi normalizzati (densita') con la densita' teorica, e CDF empiriche con la teorica
dist = [("uniforme", x_unif, 60, stats.uniform(a1, a2 - a1)),
        ("gaussiana", x_gauss, 60, stats.norm(m, sigma)),
        ("esponenziale", x_exp, 80, stats.expon(scale=mu))]
fig, ax = plt.subplots(2, 3, figsize=(12, 6))
for j, (nome, x, nb, d) in enumerate(dist):
    xs = np.linspace(-2, 8, 1000)
    ax[0, j].hist(x, bins=nb, density=True)               # MATLAB: histogram(..., 'Normalization','pdf')
    ax[0, j].plot(xs, d.pdf(xs), "r")
    ax[1, j].step(np.sort(x), np.arange(1, N + 1) / N, where="post")   # CDF empirica: i-esimo valore ordinato -> i/N
    ax[1, j].plot(xs, d.cdf(xs), "r")
    # distanza massima tra CDF empirica e teorica (statistica di Kolmogorov-Smirnov)
    print(f"{nome:13s} scarto massimo CDF empirica-teorica: {stats.kstest(x, d.cdf).statistic:.4f}")

# 3) momenti: teorici e empirici
casi = {
    "uniforme":     (x_unif,  (a1 + a2) / 2, (a1**2 + a1 * a2 + a2**2) / 3, (a2 - a1)**2 / 12),
    "gaussiana":    (x_gauss, m,             m**2 + sigma**2,               sigma**2),
    "esponenziale": (x_exp,   mu,            2 * mu**2,                     mu**2),
}
print(f"{'':13s} {'m (teo, emp)':>16s} {'M (teo, emp)':>16s} {'var (teo, emp)':>18s} {'M - m^2':>9s}")
for nome, (x, m_t, M_t, v_t) in casi.items():
    m_e, M_e, v_e = x.mean(), (x**2).mean(), x.var(ddof=1)    # ddof=1: come MATLAB var (divide per N-1)
    print(f"{nome:13s} {m_t:7.4f} {m_e:7.4f}  {M_t:7.4f} {M_e:7.4f}   {v_t:7.4f} {v_e:7.4f}   {M_e - m_e**2:8.4f}")

# 4) alfabeto M-PAM e correlazione di una sequenza a simboli indipendenti
M, d, Nsim = 4, 0.5, 2**18
A = (2 * np.arange(M) - (M - 1)) * d                 # {-1.5, -0.5, 0.5, 1.5}
x = A[rng.integers(0, M, Nsim)]                      # simboli indipendenti, equiprobabili (MATLAB: A(randi([1,M],1,N)))
print("alfabeto", A, "| frequenze relative", np.round([np.mean(x == a) for a in A], 4))
print(f"m = {x.mean():.4f}, M_x = {np.mean(x**2):.4f} (teo {d**2 * (M**2 - 1) / 3}), var = {x.var(ddof=1):.4f}")
tau_max = 20
r = np.array([x[: Nsim - k] @ x[k:] / Nsim for k in range(tau_max + 1)])   # r_x(kT) = (1/N) sum x[n] x[n+k]  ('biased')
print("r_x(0..4) =", np.round(r[:5], 4), "| max |r_x(k)| per k = 1..20:", round(np.abs(r[1:]).max(), 4))
plt.show()

Cosa si vede e valori ottenuti (seme $1$; cambiando seme le cifre cambiano di poco).

• Istogrammi: l'uniforme è un plateau di altezza circa $0{,}25$ su $[-1,3]$ (con fluttuazioni tra le barre), la gaussiana una campana centrata in $2$ con picco $\approx0{,}27$, l'esponenziale parte da $\approx1{,}0$ in $0$ e decade; le curve teoriche passano in mezzo alle barre.
• Funzioni di distribuzione: tre scalinate molto fitte, quasi indistinguibili dalle curve teoriche (rampa lineare, sigmoide, $1-e^{-a/\mu}$). Scarto massimo: $0{,}021$ (uniforme), $0{,}028$ (gaussiana), $0{,}033$ (esponenziale), tutti sotto $1{,}36/\sqrt N=0{,}043$, il limite al $5\%$ del test di Kolmogorov.
• Momenti (teorico, empirico):

Le differenze sono fluttuazioni statistiche di ordine $\sigma/\sqrt N$ (per la media gaussiana $1{,}5/\sqrt{1000}=0{,}047$: scarto di $0{,}039$). La verifica $M_x-m_x^2=\sigma_x^2$ è soddisfatta a meno del fattore $\frac{N-1}N$ tra le due definizioni di varianza (per esempio $1{,}3442$ contro $1{,}3456=1{,}3442\cdot\frac{1000}{999}$).

• M-PAM: frequenze relative $0{,}2497;\ 0{,}2505;\ 0{,}2496;\ 0{,}2502$ (teoria $0{,}25$); $m=0{,}0002$, $M_x=1{,}2498$ (teoria $1{,}25$), varianza $1{,}2499$. Correlazione: $\hat r_x(0)=1{,}2498$ e $\hat r_x(k)$ per $k=1,\dots,20$ tutti in modulo $\le0{,}005$, cioè dell'ordine dell'errore standard $M_x/\sqrt N=0{,}0024$ attorno a zero: un solo impulso in $k=0$, come per il rumore bianco.

Errori tipici.

• Usare x.var() e aspettarsi il risultato di Matlab (var divide per $N-1$) oppure, al contrario, confrontare $M_x-m_x^2$ con una varianza calcolata con denominatore diverso.
• Non normalizzare l'istogramma: le altezze dipendono da $N$ e dalla larghezza delle barre e non sono la densità.
• Confrontare la correlazione empirica non normalizzata (xcorr senza 'biased') con $M_x$.

Vedi anche: Esercizio - processi aleatori, media d'insieme, correlazione e potenza temporale in Python (laboratorio 6), Esercizio - rumore bianco gaussiano e filtro passa-basso in Python (laboratorio 6).