Forward procedure

Let ${\displaystyle \alpha _{i}(t)=P(Y_{1}=y_{1},\ldots ,Y_{t}=y_{t},X_{t}=i\mid \theta )}$, the probability of seeing the observations ${\displaystyle y_{1},y_{2},\ldots ,y_{t}}$ and being in state ${\displaystyle i}$ at time ${\displaystyle t}$. This is found recursively:

1. ${\displaystyle \alpha _{i}(1)=\pi _{i}b_{i}(y_{1}),}$
2. ${\displaystyle \alpha _{i}(t+1)=b_{i}(y_{t+1})\sum _{j=1}^{N}\alpha _{j}(t)a_{ji}.}$

Since this series converges exponentially to zero, the algorithm will numerically underflow for longer sequences.^[6] However, this can be avoided in a slightly modified algorithm by scaling ${\displaystyle \alpha }$ in the forward and ${\displaystyle \beta }$ in the backward procedure below.

Backward procedure

Let ${\displaystyle \beta _{i}(t)=P(Y_{t+1}=y_{t+1},\ldots ,Y_{T}=y_{T}\mid X_{t}=i,\theta )}$ that is the probability of the ending partial sequence ${\displaystyle y_{t+1},\ldots ,y_{T}}$ given starting state ${\displaystyle i}$ at time ${\displaystyle t}$. We calculate ${\displaystyle \beta _{i}(t)}$ as,

1. ${\displaystyle \beta _{i}(T)=1,}$
2. ${\displaystyle \beta _{i}(t)=\sum _{j=1}^{N}\beta _{j}(t+1)a_{ij}b_{j}(y_{t+1}).}$

Update

We can now calculate the temporary variables, according to Bayes' theorem:

${\displaystyle \gamma _{i}(t)=P(X_{t}=i\mid Y,\theta )={\frac {P(X_{t}=i,Y\mid \theta )}{P(Y\mid \theta )}}={\frac {\alpha _{i}(t)\beta _{i}(t)}{\sum _{j=1}^{N}\alpha _{j}(t)\beta _{j}(t)}},}$

which is the probability of being in state ${\displaystyle i}$ at time ${\displaystyle t}$ given the observed sequence ${\displaystyle Y}$ and the parameters ${\displaystyle \theta }$

${\displaystyle \xi _{ij}(t)=P(X_{t}=i,X_{t+1}=j\mid Y,\theta )={\frac {P(X_{t}=i,X_{t+1}=j,Y\mid \theta )}{P(Y\mid \theta )}}={\frac {\alpha _{i}(t)a_{ij}\beta _{j}(t+1)b_{j}(y_{t+1})}{\sum _{k=1}^{N}\sum _{w=1}^{N}\alpha _{k}(t)a_{kw}\beta _{w}(t+1)b_{w}(y_{t+1})}},}$

which is the probability of being in state ${\displaystyle i}$ and ${\displaystyle j}$ at times ${\displaystyle t}$ and ${\displaystyle t+1}$ respectively given the observed sequence ${\displaystyle Y}$ and parameters ${\displaystyle \theta }$.

The denominators of ${\displaystyle \gamma _{i}(t)}$ and ${\displaystyle \xi _{ij}(t)}$ are the same ; they represent the probability of making the observation ${\displaystyle Y}$ given the parameters ${\displaystyle \theta }$.

The parameters of the hidden Markov model ${\displaystyle \theta }$ can now be updated:

• ${\displaystyle \pi _{i}^{*}=\gamma _{i}(1),}$

which is the expected frequency spent in state ${\displaystyle i}$ at time ${\displaystyle 1}$.

• ${\displaystyle a_{ij}^{*}={\frac {\sum _{t=1}^{T-1}\xi _{ij}(t)}{\sum _{t=1}^{T-1}\gamma _{i}(t)}},}$

which is the expected number of transitions from state i to state j compared to the expected total number of transitions away from state i. To clarify, the number of transitions away from state i does not mean transitions to a different state j, but to any state including itself. This is equivalent to the number of times state i is observed in the sequence from t = 1 to t = T − 1.

• ${\displaystyle b_{i}^{*}(v_{k})={\frac {\sum _{t=1}^{T}1_{y_{t}=v_{k}}\gamma _{i}(t)}{\sum _{t=1}^{T}\gamma _{i}(t)}},}$

where

${\displaystyle 1_{y_{t}=v_{k}}={\begin{cases}1&{\text{if }}y_{t}=v_{k},\\0&{\text{otherwise}}\end{cases}}}$

is an indicator function, and ${\displaystyle b_{i}^{*}(v_{k})}$ is the expected number of times the output observations have been equal to ${\displaystyle v_{k}}$ while in state ${\displaystyle i}$ over the expected total number of times in state ${\displaystyle i}$.

These steps are now repeated iteratively until a desired level of convergence.

Note: It is possible to over-fit a particular data set. That is, ${\displaystyle P(Y\mid \theta _{\text{final}})>P(Y\mid \theta _{\text{true}})}$. The algorithm also does not guarantee a global maximum.

Multiple sequences

The algorithm described thus far assumes a single observed sequence ${\displaystyle Y=y_{1},\ldots ,y_{N}}$. However, in many situations, there are several sequences observed: ${\displaystyle Y_{1},\ldots ,Y_{R}}$. In this case, the information from all of the observed sequences must be used in the update of the parameters ${\displaystyle A}$, ${\displaystyle \pi }$, and ${\displaystyle b}$. Assuming that you have computed ${\displaystyle \gamma _{ir}(t)}$ and ${\displaystyle \xi _{ijr}(t)}$ for each sequence ${\displaystyle y_{1,r},\ldots ,y_{N_{r},r}}$, the parameters can now be updated:

• ${\displaystyle \pi _{i}^{*}={\frac {\sum _{r=1}^{R}\gamma _{ir}(1)}{R}}}$
• ${\displaystyle a_{ij}^{*}={\frac {\sum _{r=1}^{R}\sum _{t=1}^{T-1}\xi _{ijr}(t)}{\sum _{r=1}^{R}\sum _{t=1}^{T-1}\gamma _{ir}(t)}},}$
• ${\displaystyle b_{i}^{*}(v_{k})={\frac {\sum _{r=1}^{R}\sum _{t=1}^{T}1_{y_{tr}=v_{k}}\gamma _{ir}(t)}{\sum _{r=1}^{R}\sum _{t=1}^{T}\gamma _{ir}(t)}},}$

where

${\displaystyle 1_{y_{tr}=v_{k}}={\begin{cases}1&{\text{if }}y_{t,r}=v_{k},\\0&{\text{otherwise}}\end{cases}}}$

is an indicator function

Finding genes

Copy-number variation detection

Copy-number variations (CNVs) are an abundant form of genome structure variation in humans. A discrete-valued bivariate HMM (dbHMM) was used assigning chromosomal regions to seven distinct states: unaffected regions, deletions, duplications and four transition states. Solving this model using Baum-Welch demonstrated the ability to predict the location of CNV breakpoint to approximately 300 bp from micro-array experiments.^[21] This magnitude of resolution enables more precise correlations between different CNVs and across populations than previously possible, allowing the study of CNV population frequencies. It also demonstrated a direct inheritance pattern for a particular CNV.