Algoritmo de Baum-Welch 📖 Wikipedia

Na engenharia elétrica, computação estatística e bioinformática, o algoritmo de Baum–Welch é um caso especial do algoritmo de maximização da expectativa (EM) usado para encontrar os parâmetros desconhecidos de um modelo oculto de Markov (HMM). Ele faz uso do algoritmo forward-backward para calcular as estatísticas para a etapa de expectativa. O algoritmo de Baum–Welch, o método principal para inferência em modelos ocultos de Markov, é numericamente instável devido ao seu cálculo recursivo de probabilidades conjuntas. À medida que o número de variáveis cresce, essas probabilidades conjuntas tornam-se cada vez menores, levando as recursões forward a se aproximarem rapidamente de valores abaixo da precisão da máquina.^[1]

O algoritmo de Baum–Welch recebeu esse nome em homenagem aos seus inventores Leonard E. Baum e Lloyd R. Welch. O algoritmo e os Modelos Ocultos de Markov foram descritos pela primeira vez em uma série de artigos por Baum e seus colegas no Centro de Pesquisa em Comunicações do IDA, Princeton no final dos anos 1960 e início dos anos 1970.^[2] Uma das primeiras grandes aplicações dos HMMs foi no campo do processamento de fala.^[3] Na década de 1980, os HMMs estavam surgindo como uma ferramenta útil na análise de sistemas biológicos e informações, e em particular informação genética.^[4] Desde então, eles se tornaram uma ferramenta importante na modelagem probabilística de sequências genômicas.^[5]

Um modelo oculto de Markov descreve a probabilidade conjunta de uma coleção de variáveis aleatórias discretas "ocultas" e observadas. Ele se baseia na suposição de que a i-ésima variável oculta dada a (i − 1)-ésima variável oculta é independente das variáveis ocultas anteriores, e as variáveis de observação atuais dependem apenas do estado oculto atual.

O algoritmo de Baum–Welch usa o conhecido algoritmo EM para encontrar a estimativa de máxima verossimilhança dos parâmetros de um modelo oculto de Markov, dado um conjunto de vetores de características observados.

Seja ${\displaystyle X_{t}}$ uma variável aleatória discreta oculta com ${\displaystyle N}$ valores possíveis (ou seja, assumimos que existem ${\displaystyle N}$ estados no total). Assumimos que o ${\displaystyle P(X_{t}\mid X_{t-1})}$ é independente do tempo ${\displaystyle t}$, o que leva à definição da matriz de transição estocástica independente do tempo

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

A distribuição de estado inicial (ou seja, quando ${\displaystyle t=1}$) é dada por

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

As variáveis de observação ${\displaystyle Y_{t}}$ podem assumir um dos ${\displaystyle K}$ valores possíveis. Também assumimos que a observação dada o estado "oculto" é independente do tempo. A probabilidade de uma certa observação ${\displaystyle y_{i}}$ no tempo ${\displaystyle t}$ para o estado ${\displaystyle X_{t}=j}$ é dada por

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

Levando em conta todos os valores possíveis de ${\displaystyle Y_{t}}$ e ${\displaystyle X_{t}}$, obtemos a matriz ${\displaystyle N\times K}$ ${\displaystyle B=\{b_{j}(y_{i})\}}$ onde ${\displaystyle b_{j}}$ pertence a todos os estados possíveis e ${\displaystyle y_{i}}$ pertence a todas as observações.

Uma sequência de observação é dada por ${\displaystyle Y=(Y_{1}=y_{1},Y_{2}=y_{2},\ldots ,Y_{T}=y_{T})}$.

Assim, podemos descrever uma cadeia de Markov oculta por ${\displaystyle \theta =(A,B,\pi )}$. O algoritmo de Baum–Welch encontra um máximo local para ${\displaystyle \theta ^{*}=\operatorname {arg\,max} _{\theta }P(Y\mid \theta )}$ (ou seja, os parâmetros do HMM ${\displaystyle \theta }$ que maximizam a probabilidade da observação).^[6]

Defina ${\displaystyle \theta =(A,B,\pi )}$ com condições iniciais aleatórias. Elas também podem ser definidas usando informações prévias sobre os parâmetros, se disponíveis; isso pode acelerar o algoritmo e também direcioná-lo para o máximo local desejado.

Seja ${\displaystyle \alpha _{i}(t)=P(Y_{1}=y_{1},\ldots ,Y_{t}=y_{t},X_{t}=i\mid \theta )}$, a probabilidade de ver as observações ${\displaystyle y_{1},y_{2},\ldots ,y_{t}}$ e estar no estado ${\displaystyle i}$ no tempo ${\displaystyle t}$. Isso é encontrado recursivamente:

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}.}$

Como esta série converge exponencialmente para zero, o algoritmo sofrerá underflow numérico para sequências mais longas.^[7] No entanto, isso pode ser evitado em um algoritmo ligeiramente modificado escalando ${\displaystyle \alpha }$ no procedimento forward e ${\displaystyle \beta }$ no procedimento backward abaixo.

Seja ${\displaystyle \beta _{i}(t)=P(Y_{t+1}=y_{t+1},\ldots ,Y_{T}=y_{T}\mid X_{t}=i,\theta )}$ que é a probabilidade do final da sequência parcial ${\displaystyle y_{t+1},\ldots ,y_{T}}$ dado o estado inicial ${\displaystyle i}$ no tempo ${\displaystyle t}$. Calculamos ${\displaystyle \beta _{i}(t)}$ como:

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}).}$

Podemos agora calcular as variáveis temporárias, de acordo com o Teorema de Bayes:

${\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)}},}$

que é a probabilidade de estar no estado ${\displaystyle i}$ no tempo ${\displaystyle t}$ dada a sequência observada ${\displaystyle Y}$ e os parâmetros ${\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})}},}$

que é a probabilidade de estar no estado ${\displaystyle i}$ e ${\displaystyle j}$ nos tempos ${\displaystyle t}$ e ${\displaystyle t+1}$ respectivamente dada a sequência observada ${\displaystyle Y}$ e os parâmetros ${\displaystyle \theta }$.

Os denominadores de ${\displaystyle \gamma _{i}(t)}$ e ${\displaystyle \xi _{ij}(t)}$ são os mesmos; eles representam a probabilidade de fazer a observação ${\displaystyle Y}$ dados os parâmetros ${\displaystyle \theta }$.

Os parâmetros do modelo oculto de Markov ${\displaystyle \theta }$ podem agora ser atualizados:

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

que é a frequência esperada gasta no estado ${\displaystyle i}$ no tempo ${\displaystyle 1}$.

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

que é o número esperado de transições do estado i para o estado j comparado ao número total esperado de transições começando no estado i, incluindo do estado i para ele mesmo. O número de transições começando no estado i é equivalente ao número de vezes que o estado i é observado na sequência de t = 1 a 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)}},}$

onde

${\displaystyle 1_{y_{t}=v_{k}}={\begin{cases}1&{\text{se }}y_{t}=v_{k},\\0&{\text{caso contrário}}\end{cases}}}$

é uma função indicadora, e ${\displaystyle b_{i}^{*}(v_{k})}$ é o número esperado de vezes que as observações de saída foram iguais a ${\displaystyle v_{k}}$ enquanto no estado ${\displaystyle i}$ sobre o número total esperado de vezes no estado ${\displaystyle i}$.

Essas etapas são agora repetidas iterativamente até um nível desejado de convergência.

Nota: É possível ocorrer o sobreajuste (over-fit) em um conjunto de dados específico. Ou seja, ${\displaystyle P(Y\mid \theta _{\text{final}})>P(Y\mid \theta _{\text{true}})}$. O algoritmo também não garante um máximo global.

O algoritmo descrito até agora assume uma única sequência observada ${\displaystyle Y=y_{1},\ldots ,y_{T}}$. No entanto, em muitas situações, existem várias sequências observadas: ${\displaystyle Y_{1},\ldots ,Y_{R}}$. Neste caso, as informações de todas as sequências observadas devem ser usadas na atualização dos parâmetros ${\displaystyle A}$, ${\displaystyle \pi }$ e ${\displaystyle b}$. Assumindo que você calculou ${\displaystyle \gamma _{ir}(t)}$ e ${\displaystyle \xi _{ijr}(t)}$ para cada sequência ${\displaystyle y_{1,r},\ldots ,y_{N_{r},r}}$, os parâmetros podem agora ser atualizados:

• ${\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)}},}$

onde

${\displaystyle 1_{y_{tr}=v_{k}}={\begin{cases}1&{\text{se }}y_{t,r}=v_{k},\\0&{\text{caso contrário}}\end{cases}}}$

é uma função indicadora

Suponha que temos uma galinha da qual coletamos ovos ao meio-dia todos os dias. Agora, se a galinha botou ovos para coleta ou não, depende de alguns fatores desconhecidos que estão ocultos. Podemos, no entanto (por simplicidade) assumir que a galinha está sempre em um de dois estados que influenciam se a galinha bota ovos, e que este estado depende apenas do estado no dia anterior. Agora, não sabemos o estado no ponto de partida inicial, não sabemos as probabilidades de transição entre os dois estados e não sabemos a probabilidade de a galinha botar um ovo dado um estado específico.^[8][9] Para começar, primeiro deduzimos as matrizes de transição e emissão.

Então, pegamos um conjunto de observações (E = ovos, N = sem ovos): N, N, N, N, N, E, E, N, N, N

Isso nos dá um conjunto de transições observadas entre os dias: NN, NN, NN, NN, NE, EE, EN, NN, NN

O próximo passo é estimar uma nova matriz de transição. Por exemplo, a probabilidade da sequência NN e do estado ser ⁠${\displaystyle S_{1}}$⁠ e depois ⁠${\displaystyle S_{2}}$⁠ é dada pelo seguinte: ${\displaystyle P(S_{1})\cdot P(N|S_{1})\cdot P(S_{1}\rightarrow S_{2})\cdot P(N|S_{2}).}$

Sequência observada	Probabilidade mais alta de observar essa sequência se o estado for ⁠${\displaystyle S_{1}}$⁠ e depois ⁠${\displaystyle S_{2}}$⁠	Probabilidade mais alta de observar essa sequência
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

Assim, a nova estimativa para a transição de ⁠${\displaystyle S_{1}}$⁠ para ⁠${\displaystyle S_{2}}$⁠ é agora ${\displaystyle {\frac {0.22}{2.4234}}=0.0908}$ (referida como "Pseudo probabilidades" nas tabelas a seguir). Em seguida, calculamos as probabilidades de transição de ⁠${\displaystyle S_{2}}$⁠ para ⁠${\displaystyle S_{1}}$⁠, de ⁠${\displaystyle S_{2}}$⁠ para ⁠${\displaystyle S_{2}}$⁠ e de ⁠${\displaystyle S_{1}}$⁠ para ⁠${\displaystyle S_{1}}$⁠ e normalizamos cada linha da matriz de transição para que as probabilidades de transições de um determinado estado inicial somem 1. Isso nos dá a matriz de transição atualizada:

Em seguida, queremos estimar uma nova matriz de emissão,

Sequência Observada	Probabilidade mais alta de observar essa sequência se for assumido que E vem de ⁠${\displaystyle S_{1}}$⁠		Probabilidade mais alta de observar essa sequência
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

A nova estimativa para a emissão de E vindo de ⁠${\displaystyle S_{1}}$⁠ agora é ${\displaystyle {\frac {0.2394}{0.2730}}=0.8769}$.

Isso nos permite calcular a matriz de emissão conforme descrito acima no algoritmo, somando as probabilidades para as respectivas sequências observadas. Em seguida, repetimos se N veio de ⁠${\displaystyle S_{1}}$⁠ e para se N e E vieram de ⁠${\displaystyle S_{2}}$⁠ e normalizamos.

Para estimar as probabilidades iniciais, assumimos que todas as sequências começam com o estado oculto ⁠${\displaystyle S_{1}}$⁠ e calculamos a probabilidade mais alta e depois repetimos para ⁠${\displaystyle S_{2}}$⁠. Novamente, em seguida, normalizamos para fornecer um vetor inicial atualizado.

Finalmente repetimos esses passos até que as probabilidades resultantes convirjam satisfatoriamente.

Modelos Ocultos de Markov foram aplicados pela primeira vez ao reconhecimento de fala por James K. Baker em 1975.^[10] O reconhecimento de fala contínua ocorre pelos seguintes passos, modelados por um HMM. A análise de características é inicialmente realizada em características temporais e/ou espectrais do sinal de fala. Isso produz um vetor de observação. A característica é então comparada a todas as sequências das unidades de reconhecimento de fala. Essas unidades podem ser fonemas, sílabas ou unidades de palavras inteiras. Um sistema de decodificação de léxico é aplicado para restringir os caminhos investigados, de modo que apenas palavras no léxico (dicionário de palavras) do sistema sejam investigadas. Semelhante à decodificação do léxico, o caminho do sistema é ainda mais restrito pelas regras de gramática e sintaxe. Finalmente, a análise semântica é aplicada e o sistema emite o enunciado reconhecido. Uma limitação de muitas aplicações de HMM ao reconhecimento de fala é que o estado atual depende apenas do estado no passo de tempo anterior, o que não é realista para a fala, pois as dependências geralmente têm vários passos de tempo em duração.^[11] O algoritmo de Baum-Welch também tem extensas aplicações na resolução de HMMs usados no campo da síntese de fala.^[12]

O algoritmo de Baum-Welch é frequentemente usado para estimar os parâmetros de HMMs na decifração de informações ocultas ou com ruído e, consequentemente, é frequentemente usado em criptanálise. Na segurança de dados, um observador gostaria de extrair informações de um fluxo de dados sem conhecer todos os parâmetros da transmissão. Isso pode envolver engenharia reversa de um codificador de canal.^[13] HMMs e, como consequência, o algoritmo de Baum-Welch também foram usados para identificar frases faladas em chamadas VoIP criptografadas.^[14] Além disso, a criptanálise HMM é uma ferramenta importante para investigações automatizadas de dados de tempo de cache. Ela permite a descoberta automática do estado crítico do algoritmo, por exemplo, valores de chave.^[15]

O software GLIMMER (Gene Locator and Interpolated Markov ModelER) foi um programa inicial de localização de genes usado para a identificação de regiões codificantes em DNA procarioto.^[16][17] O GLIMMER usa Modelos de Markov Interpolados (IMMs) para identificar as regiões codificantes e distingui-las do DNA não-codificante. A versão mais recente (GLIMMER3) demonstrou ter maior especificidade e precisão em comparação com os seus antecessores no que diz respeito à previsão dos locais de início da tradução, demonstrando uma precisão média de 99% na localização de pontos 3' em comparação com genes confirmados em procariotos.^[18]

O servidor web GENSCAN é um localizador de genes capaz de analisar sequências eucariotas de até um milhão de pares de bases (1 Mbp) de comprimento.^[19] O GENSCAN utiliza um modelo de Markov geral não homogêneo, de período três e de quinta ordem de regiões de codificação de DNA. Além disso, esse modelo é responsável pelas diferenças na densidade e estrutura dos genes (como o comprimento dos íntrons) que ocorrem em diferentes isocoros. Enquanto a maioria dos softwares integrados de descoberta de genes (na época do lançamento do GENSCAN) presumia que as sequências de entrada contivessem exatamente um gene, o GENSCAN resolve um caso geral onde genes parciais, completos ou múltiplos (ou até mesmo nenhum gene) estão presentes.^[20] O GENSCAN demonstrou prever exatamente a localização dos éxons com 90% de precisão com 80% de especificidade em comparação a um banco de dados anotado.^[21]

As Variações do número de cópias (CNVs) são uma forma abundante de variação da estrutura do genoma em humanos. Um HMM bivariado de valor discreto (dbHMM) foi usado para atribuir regiões cromossômicas a sete estados distintos: regiões não afetadas, deleções, duplicações e quatro estados de transição. A solução deste modelo usando Baum-Welch demonstrou a capacidade de prever a localização do ponto de quebra de CNV em aproximadamente 300 pb a partir de experimentos de micro-arranjos.^[22] Essa magnitude de resolução permite correlações mais precisas entre diferentes CNVs e através de populações do que as anteriormente possíveis, permitindo o estudo das frequências populacionais das CNV. Também demonstrou um padrão de herança direta para uma CNV específica.

• Accord.NET em C#
• ghmm biblioteca em C com ligações para Python que suporta emissões discretas e contínuas.
• hmmlearn biblioteca para Python que implementa o Baum-Welch em vários HMMs de tempo discreto.
• Jajapy biblioteca para Python que implementa o Baum-Welch em vários tipos de Modelos de Markov ( HMM, MC, MDP, CTMC).
• HiddenMarkovModels.jl pacote para Julia.
• Função HMMFit no pacote RHmm para R.
• hmmtrain em MATLAB
• rustbio em Rust

• Algoritmo de Viterbi
• Modelo oculto de Markov
• Algoritmo EM
• Máxima verossimilhança
• Reconhecimento de fala
• Bioinformática

• Uma revisão abrangente de métodos HMM e software em bioinformática – Profile Hidden Markov Models
• Primeiras publicações de HMM por 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
• A Palestra Shannon de Welch, que fala sobre como o algoritmo pode ser implementado de forma eficiente:
  • Hidden Markov Models and the Baum–Welch Algorithm, IEEE Information Theory Society Newsletter, dez. de 2003.
• Uma alternativa ao algoritmo de Baum-Welch, o algoritmo de contagem de caminhos de Viterbi:
  • 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, n. 4, pp. 327–336, 2003.
• An Interactive Spreadsheet for Teaching the Forward-Backward Algorithm (planilha e artigo com passo a passo)
• Formal derivation of the Baum–Welch algorithm Arquivado em 2012-02-28 no Wayback Machine
• Implementation of the Baum–Welch algorithm