5.4 Stima dei parametri da osservazioni indipendenti

In questa sezione estendiamo i risultati della precedente al caso in cui osservano variabili indipendenti $(X_i)_{i=1}^n$, tutte gaussiane con i medesimi parametri (il termine statistico è un campione di taglia $n$, il numero di osservazioni). Per avvicinare la notazione al caso precedente, si può alternativamente descrivere la situazione dicendo che il robot stima i parametri sulla base dell’osservazione di un (singolo) vettore aleatorio gaussiano $X = (X_i)_{i=1}^n$ a valori in $\R^n$. Nella sezione successiva indichiamo come estendere al caso di $n$ osservazioni indipendenti di variabili gaussiane a loro volta vettoriali.

La situazione che stiamo descrivendo è molto comune, ad esempio se ciascuna $X_i$ rappresenta una misura della stessa quantità (rappresentata dal parametro di media) affetta da una imprecisione (un “rumore”) di una opportuna intensità (la deviazione standard): sfruttando l’indipendenza risulta che l’effetto del rumore è mitigato e si ottiene una stima più precisa.

Introduciamo quindi la seguente notazione analoga a quella della sezione precedente: i parametri di cascuna $X_i$, ossia la media $\E{X_i}$ e la varianza $\Var{X_i}$ non dipendono da $i$, e scriviamo quindi \[ m = \E{X_i} \quad v = \Var{X_i} \quad \text{per ogni $i=1, \ldots, n$.}\] Come nel caso della singola osservazione, si introducono le rispettive variabili aleatorie $M$ per la media, $V$ per la varianza (la seconda a valori positivi).

La rete bayesiana associata è rappresentata in figura e generalizza quella della sezione precedente (abbiamo introdotto la variabile congiunta $(M,V)$ per semplificare la notazione).

$Rete bayesiana tra le variabili $(M,V)$ e le $(X_i)_{i=1}^n$ per $n=6$$

Figura 5.14: Rete bayesiana tra le variabili $(M,V)$ e le $(X_i)_{i=1}^n$ per $n=6$

Poter disporre di più osservazioni indipendenti degli stessi parametri intuitivamente permette di ottenere stime più precise, un fatto che ora vediamo sia tramite l’approccio di massima verosimiglianza che quello bayesiano.

5.4.1 Stima di massima verosimiglianza

L’ipotesi di indipendenza, noti $m$ e $v$, si traduce nel fatto che la densità dell’osservazione di $X=x$, ossia la congiunta delle osservazioni $X_i=x_i$ per $i=1,\ldots, n$, è la verosimiglianza \[ \begin{split} L(m, v; x)& = p(X=x|M=m,V=v) = p(X_1 = x_1, \ldots X_n = x_n | m, v) \\ & = p(X_1 = x_1 | m,v) \cdot \ldots \cdot p(X_n = x_n | m, v) \\ & = \prod_{i=1}^n \exp\bra{ -\frac{1}{2v} (x_i -m)^2 } \frac{1}{\sqrt{ 2 \pi v}} \\ & \propto \exp\bra{ - \frac {n}{2} \bra{ \frac 1 {nv} \sum_{i=1}^n (x_i - m)^2 + \log( v )}},\end{split}\] dove nell’ultimo passaggio abbiamo omesso la costante moltiplicativa $(2 \pi)^{n/2}$ (per semplificare la notazione).

deltam <- 0.1
deltav <- 0.05

m <- seq(0, 2, by = deltam)
v <- seq(0.01, 1, by = deltav)

N_m <- length(m)
N_v <- length(v)

# osservazioni

x <- c(1, 2, 1.5, 0.5)
n <- length(x)


# creiamo una matrice con i valori
# della verosimiglianza

L <- matrix(1, nrow = N_m, ncol = N_v)

for (i in 1:N_m) {
  for (j in 1:N_v) {
    for (obs in x) {
      L[i, j] <- L[i, j] * exp(-(m[i] -
        obs)^2/(2 * v[j]))/sqrt(2 *
        pi * v[j])
    }
  }
}

# usiamo stavolta filled.contour() per
# produrre il grafico e la scala di
# valori accanto

filled.contour(m, v, L, color.palette = viridis,
  xlab = "m", ylab = "v")

Figura 5.15: verosimiglianza per $m$ e $v$, associata alle osservazioni $x = (1, 2, 1.5, 0.5)$.

Dal grafico sopra vediamo che la verosimiglianza è una funzione più interessante del caso di una singola osservazione e in particolare il massimo non è necessariamente per $v=0$. Per applicare il metodo di massima verosimiglianza, passando al logaritmo e moltiplicando per $-2/n$ e tralasciando costanti additive (che non hanno nessun ruolo utile nella procedura) si tratta quindi di determinare $m_{{\mle}}$ e $v_{{\mle}}$ che minimizzano la funzione \[ (m, v) \mapsto \frac{1}{v} \sqa{ \frac 1 {n} \sum_{i=1}^n (x_i - m)^2} + \log( v ) \] Ragionando come nel caso della singola osservazione, deriviamo rispetto ad $m$ (fissato $v$) e imponiamo che la derivata si annulli. Si trova la condizione \[ 2 \sum_{i=1}^n (x_i - m) = 0 \quad \text{da cui} \quad m = \frac{1}{n} \sum_{i=1}^n x_i \] è la media aritmetica delle osservazioni (detta anche media empirica o campionaria, in inglese sample mean), e indicata anche brevemente con \[ \bar{x} = \frac{1}{n} \sum_{i=1}^n x_i,\] Mentre se deriviamo rispetto a $v$, tenendo fisso $m$, si trova analogamente al caso della singola osservazione che \[ -\frac 1 { v^2} \frac 1 n \sum_{i=1}^n(x_i-m)^2 + \frac{1}{v} = 0, \quad \text{da cui} \quad v = \frac 1 n \sum_{i=1}^n (x_i-m)^2.\] dovendo massimizzare la funzione delle due variabili si trova quindi la coppia \[ m_{{\mle}} =\bar{x} = \frac 1 n \sum_{i=1}^n x_i \quad v_{{\mle}} = \frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2,\] dove l’ultima quantità è detta anche varianza campionaria (sample variance in inglese). Volendo specificare, questa è la versione detta anche distorta (in inglese biased) della varianza campionaria, dove la versione “corretta” o meglio non distorta (unbiased) contiene invece il fattore $n-1$ a denominatore, una differenza minima quando $n$ è grande, sulle cui ragioni non ci soffermiamo – tuttavia va tenuto presente che molte funzioni in R usano appunto la versione non distorta.

Come nel caso della singola osservazione, i calcoli sopra ci mostrano anche che, 1. se il pararametro di varianza $V=v_0 >0$ è noto (ossia $V=v_0$ è costante), la stima di massima verosimiglianza per il parametro di media è $m_{{\mle}} = \bar{x}$ (qualsiasi sia $v_0$), 2. se il parametro di media $M=m_0\in \R$ è noto (ossia $M=m_0$ è costante rispetto all’informazione a priori), allora la stima di massima verosimiglianza per il parametro di varianza è \[ v_{{\mle}} = \frac 1 n \sum_{i=1}(x_i-m_0)^2.\]

Esempio 5.3 Consideriamo il classico dataset “iris”¹³ che contiene $150$ osservazioni di esemplari diversi da $3$ specie della pianta Iris appunto. Non è necessario importarlo in R perché è sempre disponibile (nella libreria datasets precaricata vi sono alcune raccolte di dati che si usano frequentemente come esempi). Con la funzione head(), possiamo visualizzare alcuni dati per farci una idea generale. L’output del comando head(iris) è presentato sotto.

Tabella 5.1: Le prime 10 righe del dataset Iris.
Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species
5.1	3.5	1.4	0.2	setosa
4.9	3.0	1.4	0.2	setosa
4.7	3.2	1.3	0.2	setosa
4.6	3.1	1.5	0.2	setosa
5.0	3.6	1.4	0.2	setosa
5.4	3.9	1.7	0.4	setosa
4.6	3.4	1.4	0.3	setosa
5.0	3.4	1.5	0.2	setosa
4.4	2.9	1.4	0.2	setosa
4.9	3.1	1.5	0.1	setosa

Per calcolare media e varianza campionaria delle osservazioni della variabile “lunghezza del sepalo” (una parte del fiore), è sufficiente usare le funzioni mean() e var().

mean(iris$Sepal.Length)

## [1] 5.843333

var(iris$Sepal.Length)

## [1] 0.6856935

Per calcolare direttamente la deviazione standard dal campione, ossia la radice quadrata dela varianza campionaria, si può usare sd().

sd(iris$Sepal.Length)

## [1] 0.8280661

Anche prima di calcolare media e varianza, è sempre buona pratica rappresentare graficamente i dati, in questo caso ad esempio tramite un istogramma.

hist(iris$Sepal.Length, breaks = 10, xlab = "Lughezza dei sepali",
  ylab = "Frequenza", main = "", col = miei_colori[1])

Figura 5.16: Istogramma della lunghezza dei sepali osservata nel dataset Iris

5.4.2 Stima bayesiana della media, varianza nota

Vediamo ora l’approccio bayesiano, applicandolo prima alla stima della media, supponendo che la varianza $V = v_0$ sia nota (e quindi costante rispetto all’informazione a priori). Come nel caso della singola osservazione, per avere dei calcoli trattabili analiticamente, conviene supporre che $M$ a priori abbia una densità gaussiana di parametri $\mathcal{N}(m_0, \sigma_m^2)$, e avendo osservato $X=x= (x_i)_{i=1}^n$, il robot calcola poi la densità di $M$ tramite la formula di Bayes \[ \begin{split} p(M=m| X=x) & \propto p(M=m|\Omega ) p(X=x | M=m, V=v_0) \\ & \propto \exp\bra{ -\frac 1 2 \bra{ \frac{ (m-m_0)^2}{\sigma_m^2} + \sum_{i=1}^n \frac { (x_i-m)^2}{v_0} }}.\end{split}\] Si trova quindi, come nel caso $n=1$, una nuova densità gaussiana, di cui con passaggi elementari si ricavano i parametri di media e varianza (che dipendono naturalmente dall’osservazione $X=x$) \[ m_{|X=x} = (1-\alpha) \bar{x} + \alpha m_0, \quad \sigma^2_{m|X=x} = \sigma_m^2\alpha,\] dove abbiamo posto \[\alpha = \frac{1}{1+n\sigma_m^2/v_0 } \in (0,1).\] Come nel caso della singola osservazione, la nuova informazione $X=x$, ossia $X_i=x_i$ per $i=1, \ldots, n$, sposta la fiducia del robot verso la media campionaria $\bar{x}$, ossia la stima di massima verosimiglianza.

Remark. È naturale chiedersi cosa accada per $n$ grande, in particolare se $n >> v_0/\sigma_m^2$. In tal caso, si ha che $\alpha$ tende a $0$ e quindi il parametro di media della $M$ a posteriori tende alla media campionaria $\bar{x}$, ossia la stima di massima verosimiglianza. Per lo stesso motivo, anche che la varianza di $M$ tende a $0$, per cui la distribuzione di $M$ si concentra sempre più intorno al valore $\bar{x}$. Questo è un caso particolare di un teorema molto più generale, noto come legge dei grandi numeri, il quale afferma che la differenza tra la media campionaria di un gran numero di variabili aleatorie indipendenti tra loro (tutte con lo stesso valor medio e varianza) e il valor medio teorico diventa piccola con grandissima probabilità al tendere della numerosità del campione $n$ all’infinito. Ritorneremo su questo fatto nella Sezione 8.2.

5.4.3 Stima bayesiana della varianza, valor medio noto

Supponiamo ora che $M =m_0$ sia costante (rispetto all’informazione nota prima di osservare $X=x$) e introduciamo come nella sezione precedente una densità a priori per $V$ di tipo esponenziale inversa, dove $v_0>0$ è un parametro noto: \[ p( V = v|\Omega) \propto p(1/V = 1/v)\frac{1}{v^2} \propto \exp\bra{ -\frac{v_0}{2v}}\frac{1}{v^2}.\]

Usando la formula di Bayes, \[ \begin{split} p(V=v| X=x) & \propto p(V=v|\Omega ) p(X =x | M=m_0, V=v) \\ & \propto \exp\bra{ - \frac{v_0}{2v} } \frac 1 {v^2} \exp\bra{ - \frac{1}{v} \sum_{i=1}^n (x_i-m_0)^2} \frac{1}{v^{n/2}}\\ &\propto\exp\bra{ - \frac{ \bra{v_0 + \sum_{i=1}^n (x_i-m_0)^2} } {2v}} \frac{1}{v^{(4+n)/2}}.\end{split}\] Come nel caso della singola osservazione, si tratta di una densità continua molto simile a quella a priori, in cui il nuovo parametro (al posto di $v_0$) è \[ v_{|X=x} = v_0 + \sum_{i=1}^n (x_i-m_0)^2,\] ma vi è anche il termine $v^{(4+n)/2}$ a denominatore, che sempre più rilevante al crescere di $n$. Infatti, il punto di massimo della densità a posteriori (che possiamo ottenere passando al logaritmo, e imponendo la derivata nulla) è dato dall’espressione \[ \frac{ v_{|X=x}}{4+n} = \frac{ v_0 + \sum_{i=1}^n (x_i-m_0)^2}{4+n }\] che al crescere di $n$ è asintoticamente equivalente alla stima di massima verosimiglianza \[v_{{\mle}} = \frac{1}{n} \sum_{i=1}^n (x_i-m_0)^2\] (supponendo la media $m_0$ nota).

deltav <- 0.01
v <- seq(deltav, 4, by = deltav)

# densità a priori

v_0 <- 1
densita_esp_inv <- exp(-(v_0/(2 * v)))/v^2
densita_esp_inv <- densita_esp_inv/sum(densita_esp_inv *
  deltav)


# parametro m_0 e osservazione di X
m_0 <- 0
x <- c(1, 2, 1.5, 0.5)
n <- length(x)

# Usiamo la formula trovata per la
# densità a posteriori
v_x <- v_0 + sum((x - m_0)^2)
dens_post_x <- exp(-(v_x/(2 * v)))/v^((4 +
  n)/2)
dens_post_x <- dens_post_x/sum(dens_post_x *
  deltav)

# grafico e legenda

plot(v, densita_esp_inv, type = "l", xlab = "v",
  ylab = "densità", ylim = c(0, 1.5),
  col = miei_colori[1], lwd = 3)
lines(v, dens_post_x, type = "l", col = miei_colori[2],
  lwd = 3)



# Legenda

legend("topright", legend = c("a priori",
  "a posteriori"), fill = miei_colori[1:2])

$grafico della densità a priori per $V$, con parametri $v_0=1$, $m_0=0$, e della densità avendo osservato $ x=(1, 2, 1.5, 0.5)$ (la stima di massima verosimiglianza è in questo caso $v_{\mle} = 1$)$

Figura 5.17: grafico della densità a priori per $V$, con parametri $v_0=1$, $m_0=0$, e della densità avendo osservato $ x=(1, 2, 1.5, 0.5)$ (la stima di massima verosimiglianza è in questo caso $v_{\mle} = 1$)

5.4.4 Esercizi

Esercizio 5.6 Stimare tramite massima verosimiglianza (e mediante opportuni comandi R) i parametri di media e deviazione standard per la lunghezza dei petali nel dataset Iris.

Esercizio 5.7 Supponiamo di essere dei biologi che hanno già classificato la specie Iris “setosa” e determinato che la lunghezza di un petalo è $1.6 \pm 0.2$. Usando i nuovi dati raccolti sulla specie (sono le prime 50 entrate del dataset), rendere più precisa la stima della media tramite l’approccio bayesiano (supporre che $\sigma_m^2 = v^2 = (0.2)^2$ sia nota).