% \paragraph{Algoritmo di Baum}


Prendendo spunto da Rabiner \cite{Rabiner:tutorialHMM} è stato implementato un algoritmo iterativo che si rifà a una versione dell'algoritmo di Baum-Welch. Di fatto, uno dei problemi più difficili riguardanti gli HMM è quello di aggiustare i parametri del modello ($A,B,\pi$) in modo da massimizzare la probabilità della sequenza di osservazioni dato il modello. 
Non è ancora noto, però, un metodo per risolvere analiticamente questo problema, cioè non c'è garanzia di trovare i valori ottimi dei parametri. 
Tuttavia, con il metodo proposto si dovrebbe riuscire ad avere una massimizzazione almeno locale (non è infatti sempre garantita la convergenza entro i limiti desiderati).\\

Si definisce innanzitutto la variabile relativa alla passata in avanti come
\begin{equation}
\label{eq:alpha}
\alpha_t(i) = \mathcal{P}[O_1,O_2,\dots ,O_t,q_t=S_i | \lambda] 
\end{equation}
che rappresenta la probabilità che si verifichi la parziale sequenza di osservazioni $O_1,O_2,\dots ,O_t$ (fino all'istante $t$) e di trovarsi nello stato $S_i$ al tempo $t$, dato il modello $\lambda $.
Essa può essere calcolata iterativamente secondo 
\begin{equation}
\alpha_{t+1}(j) = \left[ \sum_{i=1}^{N} \alpha_t(i)a_{ij} \right] b_j(O_{t+1}) \qquad 1 \leq t \leq T-1 \qquad 1 \leq j \leq N
\end{equation}
e inizializzata con
\begin{equation}
\alpha_1(i) = \pi_i b_i (O_1) \qquad 1 \leq i \leq N \,.
\end{equation}
Si ha in particolare che la probabilità dell'intera sequenza di osservazioni $ O $ è data da 
\begin{equation}
\mathcal{P}[O | \lambda] = \sum_{i=1}^N \alpha_T(i) \,.
\end{equation}

Analogamente si definisce la variabile relativa alla passata all'indietro come
\begin{equation}
\beta_t(i) = \mathcal{P} [O_{t+1},O_{t+2},\dots ,O_T | q_t=S_i,\lambda] 
\end{equation}
che rappresenta la probabilità che si verifichi la parziale sequenza di osservazioni dall'istante $ t+1 $ fino alla fine, dato di trovarsi nello stato $S_i$ al tempo $t$ e dato il modello $\lambda $. 
Anch'essa si calcola iterativamente, secondo 
\begin{equation}
\beta_{t}(i) = \sum_{j=1}^{N} a_{ij} b_j(O_{t+1}) \beta_{t+1}(j)  \qquad t=T-1,T-2,\dots,1 \qquad 1 \leq i \leq N
\end{equation}
e inizializzata con
\begin{equation}
\beta_T(i) = 1 \qquad 1 \leq i \leq N \,.
\end{equation}

Si definisce poi una variabile $\xi_t(i,j)$ che rappresenta la probabilità di trovarsi nello stato $S_i$ al tempo $t$ e nello stato $S_j$ al tempo $t+1$, dato il modello e l'intera sequenza di osservazioni, ossia
\begin{equation}
\xi_t(i,j) = \mathcal{P}[q_t=S_i, q_{t+1}=S_j | O,\lambda ] \,. 
\end{equation}
Per come sono definite le variabili $\alpha$ e $\beta$, risulta che 
\begin{equation}
\xi_t(i,j) = \frac{\alpha_t(i) a_{ij} b_j(O_{t+1}) \beta_{t+1}(j)}{\mathcal{P}[O | \lambda]} \,.
\end{equation}

Sia infine $\gamma_t(i)$ la probabilità di trovarsi nello stato $S_i$ al tempo $ t $, date tutte le osservazioni e il modello:
\begin{equation}
\gamma_t(i) = \mathcal{P}[q_t=S_i | O,\lambda ] \,.
\end{equation}
Essa è in relazione con la precedente variabile secondo
\begin{equation}
\gamma_t(i) = \sum_{j=1}^N \xi_t(i,j) \,.
\end{equation}

Se si sommano le $\gamma$ sull'indice $t$ si ottiene ua quantità che può essere interpretata come il numero atteso di volte che lo stato $S_i$ sia visitato o, equivalentemnente, il numero atteso di transizioni che avvengono dallo stato $S_i$ (escludendo dalla somma l'istante $t = T$). 
Analogamente la somma delle $\xi$ su $t$ (tra 1 e $T-1$) può essere interpretata come il numero atteso di transizioni dallo stato $S_i$ allo stato $S_j$. Quindi
\begin{subequations}
\begin{align}
\sum_{t=1}^{T-1} \gamma_t(i) &= \text{numero atteso di transizioni da $S_i$} \\
\sum_{t=1}^{T-1} \xi_t(i,j) &=  \text{numero atteso di transizioni da $S_i$ a $S_j$} \,.
\end{align}
\end{subequations}

Si arriva così a considerare come ragionevole ""ri-stima'' dei parametri del modello le seguenti formule:
\begin{subequations}
\begin{align}
\begin{split}
\bar{\pi}_i &= \text{numero atteso di volte nello stato $S_i$ al tempo $t=1$} \\
&= \gamma_1(i) 
\end{split} \\
\begin{split}
\label{eq:ri-stima per a}
\bar{a}_{ij} &= \frac{\text{numero atteso di transizioni da $S_i$ a $S_j$}}{\text{numero atteso di transizioni da $S_i$}}  \\
&=\frac{\displaystyle\sum_{t=1}^{T-1} \xi_t(i,j)}{\displaystyle\sum_{t=1}^{T-1} \gamma_t(i)} 
\end{split} \\
\begin{split}
\bar{b}_j(k) &= \frac{\text{numero atteso di volte di essere nello stato $S_j$ e osservare $v_k$}}{\text{numero atteso di volte di essere nello stato $S_j$ }}  \\
&= \frac{\displaystyle\sum_{\underset{\text{t.c.} \, O_t=v_k}{t=1}}^{T-1} \gamma_t(j)}{\displaystyle\sum_{t=1}^{T-1} \gamma_t(i)}
\end{split}
\end{align}
\end{subequations}

E' stato dimostrato da Baum e colleghi che per il nuovo modello $ \bar{\lambda} = (\bar{A}, \bar{B}, \bar{\pi}) $ vale $ \mathcal{P}[O|\bar{\lambda}] \geq  \mathcal{P}[O|\lambda)]$ cioè con $ \bar{\lambda}$ si ha una probabilità maggiore (o al più uguale) che la sequenza sia generata secondo esso piuttosto che secondo il modello iniziale $ \lambda $. 

Nella pratica questa procedura sarà ripetuta più volte usando $ \bar{\lambda} $ al posto di 
$ \lambda $ fino ad ottenere una sorta di convergenza dei parametri. Va sottolineato che comunque questo è un algoritmo che tende a trovare un ottimo locale. E' importante quindi la scelta dei parametri iniziali: meglio questi descrivono la situazione corrente, più speranze si avrà di trovare il punto di ottimo globale.

\paragraph{Versione scalata dell'algoritmo}
\ \\

Dato che nello specifico del problema trattato la sequenza di osservazioni è molto lunga, le variabili
 $ \alpha $ e $ \beta $ tendono ad assumere valori sempre più vicini allo 0 al crescere dei campioni, fino a che la precisione del calcolatore non è più sufficiente a rappresentarli non nulli. La soluzione 
più ragionevole è quella di eseguire un riscalamento delle variabili in modo da avere valori numerici in un \emph{range} gestibile dal calcolatore. Riferendosi ad esempio ad $ \alpha_t(i) $, di fatto si moltiplica per un fattore di scala che non dipende da $i$ ma solo da $t$, e lo stesso si fa per 
$ \beta_t(i) $; alla fine tali coefficienti si compenseranno perfettamente.

Si consideri ad esempio la formula (\ref{eq:ri-stima per a}) scritta in modo esplicito in $\alpha$ e $\beta$:
\begin{equation}
\label{eq:a_bar esplicita}
\bar{a}_{ij} = \frac{\displaystyle \sum_{t=1}^{T-1} \alpha_t(i) a_{ij} b_j(O_{t+1}) \beta_{t+1}(j)}
{\displaystyle \sum_{t=1}^{T-1} \sum_{j=1}^N \alpha_t(i) a_{ij} b_j(O_{t+1}) \beta_{t+1}(j)} \,.
\end{equation}
Per quanto riguarda $ \alpha_t(j) $, si pone innanzitutto $\tilde{\alpha}_0(j) = \alpha_0(j)$, poi si definisce $c_0 = 1/\sum_{j=1}^N \tilde{\alpha}_0(j) $ e si ha la prima riscalatura $\hat{\alpha}_0(j) = c_0\tilde{\alpha}_0(j)$. 
Per $t=1,\dots,T-1$ e $j=1,\dots,N$ si calcola 
\begin{subequations}
\begin{align}
\tilde{\alpha}_t(j) &= \displaystyle \sum_{i=1}^{N} \hat{\alpha}_{t-1}(i) a_{ij}  b_j(O_t)\\
c_t &= \frac{1}{\displaystyle \sum_{j=1}^N \tilde{\alpha}_t(j)}\,.
\end{align}
\end{subequations}
e si hanno i valori scalati tramite la
\begin{equation}
\hat{\alpha}_t(j) = c_t \tilde{\alpha}_t(j) \,.
\end{equation}
Si osserva che $\hat{\alpha}_0(j) = c_0 \alpha_0(j)$ e per induzione
\begin{equation}
\label{eq:produttoria_c}
\hat{\alpha}_t(j) = c_0 c_1 \cdots c_t \alpha_t(j)= \left[ \prod_{s =1}^t c_s \right] \alpha_t(j) \,.
\end{equation}
L'aggiornamento è eseguito con la 
\begin{equation}
\hat{\alpha}_t(j) = \frac{\displaystyle \sum_{i=1}^{N} \hat{\alpha}_{t-1}(i) a_{ij}  b_j(O_t)}{\displaystyle \sum_{j=1}^{N} \sum_{i=1}^{N} \hat{\alpha}_{t-1}(i)a_{ij}  b_j(O_t)}
\end{equation}
e utilizzando la (\ref{eq:produttoria_c}) si ha anche che 
\begin{equation}
\hat{\alpha}_t(j)=  \frac{\displaystyle \sum_{i=1}^{N} \alpha_{t-1}(i)\left[ \prod_{s=1}^{t-1} c_s \right] a_{ij}  b_j(O_t)}{\displaystyle \sum_{j=1}^{N} \sum_{i=1}^{N} \alpha_{t-1}(i)\left[ \prod_{s =1}^{t-1} c_s \right] a_{ij}  b_j(O_t)} = \frac{\alpha_t(j)}{\displaystyle \sum_{j=1}^N \alpha_t(j)} \,.
\end{equation}

Per quanto riguarda $ \beta_t(i) $ si ha analogamente che
\begin{equation}
\hat{\beta}_t(i) = c_t \beta_t(i)
\end{equation}
con la particolarità che i fattori di scala sono gli stessi che per le $\alpha$.

Con queste nuove variabili la (\ref{eq:a_bar esplicita}) diventa
\begin{equation}
\bar{a}_{ij} = \frac{\displaystyle \sum_{t=1}^{T-1} \hat{\alpha}_t(i) a_{ij} b_j(O_{t+1}) \hat{\beta}_{t+1}(j)}
{\displaystyle \sum_{t=1}^{T-1} \sum_{j=1}^N \hat{\alpha}_t(i) a_{ij} b_j(O_{t+1}) \hat{\beta}_{t+1}(j)}
\end{equation}
e dato che le $\hat{\alpha}$ e le $\hat{\beta}$ si possono scrivere rispettivamente anche come
\begin{gather}
\hat{\alpha}_t(i) = \left[ \prod_{s=1}^t c_s \right] \alpha_t(i) =C_t \alpha_t(i) \\
\hat{\beta}_{t+1}(j) = \left[ \prod_{s=t+1}^T c_s \right] \beta_{t+1}(j) = D_{t+1} \beta_{t+1}(j)
\end{gather}
segue che 
\begin{equation}
\label{eq:a_bar finale}
\bar{a}_{ij} = \frac{\displaystyle \sum_{t=1}^{T-1} C_t \alpha_t(i) a_{ij} b_j(O_{t+1}) D_{t+1} \beta_{t+1}(j)}
{\displaystyle \sum_{t=1}^{T-1} \sum_{j=1}^N C_t \alpha_t(i) a_{ij} b_j(O_{t+1}) D_{t+1} \beta_{t+1}(j)} \,.
\end{equation}
Infine la quantità $C_t D_{t+1}$ risulta indipendente da $t$ infatti
\begin{equation}
C_t D_{t+1} = \prod_{s=1}^t c_s \prod_{s=t+1}^T c_s = \prod_{s=1}^T c_s = C_T
\end{equation}
e si può dunque elidere nella formula (\ref{eq:a_bar finale}) ottenendo proprio la (\ref{eq:ri-stima per a}).

Si è deciso di dare in ingresso al posto della matrice $A$, la stessa $A$ con un leggero rumore additivo al fine di non avere dei balori nulli in nessuna posizione: in pratica, si è sommato una matrice $N$x$N$ piena di $\delta = 10^{-6}$ e poi rinormalizzata per renderla ancora stocastica per righe.

Una possibile implementazione di questo algoritmo è suggerita in \cite{Stamp:implementazione_baum}, mentre  in \cite{Mann:implementazione_baum-logspace} è proposta una soluzione che utilizza non fattori di scala ma logaritmi di probabilità, sempre con lo scopo di avere qualcosa di numericamente gestibile.
 








