6.5 Stima dei parametri
In questa sezione consideriamo il problema di stimare i parametri di un processo di Markov omogeneo, sulla base di osservazioni di una traiettoria. Per quanto visto nelle sezioni precedenti, i parametri sono la matrice delle probabilità di transizione \(Q\) per una catena di Markov \((X_n)_{n}\) o delle intensità di salto \(L\) per un processo di Markov a salti \((X_t)_t\), ed eventualmente la densità discreta della marginale al tempo iniziale.
Per semplificare l’esposizione, consideriamo prima il caso di una catena di Markov \((X_n)_n\) su un insieme di stati \(E\), in cui sia noto che \(X_0 = x_0\). Richiamando per chiarezza il robot ideale, in questo caso la matrice di transizione \(Q\) non è nota al robot, ma esso viene informato che \(X\) segue un cammino \(\gamma = (x_0 \to x_1 \to \ldots \to x_n)\), ossia \(X_0 = x_0\), \(X_1=x_1\), …, \(X_n = x_n\) (brevemente scriviamo \(X = \gamma\)). L’approccio bayesiano consiste nel considerare la matrice di transizione \(Q\) come una variabile aleatoria \(\mathcal{Q}\) a valori nelle matrici quadrate \(\R^{E \times E}\) (più precisamente, sappiamo che i possibili valori di \(\mathcal{Q}\) sono matrici stocastiche). Avendo stabilito una densità a priori per \(\mathcal{Q}\), ad esempio uniforme sulle matrici stocastiche, la densità a posteriori è data dalla formula di Bayes \[ p( \mathcal{Q} = Q | X = \gamma) \propto p(\mathcal{Q} =Q ) L( \mathcal{Q} =Q ; X = \gamma),\] dove la verosimiglianza \(L\) è definita al solito come \[ L( \mathcal{Q} =Q ; X = \gamma) = P(X=\gamma| \mathcal{Q} = Q) = P(X_0 = x_0) Q_\gamma = Q_{\gamma},\] perché è noto a priori che \(X_0 = x_0\). Il peso del cammino osservato può essere riscritto raccogliendo i fattori ripetuti, ossia \[ Q_\gamma = \prod_{k=1}^n Q_{x_{k-1} \to x_{k}} = \prod_{i,j \in E} Q_{i\to j}^{\gamma_{i\to j}},\] dove \(\gamma_{i\to j}\) indica il numero di transizioni dallo stato \(i\in E\) a \(j \in E\) che avvengono nel cammino \(\gamma\). In particolare, vale \[ n = \sum_{i,j \in E} \gamma_{i\to j}.\] Se la densità a priori per \(\mathcal{Q}\) è uniforme (sull’insieme delle matrici stocastiche), ossia \(p(\mathcal{Q} = Q) \propto 1\), la densità a posteriori è \[ p(\mathcal{Q} = Q| X = \gamma) \propto L(\mathcal{Q} = Q; X=\gamma) = Q_\gamma = \prod_{i,j \in E} Q_{i\to j}^{\gamma_{i\to j}}.\] Osserviamo che la densità è un prodotto delle marginali, ma dovendo essere \(\mathcal{Q}\) una matrice stocastica le componenti relative ad una stessa riga \(i\), ad esempio \(\mathcal{Q}_{i \to j}\), \(\mathcal{Q}_{i \to k}\) non sono indipendenti (la somma deve essere \(1\)). Possiamo tuttavia affermare che variabili aleatorie associate alle righe di \(\mathcal{Q}\) sono tra loro indipendenti.
Per calcolare il punto di massimo, ossia la stima di massima verosimiglianza \(Q_{\mle}\) dobbiamo tenere conto del vincolo che la somma delle righe della matrice \(\mathcal{Q}\) valga \(1\). Il metodo generale per determinare massimi o minimi di funzioni vincolate consiste nell’introduzione di moltiplicatori di Lagrange, in modo da esprimere che nei punti critici il gradiente della funzione sia ortogonale al vincolo. Nel nostro caso, il vincolo però è così semplice che possiamo evitare l’uso di questa tecnica semplicemente esprimendo la diagonale di \(\mathcal{Q}\) in funzione delle altre entrate sulla riga: \[ Q_{i \to i} = 1 - \sum_{j \neq i} Q_{i \to j} \quad \text{per ogni $i \in E$.}\] Possiamo quindi riscrivere la verosimiglianza nel seguente modo \[ L( \mathcal{Q} = Q ; X=\gamma) = \prod_{i\in E}(1- \sum_{j \neq i} Q_{i \to j})^{\gamma_{i \to i}} \prod_{j \neq i}Q_{i\to j}^{\gamma_{i \to j}}.\] Poiché le righe sono tra loro indipendenti, possiamo ragionare separatamente per ciascuna riga \(i\), ossia determinare il massimo della funzione \[(Q_{i\to j})_{j \neq i} \mapsto (1- \sum_{j \neq i} Q_{i \to j})^{\gamma_{i \to i}} \prod_{j \neq i} Q_{i \to j}^{\gamma_{i \to j}},\] ovvero, passando al logaritmo, \[ \gamma_{i\to i} \log (1- \sum_{j \neq i} Q_{i \to j}) +\sum_{j \neq i} \gamma_{i \to j} \log(Q_{i \to j}).\] Imponendo che la derivata rispetto a ciascuna variabile \(Q_{i\to k}\) (per \(k \neq j\) si annulli, troviamo l’equazione \[0 = \frac{d}{d Q_{i\to k}} \gamma_{i\to i} \log (1- \sum_{j \neq i} Q_{i \to j}) +\sum_{j \neq i} \gamma_{i \to j} \log(Q_{i \to j}) = -\frac{\gamma_{i \to i}}{1- \sum_{j \neq i} Q_{i \to j}} + \frac{\gamma_{i \to k}}{Q_{i \to k}}.\] da cui \[ Q_{i \to k} = \gamma_{i \to k} \cdot \frac{1- \sum_{j \neq i} Q_{i \to j}}{\gamma_{i \to i}}\] Il secondo termine nel prodotto sopra non dipende da \(k\), e quindi sommando questa relazione per \(k \neq i\) troviamo che \[ \sum_{k \neq i} Q_{i \to k} = \sum_{k \neq i} \gamma_{i \to k } \frac{1- \sum_{j \neq i} Q_{i \to j}}{\gamma_{i \to i}}.\] Ricordando che \(Q_{i \to i} = 1- \sum_{j \neq i} Q_{i \to j}\) abbiamo quindi la relazione \[ 1- Q_{i \to i} = \frac{ \sum_{k \neq i} \gamma_{i \to k } }{ \gamma_{i \to i}} Q_{i \to i},\] da cui ricaviamo che \[Q _{i \to i} = \frac{\gamma_{i \to i}}{\sum_{j \in E }\gamma_{i \to j}}.\] In altre parole, abbiamo trovato che la densità discreta di probabilità \((Q_{i \to k})_{k \in E}\) è proporzionale al numero di salti osservati \((\gamma_{i \to k})_{k \in E}\), \[ Q_{i \to k} \propto \gamma_{i \to k}\] o più esplicitamente \[ Q_{i \to k} = \frac{\gamma_{i \to k}}{\sum_{j \in E } \gamma_{i \to j}}.\]
Remark. L’espressione sopra per la densità a posteriori e i calcoli per la stima di massima verosimiglianza suggerisce altre densità a priori (dette di Dirichlet) della forma \[ p(\mathcal{Q}= Q) \propto \prod_{i,j \in E} Q_{i\to j}^{\alpha_{ij}},\] per opportuni parametri \(\alpha_{ij}\ge 0\) (dove a essere precisi bisognerebbe scrivere \(Q_{i \to i} = 1- \sum_{j \neq i} Q_{i \to j}\)). Notiamo che i calcoli per la stima di massima verosimiglianza mostrano anche che la moda della densità sopra è data dalla matrice \[ Q_{i \to j} = \frac{ \alpha_{ij}}{\sum_{k \in E} \alpha_{ik}}.\] La formula di Bayes darebbe quindi come densità a posteriori \[ p(\mathcal{Q}= Q | X = \gamma ) \propto \prod_{i,j \in E} Q_{i\to j}^{\alpha_{ij} +\gamma_{i \to j}},\] e di conseguenza la stima di massima densità a posteriori \(\mathcal{Q}_{MAP}\), seguendo gli stessi calcoli della stima di massima verosimiglianza, è \[ Q_{i \to j} =\frac{ \alpha_{ij}+\gamma_{i \to j}}{ \sum_{k \in E} \alpha_{ik} + \gamma_{i \to k}}.\]
Remark. Abbiamo supposto che il cammino osservato parta al tempo \(0\) da \(x_0\). Tuttavia se iniziasse da un tempo successivo, allora si pone il problema di stimare \(X_0\). Si può assumere ad esempio che \(X\) sia stazionaria, e quindi supporre che \(\pi_0\) sia una distribuzione invariante. Il problema è che questa dipende da \(Q\) in modo tutt’altro che banale.
Nel caso di processi di Markov a salti, l’argomento è analogo ma si basa sulla formula (6.7). Per brevità non consideriamo l’approccio bayesiano ma presentiamo solo la stima di massima verosimiglianza. Si consideri un cammino \(\gamma = (x_0\to x_1 \ldots x_n)\) che rimane per un tempo \(t_1\) nello stato \(x_0\), \(t_2\) nello stato \(x_1\) ecc., e si supponga di osservare \(X = \gamma\), ossia tutta la traiettoria da \(X_0 = x_0\) fino a \(X_{t_1+\ldots +t_n} = x_n\). Allora la stima di massima verosimiglianza per la matrice di intensità di salto \(\mathcal{L}_{MLE}\) si ottiene massimizzando l’espressione \[ \prod_{k=1}^n \exp\bra{t_k L_{x_{k-1} \to x_{k-1}}} L_{x_{k-1} \to x_k} = \prod_{i \in E} \exp\bra{\gamma_{i \to i} L_{i\to i} } \prod_{ i\neq j \in E} L_{i\to j}^{\gamma_{i \to j}}\] dove stavolta si è posto \(\gamma_{i \to i}\) il tempo totale trascorso dal cammino nello stato \(i \in E\). Inoltre, poiché la somma delle righe di \(L\) è nulla, possiamo porre \[ \exp\bra{\gamma_{i \to i} L_{i\to i} } = \exp\bra{- \gamma_{i \to i} \sum_{j \neq i} L_{i\to j} }.\] Passando ai logaritmi e derivando rispetto a ciascun parametro \(L_{i \to j}\) si ottiene che \(\mathcal{L}_{MLE}\) è data dall’espressione, per \(i \neq j\), \[ L_{i \to j} = \frac{ \gamma_{i \to j}}{\gamma_{i \to i} }.\]
Remark. Negli esempi sopra si suppone di osservare completamente la catena \(X\) in un intervallo (discreto o continuo) di tempi. Più in generale ci si può chiedere cosa accada se mancano le osservazioni delle variabili \(X_k\) in alcuni tempi, oppure se si osserva solamente una funzione \(g(X_k)\) della catena invece, di \(X_k\), o più in generale una funzione \(g(X_k, Z_k)\) dove \(Z\) è un processo indipendente da \(X\). Un esempio fondamentale è dato dal caso in cui \(X\) è un segnale e \(Z\) è un “rumore” che si vorrebbe rimuovere, o filtrare. In queste situazioni si parla di modelli di Markov nascosti (in inglese Hidden Markov Models, HMM) che hanno molteplici applicazioni. Opportune modifiche degli argomenti visti sopra permettono di introdurre algoritmi specifici per stimare i parametri di un HMM, come pure stimare \(X_k\) dalle osservazioni \(g(X_k, Z_k)\) o anche effettuare previsioni.