Ragionamento probabilistico nel tempo - catene di Markov, HMM, filtraggio e smoothing

Le reti bayesiane (Reti bayesiane - semantica, indipendenze e inferenza esattaUna rete bayesiana è un grafo diretto aciclico (DAG) con una variabile per nodo e una tabella delle probabilità condizionate (CPT) $P(X_i\mid\text{Parents}(X_i))$ per nodo; rappresenta la congiunta come prodotto $P(x_1,\dots,x_n)=\prod_i P(x_i\mid\text{parents}(X_i))$, con molti meno numeri della tabella completa (rete dell'allarme: 10 contro 31). Tre strutture elementari: catena $X\to Y\to Z$ e fork $Y\leftarrow X\to Z$ (estremi indipendenti dato il centro), collider $X\to Z\leftarrow Y$ (indipendenti, ma dipendenti se si osserva $Z$ o un suo discendente). Modelli causali (genitori = cause) danno reti più piccole. Inferenza esatta per enumerazione: $\mathbf P(X\mid\mathbf e)=\alpha\sum_{\mathbf y}\mathbf P(X,\mathbf e,\mathbf y)$ calcolata come somma di prodotti, in profondità e da sinistra a destra; per l'allarme $\mathbf P(B\mid j,m)=\langle0{,}284;\ 0{,}716\rangle$. Molti prodotti si ripetono: l'eliminazione di variabili li salva come fattori e li riusa.Reti bayesiane - semantica, indipendenze e inferenza esatta →) descrivono un mondo fotografato in un istante. Qui lo stato cambia nel tempo e l'agente lo osserva solo attraverso sensori rumorosi.

Esempio guida (mondo dell'ombrello). Sei prigioniero in una cella sotterranea: non vedi il tempo, ma vedi se la guardia porta l'ombrello. Quanto è probabile che fuori piova, date le osservazioni di oggi e dei giorni precedenti?

Modello

• Tempo discreto: istanti (time slices) a distanza costante $\Delta$.
• $\mathbf X_t$: variabili di stato dell'istante $t$, non osservabili (nascoste). $\mathbf E_t$: variabili di evidenza osservabili; l'osservazione è $\mathbf E_t=\mathbf e_t$. Si scrive $e_{1:t}=e_1,\dots,e_t$.

Modello di transizione. $\mathbf P(X_t\mid X_{0:t-1})$ ha troppe variabili (crescono con $t$). Si usa l'ipotesi di Markov: lo stato dipende solo da un numero fisso di stati precedenti.

• Primo ordine: $\mathbf P(X_t\mid X_{0:t-1})=\mathbf P(X_t\mid X_{t-1})$.
• Secondo ordine: $\mathbf P(X_t\mid X_{t-2},X_{t-1})$.

Un processo di Markov del primo ordine è una catena di nodi $X_0\to X_1\to X_2\to\cdots$.

Modello del sensore. Le osservazioni dipendono solo dallo stato corrente (sensor Markov assumption): $\mathbf P(E_t\mid X_{0:t},E_{1:t-1})=\mathbf P(E_t\mid X_t)$. Dato $X_t$, $E_t$ è indipendente da tutto il resto: è un caso di catena e di fork (Reti bayesiane - semantica, indipendenze e inferenza esattaUna rete bayesiana è un grafo diretto aciclico (DAG) con una variabile per nodo e una tabella delle probabilità condizionate (CPT) $P(X_i\mid\text{Parents}(X_i))$ per nodo; rappresenta la congiunta come prodotto $P(x_1,\dots,x_n)=\prod_i P(x_i\mid\text{parents}(X_i))$, con molti meno numeri della tabella completa (rete dell'allarme: 10 contro 31). Tre strutture elementari: catena $X\to Y\to Z$ e fork $Y\leftarrow X\to Z$ (estremi indipendenti dato il centro), collider $X\to Z\leftarrow Y$ (indipendenti, ma dipendenti se si osserva $Z$ o un suo discendente). Modelli causali (genitori = cause) danno reti più piccole. Inferenza esatta per enumerazione: $\mathbf P(X\mid\mathbf e)=\alpha\sum_{\mathbf y}\mathbf P(X,\mathbf e,\mathbf y)$ calcolata come somma di prodotti, in profondità e da sinistra a destra; per l'allarme $\mathbf P(B\mid j,m)=\langle0{,}284;\ 0{,}716\rangle$. Molti prodotti si ripetono: l'eliminazione di variabili li salva come fattori e li riusa.Reti bayesiane - semantica, indipendenze e inferenza esatta →).

Definizione (HMM del primo ordine). Hidden Markov model: stati nascosti $X_t$ che formano una catena di Markov e osservazioni $E_t$ figlie dello stato. Distribuzione congiunta, con prior $\mathbf P(X_0)$: $$P(x_0,x_{1:t},e_{1:t})=P(x_0)\prod_{i=1}^t P(x_i\mid x_{i-1})\,P(e_i\mid x_i).$$ In senso stretto l'HMM ha un'unica variabile di stato discreta (le evidenze possono essere più d'una).

Mondo dell'ombrello. $X_t=\text{Rain}_t$, $E_t=\text{Umbrella}_t$. Trasformazione: $P(r_t\mid r_{t-1})=0{,}7$, $P(r_t\mid\neg r_{t-1})=0{,}3$. Sensore: $P(u_t\mid r_t)=0{,}9$, $P(u_t\mid\neg r_t)=0{,}2$. Prior $P(r_0)=0{,}5$.

Rappresentazione a matrici

Per $X_t$ con $S$ valori:

• Matrice di transizione $\mathbf T$, $S\times S$, con $T_{ij}=P(X_t=j\mid X_{t-1}=i)$ (riga = stato di partenza; ogni riga somma a $1$).
• Matrice del sensore $\mathbf O_t$, $S\times S$ diagonale, con $(\mathbf O_t)_{ii}=P(e_t\mid X_t=i)$ (probabilità di emissione). Ce n'è una per ogni valore osservato.

Ombrello (stati in ordine $[\text{pioggia},\ \text{no pioggia}]$): $$\mathbf T=\begin{pmatrix}0{,}7&0{,}3\\0{,}3&0{,}7\end{pmatrix},\quad \mathbf O^{u}=\begin{pmatrix}0{,}9&0\\0&0{,}2\end{pmatrix},\quad \mathbf O^{\neg u}=\begin{pmatrix}0{,}1&0\\0&0{,}8\end{pmatrix}.$$

I quattro compiti

Inoltre i modelli $\mathbf T$ e $\mathbf O$, se non noti, si possono apprendere dalle osservazioni.

Filtraggio

Il belief state è $\mathbf f_t=\mathbf P(X_t\mid e_{1:t})$. Si calcola ricorsivamente da $\mathbf f_{t}$ a $\mathbf f_{t+1}$ in due passi:

$$\mathbf f_{t+1}=\mathbf P(X_{t+1}\mid e_{1:t+1})\overset{\text{Bayes}}{=}\alpha\,\underbrace{\mathbf P(e_{t+1}\mid X_{t+1})}_{\text{aggiornamento}}\ \underbrace{\sum_{x_t}\mathbf P(X_{t+1}\mid x_t)\,P(x_t\mid e_{1:t})}_{\text{predizione}},\qquad \alpha=\frac1{P(e_{t+1}\mid e_{1:t})}.$$

Dimostrazione dei due passi: $P(X_{t+1}\mid e_{1:t+1})=\alpha P(e_{t+1}\mid X_{t+1},e_{1:t})P(X_{t+1}\mid e_{1:t})$ (Bayes); per l'ipotesi del sensore il primo fattore diventa $P(e_{t+1}\mid X_{t+1})$; il secondo si scrive marginalizzando $x_t$ e usando l'ipotesi di Markov $P(X_{t+1}\mid x_t,e_{1:t})=P(X_{t+1}\mid x_t)$. In forma compatta $\mathbf f_{t+1}=\text{FORWARD}(\mathbf f_t,e_{t+1})$, applicata ricorsivamente.

In matrici. Il passo di predizione è un prodotto con $\mathbf T^\top$ (si sommano le righe pesate per $\mathbf f_t$), l'aggiornamento un prodotto con $\mathbf O_{t+1}$: $$\mathbf f_{t+1}=\alpha\,\mathbf O_{t+1}\,\mathbf T^\top\,\mathbf f_t.$$

Algoritmo (forward). $\mathbf f_0=\mathbf P(X_0)$. Per ogni osservazione $e_t$: $\mathbf v=\mathbf O^{e_t}\mathbf T^\top\mathbf f_{t-1}$, $\mathbf f_t=\mathbf v/\sum\mathbf v$. Costo costante per passo: non dipende dalla lunghezza della storia.

Esempio: ombrello due giorni

Domanda: $P(r_2\mid u_1,u_2)$, con $\mathbf f_0=\langle0{,}5;\ 0{,}5\rangle$.

Giorno 1. Predizione: $\mathbf T^\top\mathbf f_0=\langle0{,}7\cdot0{,}5+0{,}3\cdot0{,}5;\ 0{,}3\cdot0{,}5+0{,}7\cdot0{,}5\rangle=\langle0{,}5;\ 0{,}5\rangle$ (con prior uniforme resta uniforme). Aggiornamento con $u_1$: $\langle0{,}9\cdot0{,}5;\ 0{,}2\cdot0{,}5\rangle=\langle0{,}45;\ 0{,}10\rangle$; $\alpha_1=1/0{,}55$: $$\mathbf f_1=\langle0{,}45;\ 0{,}10\rangle/0{,}55=\langle0{,}818;\ 0{,}182\rangle.$$

Giorno 2. Predizione: $P(r_2\mid u_1)=0{,}7\cdot0{,}818+0{,}3\cdot0{,}182=0{,}627$, $P(\neg r_2\mid u_1)=0{,}3\cdot0{,}818+0{,}7\cdot0{,}182=0{,}373$. Aggiornamento con $u_2$: $\langle0{,}9\cdot0{,}627;\ 0{,}2\cdot0{,}373\rangle=\langle0{,}5645;\ 0{,}0745\rangle$; $\alpha_2=1/0{,}639$: $$\mathbf f_2=\langle0{,}883;\ 0{,}117\rangle.$$ (Le slide scrivono $0{,}564$ e $0{,}565$ in due punti per un arrotondamento: il valore è $0{,}5645$.) Con matrici: $\mathbf f_2=\alpha\begin{pmatrix}0{,}9&0\\0&0{,}2\end{pmatrix}\begin{pmatrix}0{,}7&0{,}3\\0{,}3&0{,}7\end{pmatrix}\begin{pmatrix}0{,}818\\0{,}182\end{pmatrix}=\alpha\begin{pmatrix}0{,}5645\\0{,}0745\end{pmatrix}$.

La probabilità di pioggia sale da $0{,}5$ a $0{,}818$ al primo ombrello e a $0{,}883$ al secondo. Se al giorno 1 non si vedesse l'ombrello: $\mathbf f_1=\langle0{,}111;\ 0{,}889\rangle$ e dopo l'ombrello del giorno 2 si ottiene $0{,}7028$; se non si hanno informazioni sul giorno 1 (solo $u_2$), $0{,}818$ (Esercizio - Filtraggio nell'ombrello con pomegranate (laboratorio 6)).

Predizione

Si applica solo il passo di predizione, senza evidenza: $\mathbf P(X_{t+1}\mid e_{1:t})=\mathbf T^\top\mathbf f_t$ e in generale $\mathbf P(X_{t+k}\mid e_{1:t})=(\mathbf T^\top)^k\mathbf f_t$. Dopo $\mathbf f_2=\langle0{,}883;\ 0{,}117\rangle$: giorno 3 $\langle0{,}653;\ 0{,}347\rangle$, giorno 4 $\langle0{,}561;\ 0{,}439\rangle$, ... verso la distribuzione stazionaria $\langle0{,}5;\ 0{,}5\rangle$: l'evidenza passata viene dimenticata e la predizione tende alla distribuzione stazionaria della catena.

Catena di Markov (esempio delle slide). Voto finale annuale $A,B,C$; stesso voto dell'anno prima con probabilità $0{,}6$, voto diverso (ciascuno) $0{,}2$. $P(G_3=A\mid G_1=B)$: $$P(G_3\mid G_1=B)=\sum_{g_2}P(G_3\mid g_2)\,P(g_2\mid G_1=B)$$ (somma sul voto intermedio $g_2$; il fattore $P(G_1=B)$ del numeratore si semplifica con il denominatore), quindi per $G_3=A$ $$P(G_3=A\mid G_1=B)=\underbrace{0{,}6\cdot0{,}2}_{g_2=A}+\underbrace{0{,}2\cdot0{,}6}_{g_2=B}+\underbrace{0{,}2\cdot0{,}2}_{g_2=C}=0{,}28.$$ È l'elemento $(B,A)$ di $\mathbf T^2$: tutta la riga è $\langle0{,}28;\ 0{,}44;\ 0{,}28\rangle$.

Smoothing

Dato il passato e il presente, si rivaluta un istante $k$ passato: $\mathbf P(X_k\mid e_{1:t})$, con $1\le k<t$. Si dividono le evidenze in passate $e_{1:k}$ e future $e_{k+1:t}$: $$\mathbf P(X_k\mid e_{1:t})=\alpha\,\mathbf P(X_k\mid e_{1:k})\ \mathbf P(e_{k+1:t}\mid X_k)=\alpha\,\mathbf f_k\times\mathbf b_{k+1:t}\qquad(\times\text{ prodotto componente per componente}).$$ Il messaggio forward $\mathbf f_k$ è il filtraggio. Il messaggio backward $\mathbf b_{k+1:t}=\mathbf P(e_{k+1:t}\mid X_k)$ si calcola all'indietro: $$\mathbf b_{k+1:t}=\mathbf T\,\mathbf O_{k+1}\,\mathbf b_{k+2:t},\qquad \mathbf b_{t+1:t}=\mathbf 1=\langle1;\dots;1\rangle.$$ (Ogni componente è $\sum_{x_{k+1}}P(e_{k+1}\mid x_{k+1})\,b_{k+2:t}(x_{k+1})\,P(x_{k+1}\mid X_k)$.) L'algoritmo forward-backward fa un passaggio in avanti per i $\mathbf f$ e uno all'indietro per i $\mathbf b$, e smoothing di tutti gli istanti in $O(t)$. Il filtraggio all'ultimo istante coincide con lo smoothing.

Esempio. Osservazioni $u_1,u_2,\neg u_3$. Filtraggio: $\mathbf f_1=\langle0{,}818;\ 0{,}182\rangle$, $\mathbf f_2=\langle0{,}883;\ 0{,}117\rangle$, $\mathbf f_3=\langle0{,}191;\ 0{,}809\rangle$. All'indietro: $\mathbf b_{4:3}=\mathbf 1$; $\mathbf b_{3:3}=\mathbf T\mathbf O^{\neg u}\mathbf 1=\begin{pmatrix}0{,}7&0{,}3\\0{,}3&0{,}7\end{pmatrix}\begin{pmatrix}0{,}1\\0{,}8\end{pmatrix}=\langle0{,}31;\ 0{,}59\rangle$; $\mathbf b_{2:3}=\mathbf T\mathbf O^{u}\mathbf b_{3:3}=\mathbf T\langle0{,}9\cdot0{,}31;\ 0{,}2\cdot0{,}59\rangle=\mathbf T\langle0{,}279;\ 0{,}118\rangle=\langle0{,}2307;\ 0{,}1663\rangle$. Combinando:

• giorno 2: $\mathbf f_2\times\mathbf b_{3:3}=\langle0{,}883\cdot0{,}31;\ 0{,}117\cdot0{,}59\rangle=\langle0{,}2738;\ 0{,}0688\rangle\to\langle0{,}799;\ 0{,}201\rangle$: l'ombrello mancante al giorno 3 abbassa la pioggia al giorno 2 da $0{,}883$ a $0{,}799$;
• giorno 1: $\mathbf f_1\times\mathbf b_{2:3}=\langle0{,}818\cdot0{,}2307;\ 0{,}182\cdot0{,}1663\rangle\to\langle0{,}862;\ 0{,}138\rangle$;
• giorno 3: $\mathbf f_3\times\mathbf 1=\langle0{,}191;\ 0{,}809\rangle$ (invariato).

Con le sole osservazioni $u_1,u_2$, lo smoothing del giorno 1 dà $0{,}883$ (stesso numero di $\mathbf f_2$ per la simmetria di $\mathbf T$), e quello del giorno 0 è $0{,}653$; con la sola $u_1$, il giorno 0 vale $0{,}627$. Questi valori "all'indietro" sono quelli che la libreria pomegranate restituisce automaticamente per ogni nodo (Esercizio - Filtraggio e smoothing nell'ombrello con la nuova versione di pomegranate (materiale supplementare del laboratorio 6)). Il risultato è stato verificato per enumerazione di tutte le sequenze di stati.

Spiegazione più probabile. Si cerca la sequenza $x_{1:t}$ con $P(x_{1:t}\mid e_{1:t})$ massima: non è la sequenza degli stati più probabili singolarmente. Si risolve con l'algoritmo di Viterbi, che ha la forma del filtraggio con un massimo al posto della somma. Per $u,u,\neg u$ la sequenza più probabile è pioggia, pioggia, no pioggia.

Probabilità di una sequenza di osservazioni

Senza normalizzare, il filtraggio calcola $\mathbf s_t=\mathbf P(e_{1:t},X_t)$: $$\mathbf s_t=\mathbf O_t\mathbf T^\top\mathbf s_{t-1},\qquad P(e_{1:t})=\sum_{x_t}s_t(x_t).$$ Dimostrazione: $s_t=P(e_{1:t},X_t)=\sum_{x_{t-1}}P(e_t\mid X_t)P(X_t\mid x_{t-1})P(e_{1:t-1},x_{t-1})$ per le due ipotesi di Markov. Per l'ombrello, $P(u_1,u_2)=0{,}55\cdot0{,}639=0{,}3515$ (i due fattori di normalizzazione).

Esempio: il bambino (slide). Stato: pannolino pulito o sporco, equiprobabili all'inizio. Pulito $\to$ pulito/sporco $0{,}5/0{,}5$; sporco $\to$ pulito/sporco $0{,}8/0{,}2$. Pulito: "mama" o "papa" $0{,}4$ ciascuna, "poo-poo" o "pee-pee" $0{,}1$; sporco: $0{,}1$ e $0{,}4$. Qual è la probabilità di "mama", "mama", "pee-pee", "pee-pee"? Svolgimento in Esercizio - Filtraggio in un HMM - il bambino e il pannolino (slide del corso): $0{,}00326$.

Esempio: localizzazione di un robot

Un robot si muove a caso, con uguale probabilità, in una casella libera adiacente; i sensori (rumorosi) dicono quali dei quattro lati N, E, S, O hanno un ostacolo. Stato $X_t$: la casella, tra $S=42$ libere. Prior $P(X_0=i)=1/42$ (posizione ignota).

• Transizione: $T_{ij}=1/N(i)$ se $j$ è adiacente a $i$, altrimenti $0$ ($N(i)$ = caselle libere adiacenti). La matrice ha $42\times42=1764$ elementi.
• Sensore: $E_t$ è una sequenza di $4$ bit (es. $1010$: ostacolo a nord e sud). Ogni bit è sbagliato con probabilità $\varepsilon$ indipendente. Se la lettura $e_t$ ha $d_i$ bit sbagliati rispetto alla casella $i$: $P(e_t\mid X_t=i)=(1-\varepsilon)^{4-d_i}\,\varepsilon^{d_i}$. Con $\varepsilon=0{,}2$: tre bit giusti e uno sbagliato $0{,}8^3\cdot0{,}2=0{,}102$; la probabilità di sbagliare almeno un bit su quattro è $1-0{,}8^4=0{,}59$.
• Con $\mathbf T$ e $\mathbf O_t$ il filtraggio dà $P(X_1\mid e_1)$, $P(X_2\mid e_1,e_2)$, ... e l'errore di localizzazione scende con le osservazioni: per $\varepsilon=0{,}20$ il robot si localizza entro due caselle dopo circa $20$ osservazioni, per $\varepsilon\le0{,}10$ bastano poche osservazioni, per $\varepsilon=0{,}40$ l'errore resta grande.

Codice

import numpy as np

T = np.array([[0.7, 0.3],        # T[i, j] = P(X_t = j | X_{t-1} = i), stati: [pioggia, no pioggia]
              [0.3, 0.7]])
O = {True:  np.diag([0.9, 0.2]), # O^{u}: P(ombrello | stato) sulla diagonale
     False: np.diag([0.1, 0.8])} # O^{non u}
f0 = np.array([0.5, 0.5])        # P(X_0)

def forward(f, e):
    v = O[e] @ T.T @ f           # predizione (T^T f) e aggiornamento (O)
    return v / v.sum()           # normalizzazione alfa

def filtering(prior, evidenze):
    f, out = prior, []
    for e in evidenze:
        f = forward(f, e)
        out.append(f)
    return out                   # f_1, ..., f_t

def smoothing(prior, evidenze):
    f = filtering(prior, evidenze)
    b, out = np.ones(2), []
    for t in range(len(evidenze) - 1, -1, -1):   # da t fino a 1, messaggio all'indietro
        s = f[t] * b
        out.append(s / s.sum())                  # P(X_{t+1} | e_{1:T})
        b = T @ O[evidenze[t]] @ b               # b_{k:T}
    return out[::-1]

print(filtering(f0, [True, True]))
print(smoothing(f0, [True, True, False]))

Output: [0.818, 0.182], [0.883, 0.117] per il filtraggio; [0.862, 0.138], [0.799, 0.201], [0.191, 0.809] per lo smoothing.

Domande d'esame

1. Cosa afferma l'ipotesi di Markov del primo ordine e l'ipotesi del sensore? $P(X_t\mid X_{0:t-1})=P(X_t\mid X_{t-1})$ e $P(E_t\mid X_{0:t},E_{1:t-1})=P(E_t\mid X_t)$.
2. Nel mondo dell'ombrello con $\mathbf f_0=\langle0{,}5;\ 0{,}5\rangle$, calcolare $P(r_1\mid u_1)$. $\alpha\cdot0{,}9\cdot(0{,}7\cdot0{,}5+0{,}3\cdot0{,}5)=\alpha\cdot0{,}45$ e $\alpha\cdot0{,}2\cdot0{,}5=\alpha\cdot0{,}1$: $0{,}45/0{,}55=0{,}818$.
3. Qual è la differenza tra filtraggio e smoothing? Il filtraggio stima lo stato attuale con le evidenze fino ad ora; lo smoothing stima uno stato passato usando anche le evidenze successive (forward-backward).
4. Con i voti $A,B,C$ (stesso voto $0{,}6$, altro $0{,}2$), $P(G_3=A\mid G_1=B)$? $0{,}6\cdot0{,}2+0{,}2\cdot0{,}6+0{,}2\cdot0{,}2=0{,}28$.