3.5 Formula di Bayes per variabili aleatorie

Abbiamo visto che la densità di una variabile congiunta non è esprimibile in termini delle due densità marginali (rispetto all’informazione nota \(I\)). Un modo per aggirare questo problema è fornito dalla regola del prodotto, che garantisce (evitiamo di scrivere \(I\) per semplicità) \[ P( (X,Y) = (x,y)) = P(X=x \text{ e } Y = y ) = P(X=x) P(Y=y | X=x).\] Possiamo leggere questa identità nel seguente modo: nel caso di variabili discrete, la densità della variabile congiunta è determinata da 1. La densità della marginale \(X\) 2. La densità della marginale \(Y\), ma condizionata all’informazione \(X=x\) per ciascun \(x \in E\) (che si abbrevia dicendo semplicemente condizionata ad \(X\)).

In pratica in molti casi la legge congiunta è proprio definita tramite queste due quantità (la densità di una marginale e la densità dell’altra marginale condizionata alla prima).

Esempio 3.17 Nel caso delle due estrazioni dall’urna la densità congiunta della prima e della seconda estrazione è definita nell’ordine naturale, partendo dalla prima estrazione e poi specificando la densità della seconda condizionata alla prima (e distinguendo quindi tra estrazioni con e senza rimpiazzo). Ovviamente nulla vieta di considerare anche altre densità condizionate, ad esempio invece di aggiungere la stessa pallina si potrebbe sostituire la pallina estratta con una dell’altro colore. In tal caso, ponendo come nella sezione precedente \(X\) la variabile indicatrice della prima estrazione (\(1\) se è rossa), \(Y\) invece relativa alla seconda estrazione, si trova \[ P(Y =1 | X=1) = \frac{R-1}{N}, \quad P(Y =1 | X=0) = \frac{R+1}{N},\] che permette di definire una ulteriore densità congiunta per la variabile \((X,Y)\) (ovviamente densità congiunte diverse servono a descrivere situazioni diverse).

La formula di Bayes si può quindi riscrivere in termini di variabili aleatorie, ottenendo \[ P(X=x | Y =y ) = P(X=x) \cdot \frac{ P(Y=y | X=x)}{P(Y=y)} \propto P(X=x)L(X=x; Y=y),\] con la stessa inteprertazione del caso di sistemi di alternative (è in effetti solamente un cambio di notazione). Ricordiamo che il rapporto \[ \frac{ P(Y=y | X=x)}{P(Y=y)} = \frac{ P(X=x | Y=y)}{P(X=x)}\] non cambia se si scambiano i ruoli di \(X=x\) ed \(Y=y\), e indica quando l’osservazione di uno dei due eventi aumenti (se maggiore di \(1\)), diminuisca (se minore di \(1\)) il grado di fiducia nell’altro evento. Come nel caso delle alternative, in molti casi si è semplicemente interessati al valore più probabile della \(X\), avendo osservato \(Y=y\). Si definisce pertanto la stima di massimo a posteriori per la \(X\), avendo osservato \(y=y\), come il valore (o i valori) \[ x_{\map} \in \arg \max \cur{ P(X=x) L(X=x; Y=y) : x \in E}\] e ricordiamo che nel caso di \(X\) con densità uniforme discreta su \(E\) (rispetto all’informazione iniziale) il problema si riduce alla stima di massima verosimiglianza \[ x_{\mle} \in \arg \max \cur{L(X=x; Y=y) : x \in E}.\] In molte situazioni in cui \(E\) sia infinito, per estensione dal caso finito si considera il problema sopra (immaginando una distribuzione a priori uniforme su \(E\), che rigorosamente non esiste, ma è comunque utile). Vediamo un esempio

Esempio 3.18 Si consideri la seguente situazione: si informa il robot che un’urna contiene metà palline rosse e metà palline blu, da cui si effettuano un certo numero di estrazioni con rimpiazzo. Il robot tuttavia non viene informato del numero esatto delle estrazioni, ma si comunica al robot solamente che il numero di palline rosse estratte è \(10\). Possiamo stimare il numero di estrazione effettuate? intuitivamente, la risposta dovrebbe essere \(20\), ma vediamo come il robot può ragionare.

Definiamo una variabile aleatoria \(M\) che indica il numero di estrazioni effettuate. Non sapendo nulla su \(M\) (prima di ricevere l’informazione che \(10\) rosse sono state estratte), il robot suppone che sia uniformemente distribuita sui valori \(\cur{0,1, \ldots, \bar{m}}\), per un parametro \(\bar{m}\) molto grande (idealmente infinito), \[ P(M=m|\Omega) = \frac {1}{\bar{m}}\] per \(m \le \bar{m}\). Nel definire la densità del numero di rosse estratte sapendo \(M = m\), il robot utilizza la densità binomiale di parametri \((m,1/2)\) (perché sa che metà palline sono rosse e metà blu). Pertanto, posta \(N_R\) la variabile che indica il numero di palline rosse estratte, si ha \[ P( N_R = k | M=m) = {m\choose k} \bra{ \frac{1}{2}}^k \bra{\frac 1 2 }^{m-k} = {m \choose k} \frac {1}{2^m}.\] La formula di Bayes permette di ottenere la densità di \(M\) avendo osservato \(N_R = 10\), per \(m \le \bar{m}\), \[ P( M=m| N_R = k) = P(M=m|\Omega)\cdot \frac{P(N_R = 10| M=m)}{P(N_R = 10|\Omega)} = \frac{1}{\bar{m}} {m \choose k} \frac {1}{2^m} \cdot \frac{1}{P(N_R = 10|\Omega)}.\] Un’espressione per \(P(R_R|\Omega)\) si ottiene imponendo che l’ultimo termine sia una densità discreta (oppure tramite la formula di disintegrazione) \[P(N_R|\Omega) = \sum_{m=0}^{\bar{m}} \frac{1}{\bar{m}} {m \choose k} \frac {1}{2^m}.\] Ora non vale più la pena di proseguire con i calcoli teorici, piuttosto è meglio plottare la densità e studiare come dipenda dal parametro \(\bar{m}\).

# introduciamo due possibili parametri
# (si consiglia di ripetere con altri
# valori)

possibili_bar_m <- c(25, 50)

# Inizializziamo una matrice che
# conterrà le densità al variare del
# parametro bar_m (utile il barplot).

dens_matrice_M_NR_10 <- matrix(0, nrow = 2,
  ncol = 50)


# dovendo ripetere gli stessi calcoli
# al variare di bar_m, utilizziamo un
# ciclo for

for (iter in 1:2) {
  bar_m <- possibili_bar_m[iter]
  m <- 0:bar_m

  # scriviamo la verosimiglianza

  likelihood <- dbinom(10, m, 1/2)

  # la densità discreta di M avendo
  # osservato N_R = 10 si ottiene
  # moltiplicando la densità a priori
  # (che è costante uguale a 1/bar_m,)
  # e normalizzando in modo che sommi
  # ad 1 -- quindi non serve neppure
  # moltiplicare per 1/bar_m, ma lo
  # facciamo per chiarezza

  dens_M_NR_10 <- likelihood/bar_m
  dens_M_NR_10 <- dens_M_NR_10/sum(dens_M_NR_10)
  dens_matrice_M_NR_10[iter, ] <- dens_M_NR_10[1:50]
}

colori <- miei_colori[1:2]

barplot(dens_matrice_M_NR_10, beside = TRUE,
  col = colori, names.arg = as.character(0:49),
  xlab = "numero estrazioni", ylab = "densità discreta")


# aggiungiamo una legenda

legend("topright", legend = c("bar_m = 25",
  "bar_m = 50"), fill = colori)

Figura 3.9: grafico della densità discreta di \(M\) sapendo \(N_R=10\), al variare del parametro \(\bar{m}\)

Osserviamo che se \(\bar{m}\) è piccolo, la densità di \(M\) varia abbastanza, tuttavia al crescere di \(\bar{m}\) si stabilizza in un opportuno profilo. Pertanto il risultato non dipende essenzialmente da \(\bar{m}\) (che si può anche mandare all’infinito, per essere formali).

Inoltre, la stima di massima verosimiglianza \(m_{\mle}\) si può ottenere dai calcoli svolti sopra.

# la funzione which.max restituisce
# l'indice corrispondente al punto di
# massimo di un vettore

indice_max <- which.max(dens_M_NR_10)

# di conseguenza m_max si ottiene dalla
# componente di indice_max del vettore
# m

(m_max <- m[indice_max])

## [1] 19

# la risposta potrebbe stupirci, ma
# confrontiamo le probabilità a
# posteriori

dens_M_NR_10[indice_max]  # per la probabilità che M=19

## [1] 0.08809918

dens_M_NR_10[indice_max + 1]  # per la probabilità che M=20

## [1] 0.08809918

Questa costruzione, che è essenzialmente la regola del prodotto e la conseguente formula di Bayes, si estende in diversi modi. Il primo è di considerare variabili continue. In tal caso, per definire la densità continua della congiunta, si deve mostrare una versione della regola del prodotto per le densità, precisamente, per \(X \in \R^d\), \(Y \in \R^k\), \[ p((X,Y)=(x,y) ) = p(X=x, Y=y) = p(X=x) p(Y=y | X=x),\] dove il termine \(p(Y=y|X=x)\) è detto anche densità condizionale di \(Y\) rispetto ad \(X\) ed è la densità della variabile \(Y\) sapendo, oltre all’informazione \(I\) qui non scritta, anche che \(X\) assume il valore \(x\). Il modo più semplice per pensare a questa formula è che sia una definizione della densità della congiunta partendo dalla densità continua di \(X\) e dalla densità di \(Y\) sapendo \(X=x\). Si può anche leggere al contrario, ossia conoscendo la densità della variabile congiunta, si ricava una versione della fomula di Kolmogorov per densità continue \[ p(Y =y | X=x) = \frac{ p(X=x, Y=y)}{p(X=x)}.\]

La costruzione è perfettamente simmetrica, e invertendone l’ordine si ottiene la formula di Bayes per densità continue: \[ p(X=x | Y =y ) = p(X=x) \cdot \frac{ p(Y=y | X=x)}{p(Y=y)} \propto p(X=x) L(X=x; Y=y),\] con la stessa intepretazione che nel caso discreto (eccetto che stiamo trattando densità continue e non probabilità). Per determinare il denominatore, possiamo sempre ricordare che l’integrale della densità a sinistra (come funzione di \(x\)) deve valere \(1\), da cui si trova una formula di disintegrazione per densità continue: \[ p(Y=y) = \int_{\R^d} p(Y=y | X=x) p(X=x)d x = \int_{\R^d} L(X=x; Y=y) p(X=x) dx.\] Analogamente a quanto accadeva nel caso discreto, il rapporto \[ \frac{ p(Y=y | X=x)}{p(Y=y)} = \frac{ p(X=x | Y=y)}{p(X=x)} \] non cambia se si scambiano gli eventi \(Y=y\), \(X=x\). Come nel caso discreto, per determinare l’alternativa \(\cur{X=x}\) più probabile, avendo osservato \(\cur{Y=y}\), è sufficiente trovare il (o i) punti di massimo della funzione \[ x \mapsto p(X=x) L(X=x; Y=y),\] determinando così la stima del massimo a posteriori per \(X\) avendo osservato \(Y=y\).

Remark (massima verosimiglianza nel caso continuo). Nel caso in cui \(X\) sia una variabile uniforme (ad esempio, su un intervallo \([a,b]\)), come nel caso discreto il problema di determinare i punti di massimo per la distribuzione condizionata all’osservazione di \(Y\), si riduce a \[ x_{\mle} \in \operatorname{arg}\cur{ \max p(Y=y | X=x): x \in [a,b]}.\] Tale \(x_{\mle}\) è il caso continuo della stima di massima verosimiglianza della variabile \(X\), avendo osservato \(Y=y\). Spesso, si massimizza direttamente al variare di \(x \in \R\) o su un intervallo illimitato, come una semiretta (anche se propriamente non esiste una densità continua uniforme su intervalli illimitati).

Esempio 3.19 Vediamo un esempio di applicazione della formula di Bayes e della stima massima verosimiglianza nel caso continuo. Si vuole modellizare la durata della carica di un dispositivo (ad esempio, uno smartphone, o un drone) tramite una variabile aleatoria \(T\) avente densità continua esponenziale di un certo parametro \(\lambda\). Questo semplice modello contiene il solo parametro \(\lambda\), appunto, che inizialmente possiamo supporre uniformemente distribuito nei valori \([0,1]\) (ad esempio, misurando in giorni, e ricordando che tanto più piccolo è \(\lambda\), maggiori sono i valori assunti da \(T\)). Il robot considera quindi una variabile \(\Lambda\) uniforme continua su \([0,1]\) (rispetto alla informazione iniziale \(\Omega\)). Condizionata a \(\Lambda= \lambda\), \(T\) ha una densità esponenziale di parametro \(\lambda\), \[p(T = t | \Lambda = \lambda) = \lambda e^{-\lambda t } \quad \text{per $t \ge 0$.}\] Avendo osservato che (per un dispositivo) la sua durata è \(T=10\), possiamo ottenere la densità di \(\Lambda\) aggiornata a questa informazione: per \(\lambda \in [0,1]\), \[ p(\Lambda = \lambda| T = 10) = p(\Lambda = \lambda|\Omega) p(T=10 | \Lambda = \lambda) \cdot \frac{1}{p(T=10|\Omega)} = p(T=10 | \Lambda = \lambda) \cdot \frac{1}{p(T=10|\Omega)},\] ricordando l’ipotesi che \(\Lambda\) sia uniforme (rispetto a \(\Omega\)). Per determinare il denominatore esplicitamente, basta imporre che il membro a destra sia una densità continua (rispetto a \(\lambda\)) e quindi integrare ad \(1\). Si trova la formula \[ p(T=10|\Omega) = \int_0^1 \lambda e^{- 10 \lambda } d \lambda.\] Non vale la pena di cercare una espressione in termini di funzioni elementari, ma è più utile tracciarne un grafico approssimato.

delta_lambda <- 0.01
lambda <- seq(0, 1, by = delta_lambda)

# scriviamo la verosimiglianza, ossia
# P(T=10|Lambda=lambda)

likelihood <- lambda * exp(-lambda * 10)

# otteniamo la distribuzione di Lambda
# moltiplicando la densità a priori
# (che vale 1 in questo caso) per la
# verosimiglianza e normalizzando
# dividendo per l'integrale
# approssimato tramite somme di Riemann

dens_Lambda_T_10 <- likelihood
dens_Lambda_T_10 <- dens_Lambda_T_10/(sum(dens_Lambda_T_10) *
  delta_lambda)

# plottiamo sia la distribuzione a
# priori (uniforme) che quella
# condizionata

plot(NULL, xlim = c(0, 1), ylim = c(0, 4),
  xlab = "lambda", ylab = "densità")

lines(lambda, dunif(lambda), type = "l",
  lwd = 3, col = miei_colori[1])

lines(lambda, dens_Lambda_T_10, type = "l",
  lwd = 3, col = miei_colori[2])



# aggiungiamo una legenda

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

Figura 3.10: densità per \(\Lambda\): a priori e avendo osservato \(T=10\)

Se si è interessati alla stima di massima verosimiglianza, si determinare \[ \lambda_{\mle} \in \arg \max \cur{ \lambda e^{- 10 \lambda } : \lambda \in [0,1]},\] e in questo caso il problema ha una semplice soluzione analitica. Basta derivare e imporre che la derivata sia nulla, trovando \[ e^{- 10 \lambda } - 10 \lambda e^{-10 \lambda} = 0 \quad \text{ossia} \quad \lambda = 1/10.\] (andrebbe anche verificato che il massimo non sia raggiunto agli estremi dell’intervallo, ma dal grafico sopra è evidente). Per modelli più complicati, si deve ricorrere a metodi numerici per determinare la stima di massima verosimiglianza. Possiamo anche accontentarci del risultato approssimato ottenuto dai valori plottati sopra.

indice_max <- which.max(dens_Lambda_T_10)

# per ottenere la stima di massima
# verosimiglianza basta ottere la
# componente dal vettore lambda

lambda[indice_max]

## [1] 0.1

consideri una variabile aleatoria \(\Lambda\) a valori in \([0, 10]\) uniforme (rispetto ad una informazione iniziale \(\Omega\))

Vale la pena di menzionare anche il caso misto, in cui \(X \in E\) ha densità discreta mentre \(Y\), condizionata ad \(\cur{X=x}\), è continua, per ogni \(x \in E\). In tal caso il problema è che in generale la variabile congiunta non ha densità né continua né discreta, tuttavia la formula di Bayes rimane valida, nella forma \[ P(X = x | Y=y ) = P(X =x)\cdot \frac{ p(Y= y| X=x)}{p(Y=y)} \propto P(X=x) L(X=x; Y=y)\] e al solito imponendo che sia una densità discreta e quindi sommi ad \(1\) si trova per il denominatore l’espressione \[ p(Y=y) = \sum_{x \in E} p(Y= y| X=x)P(X =x).\] Vale anche la formula nel caso simmetrico (ossia \(X\in \R^d\) è continua mentre \(Y\), condizionata ad \(\cur{X=x}\) è discreta): \[ p(X=x | Y=y) = p(X=x) \cdot \frac{ P(Y= y| X=x)}{P(Y=y)} \propto p(X=x) L(X=x; Y=y)\] con \[ P(Y=y) = \int_{\R^d} P(Y= y| X=x)p(X =x) dx.\] Anche in questo caso è utile studiare la stima di massimo a posteriori o di massima verosimiglianza per \(X\) avendo osservato \(Y=y\) (non ripetiamo la definizione, che è praticamente la stessa, con le dovute modifiche).

Esempio 3.20 Per fare un esempio di questo caso “misto”, si supponga che il robot sia stato informato che un’urna contiene palline rosse oppure blu, ma non del numero totale né della frazione di palline rosse sul totale. Successivamente viene informato che, avendo effettuato \(10\) estrazioni con rimpiazzo, sono state osservate \(3\) palline rosse. Come stimare la frazione di palline rosse sul totale? Intuitivamente capiamo che la risposta è \(3/10\), ma vediamo come ragiona il robot.

Il robot introduce una variabile aleatoria \(F_R \in [0,1]\) per indicare la frazione di palline rosse. Rispetto alla informazione iniziale (prima di sapere delle \(10\) estrazioni), suppone che abbia densità uniforme continua. Introduce poi la variabile \(N_R\) che indica il numero di palline rosse estratte nelle \(10\) estrazioni. Sapendo \(F_R = r\), \(N_R\) ha densità discreta Binomiale di parametri \(n =10\) (numero di estrazioni) e \(p = r\) (probabilità di estrarre una rossa). Pertanto, per \(k \in \cur{0, \ldots, 10}\), \[ P(N_R = k | F_R = r) = { 10 \choose k} r^k (1-r)^{10-k}.\] Avendo osservato \(N_R =3\), la formula di Bayes (nel caso “misto” discreto e continuo) diventa, per \(r \in [0,1]\), \[ p( F_R = r | N_R = 3) = p(F_R = r|\Omega) \frac{ P(N_R=3 | F_R = r)}{P(N_R = 3|\Omega)} = { 10 \choose 3} r^3 (1-r)^{7} \cdot \frac{1}{P(N_R = 3 | \Omega)},\] icordando che \(p(F_R = r| \Omega) = 1\) (la densità a priori è uniforme). Troviamo anche una espressione per il denominatore – anche se abbiamo capito che per molti aspetti non serve – imponendo che il membro a destra sia una densità continua \[ P(N_R = 3 | \Omega) = \int_0^1 { 10 \choose 3} r^3 (1-r)^{7} d r.\] Si tratta di integrare un polinomio, quindi con un po’ di lavoro (oppure l’uso di un integratore simbolico) si può anche ottenere una formula esplicita. Tuttavia possiamo anche limitarci allo studio numerico.

delta_fr <- 0.01
fr <- seq(0, 1, by = delta_fr)

# scriviamo la verosimiglianza

likelihood <- dbinom(3, 10, fr)

# otteniamo la distribuzione di F_R
# moltiplicando la densità a priori
# (che vale 1 in questo caso) per la
# verosimiglianza e normalizzando
# dividendo per l'integrale
# approssimato tramite somme di Riemann

dens_FR_NR_3 <- likelihood
dens_FR_NR_3 <- dens_FR_NR_3/(sum(dens_FR_NR_3) *
  delta_fr)

# plottiamo sia la distribuzione a
# priori (uniforme) che quella
# condizionata

plot(NULL, xlim = c(0, 1), ylim = c(0, 3),
  xlab = "frazione di palline rosse", ylab = "densità")

lines(fr, dunif(fr), type = "l", col = miei_colori[1],
  lwd = 3)


lines(fr, dens_FR_NR_3, type = "l", col = miei_colori[2],
  lwd = 3)


# aggiungiamo una legenda

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

Figura 3.11: densità continua della frazione di palline rosse a priori e avendo osservato in \(10\) estrazioni con rimpiazzo \(3\) palline rosse

Osserviamo in particolare come la densità si sia accumulata intorno al punto di massima verosimiglianza, che possiamo calcolare numericamente al solito modo:

fr[which.max(dens_FR_NR_3)]

## [1] 0.3

che conferma la nostra intuizione. Osserviamo anche che la densità vale \(0\) nei casi estremi \(F_R=0\) oppure \(F_R = 1\): questo riflette il fatto che il robot abbia osservato almeno una pallina rossa e almeno una blu dentro l’urna.

3.5.1 Esercizi

Esercizio 3.10 Ripetere gli esempi con densità non uniformi (ad esempio partendo dalle densità a posteriori ottenute)