Baum–Welch algorithm

Template:Short description In electrical engineering, statistical computing and bioinformatics, the Baum–Welch algorithm is a special case of the expectation–maximization algorithm used to find the unknown parameters of a hidden Markov model (HMM). It makes use of the forward-backward algorithm to compute the statistics for the expectation step. The Baum–Welch algorithm, the primary method for inference in hidden Markov models, is numerically unstable due to its recursive calculation of joint probabilities. As the number of variables grows, these joint probabilities become increasingly small, leading to the forward recursions rapidly approaching values below machine precision.<ref>{{#invoke:citation/CS1|citation |CitationClass=web }}</ref>

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.<ref>{{#invoke:citation/CS1|citation |CitationClass=web }}</ref> One of the first major applications of HMMs was to the field of speech processing.<ref>Template:Cite journal</ref> In the 1980s, HMMs were emerging as a useful tool in the analysis of biological systems and information, and in particular genetic information.<ref>Template:Cite journal</ref> They have since become an important tool in the probabilistic modeling of genomic sequences.<ref name="Durbin1998">Template:Cite book</ref>

Description

A hidden Markov model describes the joint probability of a collection of "hidden" and observed discrete random variables. It relies on the assumption that the i-th hidden variable given the (i − 1)-th hidden variable is independent of previous hidden variables, and the current observation variables depend only on the current hidden state.

The Baum–Welch algorithm uses the well known EM algorithm to find the maximum likelihood estimate of the parameters of a hidden Markov model given a set of observed feature vectors.

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

<math>A=\{a_{ij}\}=P(X_t=j\mid X_{t-1}=i).</math>

The initial state distribution (i.e. when <math>t=1</math>) is given by

<math>\pi_i=P(X_1 = i).</math>

The observation variables <math>Y_t</math> can take one of <math>K</math> possible values. We also assume the observation given the "hidden" state is time independent. The probability of a certain observation <math>y_i</math> at time <math>t</math> for state <math>X_t = j</math> is given by

<math>b_j(y_i)=P(Y_t=y_i\mid X_t=j).</math>

Taking into account all the possible values of <math>Y_t</math> and <math>X_t</math>, we obtain the <math>N \times K</math> matrix <math>B=\{b_j(y_i)\}</math> where <math>b_j</math> belongs to all the possible states and <math>y_i</math> belongs to all the observations.

An observation sequence is given by <math>Y= (Y_1=y_1,Y_2=y_2,\ldots,Y_T=y_T)</math>.

Thus we can describe a hidden Markov chain by <math>\theta = (A,B,\pi)</math>. The Baum–Welch algorithm finds a local maximum for <math>\theta^* = \operatorname{arg\,max}_\theta P(Y\mid\theta)</math> (i.e. the HMM parameters <math>\theta</math> that maximize the probability of the observation).<ref>Template:Cite book</ref>

Algorithm

Set <math>\theta = (A, B, \pi)</math> 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 <math>\alpha_i(t)=P(Y_1=y_1,\ldots,Y_t=y_t,X_t=i\mid\theta)</math>, the probability of seeing the observations <math>y_1,y_2,\ldots,y_t</math> and being in state <math>i</math> at time <math>t</math>. This is found recursively:

1. <math>\alpha_i(1)=\pi_i b_i(y_1),</math>
2. <math>\alpha_i(t+1)=b_i(y_{t+1}) \sum_{j=1}^N \alpha_j(t) a_{ji}.</math>

Since this series converges exponentially to zero, the algorithm will numerically underflow for longer sequences.<ref>{{#invoke:citation/CS1|citation |CitationClass=web }}</ref> However, this can be avoided in a slightly modified algorithm by scaling <math>\alpha</math> in the forward and <math>\beta</math> in the backward procedure below.

Backward procedure

Let <math>\beta_i(t)=P(Y_{t+1}=y_{t+1},\ldots,Y_T=y_{T}\mid X_t=i,\theta)</math> that is the probability of the ending partial sequence <math>y_{t+1},\ldots,y_T</math> given starting state <math>i</math> at time <math>t</math>. We calculate <math>\beta_i(t)</math> as,

1. <math>\beta_i(T)=1,</math>
2. <math>\beta_i(t)=\sum_{j=1}^N \beta_j(t+1) a_{ij} b_j(y_{t+1}).</math>

Update

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

<math>\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)},</math>

which is the probability of being in state <math>i</math> at time <math>t</math> given the observed sequence <math>Y</math> and the parameters <math>\theta</math>

<math>\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}) }, </math>

which is the probability of being in state <math>i</math> and <math>j</math> at times <math>t</math> and <math>t+1</math> respectively given the observed sequence <math>Y</math> and parameters <math>\theta</math>.

The denominators of <math>\gamma_i(t)</math> and <math>\xi_{ij}(t)</math> are the same ; they represent the probability of making the observation <math>Y</math> given the parameters <math>\theta</math>.

The parameters of the hidden Markov model <math>\theta</math> can now be updated:

• <math>\pi_i^* = \gamma_i(1),</math>

which is the expected frequency spent in state <math>i</math> at time <math>1</math>.

• <math>a_{ij}^*=\frac{\sum^{T-1}_{t=1}\xi_{ij}(t)}{\sum^{T-1}_{t=1}\gamma_i(t)},</math>

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.

• <math>b_i^*(v_k)=\frac{\sum^T_{t=1} 1_{y_t=v_k} \gamma_i(t)}{\sum^T_{t=1} \gamma_i(t)},</math>

where

<math>

1_{y_t=v_k}= \begin{cases} 1 & \text{if } y_t=v_k,\\ 0 & \text{otherwise} \end{cases} </math> is an indicator function, and <math>b_i^*(v_k)</math> is the expected number of times the output observations have been equal to <math>v_k</math> while in state <math>i</math> over the expected total number of times in state <math>i</math>.

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, <math>P(Y\mid\theta_\text{final}) > P(Y \mid \theta_\text{true}) </math>. The algorithm also does not guarantee a global maximum.

Multiple sequences

The algorithm described thus far assumes a single observed sequence <math>Y = y_1, \ldots, y_T</math>. However, in many situations, there are several sequences observed: <math>Y_1, \ldots, Y_R</math>. In this case, the information from all of the observed sequences must be used in the update of the parameters <math>A</math>, <math>\pi</math>, and <math>b</math>. Assuming that you have computed <math>\gamma_{ir}(t)</math> and <math>\xi_{ijr}(t)</math> for each sequence <math>y_{1,r},\ldots,y_{N_r,r}</math>, the parameters can now be updated:

• <math>\pi_i^* = \frac{\sum_{r=1}^{R}\gamma_{ir}(1)}{R}</math>
• <math>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)},</math>
• <math>b_i^*(v_k)=\frac{\sum_{r=1}^{R} \sum^T_{t=1} 1_{y_{tr}=v_k} \gamma_{ir}(t)}{\sum_{r=1}^{R} \sum^T_{t=1} \gamma_{ir}(t)},</math>

where

<math>

1_{y_{tr}=v_k}= \begin{cases} 1 & \text{if } y_{t,r}=v_k,\\ 0 & \text{otherwise} \end{cases} </math> 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.<ref>{{#invoke:citation/CS1|citation |CitationClass=web }}</ref><ref>{{#invoke:citation/CS1|citation |CitationClass=web }}</ref> To start we first guess the transition and emission matrices.

Transition State 1 State 2 State 1 0.5 0.5 State 2 0.3 0.7

Emission No Eggs Eggs State 1 0.3 0.7 State 2 0.8 0.2

Initial State 1 0.2 State 2 0.8

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 Template:Tmath then Template:Tmath is given by the following, <math>P(S_1) \cdot P(N|S_1) \cdot P(S_1 \rightarrow S_2) \cdot P(N|S_2).</math>

Observed sequence	Highest probability of observing that sequence if state is Template:Tmath then Template:Tmath	Highest Probability of observing that sequence
NN	0.024 = 0.2 × 0.3 × 0.5 × 0.8	0.3584	Template:Tmath,Template:Tmath
NN	0.024 = 0.2 × 0.3 × 0.5 × 0.8	0.3584	Template:Tmath,Template:Tmath
NN	0.024 = 0.2 × 0.3 × 0.5 × 0.8	0.3584	Template:Tmath,Template:Tmath
NN	0.024 = 0.2 × 0.3 × 0.5 × 0.8	0.3584	Template:Tmath,Template:Tmath
NE	0.006 = 0.2 × 0.3 × 0.5 × 0.2	0.1344	Template:Tmath,Template:Tmath
EE	0.014 = 0.2 × 0.7 × 0.5 × 0.2	0.0490	Template:Tmath,Template:Tmath
EN	0.056 = 0.2 × 0.7 × 0.5 × 0.8	0.0896	Template:Tmath,Template:Tmath
NN	0.024 = 0.2 × 0.3 × 0.5 × 0.8	0.3584	Template:Tmath,Template:Tmath
NN	0.024 = 0.2 × 0.3 × 0.5 × 0.8	0.3584	Template:Tmath,Template:Tmath
Total	0.22	2.4234

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

Old Transition Matrix State 1 State 2 State 1 0.5 0.5 State 2 0.3 0.7

New Transition Matrix (Pseudo Probabilities) State 1 State 2 State 1 0.0598 0.0908 State 2 0.2179 0.9705

New Transition Matrix (After Normalization) State 1 State 2 State 1 0.3973 0.6027 State 2 0.1833 0.8167

Next, we want to estimate a new emission matrix,

Observed Sequence	Highest probability of observing that sequence if E is assumed to come from Template:Tmath		Highest Probability of observing that sequence
NE	0.1344	Template:Tmath,Template:Tmath	0.1344	Template:Tmath,Template:Tmath
EE	0.0490	Template:Tmath,Template:Tmath	0.0490	Template:Tmath,Template:Tmath
EN	0.0560	Template:Tmath,Template:Tmath	0.0896	Template:Tmath,Template:Tmath
Total	0.2394		0.2730

The new estimate for the E coming from Template:Tmath emission is now <math>\frac{0.2394}{0.2730}=0.8769</math>.

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 Template:Tmath and for if N and E came from Template:Tmath and normalize.

Old Emission Matrix No Eggs Eggs State 1 0.3 0.7 State 2 0.8 0.2

New Emission Matrix (Estimates) No Eggs Eggs State 1 0.0404 0.8769 State 2 1.0000 0.7385

New Emission Matrix (After Normalization) No Eggs Eggs State 1 0.0441 0.9559 State 2 0.5752 0.4248

To estimate the initial probabilities we assume all sequences start with the hidden state Template:Tmath and calculate the highest probability and then repeat for Template:Tmath. 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 James K. Baker in 1975.<ref>Template:Cite journal</ref> Continuous speech recognition occurs by the following steps, modeled by a HMM. Feature analysis is first undertaken on temporal and/or spectral features of the speech signal. This produces an observation vector. The feature is then compared to all sequences of the speech recognition units. These units could be phonemes, syllables, or whole-word units. A lexicon decoding system is applied to constrain the paths investigated, so only words in the system's lexicon (word dictionary) are investigated. Similar to the lexicon decoding, the system path is further constrained by the rules of grammar and syntax. Finally, semantic analysis is applied and the system outputs the recognized utterance. A limitation of many HMM applications to speech recognition is that the current state only depends on the state at the previous time-step, which is unrealistic for speech as dependencies are often several time-steps in duration.<ref>Template:Cite journal</ref> The Baum–Welch algorithm also has extensive applications in solving HMMs used in the field of speech synthesis.<ref>Template:Cite journal</ref>

Cryptanalysis

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.<ref>Template:Cite journal</ref> HMMs and as a consequence the Baum–Welch algorithm have also been used to identify spoken phrases in encrypted VoIP calls.<ref>Template:Cite journal</ref> 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.<ref>Template:Cite book</ref>

Applications in bioinformatics

Finding genes

Prokaryotic

The GLIMMER (Gene Locator and Interpolated Markov ModelER) software was an early gene-finding program used for the identification of coding regions in prokaryotic DNA.<ref name="GLIMMER paper">Template:Cite journal</ref><ref name="GLIMMER web">{{#invoke:citation/CS1|citation |CitationClass=web }}</ref> GLIMMER uses Interpolated Markov Models (IMMs) to identify the coding regions and distinguish them from the noncoding DNA. The latest release (GLIMMER3) has been shown to have increased specificity and accuracy compared with its predecessors with regard to predicting translation initiation sites, demonstrating an average 99% accuracy in locating 3' locations compared to confirmed genes in prokaryotes.<ref>Template:Cite journal</ref>

Eukaryotic

The GENSCAN webserver is a gene locator capable of analyzing eukaryotic sequences up to one million base-pairs (1 Mbp) long.<ref>{{#invoke:citation/CS1|citation |CitationClass=web }}</ref> GENSCAN utilizes a general inhomogeneous, three periodic, fifth order Markov model of DNA coding regions. Additionally, this model accounts for differences in gene density and structure (such as intron lengths) that occur in different isochores. While most integrated gene-finding software (at the time of GENSCANs release) assumed input sequences contained exactly one gene, GENSCAN solves a general case where partial, complete, or multiple genes (or even no gene at all) is present.<ref>Template:Cite journal</ref> GENSCAN was shown to exactly predict exon location with 90% accuracy with 80% specificity compared to an annotated database.<ref>Template:Cite journal</ref>

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.<ref>Template:Cite journal</ref> 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.

Implementations

• Accord.NET in C#
• ghmm C library with Python bindings that supports both discrete and continuous emissions.
• hmmlearn Python library that implements Baum-Welch on various discrete-time HMMs
• 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

External links

• A comprehensive review of HMM methods and software in bioinformatics – Profile Hidden Markov Models
• Early HMM publications by Baum:
  • A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains
  • An inequality with applications to statistical estimation for probabilistic functions of Markov processes and to a model for ecology
  • Statistical Inference for Probabilistic Functions of Finite State Markov Chains
• The Shannon Lecture by Welch, which speaks to how the algorithm can be implemented efficiently:
  • Hidden Markov Models and the Baum–Welch Algorithm, IEEE Information Theory Society Newsletter, Dec. 2003.
• An alternative to the Baum–Welch algorithm, the Viterbi Path Counting algorithm:
  • Davis, Richard I. A.; Lovell, Brian C.; "Comparing and evaluating HMM ensemble training algorithms using train and test and condition number criteria", Pattern Analysis and Applications, vol. 6, no. 4, pp. 327–336, 2003.
• An Interactive Spreadsheet for Teaching the Forward-Backward Algorithm (spreadsheet and article with step-by-step walkthrough)
• Formal derivation of the Baum–Welch algorithm Template:Webarchive
• Implementation of the Baum–Welch algorithm