5.3 Stima dei parametri da una singola osservazione

Data una variabile aleatoria $X \in \R^d$ con densità gaussiana $\mathcal{N}(m, \Sigma)$, come stimare i parametri basandosi sull’osservazione di $X$? Come abbiamo visto nella Sezione 3.5, tale problema può essere naturalmente studiato dal punto di vista bayesiano, supponendo quindi che i parametri $m$, $\Sigma$ siano delle variabili aleatorie con opportune densità a priori, e determinando quindi la densità avendo osservato $X=x$ tramite la formula di Bayes. A questo metodo si affianca l’alternativa più diretta, ma meno informativa, di limitarsi ad una stima di massima verosimiglianza.

Per illustrare il metodo, iniziamo con lo studio in questa sezione del caso di una variabile aleatoria $X \in \R$ di cui bisogna stimare i parametri $m$, $\sigma^2$ sulla base della sola osservazione di $X$. Nella sezione successiva, generalizzeremo al caso in cui vi siano più osservazioni e argomenteremo che la stima diverrà allora più precisa.

Supponiamo in questa sezione che il robot abbia modellizato una quantità aleatoria reale $X \in \R$ come una variabile con densità gaussiana $\mathcal{N}(m,\sigma^2)$. I parametri non sono noti a priori, e vengono stimati sulla base di una osservazione $X=x$. Pertanto seguendo il metodo bayesiano il robot introduce le rispettive variabili aleatorie $M$ per la media e $V$ per la varianza (questa a valori positivi).

La verosimiglianza, ossia la densità di $X$ supponendo note la media $M=m$ e la varianza $V = v $ (usiamo la lettera $v$ al posto di $\sigma^2$, per alleggerire la notazione) è quindi \[ L( m,v; x) = p(X =x | M=m V=v) = \exp\bra{ - \frac 1 {2 v} (x-m)^2 }\frac{1}{\sqrt{ 2 \pi v}}.\] È importante specificare tutta la densità (anche se il termine $\sqrt{2 \pi}$ non sarebbe rilevante), perché interessa la dipendenza dai parametri $m$ e $v$. Possiamo rappresentare graficamente la funzione dei due parametri $(m,v)$ con una heatmap.

deltam <- 0.1
deltav <- 0.05

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

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

# osservazione

x <- 1


# creiamo una matrice con i valori
# della verosimiglianza

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

for (i in 1:N_m) {
  for (j in 1:N_v) {
    L[i, j] <- exp(-(m[i] - x)^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.10: verosimiglianza per $M$ e $V$ avendo osservato $X=1$.

5.3.1 Stima di massima verosimiglianza

Dal grafico sopra è evidente che la verosimiglianza $L$ è massima per $m=x$ (in tal caso era $X=1$) e $v=0$, un fatto che ora giustifichiamo analiticamente, determinando la stima di massima verosimiglianza per i parametri. Passando al logaritmo e moltiplicando per $-2$, invece di massimizzare $L$ basta minimizzare la funzione \[ (m,v) \mapsto - 2 \log L(m,v; x) = \frac 1 { v} (x-m)^2 + \log( 2 \pi v).\] Derivando rispetto a $m$ e imponendo che la derivata si annulli si ottiene, per ogni $v>0$, \[ \frac 2 { v} (x-m) =0 \quad \text{da cui $m=x$,} \] mentre se deriviamo rispetto a $v$, tenendo fisso $m$, si trova \[ -\frac 1 { v^2} (x-m)^2 + \frac{1}{v} = 0, \quad \text{da cui $v = (x-m)^2$. }\] dovendo massimizzare congiuntamente si trova quindi la coppia $m_{{\mle}} =x$, $v_{{\mle}} = 0$ (anche se, ad essere rigorosi, per $v=0$ non è ben definita la verosimiglianza).

Con calcoli analoghi a quanto fatto sopra si ottiene anche che 1. se il parametro di varianza $V=v_0 >0$ è noto (ossia $V=v_0$ è costante rispetto all’informazione a priori), allora la stima di massima verosimiglianza per il parametro di media è $m_{\mle} = x$ (qualsiasi sia $v_0$). 2. se il parametro di media $M=m_0 \in \R$ è noto (ossia $M=m_0$ è costante a priori), allora la stima di massima verosimiglianza per il parametro di varianza è $v_{{\mle}} = (x-m_0)^2$.

5.3.2 Stima bayesiana per la media, varianza nota

Per applicare il metodo bayesiano, il robot deve proporre delle densità a priori per media $M$ e varianza $V$ che rappresentano l’informazione di cui dispone prima di osservare $X$, e applicare la formula di Bayes per ottenere le densità condizionate a $X=x$. Ovviamente le densità a priori dipendono dalla natura dell’informazione iniziale di cui dispone. La semplice rete bayesiana che rappresenta il problema è rappresentata in figura.

Figura 5.11: Rete bayesiana tra le variabili $M$, $V$ e $X$

Affrontiamo prima due casi estremi, in cui i risultati si possono ottenere analiticamente, se si scelgono delle densità a priori particolari: consideriamo prima il caso in cui $V = v_0$ sia nota (e quindi sia una variabile costante) e successivamente quello in cui $M = m_0$ sia invece nota. In questi casi si può rimuovere il nodo corrispondente alle variabili costanti dalla rete bayesiana.

Supponendo che $V=v_0$ sia costante, i calcoli risultano particolarmente semplici se si suppone che $M$ sia una variabile gaussiana con parametri (noti) $\mathcal{N}(m_0, \sigma_m^2)$, ossia \[ p(M= m| \Omega) \propto \exp\bra{-\frac 1 {2\sigma_m^2} (m-m_0)^2}.\] Questa densità codifica una informazione nota al robot riguardante il parametro di media $m$: esso è localizzato intorno al parametro $m_0$, con una dispersione (deviazione standard) $\sigma_m$. Per rappresentare maggiore incertezza su $M$ è sufficiente far crescere $\sigma_m$, e nel limite $\sigma_m \to \infty$ vedremo che si ricade nella stima di massima verosimiglianza.

Avendo osservato $X=x$, dalla formula di Bayes segue che la densità per $M$ è \[ \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} + \frac { (x-m)^2}{v_0} }}.\end{split}\] Si tratta evidentemente di una nuova densità gaussiana essendo l’esponenziale di un polinomio di secondo grado nella variabile $m$. Con semplici passaggi algebrici si ricavano i parametri di media e varianza, che dipendono naturalmente dall’osservazione $x$, \[ m_{|X=x} = (1-\alpha) x + \alpha m_0, \quad \sigma^2_{m|X=x} = \sigma_m^2\alpha,\] dove abbiamo posto, per semplicità, \[\alpha = \frac{1}{1+\sigma_m^2/v_0 } \in (0,1).\] Osserviamo quindi che la nuova informazione $X=x$ “sposta” la fiducia del robot verso il valore osservato $x$, ma non completamente (come invece accade con la massima verosimiglianza). Molto dipende dal valore di $\alpha$, ossia dal rapporto tra le varianze $\sigma_m^2/v_0$ (entrambe note a priori). I casi limite sono particolarmente rilevanti: 1. se $v_0<< \sigma_m^2$, allora $\alpha \sim 0$ e ci si avvicina alla stima di massima verosimiglianza (il che non stupisce perché $\sigma_m$ grande significa che la densità a priori per $M$ era poco informativa) 2. se $v_0>> \sigma_m^2$, allora $\alpha \sim 1$ e la media $m_{|X=x}$ rimane praticamente invariata, essenzialmente perché l’informazione iniziale era molto precisa su $M$, e una singola osservazione non la modifica molto.

5.3.3 Stima bayesiana per la varianza, media nota

Supponiamo ora che $M =m_0$ sia costante rispetto all’informazione iniziale disponibile al robot e introduciamo una densità a priori per la varianza $V$. Uno dei problemi principali è che $V$ deve essere non-negativa, ma vi sono diverse scelte comode per ottenere risultati analitici. Una di queste è una densità continua del tipo esponenziale inversa, ossia $1/V$ è esponenziale con parametro $\lambda = v_0/2$ (dove $v_0$ è un parametro che riteniamo noto mentre il coefficiente $1/2$ è solo per semplificare i calcoli). Dalla formula di cambio di variabile si trova che \[ 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}.\]

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

v_0 <- 1

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

plot(v, densita_esp_inv, type = "l", xlab = "v",
  ylab = "densità", col = miei_colori[1],
  lwd = 3)

Figura 5.12: grafico della densità esponenziale inversa con $v_0=1$

Con questa scelta, avendo osservato $X=x$, la formula di Bayes implica che la densità per $V$ è \[ \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 {2v} (x-m_0)^2 } \frac{1}{\sqrt{v}}\\ &\propto\exp\bra{ - \frac{ \bra{v_0 + (x-m_0)^2} } {2v}} \frac{1}{v^{5/2}}\end{split}\] 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 + (x-m_0)^2.\] A prima vista potrebbe quindi sembrare che la varianza cresca sempre, ma bisogna anche notare che il termine $v^2$ al denominatore è sostituito con $v^{5/2}$. Questo ha l’effetto di trasformare la densità circa in modo che, se $(x-m_0)^2$ è minore di $v_0$, allora la densità di $V$ si concentra verso valori più vicini a $0$, viceversa se $(x-m_0)^2$ è maggiore, allora la densità di $V$ si concentra verso valori maggiori.

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_1 <- 0
x_2 <- 2

## Usiamo direttamente la formula
## trovata sopra
v_x_1 <- v_0 + (x_1 - m_0)^2
dens_post_x_1 <- exp(-(v_x_1/(2 * v)))/v^{
  5/2
}
dens_post_x_1 <- dens_post_x_1/sum(dens_post_x_1 *
  deltav)

v_x_2 <- v_0 + (x_2 - m_0)^2
dens_post_x_2 <- exp(-(v_x_2/(2 * v)))/v^{
  5/2
}
dens_post_x_2 <- dens_post_x_2/sum(dens_post_x_2 *
  deltav)

# grafico e legenda

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

# Legenda

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

Figura 5.13: grafico della densità a priori per $V$, con parametri $v_0=1$, $m_0=0$, e della densità avendo osservato $X=0$ (in rosso) oppure $X=2$ (in blu)

Lo studio del caso generale dal punto di vista bayesiano, ossia della variabile congiunta $(M,V)$ (senza supporre che una sia costante nota a priori), è più complicato analiticamente. Tuttavia possiamo renderci conto di come l’informazione a priori sia rilevante osservando un risultato numerico ottenuto assumendo che $M$ e $V$ siano a priori indipendenti, la prima uniformemente distribuita su $[2,4]$, la seconda sull’intervallo $[1, 3]$, e si osserva $X=1$. La stima di massima verosimiglianza (calcolata in precedenza su tutti i valori di $m$ e $v$ sulla stessa osservazione $X=1$) era $m_{{\mle}}=1$, $v_{{\mle}} = 0$, ora l’informazione a priori induce una densità congiunta (le cui marginali non sono indipendenti) con un punto di massimo in $m=2$, $v=3$. Possiamo intepretare parzialmente questo fatto notando che essendo media a priori tra $2$ e $4$, l’osservazione $X=1$ sposta la media verso il valore più basso possibile $m=2$.

deltam <- 0.1
deltav <- 0.05

m <- seq(2, 4, by = deltam)
v <- seq(1, 3, by = deltav)

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

# osservazione

x <- 0


# creiamo una matrice con i valori
# della verosimiglianza

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

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

posteriori <- L/sum(L * deltam * deltav)


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

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