Baum–Welch algoritm

Inom elektroteknik , statistisk beräkning och bioinformatik är Baum -Welch-algoritmen ett specialfall av algoritmen för förväntan-maximering som används för att hitta de okända parametrarna för en dold Markov-modell (HMM). Den använder sig av framåt-bakåt-algoritmen för att beräkna statistiken för förväntningssteget.

Historia

Baum-Welch-algoritmen fick sitt namn efter dess uppfinnare Leonard E. Baum och Lloyd R. Welch . Algoritmen och Hidden Markov-modellerna beskrevs först i en serie artiklar av Baum och hans kamrater vid IDA Center for Communications Research, Princeton i slutet av 1960-talet och början av 1970-talet. En av de första stora tillämpningarna av HMM var till området för talbehandling . På 1980-talet växte HMMs fram som ett användbart verktyg för analys av biologiska system och information, och i synnerhet genetisk information . De har sedan dess blivit ett viktigt verktyg i probabilistisk modellering av genomiska sekvenser.

Beskrivning

En dold Markov-modell beskriver den gemensamma sannolikheten för en samling " dolda " och observerade diskreta slumpvariabler. Den förlitar sig på antagandet att den i -te dolda variabeln givet ( i − 1)-th dolda variabeln är oberoende av tidigare dolda variabler, och de aktuella observationsvariablerna beror endast på det aktuella dolda tillståndet.

Baum–Welch-algoritmen använder den välkända EM-algoritmen för att hitta den maximala sannolikhetsuppskattningen av parametrarna för en dold Markov-modell givet en uppsättning observerade funktionsvektorer.

Låt $X_{t}$ vara en diskret dold slumpvariabel med $N$ möjliga värden (dvs. Vi antar att det finns $N$ tillstånd totalt). Vi antar att $P(X_{t}\mid X_{t-1})$ är oberoende av tiden $t$ , vilket leder till definitionen av tidsoberoende stokastisk övergångsmatris

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

Den initiala tillståndsfördelningen (dvs när $t=1$ ) ges av

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

Observationsvariablerna $Y_{t}$ kan ta ett av $K$ möjliga värden. Vi antar också att observationen givet det "dolda" tillståndet är tidsoberoende. Sannolikheten för en viss observation $y_{i}$ vid tidpunkten $t$ för tillstånd $X_{t}=j$ ges av

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

Med hänsyn till alla möjliga värden för $Y_{t}$ och $X_{t}$ får vi $N\times K$ matrisen $B=\{b_{j}(y_{i})\}$ där $b_{j}$ tillhör alla möjliga tillstånd och $y_{ i}$ tillhör alla observationer.

En observationssekvens ges av $Y=(Y_{1}=y_{1},Y_{2}= y_{2},\ldots ,Y_{T}=y_{T})$ .

Således kan vi beskriva en dold Markov-kedja med $\theta =(A,B,\pi )$ . Baum–Welch-algoritmen hittar ett lokalt maximum för $\theta ^{*}=\operatörsnamn {arg\,max} _{\theta } P(Y\mid \theta )$ (dvs. HMM-parametrarna $\theta$ som maximerar sannolikheten för observationen).

Algoritm

Ange $\theta =(A,B,\pi )$ med slumpmässiga initialvillkor. De kan också ställas in med hjälp av tidigare information om parametrarna om den är tillgänglig; detta kan påskynda algoritmen och även styra den mot det önskade lokala maximumet.

Vidarebefordran procedur

Låt ${\displaystyle \alpha _{i}(t)=P(Y_{1}= y_{1},\ldots ,Y_{t}=y_{t},X_{t}=i\mid \theta )} , sannolikheten att$ se observationerna $y_{1},y_{2},\ldots ,y_{t}$ och är i tillståndet $i$ vid tidpunkten $t$ . Detta hittas rekursivt:

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

Eftersom denna serie konvergerar exponentiellt till noll, kommer algoritmen att sjunka numeriskt för längre sekvenser. Detta kan dock undvikas i en något modifierad algoritm genom att skala $\alpha$ framåt och $\beta$ i bakåtproceduren nedan.

Bakåtprocedur

Låt ${\displaystyle \beta _{i}(t)=P(Y_ {t+1}=y_{t+1},\ldots ,Y_{T}=y_{T}\mid X_{t}=i,\theta )} som är sannolikheten för den slutliga partiella$ sekvensen $y_{t+1},\ldots ,y_{T}$ givet starttillstånd $i$ vid tidpunkten $t$ . Vi beräknar $\beta _{i}(t)$ som,

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

Uppdatering

Vi kan nu beräkna de temporära variablerna, enligt Bayes teorem:

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

vilket är sannolikheten att vara i tillstånd $i$ vid tidpunkten $t$ givet den observerade sekvensen $Y$ och parametrarna $\theta$

\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}\summa _{w=1}^ {N}\alpha _{k}(t)a_{kw}\beta _{w}(t+1)b_{w}(y_{t+1})}},

vilket är sannolikheten att vara i tillståndet $i$ och $j$ vid tidpunkterna $t$ respektive $t+1$ givet den observerade sekvensen $Y$ och parametrarna $\theta$ .

Nämnarna för $\gamma _{i}(t)$ och $\xi _{ij}(t)$ är desamma ; de representerar sannolikheten att göra observationen $Y$ givet parametrarna $\theta$ .

Parametrarna för den dolda Markov-modellen $\theta$ kan nu uppdateras:

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

vilket är den förväntade frekvensen som spenderas i tillstånd $i$ vid tidpunkten $1$ .

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

vilket är det förväntade antalet övergångar från tillstånd i till tillstånd j jämfört med det förväntade totala antalet övergångar bort från tillstånd i . För att förtydliga, betyder antalet övergångar bort från tillstånd i inte övergångar till ett annat tillstånd j , utan till vilket tillstånd som helst inklusive sig själv. Detta motsvarar antalet gånger tillstånd i observeras i sekvensen från t = 1 till t = T − 1.

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

var

1_{y_{t}=v_{k}}={\begin{cases}1&{\text{if }}y_{t }=v_{k},\\0&{\text{annars}}\end{case}}

är en indikatorfunktion, och $b_{i}^{*}(v_{k})$ är det förväntade antalet gånger utdataobservationerna har varit lika med $v_ {k}$ i tillstånd $i$ över det förväntade totala antalet gånger i tillstånd $i$ .

Dessa steg upprepas nu iterativt tills en önskad nivå av konvergens.

Obs: Det är möjligt att överanpassa en viss datamängd. Det vill säga $P(Y\mid \theta _{\text{final}})>P(Y\mid \theta _{\text{ sant}})$ . Algoritmen inte heller ett globalt maximum.

Flera sekvenser

Algoritmen som hittills beskrivits antar en enda observerad sekvens $Y=y_{1},\ldots ,y_{N}$ . Men i många situationer finns det flera sekvenser observerade: $Y_{1},\ldots ,Y_{R}$ . I det här fallet måste informationen från alla observerade sekvenser användas i uppdateringen av parametrarna $A$ , $\pi$ och $b$ . Förutsatt att du har beräknat $\gamma _{ir}(t)$ och $\xi _{ijr}(t)$ för varje sekvens $y_{1,r},\ldots ,y_{N_{r},r}$ , parametrarna kan nu uppdateras:

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

var

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

är en indikatorfunktion

Exempel

Anta att vi har en kyckling som vi plockar ägg från vid middagstid varje dag. Huruvida kycklingen har lagt ägg för insamling eller inte beror på några okända faktorer som är dolda. Vi kan dock (för enkelhetens skull) anta att kycklingen alltid befinner sig i ett av två tillstånd som påverkar om hönan lägger ägg, och att detta tillstånd bara beror på tillståndet föregående dag. Nu vet vi inte tillståndet vid den initiala utgångspunkten, vi vet inte övergångssannolikheterna mellan de två tillstånden och vi vet inte sannolikheten att kycklingen lägger ett ägg givet ett visst tillstånd. Till att börja med gissar vi först övergångs- och emissionsmatriserna.

Övergång
	Tillstånd 1	Tillstånd 2
Tillstånd 1	0,5	0,5
Tillstånd 2	0,3	0,7

Utsläpp
	Inga ägg	Ägg
Tillstånd 1	0,3	0,7
Tillstånd 2	0,8	0,2

Första
Tillstånd 1	0,2
Tillstånd 2	0,8

Vi tar sedan en uppsättning observationer (E = ägg, N = inga ägg): N, N, N, N, N, E, E, N, N, N

Detta ger oss en uppsättning observerade övergångar mellan dagar: NN, NN, NN, NN, NE, EE, EN, NN, NN

Nästa steg är att uppskatta en ny övergångsmatris. Till exempel, sannolikheten för att sekvensen NN och tillståndet är $S_{1}$ sedan $S_{2}$ ges av följande, $P(S_{1})*P(N|S_{1})*P(S_{1}\högerpil S_{2})*P(N|S_{2}).$

Observerad sekvens	Högsta sannolikhet att observera den sekvensen om tillståndet är $S_{1}$ sedan $S_{2}$	Högsta sannolikhet att observera den sekvensen
NN	0,024 = 0,2 * 0,3 * 0,5 * 0,8	0,3584	$S_{2}$ , $S_{2}$
NN	0,024 = 0,2 * 0,3 * 0,5 * 0,8	0,3584	$S_{2}$ , $S_{2}$
NN	0,024 = 0,2 * 0,3 * 0,5 * 0,8	0,3584	$S_{2}$ , $S_{2}$
NN	0,024 = 0,2 * 0,3 * 0,5 * 0,8	0,3584	$S_{2}$ , $S_{2}$
NE	0,006 = 0,2 * 0,3 * 0,5 * 0,2	0,1344	$S_{2}$ , $S_{1}$
EE	0,014 = 0,2 * 0,7 * 0,5 * 0,2	0,0490	$S_{1}$ , $S_{1}$
SV	0,056 = 0,2 * 0,7 * 0,5 * 0,8	0,0896	$S_{1}$ , $S_{2}$
NN	0,024 = 0,2 * 0,3 * 0,5 * 0,8	0,3584	$S_{2}$ , $S_{2}$
NN	0,024 = 0,2 * 0,3 * 0,5 * 0,8	0,3584	$S_{2}$ , $S_{2}$
Total	0,22	2,4234

Den nya uppskattningen för övergången $S_{1}$ till $S_{2}$ är nu ${\frac {0,22}{2,4234}}=0,0908$ (kallas "Pseudo sannolikheter" i följande tabeller). Vi beräknar sedan $S_{2}$ till $S_{1}$ , $S_{2}$ till $S_{2}$ och $S_{1}$ till $S_{1}$ övergångssannolikheter och normalisera så att de adderas till 1. Detta ger oss den uppdaterade övergångsmatrisen:

Gammal övergångsmatris
	Tillstånd 1	Tillstånd 2
Tillstånd 1	0,5	0,5
Tillstånd 2	0,3	0,7

Ny övergångsmatris (pseudo sannolikheter)
	Tillstånd 1	Tillstånd 2
Tillstånd 1	0,0598	0,0908
Tillstånd 2	0,2179	0,9705

Ny övergångsmatris (efter normalisering)
	Tillstånd 1	Tillstånd 2
Tillstånd 1	0,3973	0,6027
Tillstånd 2	0,1833	0,8167

Därefter vill vi uppskatta en ny emissionsmatris,

Observerad sekvens	Högsta sannolikhet att observera den sekvensen om E antas komma från $S_{1}$		Högsta sannolikhet att observera den sekvensen
NE	0,1344	$S_{2}$ , $S_{1}$	0,1344	$S_{2}$ , $S_{1}$
EE	0,0490	$S_{1}$ , $S_{1}$	0,0490	$S_{1}$ , $S_{1}$
SV	0,0560	$S_{1}$ , $S_{2}$	0,0896	$S_{2}$ , $S_{2}$
Total	0,2394		0,2730

Den nya uppskattningen för E som kommer från $S_{1}$ emission är nu ${\frac {0,2394}{0,2730}}=0,8769$ .

Detta gör att vi kan beräkna emissionsmatrisen enligt beskrivningen ovan i algoritmen, genom att lägga ihop sannolikheterna för respektive observerade sekvenser. Vi upprepar sedan för om N kom från $S_{1}$ och för om N och E kom från $S_{2}$ och normaliserar.

Gammal emissionsmatris
	Inga ägg	Ägg
Tillstånd 1	0,3	0,7
Tillstånd 2	0,8	0,2

Ny utsläppsmatris (uppskattningar)
	Inga ägg	Ägg
Tillstånd 1	0,0404	0,8769
Tillstånd 2	1 0000	0,7385

Ny emissionsmatris (efter normalisering)
	Inga ägg	Ägg
Tillstånd 1	0,0441	0,9559
Tillstånd 2	0,5752	0,4248

För att uppskatta de initiala sannolikheterna antar vi att alla sekvenser börjar med det dolda tillståndet $S_{1}$ och beräknar den högsta sannolikheten och upprepar sedan för $S_{2}$ . Återigen normaliserar vi sedan för att ge en uppdaterad initial vektor.

Slutligen upprepar vi dessa steg tills de resulterande sannolikheterna konvergerar på ett tillfredsställande sätt.

Ansökningar

Taligenkänning

Dolda Markov-modeller användes först för taligenkänning av James K. Baker 1975. Kontinuerlig taligenkänning sker genom följande steg, modellerade av en HMM. Särdragsanalys utförs först på temporala och/eller spektrala egenskaper hos talsignalen. Detta producerar en observationsvektor. Funktionen jämförs sedan med alla sekvenser av taligenkänningsenheterna. Dessa enheter kan vara fonem , stavelser eller helordsenheter. Ett lexikonavkodningssystem används för att begränsa de sökvägar som undersöks, så endast ord i systemets lexikon (ordlexikon) undersöks. I likhet med lexikonavkodningen är systemvägen ytterligare begränsad av reglerna för grammatik och syntax. Slutligen tillämpas semantisk analys och systemet matar ut det igenkända yttrandet. En begränsning för många HMM-applikationer för taligenkänning är att det aktuella tillståndet endast beror på tillståndet vid det föregående tidssteget, vilket är orealistiskt för tal eftersom beroenden ofta är flera tidssteg i varaktighet. Baum–Welch-algoritmen har också omfattande tillämpningar för att lösa HMMs som används inom området för talsyntes.

Kryptanalys

Baum–Welch-algoritmen används ofta för att uppskatta parametrarna för HMMs vid dechiffrering av dold eller brusig information och används därför ofta i kryptoanalys . Inom datasäkerhet skulle en observatör vilja extrahera information från en dataström utan att känna till alla parametrar för överföringen. Detta kan innebära reverse engineering av en kanalkodare . HMMs och som en konsekvens av Baum-Welch-algoritmen har också använts för att identifiera talade fraser i krypterade VoIP-samtal. Dessutom är HMM-krypteringsanalys ett viktigt verktyg för automatiserade undersökningar av cache-timingsdata. Det möjliggör automatisk upptäckt av kritiska algoritmtillstånd, till exempel nyckelvärden.

Tillämpningar inom bioinformatik

Att hitta gener

Prokaryot

Programvaran GLIMMER (Gene Locator and Interpolated Markov ModelER) var ett tidigt gensökningsprogram som användes för identifiering av kodande regioner i prokaryot DNA. GLIMMER använder Interpolated Markov Models (IMM) för att identifiera de kodande regionerna och skilja dem från det icke-kodande DNA:t . Den senaste utgåvan (GLIMMER3) har visat sig ha ökad specificitet och noggrannhet jämfört med sina föregångare när det gäller att förutsäga translationsinitieringsställen, vilket visar en genomsnittlig 99% noggrannhet i lokalisering av 3'-platser jämfört med bekräftade gener i prokaryoter.

Eukaryot

GENSCAN - webbservern är en genlokaliserare som kan analysera eukaryota sekvenser upp till en miljon baspar (1 Mbp) långa. GENSCAN använder en allmän inhomogen, tre periodiska, femte ordningens Markov-modell av DNA-kodande regioner. Dessutom står denna modell för skillnader i gendensitet och struktur (som intronlängder) som förekommer i olika isokorer . Medan de flesta integrerade gensökningsprogramvaran (vid tidpunkten för GENSCANs release) antog att indatasekvenser innehöll exakt en gen, löser GENSCAN ett allmänt fall där partiella, fullständiga eller flera gener (eller till och med ingen gen alls) är närvarande. GENSCAN visades exakt förutsäga exonlokalisering med 90 % noggrannhet med 80 % specificitet jämfört med en kommenterad databas.

Detektering av kopia-nummervariationer

Kopieringsnummervariationer (CNV) är en riklig form av genomstrukturvariation hos människor. En diskret värderad bivariat HMM (dbHMM) användes för att tilldela kromosomala regioner till sju distinkta tillstånd: opåverkade regioner, deletioner, duplikationer och fyra övergångstillstånd. Att lösa denna modell med hjälp av Baum-Welch visade förmågan att förutsäga platsen för CNV-brytpunkten till cirka 300 bp från mikroarrayexperiment . Denna upplösningsstorlek möjliggör mer exakta korrelationer mellan olika CNV och över populationer än vad som tidigare varit möjligt, vilket gör det möjligt att studera CNV-populationsfrekvenser. Det visade också ett direkt arvsmönster för en viss CNV .

Genomföranden

Accord.NET i C#
ghmm C-bibliotek med Python- bindningar som stöder både diskreta och kontinuerliga emissioner.
Jajapy Python- bibliotek som implementerar Baum-Welch på olika typer av Markov-modeller ( HMM , MC , MDP , CTMC , GOHMM och MGOHMM).
HMMBase- paket för Julia .
HMMFit-funktion i RHmm -paketet för R .
hmmträna i MATLAB
rustbio i Rust

Se även

externa länkar

En omfattande genomgång av HMM-metoder och mjukvara inom bioinformatik – Profile Hidden Markov Models
Tidiga HMM-publikationer av Baum:
The Shannon Lecture av Welch, som talar om hur algoritmen kan implementeras effektivt:
- Hidden Markov Models and the Baum–Welch Algorithm , IEEE Information Theory Society Newsletter, dec. 2003.
Ett alternativ till Baum–Welch-algoritmen, Viterbi Path Counting-algoritmen:
- Davis, Richard IA; Lovell, Brian C.; "Jämföra och utvärdera HMM-ensembleträningsalgoritmer med hjälp av tåg- och test- och tillståndsnummerkriterier", Pattern Analysis and Applications, vol. 6, nr. 4, s. 327–336, 2003.
Ett interaktivt kalkylblad för att lära ut framåt-bakåtalgoritmen (kalkylblad och artikel med steg-för-steg-genomgång)
Formell härledning av Baum–Welch-algoritmen
Implementering av Baum–Welch-algoritmen