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

In questa pagina 4

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

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)$ ).
Funzione di distribuzione. Rappresentare la funzione di distribuzione empirica $F_x(a)=P[x\le a]$ .
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$ .
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 continue, Distribuzioni uniforme continua ed esponenziale, Distribuzione gaussiana (normale), Funzione di distribuzione, Valore atteso, Varianza e momenti, Variabili aleatorie discrete e densità discreta, Processi aleatori - definizioni, media e autocorrelazione, 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):

	densità $f_x(a)$	$F_x(a)$	$m_x$	$M_x$	$\sigma_x^2$
uniforme $U[a_1,a_2]$	$\frac1{a_2-a_1}$ su $[a_1,a_2]$	$\frac{a-a_1}{a_2-a_1}$ su $[a_1,a_2]$	$\frac{a_1+a_2}2$	$\frac{a_1^2+a_1a_2+a_2^2}3$	$\frac{(a_2-a_1)^2}{12}$
gaussiana $N(m,\sigma^2)$	$\frac1{\sqrt{2\pi}\,\sigma}e^{-\frac{(a-m)^2}{2\sigma^2}}$	$\Phi\!\left(\frac{a-m}\sigma\right)$	$m$	$m^2+\sigma^2$	$\sigma^2$
esponenziale (media $\mu$ )	$\frac1\mu e^{-a/\mu}1(a)$	$\left(1-e^{-a/\mu}\right)1(a)$	$\mu$	$2\mu^2$	$\mu^2$

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 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$ .

Grafico interattivo: Le tre densità teoriche del laboratorio: uniforme su [−1, 3] (altezza 1/4), gaussiana N(2, 1,5²) (picco 0,266) ed esponenziale di media 0,99 (valore 1,01 in 0)

Grafico interattivo: Funzioni di distribuzione teoriche dell'uniforme su [−1, 3] (rampa da 0 a 1) e dell'esponenziale di media 0,99 (1 − e^{−a/0,99})

(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): 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 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 distribuzione, 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 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 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 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$ .

python

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):

	$m_x$	$M_x$	$\sigma_x^2$	$M_x-m_x^2$ (empirico)
uniforme	$1{,}0000$ ; $1{,}0112$	$2{,}3333$ ; $2{,}3668$	$1{,}3333$ ; $1{,}3456$	$1{,}3442$
gaussiana	$2{,}0000$ ; $2{,}0394$	$6{,}2500$ ; $6{,}5428$	$2{,}2500$ ; $2{,}3860$	$2{,}3836$
esponenziale	$0{,}9900$ ; $0{,}9661$	$1{,}9602$ ; $1{,}8739$	$0{,}9801$ ; $0{,}9415$	$0{,}9406$

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).

Versione ripasso

Generazione. Uniforme: $a_1+(a_2-a_1)\,u$ con rng.random. Gaussiana: $m+\sigma z$ con standard_normal. Esponenziale: rng.exponential(mu), con media $\mu$ .
Teoria. Uniforme su $[-1,3]$ : $m=1$ , $M=\frac73$ , $\sigma^2=\frac43$ . Gaussiana $N(2,1{,}5^2)$ : $m=2$ , $M=6{,}25$ , $\sigma^2=2{,}25$ . Esponenziale con $\mu=0{,}99$ : $m=0{,}99$ , $M=2\mu^2=1{,}9602$ , $\sigma^2=\mu^2$ . Vale sempre $M=\sigma^2+m^2$ .
Istogrammi. Con density=True l'altezza delle barre approssima $f_x(a)$ , perché $P[a\in I]\approx f_x\,\Delta$ . Con $1000$ campioni le barre sono frastagliate, con errore relativo circa $1/\sqrt{k_j}$ .
Funzione di distribuzione empirica. La scalinata passa a $\frac iN$ in $x_{(i)}$ . Lo scarto massimo dalla teoria è $0{,}021$ , $0{,}028$ , $0{,}033$ per uniforme, gaussiana, esponenziale, sotto il limite $1{,}36/\sqrt N=0{,}043$ del test di Kolmogorov al $5\%$ .
Momenti empirici. x.var(ddof=1) è il comando var di Matlab, che divide per $N-1$ . Le differenze dai teorici sono dell'ordine di $\sigma/\sqrt N$ , come per la media gaussiana con scarto $0{,}039$ .
M-PAM, $M=4$ , $d=0{,}5$ . Alfabeto $\{\pm0{,}5,\pm1{,}5\}$ , frequenze relative $\approx\frac14$ . Teoria: $m=0$ , $M_x=\frac{d^2(M^2-1)}3=1{,}25$ ; simulazione $1{,}2498$ .
Correlazione. $\hat r_x(0)=1{,}2498$ , e $\hat r_x(k)$ per $k=1,\dots,20$ sono in modulo $\le0{,}005$ , cioè dell'ordine di $M_x/\sqrt N=0{,}0024$ : rumore bianco a tempo discreto. Con la divisione per $N$ (stima biased) i valori sono confrontabili con $M_x$ .