2.7 Statistica bayesiana

Dato un sistema di alternative $(A_i)_{i=1}^n$ (rispetto ad una informazione $I$ che qui sottointendiamo per alleggerire la notazione) e una qualsiasi affermazione $B$, possiamo applicare la formula di Bayes a ciascuna alternativa e ottenere \[\begin{equation} P(A_i | B ) = P(A_i) P(B | A_i ) \cdot \frac{1}{P(B)}, \tag{2.4} \end{equation}\] dove abbiamo messo in evidenza il denominatore $P(B)$ per due ragioni. Da un lato, possiamo calcolarlo sempre tramite lo stesso sistema di alternative e la formula di disintegrazione, ossia \[ P(B) = \sum_{i=1}^n P(A_i) P(B | A_i ).\] Questa identità permette di ottenere una formula esplicita per la densità discreta del sistema di alternative, rispetto alla nuova informazione che include $B$ (oltre ad $I$).

D’altra parte, non è neppure necessario usare la formula di disintegrazione per calcolare $P(B)$ perché, osservando la formula (2.4), il membro a sinistra definisce la densità discreta associata al sistema di alternative $(A_i)_{i=1}^n$ rispetto alla nuova informazione che include $B$ oltre ad $I$. Perciò, possiamo sempre indicarlo a meno di una costante moltiplicativa comune (in questo caso appunto $P(B)$): \[\begin{equation} P(A_i | B ) \propto P(A_i) P(B | A_i ) = P(A_i)L(A_i;B), \tag{2.5} \end{equation}\] avendo usato la verosimiglianza nella seconda espressione. La probabilità di ciascuna alternativa rispetto alla nuova informazione $B$ è dunque ottenuta moltiplicando la probabilità iniziale per la verosimiglianza associata. Ovviamente imponendo che la somma delle probabilità sia $1$ si trova la stessa formula per $P(B)$, o equivalentemente usando la notazione per la verosimiglianza \[ P(B) = \sum_{i=1}^n P(A_i) L(A_i; B).\]

È quindi utile riconoscere che il denominatore $P(B)$ ha un ruolo secondario, e volendo si può tenere a mente solo la versione (2.5) della formula. Un altra ragione per non concentrarsi troppo sul denominatore comune $P(B)$ è perché spesso non si richiede di determinare quanto sia probabile ciascuna alternativa $A_i$, rispetto alla nuova informazione $B$, ma semplicemente ci si limita a determinare quale sia l’alternativa più probabile della nuova densità discreta $(P(A_i|B))_{i=1}^n$. Abbiamo introdotto in precedenza il termine statistico moda, ma nel contesto della formula di Bayes è detta stima del massimo a posteriori (in inglese maximum a posteriori probability estimate, MAP). Per calcolare tale $i_{\map} \in \cur{1, \ldots n}$, possiamo quindi tralasciare costanti moltiplicative (positive) e allora vale \[ i_{\map} \in \operatorname{arg} \max\cur{ P(A_i) P(B | A_i )\, : \, i \in \cur{1, \ldots, n} }, \] oppure, usando la notazione della verosimiglianza, \[ i_{\map} \in \operatorname{arg} \max\cur{ P(A_i) L( A_i ;B )\, : \, i \in \cur{1, \ldots, n} }, \] dove ricordiamo che il simbolo di appartenenza ($\in$) è perché potrebbero esserci più indici che raggiungono il massimo.

Remark (massima verosimiglianza). Nel caso speciale in cui la densità discreta a priori sia uniforme, ossia $P(A_i) \propto 1$, si può ridurre ulteriormente il problema di massimizzare $P(A_i) P(B | A_i )$ alla semplice massimizzazione della verosimiglianza $P(B| A_i)=L(A_i;B)$. Tale approccio, molto utilizzato in pratica perché spesso non si è in grado di precisare densità a priori non uniformi, è detto stima di massima verosimiglianza (in inglese maximum likelihood estimation, MLE): dato un sistema di alternative $(A_i)_{i=1}^n$ e una affermazione $B$, avendo osservato $B$ si determina $i_{\mle}$ tale che \[ L(A_{i_{\mle}}; B) = \max_{i=1, \ldots, n}L(A_i; B).\] Questo metodo si può presentare anche senza la formula di Bayes, ma alla luce di questa ne abbiamo una giustificazione in termini del calcolo delle probabilità, e anche una sua estensione al caso in cui le probabilità a priori non siano uniformi.

Remark (passaggio al logaritmo). Ricordando che è possibile anche passare al logaritmo per determinare la moda, si può quindi determinare \[ i_{\map} \in \operatorname{arg} \max\cur{ \log( P(A_i)) + \log( L( A_i ;B ))\, : \, i \in \cur{1, \ldots, n} }, \] che nel caso di densità a priori uniforme si ridure a massimizzare la $\log$-verosimiglianza ($\log$-likelihood in inglese) $\log( L(A_i; B))$ sulle possibili alternative $(A_i)_{i\in I}$.

Vediamo in pratica un esempio dal modello delle estrazioni senza rimpiazzo dall’urna. Questa volta, il robot non è inizialmente informato sul numero di palline rosse contenute, ma solamente sul numero totale $N=3$. Dopo un’estrazione, si aggiorna il robot con l’informazione $R^1$, ossia che una pallina rossa è stata estratta. Cosa può dedurre circa il contenuto dell’urna? Usando il calcolo delle probabilità sviluppato finora e in particolare la formula di Bayes, diamo una risposta (questo modo di procedere è detto anche statistica Bayesiana).

Introduciamo un sistema di alternative relativo al contenuto dell’urna: \[ A_i = \text{l'urna contiene $R=i$ palline rosse}, \] per $i \in \cur{0,1, 2, 3}$. Prima di poter applicare la formula di Bayes, il robot deve stabilire la densità discreta associata al sistema rispetto all’informazione inizale (la probabilità a priori). È chiaro che potrebbero esserci molteplici scelte, ma non avendo ragioni per favorire una alternativa rispetto all’altra, questa è proprio una situazione in cui usare la densità discreta uniforme. Pertanto, indicando con $\Omega$ l’informazione iniziale, per ciascun $i\in \cur{0,1,2,3}$, il robot pone \[ P(A_i | \Omega) = \frac 1 4.\] Successivamente, il robot calcola la probabilità condizionata (il problema conoscendo il numero di palline rosse e blu è già stato affrontato nella sezione precedente) \[ P( R^1 | A_i) = \frac{i}{3},\] e conclude, usando la formula di Bayes \[ P(A_i | R^1 ) = \frac{1}{4} \cdot \frac{i}{3} \cdot \frac{1}{P(R^1|\Omega)} \propto i.\] Anche senza calcolare $P(R^1|\Omega)$, si vede subito che l’alternativa più probabile è $i_{\max} = 3$, ossia l’urna contiene solo palline rosse – questa è anche la risposta mediante massima verosimglianza. Osserviamo anche un fatto banale, ma che ci rassicura: vale $P(A_0 | R^1) = 0$, perché il robot è ora certo che vi sia almeno una pallina rossa.

Possiamo confrontare la densità discreta con un grafico a barre:

# Calcoliamo vettori delle densità a
# priori e dopo aver osservato R^1

(dens_apriori <- rep(1/4, 4))

## [1] 0.25 0.25 0.25 0.25

dens_R1 <- 0:3
(dens_R1 <- dens_R1/sum(dens_R1))

## [1] 0.0000000 0.1666667 0.3333333 0.5000000

# Per il grafico a barre vanno inseriti
# in una matrice

dens_matrice <- matrix(c(dens_apriori, dens_R1),
  nrow = 2, byrow = TRUE)

# Alcuni parametri per il grafico e il
# comando barplot()

alternative <- as.character(0:3)
colori <- miei_colori[1:2]

barplot(dens_matrice, beside = TRUE, col = colori,
  names.arg = alternative, ylab = "probabilità",
  xlab = "alternativa")

# Aggiungiamo infine una legenda

legend("top", fill = colori, legend = c("A priori",
  "Sapendo R^1"), cex = 0.8)

Figura 2.14: Densità discreta delle alternative a priori e dopo una estrazione.

Remark (decisioni e test statistici). Dopo aver osservato $R^1$ e determinato l’alternativa più probabile, in questo caso $A_3$, ossia l’urna contiene solo palline rosse, si può affermare con abbastanza sicurezza che la realtà sia proprio questa? L’approccio di massima verosimiglianza darebbe appunto questa indicazione, ma è ovvio in questo caso che, essendo la probabilità $P(A_3|R^1) = 1/2$, l’incertezza è ancora grande.

Nella pratica, il calcolo delle probabilità serve spesso appunto a valutare l’incertezza sugli effetti di intraprendere determinate iniziative e indicare quindi quali decisioni mettere in pratica. Almeno una volta nella vita avremo guardato le previsioni del tempo e sulla base della probabilità che piova abbiamo stabilito se organizzare un’uscita nel fine settimana. Una decisione in ogni caso realistico è basata su molti altri aspetti, oltre al fatto che l’alternativa sia la più probabile: esiste una intera teoria della decisione che si occupa di questo problema, specialmente dal punto di vista economico/utilitario.

Un punto di vista analogo, ma tradizionalmente legato al metodo scientifico di alcune discipline, è fornito dalla teoria dei test statistici: le alternative $A_i$ sono pensate come “ipotesi” e sulla base delle osservazioni si decide quali siano confutate/falsificate e pertanto vadano scartate (il termite tecnico è rifiutate). Le rimanenti ipotesi sono quindi “accettate”, ma sempre con il beneficio del dubbio, nello stesso senso in cui una teoria scientifica rimane valida fintanto che non si trova un esperimento che ci indichi che qualcosa non va e pertanto debba essere modificata.

Tecnicamente, nei test statistici le alternative vengono collezionate in due sottoinsiemi disgiunti, detti rispettivamente l’ipotesi nulla (e indicato con $\mathcal{H}_0$) e l’alternativa (indicato con $\mathcal{H}_1$). L’attenzione principale è rivolta a rifiutare l’ipotesi nulla, ossia scartare tutte le alternative $A_i \in \mathcal{H}_0$ (questa è in un certo senso la decisione da prendere o meno sulla base della evidenza $B$). Se si osserva che $B$ vale, si decide di rifiutare l’ipotesi nulla se la probabilità che $B$ sia vero, condizionata a una qualsiasi $A_i$ tra quelle dell’ipotesi nulla $\mathcal{H}_0$ è troppo bassa. Si stabilisce quindi una soglia (detto tecnicamente livello di significatività del test, spesso indicato con $\alpha \in (0,1)$) sotto la quale si ritiene troppo poco probabile che $B$ possa essere vero se vale una qualsiasi $A_i$ dell’ipotesi nulla, e quindi $\mathcal{H}_0$ va scartato, a favore di $\mathcal{H}_1$. La quantità centrale della teoria dei test statistici è quindi il valore $p$ (in inglese $p$-value) introdotto da Fisher, definito (con la nostra notazione), come \[ p = \max_{A_i \in \mathcal{H}_0} P(B | A_i) = \max_{A_i \in \mathcal{H}_0} L(A_i; B).\] Se $p$ è minore del livello di significatività $\alpha$, significa che assumendo una qualsiasi $A_i$ tra quelli dell’ipotesi nulla $\mathcal{H}_0$, la probabilità che $B$ sia vero è comunque minore di $\alpha$, e quindi se si osserva $B$ si decide di rifiutare $\mathcal{H}_0$. Più piccolo è il valore $p$, minore è la probabilità che $B$ sia vero rispetto a l’ipotesi $\mathcal{H}_0$, e quindi saremo più sicuri nel rifiutarla.

Osserviamo però che con questo metodo non ci si domanda quanto $B$ fosse probabile se si assume una $A_i$ tra quelle delle $\mathcal{H}_1$: si potrebbe criticare quindi che l’ipotesi venga scartata semplicemente perché $B$ è sempre poco probabile (ma comunque la decisione interviene solo se lo si è osservato nella realtà). Una ulteriore critica è che, per determinare il valore $p$, bisognerebbe considerare la probabilità $P(A_i |B)$, e non la verosimiglianza $P(B|A_i) = L(A_i; B)$, perché non necessariamente tutte le ipotesi a priori sono ugualmente probabili. Ovviamente vi sono ulteriori tecniche per mitigare i possibili usi errati della teoria dei test statistici che sorgono dalle osservazioni critiche sopra. Quello dei test rimane sicuramente uno strumento importante, ma spesso poco compreso (e a volte usato erroneamente).

Torniamo all’esempio dell’urna, e chiediamoci cosa deduce il robot se viene informato che alla seconda estrazione (senza rimpiazzo) la pallina estratta è blu. Si tratta allora di condizionare anche rispetto a $B^2$. Notiamo che non serve tornare alla densità iniziale (uniforme), ma basta usare come nuova densità a priori la densità discreta rispetto all’informazione $R^1$, e si trova \[ P(A_i | R^1,B^2 ) = P(A_i|R^1) P(B^2 | R^1,A_i) \cdot \frac{1}{P(B^2|R^1)}.\] Se vale $A_i$, allora sapendo $R^1$ ci sono due palline nell’urna, di cui $3-i$ blu (purché $i \neq 0$), quindi \[P(B^2 |R^1, A_i) = \frac{3-i}{2}.\] Si ottiene allora che \[ P(A_i|R^1, B^2) \propto i (3-i)\] e quindi la probabilità rilusta uniforme, ma solamente sulle alternative $A_1$, $A_2$ (d’altra parte $A_0$ e $A_3$ sono diventate trascurabili). Possiamo anche confrontarle in un grafico a barre.

# La 'nuova' densità a priori è quella
# ottenuta prima, avendo osservato R^1,
# perciò basta definire la densità
# avendo osservato anche B^2

dens_B2_R1 <- dens_R1 * 3:0
(dens_B2_R1 <- dens_B2_R1/sum(dens_B2_R1))

## [1] 0.0 0.5 0.5 0.0

# Per il grafico a barre le inseriamo
# tutte in una matrice

dens_matrice <- matrix(c(dens_apriori, dens_R1,
  dens_B2_R1), nrow = 3, byrow = TRUE)

# Alcuni parametri per il grafico e il
# comando barplot()

alternative <- as.character(0:3)
colori <- miei_colori[1:3]
barplot(dens_matrice, beside = TRUE, col = colori,
  names.arg = alternative, ylab = "probabilità",
  xlab = "alternativa")

# Legenda

legend("topleft", fill = colori, legend = c("A priori",
  "Sapendo R^1", "Sapendo R^1 e B^2"),
  cex = 0.8)

Figura 2.15: Densità discreta delle alternative fino a due estrazioni.

È ovvio ma interessante comunque osservare come l’alternativa $A_3$, che dopo la prima estrazione era la più probabile, ora è trascurabile. Inoltre il robot risulta completamente indeciso tra l’alternativa $A_1$ e l’alternativa $A_2$. Aggiungere informazione in questo caso ha aumentato la sua incertezza, e non rimane che effettuare l’ultima estrazione per scoprire quale sia il colore della terza pallina!

2.7.1 Esercizi

Esercizio 2.13 Cosa accade nell’esempio sopra se il robot viene informato che le prime due estrazioni sono invece $R^1$, $R^2$? calcolare le probabilità e confrontare visivamente le densità discrete ottenute mediante grafici a barre.

Esercizio 2.14 Calcolare e plottare una variante dell’esempio sopra in cui $N=6$, inizialmente il robot non sa quante palline rosse vi siano e si informa poi che nelle prime $3$ estrazioni le palline estratte sono tutte dello stesso colore (ma non si specifica quale).