Baum–Welch algorithm

In

History

The Baum–Welch algorithm was named after its inventors Leonard E. Baum and Lloyd R. Welch. The algorithm and the Hidden Markov models were first described in a series of articles by Baum and his peers at the IDA Center for Communications Research, Princeton in the late 1960s and early 1970s.^[1] One of the first major applications of HMMs was to the field of speech processing.^[2] In the 1980s, HMMs were emerging as a useful tool in the analysis of biological systems and information, and in particular genetic information.^[3] They have since become an important tool in the probabilistic modeling of genomic sequences.^[4]

Description

A

The Baum–Welch algorithm uses the well known EM algorithm to find the

Let ${\displaystyle X_{t}}$ be a discrete hidden random variable with ${\displaystyle N}$ possible values (i.e. We assume there are ${\displaystyle N}$ states in total). We assume the ${\displaystyle P(X_{t}\mid X_{t-1})}$ is independent of time ${\displaystyle t}$, which leads to the definition of the time-independent stochastic transition matrix

${\displaystyle A=\{a_{ij}\}=P(X_{t}=j\mid X_{t-1}=i).}$

The initial state distribution (i.e. when ${\displaystyle t=1}$) is given by

${\displaystyle \pi _{i}=P(X_{1}=i).}$

The observation variables ${\displaystyle Y_{t}}$ can take one of ${\displaystyle K}$ possible values. We also assume the observation given the "hidden" state is time independent. The probability of a certain observation ${\displaystyle y_{i}}$ at time ${\displaystyle t}$ for state ${\displaystyle X_{t}=j}$ is given by

${\displaystyle b_{j}(y_{i})=P(Y_{t}=y_{i}\mid X_{t}=j).}$

Taking into account all the possible values of ${\displaystyle Y_{t}}$ and ${\displaystyle X_{t}}$, we obtain the ${\displaystyle N\times K}$ matrix ${\displaystyle B=\{b_{j}(y_{i})\}}$ where ${\displaystyle b_{j}}$ belongs to all the possible states and ${\displaystyle y_{i}}$ belongs to all the observations.

An observation sequence is given by ${\displaystyle Y=(Y_{1}=y_{1},Y_{2}=y_{2},\ldots ,Y_{T}=y_{T})}$.

Thus we can describe a hidden Markov chain by ${\displaystyle \theta =(A,B,\pi )}$. The Baum–Welch algorithm finds a local maximum for ${\displaystyle \theta ^{*}=\operatorname {arg\,max} _{\theta }P(Y\mid \theta )}$ (i.e. the HMM parameters ${\displaystyle \theta }$ that maximize the probability of the observation).^[5]

Algorithm

Set ${\displaystyle \theta =(A,B,\pi )}$ with random initial conditions. They can also be set using prior information about the parameters if it is available; this can speed up the algorithm and also steer it toward the desired local maximum.

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

Example

Suppose we have a chicken from which we collect eggs at noon every day. Now whether or not the chicken has laid eggs for collection depends on some unknown factors that are hidden. We can however (for simplicity) assume that the chicken is always in one of two states that influence whether the chicken lays eggs, and that this state only depends on the state on the previous day. Now we don't know the state at the initial starting point, we don't know the transition probabilities between the two states and we don't know the probability that the chicken lays an egg given a particular state.^[7][8] To start we first guess the transition and emission matrices.

We then take a set of observations (E = eggs, N = no eggs): N, N, N, N, N, E, E, N, N, N

This gives us a set of observed transitions between days: NN, NN, NN, NN, NE, EE, EN, NN, NN

The next step is to estimate a new transition matrix. For example, the probability of the sequence NN and the state being ${\displaystyle S_{1}}$ then ${\displaystyle S_{2}}$ is given by the following, ${\displaystyle P(S_{1})\cdot P(N|S_{1})\cdot P(S_{1}\rightarrow S_{2})\cdot P(N|S_{2}).}$

Observed sequence	Highest probability of observing that sequence if state is ${\displaystyle S_{1}}$ then ${\displaystyle S_{2}}$	Highest Probability of observing that sequence
NN	0.024 = 0.2 × 0.3 × 0.5 × 0.8	0.3584	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
NN	0.024 = 0.2 × 0.3 × 0.5 × 0.8	0.3584	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
NN	0.024 = 0.2 × 0.3 × 0.5 × 0.8	0.3584	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
NN	0.024 = 0.2 × 0.3 × 0.5 × 0.8	0.3584	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
NE	0.006 = 0.2 × 0.3 × 0.5 × 0.2	0.1344	${\displaystyle S_{2}}$,${\displaystyle S_{1}}$
EE	0.014 = 0.2 × 0.7 × 0.5 × 0.2	0.0490	${\displaystyle S_{1}}$,${\displaystyle S_{1}}$
EN	0.056 = 0.2 × 0.7 × 0.5 × 0.8	0.0896	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
NN	0.024 = 0.2 × 0.3 × 0.5 × 0.8	0.3584	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
NN	0.024 = 0.2 × 0.3 × 0.5 × 0.8	0.3584	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
Total	0.22	2.4234

Thus the new estimate for the ${\displaystyle S_{1}}$ to ${\displaystyle S_{2}}$ transition is now ${\displaystyle {\frac {0.22}{2.4234}}=0.0908}$ (referred to as "Pseudo probabilities" in the following tables). We then calculate the ${\displaystyle S_{2}}$ to ${\displaystyle S_{1}}$, ${\displaystyle S_{2}}$ to ${\displaystyle S_{2}}$ and ${\displaystyle S_{1}}$ to ${\displaystyle S_{1}}$ transition probabilities and normalize so they add to 1. This gives us the updated transition matrix:

Next, we want to estimate a new emission matrix,

Observed Sequence	Highest probability of observing that sequence if E is assumed to come from ${\displaystyle S_{1}}$		Highest Probability of observing that sequence
NE	0.1344	${\displaystyle S_{2}}$,${\displaystyle S_{1}}$	0.1344	${\displaystyle S_{2}}$,${\displaystyle S_{1}}$
EE	0.0490	${\displaystyle S_{1}}$,${\displaystyle S_{1}}$	0.0490	${\displaystyle S_{1}}$,${\displaystyle S_{1}}$
EN	0.0560	${\displaystyle S_{1}}$,${\displaystyle S_{2}}$	0.0896	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
Total	0.2394		0.2730

The new estimate for the E coming from ${\displaystyle S_{1}}$ emission is now ${\displaystyle {\frac {0.2394}{0.2730}}=0.8769}$.

This allows us to calculate the emission matrix as described above in the algorithm, by adding up the probabilities for the respective observed sequences. We then repeat for if N came from ${\displaystyle S_{1}}$ and for if N and E came from ${\displaystyle S_{2}}$ and normalize.

To estimate the initial probabilities we assume all sequences start with the hidden state ${\displaystyle S_{1}}$ and calculate the highest probability and then repeat for ${\displaystyle S_{2}}$. Again we then normalize to give an updated initial vector.

Finally we repeat these steps until the resulting probabilities converge satisfactorily.

Applications

Speech recognition

Hidden Markov Models were first applied to speech recognition by

The Baum–Welch algorithm is often used to estimate the parameters of HMMs in deciphering hidden or noisy information and consequently is often used in cryptanalysis. In data security an observer would like to extract information from a data stream without knowing all the parameters of the transmission. This can involve reverse engineering a channel encoder.^[12] HMMs and as a consequence the Baum–Welch algorithm have also been used to identify spoken phrases in encrypted VoIP calls.^[13] In addition HMM cryptanalysis is an important tool for automated investigations of cache-timing data. It allows for the automatic discovery of critical algorithm state, for example key values.^[14]

Applications in bioinformatics

Finding genes

Prokaryotic

The

• Accord.NET in C#
• ghmm C library with Python bindings that supports both discrete and continuous emissions.
• Jajapy Python library that implements Baum-Welch on various kind of Markov Models ( HMM, MC, MDP, CTMC).
• HiddenMarkovModels.jl package for Julia.
• HMMFit function in the RHmm package for R.
• hmmtrain in MATLAB
• rustbio in Rust

See also

• Viterbi algorithm
• Hidden Markov model
• EM algorithm
• Maximum likelihood
• Speech recognition
• Bioinformatics
• Cryptanalysis

References