Dolda Markovmodeller

(1)

Institutionen för naturvetenskap och teknik

Självständigt arbete

Dolda Markovmodeller

(2)

Örebro universitet

Institutionen för naturvetenskap och teknik

Självständigt arbete för kandidatexamen i matematik, 15 hp

Självständigt arbete

Dolda Markovmodeller

Anton Wirén 21 juni 2018

Handledare: Mårten Gulliksson Examinator: Niklas Eriksen

(3)

Sammanfattning

Denna uppsats bygger på de tre klassiska problemen för en dold Markov-modell. Alla tre problemen kommer att matematiskt beskrivas på ett full-ständigt sätt. Bland annat kommer problemet med hur man tränar en mo-dell alltså lösas analytiskt med Lagrange multiplikator samt motivera varför Expectation-Maximization fungerar. En väsentlig del kommer också vara att introducera beräkningseffektiva algoritmer för att göra det möjligt att lösa dessa problem.

(4)

(5)

Innehåll

1 Introduktion 5

1.1 Grundläggande definitioner för sannolikheter . . . 5

1.2 Markovkedja . . . 6

1.3 Dold Markovmodell . . . 8

1.3.1 Problemformulering för dolda Markovmodeller . . . 10

2 Sannolikheten för en observationssekvens 11 2.1 Den direkta metoden . . . 11

2.2 Framåtalgoritmen . . . 13

2.2.1 Den skalade framåtalgoritmen . . . 16

2.3 Bakåtalgoritmen . . . 18

2.3.1 Den skalade bakåtalgoritmen . . . 21

3 Den mest sannolika tillståndssekvensen 22 3.1 Metod 1 - Lokal avkodning . . . 22

3.2 Metod 2 - Global avkodning (Viterbialgoritmen) . . . 24

4 Att träna en modell 28 4.1 Expectation-Maximization . . . 28

4.2 Härledning av Expectation-Maximization . . . 29

4.3 Härledning av Baum-Welch algoritm . . . 31

4.3.1 Optimum för πi . . . 35

4.3.2 Optimum för aij . . . 36

4.3.3 Optimum för bij . . . 37

4.3.4 Att tolka resultaten . . . 38

(6)

(7)

Kapitel 1

Introduktion

Dolda Markovkedjor introducerades av L.E Baum och Petrie (1966) och har sedan dess använts inom bland annat för taligenkänning, signalbehandling, bioinformatik och ekonomi. Med dolda Markovmodeller kan observerade data användas för att identifiera den underliggande dolda datan som antas komma från en Markovkedja. Huvuddelen av denna uppsats bygger på [3]. I kapitel 4 har även källorna [1], [2], [4] och [5] använts som komplement till [3].

1.1 Grundläggande definitioner för sannolikheter

Innan vi går in på ämnet dolda Markovmodeller ges här några elementära begrepp för sannolikheter.

Definition 1.1.1. Låt A₁, . . . , An vara händelser. Låt P (A1, . . . , An) = P n \ i=1 Ai !

och utläses sannolikheten för A1 och A2 . . . och An.

Skrivsättet P (A1, . . . , An) kommer alltid att användas när vi beskriver sannolikheter.

Definition 1.1.2 (Betingad sannolikhet). Låt A₁ och A₂ vara händelser. Sannolikheten för A1 givet A2 ges av

P (A1 | A2) =

P (A1, A2) P (A2)

.

Definition 1.1.3 (Radstokastisk matris). Låt C = {c_ij} vara en n × m matris. Vi säger att C är en radstokastisk matris om alla element är icke negativa och varje rad har summa 1. Alltså om

(8)

och _m X j=1 cij = 1, i = 1, . . . , n (1.2) är uppfyllt.

1.2 Markovkedja

Definition 1.2.1(Markovegenskapen). En stokastisk process {x_t_{: t ∈ Z}+} med tillståndsrum S har Markovegenskapen om

P (xt+1= j | xt= i, xt−1= k, . . . , x1 = m) = P (xt+1 = j | xt= i) (1.3) för alla tillstånd i, j, k, . . . , m ∈ S och alla t ≥ 1, och brukar då kallas för en Markovkedja.

Den stokastiska variabeln xtbeskriver alltså vilket tillstånd man befinner sig i kedjan vid tidpunkt t och är det enda som behövs för att kunna veta vilka olika tillstånd xt+1 kan anta med en viss sannolikhet. Man brukar säga att en Markovkedja inte har något minne eftersom det endast beror på var man befinner sig i nuet och inte på historien. För ett konkret exempel på en Markovkedja, antag att man singlar slant och vinner 1 kr om myntet visar krona eller förlorar 1 kr om myntet visar klave. Om vi vet att man efter t kast av myntet har 10 kr så kommer man efter nästa kast antingen att ha 9 kr eller 11 kr oberoende av hur de t − 1 första kasten sett ut. Detta är alltså en Markovkedja. I fall då det inte är självkart att vi kan beskriva en stokastisk process som en Markovkedja kan man såklart anta att (1.3) är uppfylld. Men för ett exempel på en stokastisk process som inte kan beskrivas av en Markovkedja, antag att vi har en låda med ett känt antal röda, gröna och blåa bollar där man i varje tidpunkt drar en boll slumpmässigt utan återläggning. Om man vid tidpunkt t vill beräkna sannolikheten att vid t+1 dra en viss typ av färg på bollen, beror detta på alla tidigare t dragningar. Detta är inte en Markovkedja. Generellt är allt som dras utan återläggning inte en Markovkedja, eftersom det skapas ett “minne” för varje gång man drar något.

Definition 1.2.2(Homogen Markovkedja). En Markovkedja med tillstånds-rum S sägs ha konstanta övergångssannolikheter om för alla i, j ∈ S är san-nolikheten att gå från tillstånd i till tillstånd j på en steglängd oberoende av när det inträffar, alltså när P (xt+1 = j | xt = i) är oberoende av t. Vi säger då att Markovkedjan är homogen.

Vi kommer alltid anta att Markovkedjan är homogen och använda no-tationen aij = P (xt+1 = j | xt = i) för övergångssannolikheterna. Med en Markovkedja brukar man också definiera en vektor πi = P (x1 = si) som

(9)

s1 s2

π1 = 0.5 π2= 0.5

0.4 0.6

0.3 0.7

Figur 1.1: Grafisk beskrivning av en Markovkedja för de två tillstånden s1 och s2.

används för att generera det första tillståndet vid tidpunkt t = 1. En full-ständig beskrivning av en Markovkedja kan göras med en 3-tupel (S,π,A) där

• S = {s1, . . . , sN}är en mängd av N tillstånd (states).

• π = {πi} är den initiala tillståndsfödelningsvektorn där varje element beskriver sannolikheten att vid tidpunkt t = 1 vara i ett visst tillstånd, πi = P (x1 = si), i ∈ S med πi∈ [0, 1] och PNi=1πi= 1.

• A = {aij}är en övergångssannolikhetsmatris där varje element beskri-ver sannolikheten att gå från tillstånd i till tillstånd j på en tidsenhet, aij = P (xt+1 = sj | xt = si), i, j ∈ S, t ∈ Z+ med aij ∈ [0, 1] och PN

j=1aij = 1för alla i ∈ S.

Notera att både π och A är radstokastisk.

Exempel 1.2.1. Modellen i figur 1.1 beskriver en stokastisk process för två tillstånd s1 och s2, alltså S = {s1, s2} så att N = 2. Sannolikheten att vid t = 1 befinna sig i tillstånd s1 är π1 = 0.5 och densamma för s2. Om man befinner sig i s1 vid tidpunkt t ≥ 1 så kommer man vid t + 1 att vara i s1 med sannolikheten 0.6 eller s2 med 0.4. Om man vid tidpunkt t ≥ 1 befinner sig i s2 kommer man vid t + 1 att vara i s1 med sannolikheten 0.3 eller s2 med sannolikheten 0.7. Vi beskriver grafen på matrisform som

π = 0.5 0.5 A = 0.6 0.4 0.3 0.7

Lemma 1.2.1. Låt X = (x1, . . . , xT) vara en tillståndssekvens från en Mar-kovkedja där x_t∈ S, 1 ≤ t ≤ T . Då är P (X | π, A) = πx1 T Y t=2 axT −1xT. (1.4)

Exempel 1.2.2. Antag vi har tillståndssekvensen X = (x₁ = s2, x2 = s1, x3 = s1) som vi förkortar till X = (s2, s1, s1), som genererats från modellen från figur 1.1. Sannolikheten för denna sekvens är 0.09 eftersom P (X | π, A) = π2· a21· a11= 0.5 · 0.3 · 0.6 = 0.09.

(10)

1.3 Dold Markovmodell

Med en dold Markovmodell, eller HMM (Hidden Markov Model) antas den stokastiska processen vara en Markovkedja men att tillstånden inte direkt kan observeras, är dolda. Istället emitteras en observerad symbol från ett s.k vokabulär V = {v1, . . . , vm} med en given sannolikhet vid varje tidpunkt. Sannolikheten för att en viss symbol ska emitteras beror på vilket tillstånd Markovkedjan befinner sig i och som vi betecknar med bij = P (yt= vj | xt= si) och är element från emissionsmatrisen B. Notationen för Markovkedjor fortsätter att gälla även för en HMM. Formellt kan man representera en dold Markovmodell som en 5-tupel (S,V,π,A,B) där

• S = {s1, . . . , sN}är en mängd av N tillstånd (states).

• V = {v1, . . . , vM}är en mängd av M symboler från ett vokabulär. • π = {πi} är den initiala tillståndsfödelningsvektorn där varje element

beskriver sannolikheten att starta i ett visst tillstånd, πi = P (x1 = si), i ∈ S med πi ∈ [0, 1]och PNi=1πi = 1.

• A = {aij}är en övergångssannolikhetsmatris där varje element beskri-ver sannolikheten att gå från tillstånd i till tillstånd j på en tidsenhet, aij = P (xt+1 = sj | xt = si), i, j ∈ S, t ∈ Z+ med aij ∈ [0, 1] och PN

j=1aij = 1för alla i ∈ S.

• B = {bij}är en emissionssannolikhetsmatris där varje element beskri-ver sannolikheten att emittera en viss symbol från vokabuläret givet ett visst tillstånd, bij = P (yt = vj | xt = si), i ∈ S, j ∈ V, t ∈ Z+ med bij ∈ [0, 1] och PMj=1bij = 1.

Som vi kan se är π, A och B radstokastisk. För att underlätta notationen kommer vi att använda λ = (π, A, B) när vi beskriver en modell på kompakt form. En tillståndssekvens skrivs som X = (x1, . . . , xT) där xt ∈ S för t = 1, . . . , T och beskriver vilket tillstånd man varit i vid tidpunkt t. Notationen för en observationssekvens är Y = (y1, . . . , yT) där yt∈ V för t = 1 . . . , T .

(11)

x1 x2 _{. . .} xT x0 y1 y2 yT A A A π B B B

Figur 1.2: Vi generar det första tillståndet x1 med π. Från det första till-ståndet x1 genereras den första observationen y1 med B. Andra tillståndet x2 genereras med A och väl i det andra tillståndet genereras den andra ob-servationen y2. Denna process upprepas tills längden på sekvenserna är T .

I figuren 1.2 ser vi hur en tillståndssekvens X och en observationssekvens Y genereras. Dock kommer tillståndssekvensen alltid att vara dold för oss och den enda bevis vi har är just observationssekvensen som alltid kommer att vara synlig. Som exempel, antag vi har två mynt, mynt 1 och mynt 2, som har utfallen krona eller klave. Antag sannolikheten att få klave är 0.2 och 0.5 för mynt 1 respektive mynt 2. Ett mynt väljs med återläggning vid varje tidpunkt och som singlas. Vi får inte veta vilket mynt som använts vid vilken tidpunkt, men däremot får vi veta att utfallen blev klave, krona, klave. Utfallen skulle då bilda en observationssekvens, och de mynt som använts vid varje tidpunkt en tillståndssekvens som för oss är dold.

Exempel 1.3.1. I figuren 1.3 har modellen från exempel 1.2.1 utvidgats till en HMM. Till exempel kan vi ser att sannolikheten att emittera symbolen v2 givet att vi befinner oss i tillstånd s1 är 0.3. Denna modell skrivs som λ = (π, A, B)där π = 0.5 0.5 A = 0.6 0.4 0.3 0.7 B = 0.5 0.3 0.2 0.2 0.4 0.4

(12)

s1 s2 π1 = 0.5 π2= 0.5 v1 v2 v3 v1 v2 v3 0.4 0.6 0.3 0.7 0.5 0.3 0.2 0.2 0.4 0.4

Figur 1.3: Grafisk beskrivning av en dold Markovmodell med två tillstånd s1 och s2 samt emissionssymbolerna v1, v2, v3.

Lemma 1.3.1. Låt X = (x1, . . . , xT) vara en tillståndssekvens och låt Y = (y1, . . . , yT) vara en observationssekvens där xt ∈ S och yt ∈ V, 1 ≤ t ≤ T . Då är P (Y | X, λ) = T Y t=1 bxtyt. (1.5)

Exempel 1.3.2. Från modellen i figur 1.3, antag vi vill beräkna sannolikhe-ten att observera symbolerna v1, v3, v1givet tillståndssekvensen s2, s1, s1. Vi använder (1.5) som då ger att P (Y | X, λ) = b21·b13·b11= 0.2·0.2·0.5 = 0.02.

1.3.1 Problemformulering för dolda Markovmodeller

Men en dold Markovmodell finns det tre centrala problem som man behöver lösa för att kunna tillämpa modellen.

Sannolikheten för en observationssekvens (Problem 1)

Beräkna sannolikheten för en observationssekvens givet en modell. Alltså P (Y | λ).

Den mest sannolika tillståndssekvens (Problem 2)

Givet en observationssekvens och en modell, finn den mest sannolika (dolda) tillståndssekvens.

Att träna en modell (Problem 3)

Givet en eller flera observationssekvens, hitta den modell som bäst beskriver observationerna. Matematiskt arg max

λ

(13)

Kapitel 2

Sannolikheten för en

observationssekvens

2.1 Den direkta metoden

Antag man nu vill ta reda på sannolikheten för en observerad sekvens gi-vet en modell. Denna beräkning är användbar då man exempelvis har flera modeller och vill ta reda på vilken modell som bäst beskriver observations-sekvensen. Här antar man att en observationssekvens Y och en modell λ är given, men att tillståndssekvensen X är dold. För att nu kunna beräkna P (Y | λ) behöver vi uttrycka den i termer av de olika parameterna så att man får en sluten formel. Först skriver vi om

P (Y | λ) =X X

P (X, Y | λ) (2.1)

där man summerar över alla tillståndsekvenser X. Alla tillståndssekvenser innebär alltså alla möjliga NT _{kombinationer av tillståndssekvenser där N} är antalet tillstånd i modellen och T är längden av sekvensen. Exempelvis om längden på tillståndssekvensen är T = 3 och antalet tillstånd är N = 2 innebär detta att första sekvensen är s1, s1, s1, den andra s1, s1, s2, den trejdje s1, s2, s1 osv. ända tills den sista då s2, s2, s2. I detta fall totalt 23 tillståndssekvenser.

Sats 2.1.1 (Kedjeregeln för sannolikheter). Låt A₁, . . . , An vara händelser. Den sammansatta sannolikheten för alla händelserna kan då skrivas som

P n \ i=1 Ai ! = n Y i=1 P  Ai i−1 \ j=1 Aj  . (2.2)

Bevis. Enligt definitionen för betingade sannolikheter är

(14)

och om vi upprepar denna process för varje slutterm kommer det att ge

(2.2).

I nästa steg behöver vi formulera om uttrycket P (X, Y | λ) = P (X | λ)P (Y | X, λ)

som vi kan göra med hjälp av kedjeregeln för sannolikheter. Om vi kombi-nerar resultaten från (1.4) och (1.5) ger det att

P (X | λ)P (Y | X, λ) = πx1bx1y1 T Y t=2

axt−1xtbxtyt. (2.3)

Nu har vi allt som behövs för att kunna uttrycka P (Y | λ) som en formel i termer av parameter och kan sammanfatta att

P (Y | λ) =X X P (X, Y | λ) =X X P (X | λ)P (Y | X, λ) = X x1,...,xT πx1bx1y1 T Y t=2 axt−1xtbxtyt (2.4)

som vi kallar för den direkta metoden och är lösningen till problem 1. För att förstå hur man kan använda denna formel är det lämpligt med ett exempel. Exempel 2.1.1. Låt oss använda modellen från exempel 1.3.1 och göra ett mer illustrativt exempel. Antag vi har två olika lådor, låda 1 och låda 2, som innehåller röda, gröna och blåa bollar.

π = 0.5 0.5 A = 0.6 0.4 0.3 0.7 B = 0.5 0.3 0.2 0.2 0.4 0.4 Lådorna väljs stokastiskt enligt π och A där låda 1 representerar de första tillståndet och låda 2 de andra, formellt S = {s1 = Låda 1, s2 = Låda 2}. Bollarna väljs sedan stokastiskt med återläggning enligt B, där röd repre-senterar den första observationsklassen, grön den andra och blå den tredje, dvs V = {v1 = röd, v2 = grön, v3 = blå}. Om man utifrån den givna mo-dellen vill ta reda på hur sannolikt det är att först generera en röd boll, sen en blå boll och sist en röd boll använder man (2.4). Observationssekvensen för detta problem formuleras som Y = (v1, v3, v1) och eftersom längden på sekvensen är 3 och antalet tillstånd 2, har vi 23 _{beräkningar att genomföra.} En av dessa beräkningar är exempelvis

P (X = (s2, s1, s1), Y = (v1, v3, v1) | λ) = π2· b21· a21· b13· a11· b11 = 0.5 · 0.2 · 0.3 · 0.2 · 0.6 · 0.5 = 0.0018

(15)

som innebär att sannolikheten för tillståndssekvensen Låda 2, Låda 1, Låda 1 då vi observerat röd, blå, röd givet modellen är 0.0018. Fortsätter man att beräkna på samma sätt som ovan för alla tillståndsekvenser och summe-rar resultaten får man att P (Y | λ) = 0.0334 , dvs att sannolikheten att observera först röd boll, sen blå boll, och sist röd boll givet modellen λ är 0.0334.

Tabell 2.1: Sannolikheten att observera sekvensen röd, blå, röd givet en viss tillståndssekvens och λ. Summeras alla värden i tabellen får man att P (Y | λ) = 0.0334. Tillståndssekvens P (X, Y | λ) s1, s1, s1 0.009 s1, s1, s2 0.0024 s1, s2, s1 0.006 s1, s2, s2 0.0056 s2, s1, s1 0.0018 s2, s1, s2 0.00048 s2, s2, s1 0.0042 s2, s2, s2 0.00392

Att använda den direkta metoden (2.4) för att beräkna P (Y | λ) är det enklaste sättet att attackera detta problem eftersom vi kan använda oss av en sluten formel. Dock kommer denna metod att vara mycket beräknings-intensiv för längre sekvenser. Metoden kräver totalt (2T − 1)NT multipli-kationer och NT _{− 1} _{additioner [3]. I praktiken innebär detta att (2.4) inte} bör användas för långa observationssevenser ty den exponentiella ökningen av beräkningar då T ökar. För att belysa detta, antag vi har en modell med N = 2 tillstånd och T = 100 observationer. Antalet multiplikationer är då 2.52 · 1032 och antalet additioner 1.27 · 1030 vilket såklart inte är hanterbart.

2.2 Framåtalgoritmen

För att minimera antalet beräkningar för att ta fram P (Y | λ) kan man använda sig av den s.k framåtalgoritmen som använder rekursion för att på det mest effektiva sätt att beräkna P (Y | λ). För att göra detta behöver vi definiera en framåtvariabel som

αt(i) = P (y1, . . . , yt, xt= si| λ). (2.5) Det innebär att αt(i) är sannolikheten att observera delsekvensen y1, . . . , yt där man befinner sig i tillstånd si vid tidpunkt t givet λ. Att ta fram denna sannolikheter är inte svårt, men för att göra det effektivt behövs en algoritm.

(16)

s1 αt(1) s2 αt(2) ... sN αt(N ) sj yt+1 a 1j a_2j aN j bjyt+1

Figur 2.1: Visualisering av rekursionssteget. Anta vi vet αt(1), . . . , αt(N ). Sedan beräknas PN

i=1αt(i)aij för fixt j och vi multiplicerar resultatet av beräkningen med bjyt+1 som kommer att ge αt+1(j) enligt definition.

Att beräkna αt(i) för t = 1 kan göras på följande sätt,

α1(i) = P (y1, x1= si | λ) = P (x1 = si)P (y1, x1 = si) = πibiy1 (2.6)

där vi använt kedjeregeln och det är underförstått att t.ex P (x1 = si| λ) = P (x1 = si). Härnäst vill vi beräkna at(i) för 2 ≤ t ≤ T , vilket kan göras rekursivt. Med kedjeregeln gäller

Att vi rekursivt kan beräkna αt(i) är anledningen till varför algoritmen är så effektiv. Man fortsätter denna rekursion ända tills αT(i) är beräknad

(17)

och som vi kan använda för att ta fram P (Y | λ) på följande sätt, P (Y | λ) = N X i=1 P (Y, xT = si | λ) = N X i=1 P (y1, . . . , yT, xT = si| λ) = N X i=1 αT(i). (2.8) Om vi sammanfattar dessa resultat får man den s.k framåtalgoritmen som lyder: 1. Initialisering α1(i) = πibiy1, i = 1 . . . , N 2. Rekursion αt+1(j) = bjyt+1 N X i=1 αt(i)aij, j = 1, . . . , N t = 1, . . . , T − 1 3. Slutfas P (Y | λ) = N X i=1 αT(i)

I första steget beräknas α1(i) genom att multiplicera övergångssannolik-heten för det initiala tillståndet med observationssannolikövergångssannolik-heten för den givna observationen och som görs för alla tillstånd. I andra steget utnyttjar man att αt+1(j) kan beräknas rekursivt genom att lista alla möjliga övergångar från tid t till t + 1. Detta görs genom att multiplicerar αt(i)med övergångs-sannolikheterna aij för alla tillstånd i, då vi håller j fixt, och summerar dessa sannolikheter, se Figur 2.1. Sedan multipliceras resultatet med observations-sannolikheten som ger αt+1(j). Processen upprepas för alla j och alla t upp till och med T − 1. Till sist summeras alla αT(i)för alla N tillstånd och som kommer att ge P (Y | λ).

Exempel 2.2.1. Vi kan givetvis lösa problemet från exempel 2.1.1 med framåtalgoritmen för att beräkna P (Y | λ) där vi använder samma modell och samma observationssekvens, alltså Y = (v1, v3, v1). Med framåtalgorit-men beräknas αt(i)för alla tidpunkter och alla tillstånd.

För t = 1: α1(1) = π1b11= 0.5 · 0.5 = 0.25 α1(2) = π2b21= 0.5 · 0.2 = 0.1 För t = 2: α2(1) = b13(α1(1)a11+ α1(2)a21) = 0.2(0.25 · 0.6 + 0.1 · 0.3) = 0.036 α2(2) = b23(α1(1)a12+ α1(2)a22) = 0.4(0.25 · 0.4 + 0.1 · 0.7) = 0.068

(18)

För t = 3:

α3(1) = b11(α2(1)a11+ α2(2)a21) = 0.5(0.036 · 0.6 + 0.068 · 0.3) = 0.021 α3(2) = b21(α2(1)a12+ α2(2)a22) = 0.2(0.036 · 0.4 + 0.068 · 0.7) = 0.0124. Genom att använda slutfassteget i framåtalgoritmen får man att P (Y | λ) = α3(1) + α3(2) = 0.021 + 0.0124 = 0.0334vilket stämmer överens med resul-tatet från exempel 2.1.1.

Framåtalgoritmen kräver totalt N(N +1)(T −1)+N multiplikationer och N (N − 1)(T − 1)antal additioner [3] vilket är en kraftig minskning jämfört med den direkta metoden (2.4), speciellt då T är stort. Antalet beräkningar ökar linjärt med avseende på T till skillnad från den direkta metoden som ökade exponentiellt. Om vi använder samma antal tillstånd och observationer som för den direkta metoden, N = 2 och T = 100, är antalet multiplikationer 596och antalet additioner 198 vilket nu är hanterbart.

2.2.1 Den skalade framåtalgoritmen

När man beräknar αt(i)kommer den innehålla multiplikationer av sannolik-heter som i de flesta fall är strikt mindre än 1. Detta leder till att αt(i) går exponentiellt mot noll då t ökar. Vid tillräckligt stort t kommer αt(i) vara mindre än maskinnoggranheten och vi får noggrannhetsförlust (underflow), dvs att talen är så pass små att beräkningarna inte blir exakta. Storleken på tnär αt(i)får noggrannhetsförlust varierar beroende på modell, så det finns ingen generell regel. Men för att undvika detta normaliserar vi varje αt(i) genom att dividera summan över j av αt(j) och som vi kallar ˆαt(i), alltså

ˆ αt(i) = αt(i) PN j=1αt(j) (2.9) och som gör att storleken på ˆαt(i) endast beror på den relativa storleken av αt(i)jämfört med summan för alla tillstånden vid samma tidpunkt. Det-ta innebär att storleken på ˆαt(i) inte minskar då t ökar, storleken är alltså oberoende av tiden. Men eftersom αt(i) kan vara mindre än maskinnoggran-heten vid stora t kan man inte direkt använda (2.9). Alltså även om ˆαt(i) är ett hanterbart tal behöver inte täljeren eller nämnaren vara detta. För att kunna använda (2.9) måste vi hela tiden öka storleken på αt(i)vid varje tidpunkt. Detta är precis vad den skalade framåtalgoritmen gör. Men när man ändrar skalan på alla αt(i) kommer man inte kunna beräkna P (Y | λ) på samma sätt som innan, utan man får istället använda sig av ett indirekt sätt med hjälp av skalärer. Vi ger algoritmen först och sen motivera varför den fungerar.

1. Initialisering

(19)

˜ α1(i) = α1(i), i = 1 . . . , N c1 = 1 PN i=1α1(i) ˆ α1(i) = c1α1(i) 2. Rekursion ˜ αt(j) = bjyt N X i=1 ˆ αt−1(i)aij, j = 1, . . . , N t = 2, . . . , T ct= 1 PN i=1α˜t(i) , t = 2, . . . , T ˆ αt(i) = ctα˜t(i), i = 1, . . . , N t = 2, . . . , T 3. Slutfas ln P (Y | λ) = − T X t=1 ln ct

Anledningen till varför man inte kan beräkna P (Y | λ) som i den vanliga framåtalgoritmen är att man vid varje tidpunkt skalar αt(i)som gör att den tappar sin definition. Man får alltså istället använda sig av alla skalärer ct som vi kan se i det sista steget. Det är såklart väsentligt att visa varför man kan göra på detta sätt. Först behöver vi visa att

ˆ

αt(i) = ctα˜t(i) =

αt(i) PN

j=1αt(j) vilket vi gör med induktion. Vi antar att

ˆ αt(i) = t Y k=1 ck ! αt(i). (2.10)

För fallet t = 1 är ˆα1(i) = c1α1(i) enligt det första steget i algoritmen. För t + 1blir ˆ αt+1(i) = ct+1α˜t+1(i) = ct+1bjyt+1 N X i=1 ˆ αt(i)aij = t+1 Y k=1 ck ! bjyt+1 N X i=1 αt(i)aij = t+1 Y k=1 ck ! αt+1(i)

(20)

som visar att induktionsantagendet stämmer för alla t. Med hjälp av (2.10) och definitionen för ˜αt(i) samt ˆαt(i) kan man få fram att

ˆ αt(i) = αt(i) PN j=1αt(j) (2.11)

vilket innebär att ˆαt(i)behåller sin definition även vid skalning för alla t. En konsekvens av (2.11) är att N X i=1 ˆ αT(i) = 1 (2.12) och från (2.10) får vi att N X i=1 ˆ αT(i) = T Y k=1 ck ! _N X i=1 αT(i) = T Y k=1 ck ! P (Y | λ). (2.13) Om man kombinerar resultaten (2.12) och (2.13) ger det att

P (Y | λ) = 1 T Y k=1 ck . (2.14)

Dock kan denna form vara problemetisk vid stora t eftersom att alla ct> 1 som kan leda till noggrannhetsförlust. För att undvika detta beräknas istället den logaritmerade versionen

ln P (Y | λ) = − T X k=1

ln ck (2.15)

vilket motiverar det sista steget i den skalade framåtalgoritmen.

2.3 Bakåtalgoritmen

På samma sätt som vi beräknar en framåtvariabel rekursivt kan man även göra det med en bakåtvariabel, fast bakåt i tiden, och som vi definierar

βt(i) = P (yt+1, . . . , yT | xt= si, λ) (2.16) som betyder sannolikheten att observera delsekvensen yt+1, . . . , yT givet att man i tidpunkt t varit i tillstånd sioch givet en modell. På liknande sätt som framåtalgoritmen vill vi ta fram alla βt(i) på ett effektivt sätt. För t = T väljer man att låta βT(i) = 1 godtyckligt. För t ≤ T − 1 kommer vi att

(21)

rekursivt beräkna bakåtvariablerna. Med hjälp av kedjeregeln kan man visa att βt(i) = P (yt+1, . . . , yT | xt= si, λ) = N X j=1 P (yt+1, . . . , yT, xt+1= sj | xt= si, λ) = N X j=1 P (yt+1, . . . , yT | xt= si, xt+1 = sj, λ)P (xt+1= sj | xt= si, λ) = N X j=1 P (yt+2, . . . , yT | yt+1, xt= si, xt+1= sj, λ)P (yt+1| xt= si, xt+1= sj, λ)aij = N X j=1 P (yt+2, . . . , yT | xt= si, xt+1 = sj, λ)P (yt+1| xt+1= sj, λ)aij = N X j=1 βt+1(j)bjyt+1aij (2.17)

enligt definition. Notera att rekursionen går bakåt i tiden. Denna rekursion fortsätter tills man når t = 1 så att β1(i)för alla tillstånd är beräknad. Dessa kan användas för att ta fram P (Y | λ). Vi använder resultaten från (2.17) för att visa att

För att koppla samman alla dessa resultat kan vi skriva den fullständiga bakåtalgoritmen som 1. Initialisering βT(i) = 1, i = 1, . . . , N 2. Rekursion βt(i) = N X j=1 βt+1(i)bjyt+1aij, i = 1, . . . , N t = T − 1, . . . , 1

(22)

s1 βt+1(1) s2 βt+1(2) ... sN βt+1(N ) si yt+1 ai1 ai2 a_iN

Figur 2.2: Visualisering av rekursionssteget i bakåtalgoritmen. 3. Slutfas P (Y | λ) = N X j=1 β1(j)bjy1πj

I första steget sätter man βT(i) = 1för alla i. Sen beräknar man rekursivt bakåt genom att lista alla övergångar från t + 1 till t för alla kombinatio-ner av övergångar. Man gör detta genom att multiplicera βt+1(i) med över-gångssannolikheten aij samt observationssannolikheten bjyt+1 för alla

till-stånd och som görs för alla tidpunkter t = T −1, . . . , 1. I slutfasen summerar vi β1(j)πjbjy1 för alla tillstånd som ger P (Y | λ).

Exempel 2.3.1. Antag vi vill lösa problemet från exempel 2.1.1 med bak-åtalgoritmen. Först beräknar vi alla bakåtvariabler βt(i) för alla tidpunkter och alla tillstånd.

För t = 3: β3(1) = 1 β3(2) = 1 För t = 2: β2(1) = β3(1)a11b11+ β3(2)a12b21= 1 · 0.6 · 0.5 + 1 · 0.4 · 0.2 = 0.38 β2(2) = β3(1)a21b11+ β3(2)a22b21= 1 · 0.3 · 0.5 + 1 · 0.7 · 0.2 = 0.29 För t = 1: β1(1) = β2(1)a11b13+ β2(2)a12b23= 0.38 · 0.6 · 0.2 + 0.29 · 0.4 · 0.4 = 0.092 β1(2) = β2(1)a21b13+ β2(2)a22b23= 0.38 · 0.3 · 0.2 + 0.29 · 0.7 · 0.4 = 0.104.

(23)

Till sist beräknar vi sannolikheten för observationssekvensen givet modellen som

P (Y | λ) = β1(1)b11π1+ β1(2)b21π2

= 0.029 · 0.5 · 0.5 + 0.104 · 0.2 · 0.5 = 0.0334.

Antalet beräkningar för bakåtalgoritmen är densamma som för framå-talgoritmen, så det är ingen fördel att använda den ena eller den andra. Det är även möjligt att kombinera framåt och bakåtvariabler för att ta fram sannolikheten för en observationssekvens. Eftersom

P (Y, xt= si| λ) = P (y1, . . . , yT, xt= si | λ) = P (y1, . . . , yt, xt= si| λ)P (yt+1, . . . , yT | y1, . . . , yt, xt= si, λ) = P (y1, . . . , yt, xt= si| λ)P (yt+1, . . . , yT | xt= si, λ) = αt(i)βt(i) (2.18) och P (Y | λ) = N X i=1 P (Y, xt= si | λ) = N X i=1 αt(i)βt(i) (2.19) innebär det att vi har flera sätt att ta fram sannolikheten för en observa-tionssekvens. Denna relation kommer att kunna utnyttjas vid senare tillfälle. 2.3.1 Den skalade bakåtalgoritmen

Precis som med framåtalgoritmen kommer bakåtalgoritmen också få nog-grannhetsförlust vid stora t för βt(i)variablerna. För att undvika detta kan vi använda skalärerna ct från den skalade framåtalgoritmen för att skalera βt(i), alltså ˆβt(i) = ctβt(i). Detta innebär att den skalade bakåtalgoritmen blir 1. Initialisering βT(i) = 1, i, . . . , N ˆ βT(i) = cTβT(i), i = 1, . . . , N 2. Rekursion ˜ βt(i) = N X j=1 ˆ βt+1(j)aijbjyt+1, i = 1, . . . , N ˆ βt(i) = ctβ˜t(i), i = 1, . . . , N.

(24)

Kapitel 3

Den mest sannolika

tillståndssekvensen

Antag man är ute efter att hitta den dolda tillståndssekvens som är den mest troliga givet en observationssekvens och en modell. Vi kan se detta som ett avkodningsproblem där man är ute efter att avkoda hur man förflyttat sig mellan de olika dolda tillstånden. För detta problem finns det ingen entydig lösning utan man kan göra det på olika sätt.

3.1 Metod 1 - Lokal avkodning

Ett sätt är hitta den mest sannolika tillståndssekvens är att vid varje tid-punkt hitta de tillstånd som är det mest troliga. Man vill alltså maximera de förväntade korrekta tillstånd genom att beräkna sannolikheten att vara i ett tillstånd vid en viss tidpunkt och välja den som är det mest sannolika. Gör man det för alla tidpunkter kommer vi att få den tillståndssekvens som maximerar antalet förväntade korrekta tillstånd. Att beräkna sannolikheten att vara i tillstånd si givet en observationssekvens samt en modell kan gö-ras på ett effektivt sätt med hjälp av framåt och bakåtvariabler genom att använda definitionen av betingad sannolikhet samt (2.18) och (2.19) för att uttrycka att P (xt= si | Y, λ) = P (xt= si, Y | λ) P (Y | λ) = αt(i)βt(i) PN

i=1αt(i)βt(i)

. (3.1) Vi kan även visa att (se [3] s.46)

αt(i)βt(i) PN

i=1αt(i)βt(i)

= αˆt(i) ˆβt(i) PN

(25)

för att beräkna P (xt = si | Y, λ) på ett effektivt sätt vid stora t. Vi finner det mest sannolika tillstånd för varje tidpunkt som

x∗_t = arg max 1≤i≤N

P (xt= si | Y, λ) t = 1, . . . , T (3.2) och betyder att vi väljer det tillstånd vid varje tidpunkt som är mest troligt och blir den optimala tillståndsekvens som maximerar antalet förväntade korrekta tillstånd. Vi kan se det som att denna metod lokalt avkodar till-stånden.

Exempel 3.1.1. Antag vi har modellen från exempel 2.1.1 med samma observationssekvens och vill hitta den tillståndssekvens som maximerar de förväntade korrekta tillstånden. Från exempel 2.2.1 fick vi att

α1(1) = 0.25, α2(1) = 0.036, α3(1) = 0.021 α1(2) = 0.1, α2(2) = 0.068, α3(2) = 0.0124 och från exemplet 2.3.1 beräknades att

β1(1) = 0.092, β2(1) = 0.38, β3(1) = 1 β1(2) = 0.104, β2(2) = 0.29, β3(2) = 1.

Vi kan nu använda oss av (3.1) för att beräkna sannolikheten att vara i ett viss tillstånd vid varje tidpunkt,

P (x1 = s1 | Y, λ) = 0.25 · 0.092 0.25 · 0.092 + 0.1 · 0.104 ≈ 0.69 P (x2 = s1 | Y, λ) = 0.036 · 0.38 0.0334 ≈ 0.41 P (x3 = s1 | Y, λ) = 0.021 · 1 0.0334 ≈ 0.63

och kan göra liknande beräkningar för tillstånd s2 för alla tidpunkter.

Tidpunkt 1 2 3

P(Låda 1) 0.69 0.41 0.63 P(Låda 2) 0.31 0.59 0.37

Utifrån denna tabell kan vi få fram enligt (3.2) att den tillståndssekvens som maximerar de förväntade antal korrekta tillstånd är Låda 1, Låda 2, Låda 1. Man kan givetvis lösa detta problem utan framåt och bakåtvariabler, men då denna metod inte är beräkningseffektiv bör den undvikas, speciellt vid långa sekvenser.

(26)

Ett problem som kan uppstå med denna metod är att den kan ge en icke giltig tillståndssekvens eftersom den inte tar hänsyn till övergångssan-nolikheterna. Antag tillstånd si är den mest sannolika vid tidpunkt t och att tillstånd sj är mest sannolik vid t + 1, men att aij = 0. Denna övergång är alltså inte möjlig vilket innebär att den optimala tillståndsekvensen inte är giltig. Vi behöver alltså en ny metod som tar hänsyn till övergångssannolik-heterna så att man alltid får en giltig tillståndssekvens. En naturlig lösning till detta problem är att ändra sättet vi hittar den optimala tillståndssekven-sen genom att parvis hitta de tillstånd xt, xt+1som är mest sannolika för alla tidpunkter. Även fast den metoden kan fungera brukar man istället med en HMM hitta den hela sekvens x1, . . . , xT som är den mest sannolika givet en modell som vi kan göra med Viterbialgoritmen.

3.2 Metod 2 - Global avkodning

(Viterbialgorit-men)

Med Viterbialgoritmen avkodas hela sekvensen på en gång till skillnad från (3.2) som individuellt avkodar varje tillstånd vid varje tidpunkt. Man kan säga att Viterbialgoritmen globalt avkodar den dolda tillståndssekvensen. Som vi såg från exempel 2.1.1 i tabellen 2.1 är tillståndssekvensen s1, s1, s1. Det är just denna vi vill finna med Viterbialgoritmen. Vi definierar δt(j)som den mest sannolika tillståndssekvens med längd t som slutar i tillstånd sj och emitterat de t första observationerna, alltså

δt(j) = max x1,...,xt−1

P (y1, . . . , yt, x1, . . . , xt= sj | λ). (3.3) Men för att beräkna δt(i) på ett effektivt sätt behöver vi även här använda oss av rekursion. Med kedjeregeln kan vi visa att

δt(j) = max x1,...,xt−1 P (y1, . . . , yt, x1, . . . , xt= sj | λ) = max x1,...,xt−1 P (yt, xt= sj | y1, . . . , yt−1, x1, . . . , xt−1, λ) P (y1, . . . , yt−1, x1, . . . , xt−1| λ) = max x1,...,xt−1=si P (yt, xt= sj | xt−1= si) P (y1, . . . , yt−1, x1, . . . , xt−1= si | λ) = max x1,...,xt−1=si P (y1, xt= sj, xt−1= si) P (xt−1=sj) P (y1, . . . , yt−1, x1, . . . , xt−1= si | λ) = max x1,...,xt−1=si P (yt| xt= sj)P (xt= sj | xt−1 = si)

(27)

P (y1, . . . , yt−1, x1, . . . , xt−1= si | λ) = max x1,...,xt−1=si bjytaijP (y1, . . . , yt−1, x1, . . . , xt−1= si | λ) = bjyt max xt−1=si aij max x1,...,xt−2 P (y1, . . . , yt−1, x1, . . . , xt−1= si| λ) = bjyt max 1≤i≤N δt(i)aij (3.4)

Kopplingen mellan Viterbialgoritmen och framåtalgoritmen är att istället för att summera över tidigare tillstånd, som görs i framåtalgoritmen, används maximering i rekursionsteget samt slutfasen. Men eftersom man vill ta fram den mest sannolika tillståndssekvensen behöver vi också hålla koll på vilka argumenet som maximerar δt(i) för alla tillstånd vid alla tidpunkter. Detta görs genom att spara argumentet i en matris ψ och används i slutet av algo-ritmen för att ta fram den optimala tillståndssekvensen, ett s.k backtracking steg. 1. Initialisering δ1(i) = πibiy1, i = 1, . . . , N ψ1(i) = 0, i = 1, . . . , N 2. Rekursion δt(j) = bjyt max 1≤i≤Nδt−1(i)aij, j = 1, . . . , N, t = 2, . . . , T ψt(j) = arg max 1≤i≤N δt−1(i)aij, j = 1, . . . , N, t = 2, . . . , T 3. Slutfas P (Y, X∗ | λ) = max 1≤i≤NδT(i) x∗_T = arg max 1≤i≤NδT(i) 4. Backtracking x∗_t = ψt+1(x∗t+1), t = T − 1 . . . 1

För att ge en mer intuitiv förklaring vad denna algoritm faktiskt gör använder vi ett exempel.

Exempel 3.2.1. Antag vi vill finna den mest sannolika sekvens då vi har modellen från 2.1.1 med samma observationssekvens, dvs Y = (y1 = v1, y2 = v3, y3 = v1), och som vi avkodar med Viterbi algoritmen. Vi har alltså N = 2 och T = 3.

(28)

För t = 1: δ1(1) = π1b11= 0.5 · 0.5 = 0.25 ψ1(1) = 0 δ1(2) = π2b21= 0.5 · 0.2 = 0.1 ψ1(2) = 0 För t = 2:

δ2(1) = b13max{δ1(1)a11, δ1(2)a21} = 0.2 · max{0.15, 0.03} = 0.03 ψ2(1) = 1, eftersom 0.15>0.03.

δ2(2) = b23max{δ1(1)a12, δ1(2)a22} = 0.4 · max{0.1, 0.07} = 0.04 ψ2(2) = 1

För t = 3:

δ3(1) = b11max{δ2(1)a11, δ2(2)a21} = 0.5 · max{0.018, 0.012} = 0.009 ψ3(1) = 1

δ3(2) = b21max{δ2(1)a12, δ2(2)a22} = 0.2 · max{0.012, 0.021} = 0.0056 ψ3(2) = 2

Eftersom max{δ3(1), δ3(2)} = δ3(1) = 0.009 vet vi att den mest troliga tillståndsekvens har sannolikheten 0.009 och den optimala sekvensen får vi ut genom backtracking-steget,

x∗₃ = arg max{δ3(1), δ3(2)} = 1 x∗₂ = ψ3(x∗3) = 1 x∗₁ = ψ2(x∗2) = 1 Vi har alltså att X∗ _{= (s}

1, s1, s1) vilket stämmer överens med vad vi såg i tabellen 2.1 från exemplet 2.1.1.

(29)

x0 s1 δ1(1) = 0.25 s2 δ1(2) = 0.1 t = 1 s1 δ2(1) = 0.03 s2 δ2(2) = 0.04 t = 2 s1 δ3(1) = 0.009 s2 δ3(2) = 0.0056 t = 3

Figur 3.1: Visuallisering av Viterbialgoritmen där vi listar alla övergångar vid varje tidpunkt. Bågarna visar den mest sannolika sekvens som slutar i ett visst tillstånd. I fallet då det finns en unik optimal tillståndssekvens kommer endast en båge att riktas mot ett tillstånd vid varje tidpunkt. Genom att ta den mest sannolika tillståndsekvens som slutar i ett visst tillstånd vid t = 3, alltså max{δ3(1), δ3(2)} = δ3(1), och följa baklänges (backtracking) från tillstånd s1vid t = 3 på de heldragna linjerna finner vi att den optimala tillståndssekvensen är s1, s1, s1. De gråa linjerna visar backtracking-steget.

En fördel med Viterbialgoritmen är just att den tar hänsyn till hela till-ståndssekvensen och inte enbart vid varje tidpunkt för sig vilket medför att den kommer att vara mer stabil i avkodningen. Alltså om någon osannolik händelse skulle inträffa vid en tidpunkt givet ett tillstånd kommer Viterbi-algoritmen att göra en mer sammanvägd bedömning för att avkoda detta tillstånd jämfört med (3.2) som mer frekvent byter tillstånd i avkodningen.

(30)

Kapitel 4

Att träna en modell

I detta problem vill man finna den modell λ som maximerar P (Y | λ). Ett annat sätt att se på det är att vi vill hitta en modell som bäst beskriver den observerade datan. Här tillåter vi alltså att parameterna får variera för att hitta denna modell, men med krav att π, A och B fortfarande ska vara radstokastisk. På grund av att tillståndssekvensen X är dold finns det ingen analytisk direkt metod för att hitta den optimala modellen [3]. Man får istället använda sig av en iterativ metod som kallas för Expectation-Maximization (EM). I fallet med en HMM brukar man kalla denna algoritm för Baum-Welch som är ett specialfall av EM.

4.1 Expectation-Maximization

EM algoritmen är en iterativ metod för att beräkna maximum-likelihood skattningen för en model vars datamaterial inte är fullständigt eller har dol-da variabler. Varje iteration i EM algoritmen består av två steg: expectation step (E-steg) och maximization step (M-steg). I E-steget beräknas vänte-värdet givet observerad data och det nuvarande skattning av parameterna i modellen. I M-steget maximeras likelihood-funktionen genom att använda skattningen som beräknats i E-steget. Dessa två steg upprepas ända tills kon-vergens är uppnådd. EM algoritmen garanterar att öka likelihood-funktionen för varje iteration så att den konvergerar mot ett lokalt maximum eller en sadelpunkt i specialfall. Generellt sett är det svårt att hitta ett globalt maxi-mum och det enda man kan göra är att ändra de initiala värdena på modellen i förhoppning att algoritmen ska konvergera till en bättre modell. Det är dock lämpligt att välja de initiala värdena i modellen på ett sätt som man tror är närmast den optimala för att öka chanserna att hitta just denna modell. Men om man inte har någon kvalificerad gissning går det bra att slumpa alla värdena med enda krav att den ska vara radstokastisk.

Låt Y vara en mängd observerad data och X dold data. Låt λ vara en modell med okända parametrar och λn vara skattningen av λ i iteration n.

(31)

EM algoritmen lyder då 1. E-steg: Q(λ, λn) =

X X

P (X | Y, λn) ln P (X, Y | λ) 2. M-steg: λn+1= arg max

λ

Q(λ, λn)

och är en iterativ process som används i syfte att hitta den modell λ∗ _som är ett lokalt maximum för P (Y | λ). Algoritmen bygger på att introducera en hjälpfunktion Q(λ, λn) som uttnyttjar att för alla λ som ökar Q(λ, λn) medför att P (Y | λ) också ökar och som upprepas tills konvergens.

Det finns nu två frågor som ska besvaras. Den ena är varför EM algo-ritmen fungerar och den andra hur den bästa skattningen tas fram i varje iteration i fallet med en HMM för att kunna träna en modell.

4.2 Härledning av Expectation-Maximization

Den modell λ som maximerar P (Y | λ) kallas för maximum-likelihood skatt-ningen av λ. Men att finna denna är svårt, så vi nöjer oss med en modell λ som lokalt maximerar P (Y | λ). För att hitta denna modell behöver vi introducera en logaritmerad likelihood-funktion definierad som

L(λ) = ln P (Y | λ). (4.1)

Denna likelihood-funktion är alltså en funktion av parameterna λ, och ef-tersom ln(x) är en strikt växande funktion kommer samma värde som max-imerar L(λ) också att maximera P (Y | λ). Antag att den aktuella modellen är λnoch man nu vill hitta en bättre modell som gör att L(λ) ökar. Vi skul-le alltså kunna uttrycka det som L(λ) > L(λn), eller ekvivalent att vi ska maximera differensen av

L(λ) − L(λn) = ln X

X

P (X | λ)P (Y | X, λ) − ln P (Y | λn). (4.3) För att gå vidare härifrån behöver vi använda oss av några satser.

Sats 4.2.1 (Jensens olikhet [1]). Låt f vara en reell kontinuerlig konvex funktion definierad på intervallet I. Om x1, . . . , xn ∈ I och φ1, . . . , φn ≥ 0 samtPn i=1φi = 1, då gäller f n X i=1 φixi ! ≤ n X i=1 φif (xi) och med omvänd olikhet om f är en konkav funktion.

(32)

Sats 4.2.2. [1] Om f(x) har en andraderivata på [a, b] och f00(x) ≤ 0 på [a, b], då är f (x) konkav på [a, b].

Låt f(x) = ln(x) vilket ger att f00_{(x) = −}1

x2 < 0 för x ∈ (0, ∞) och vi

har att ln(x) är en strikt konkav funktion på (0, ∞). Använder man Jensens olikhet blir alltså

ln n X i=1 φixi ≥ n X i=1 φiln(xi) för φ ≥ 0 och Pn

L(λ) ≥ L(λn) + ∆(λ | λn) , `(λ | λn). (4.4) I denna relation har man alltså att funktionen `(λ | λn)är uppåt begränsad

(33)

som innebär att om λ = λn kommer funktionerna `(λ | λn) och L(λ) skära varandra. Använder vi de båda resultaten från (4.4) och (4.5) kan man konstatera att för alla λ som ökar `(λ | λn) kommer också att öka L(λ), alltså

Q(¯λ, λn) > Q(λ, λn) −→ L(¯λ) > L(λ).

Syftet är att hitta den modell som maximerar L(λ). Men då L(λ) inte kan optimeras på ett ett direkt sätt [3] får man använda sig av en indirekt metod. Vad man gör är att istället hitta den modell som maximerar `(λ | λn)och som blir den nya modellen λn+1och kommer att ge den största möjliga ökningen av L(λ). Det som nu återstår är att ta arg max `(λ, λn) med avseende på λ och förenkla bort alla konstanter. Vi hittar då maximum av `(λ, λn) som

λn+1 = arg max λ `(λ | λn) = arg max λ ( L(λn) + X X P (X | Y, λn) ln P (X | λ)P (Y | X, λ) P (X | Y, λn)P (Y | λn) ) = arg max λ ( X X P (X | Y, λn) ln (P (X | λ)P (Y | X, λ)) ) = arg max λ ( X X P (X | Y, λn) ln P (X, Y, λ)P (X | λ) P (X, λ)P (λ) ) = arg max λ ( X X P (X | Y, λn) ln P (X, Y | λ) ) = arg max λ {Q(λ, λ_n)} (4.6)

vilket vi känner igen från M-step i algoritmbeskrivningen av EM.

4.3 Härledning av Baum-Welch algoritm

Vi ska nu använda oss av EM algoritmen för att finna en iterativ metod för hur man med en HMM kan träna en modell så att den bäst beskriver

(34)

observerade data. Detta specialfall av EM algoritmen brukar man kalla för Baum-Welch algoritm. Men att enbart använda sig av en observationssekvens kommer leda till att modellen övertränas till just detta urval, alltså att den skattade modellen troligtvis inte kommer att vara nära den sanna. Det man istället bör göra är att använda flera observationssekvenser samtidigt för att öka träffsäkerheten av den skattade modell. Vi definierar

Y = {Y(1), . . . , Y(K)} (4.7) där

Y(k)= (y(k)₁ , . . . , y_T(k)

k), k = 1, . . . , K (4.8)

är individuella observationssekvenser och som antas vara genererade från samma modell. Notera att varje observationssekvens kan vara av olika längd. Dock skulle dessa kunna vara korrelerade med varandra så för att göra det allmänt använder vi kedjeregeln för att skriva [2]:

           P (Y | λ) = P (Y(1) | λ)P (Y(2) _{| Y}(1)_{, λ) . . . P (Y}(K) _{| Y}(K−1)_{, . . . , Y}(1)_{, λ)} P (Y | λ) = P (Y(2) | λ)P (Y(3) _{| Y}(2)_{, λ) . . . P (Y}(1)_{| Y}(K)_{, . . . , Y}(2)_{, λ)} ... P (Y | λ) = P (Y(K) _{| λ)P (Y}(1) _{| Y}(K)_{, λ) . . . P (Y}(K−1)_{| Y}(1)_{, . . . , Y}(K−2)_{, Y}(K)_{, λ).} Om vi använder uttrycken ovan kan vi formulera de multipla

observations-sekvenserna givet en modell som

P (Y | λ) = K X k=1 wkP (Y(k)| λ) (4.9) där            w1= _K1P (Y(2) | Y(1), λ) . . . P (Y(K) | Y(K−1), . . . , Y(1), λ) w2= _K1P (Y(3) | Y(2), λ) . . . P (Y(1)| Y(K), . . . , Y(1), λ) ... wK = _K1P (Y(1) | Y(K), λ) . . . P (Y(K−1)| Y(K), Y(K−2). . . , Y(1), λ) (4.10) är vikter som uppfyller wk ∈ [0, 1] och PK_k=1wk = 1. Om vi nu använder hjälpfunktionen från EM algoritmen kan vi nu uttrycka att

Q(λ, λn) = K X k=1 wkQk(λ, λn) = K X k=1 wk X X P (X | Y(k), λn) ln(X, Y(k)| λ) (4.11)

(35)

vilket är den funktion som vi vill maximera med avseende på λ. För att göra detta behöver vi uttrycka (4.11) i termer av parameterna i en HMM och som kan maximeras med hjälp av Lagrangemultiplikator. Men innan vi tar oss an detta används att P (X | Y(k)_{, λ}

vilket vi kan göra eftersom P (X | λn) är oberoende av λ och inte kommer påverka den optimala lösningen, med andra ord P (X | λn) är en konstant. Som vi sett tidigare från (2.4) är

P (X, Y(k) | λ) = πx1b_x 1y1(k) Tk Y t=2 axt−1xtb_x ty(k)t (4.12)

och om vi logaritmerar och använder räknerregler fås

ln P (X, Y(k)| λ) = ln π_x₁+ Tk X t=2 ln ax1x2+ Tk X t=1 ln b x1y₁(k). (4.13)

Använder vi detta resultat blir ˆ Q(λ, λn) = K X k=1 wk X X ln πx1P (X, Y (k)_{| λ} n) + K X k=1 wk X X Tk X t=2 ln ax1x2P (X, Y (k)_{| λ} n) + K X k=1 wk X X Tk X t=1 ln b_x 1y(k)1 P (X, Y(k)| λn) (4.14) Nu när vi har parameterna uppdelade i tre delar kan vi ändra på hur vi summerar över alla tillståndssekvenser. Exempelvis kan vi skriva den första delen i (4.14) som K X k=1 wk X X ln πx1P (X, Y (k)_{| λ} n) = K X k=1 wk N X i=1 ln πiP (X, Y(k)| λn). (4.15) Anledningen till detta är att X betyder alla sätt man kan vara vid tidpunkt t = 1 eftersom parametern är πi. Detta kan istället beskrivas med summan

(36)

över alla tillstånd. Samma resonemang kan användas för den andra delen i (4.14) där aij är parameter. I detta fall betyder X alla möjliga övergångar från tillstånd i till tillstånd j. Detta kan beskrivas med en dubbelsumma över i och j för alla tillstånd. För den sista delen där bij är parameter kan man istället använda en summa över i för alla tillstånd. Matematiskt kan vi alltså skriva ˆ Q(λ, λn) = K X k=1 wk N X i=1 ln πiP (X, Y(k)| λn) + K X k=1 wk Tk X t=2 N X i=1 N X j=1 ln aijP (X, Y(k)| λn) + K X k=1 wk Tk X t=1 N X i=1 ln b iy(k)₁ P (X, Y (k)_{| λ} n).

Nu är ˆQ(λ, λn) formulerad i de olika HMM parameterna och vi kan an-vända Lagrangemultiplikator i syfte att hitta optimum.

Sats 4.3.1 (Lagrangemultiplikator [5]). Låt U ⊂ Rn och låt f, g₁, . . . , gm : U → R vara kontinuerlig. Låt x∗ vara en punkt som lokalt maximerar f med bivillkor gi = 0, i = 1, . . . , m. Antag f, g1, . . . , gm är kontinuerligt deri-verbar i punkten x∗ och att alla ∇g_i är linjärt oberoende. Då existerar det θ1, . . . , θm∈ R sådan att x∗ är en stationär punkt för Lagrangefunktionen

L(x1, . . . , xn, θ1, . . . , θM) = f (x1, . . . , xn) + m X

i=1

θigi(x1, . . . , xn) (4.16) Denna sats säger är att om vi löser ∇L(x1, . . . , xn, θ1, . . . , θm) = 0 kom-mer detta ge oss extrempunkter till problemet f(x) under bivillkoren gi(x) = 0, i = 1, . . . , m. Vi kan nu använda detta resultat för att hitta en optimal lösning till det ursprungliga problemet ˆQ(λ, λn) där vi har enligt definition att bivillkoren är 1 − N X i=1 πi= 0 1 − N X j=1 aij = 0, i = 1, . . . , N 1 − M X j=1 bij = 0, i = 1, . . . , N.

(37)

till detta problem som blir L(λ, λn) = ˆQ(λ, λn) + θπ 1 − N X i=1 πi ! + N X i=1 θai  1 − N X j=1 aij  + N X j=1 θbj 1 − M X i=1 bij ! (4.17) där θπ, θa1, . . . , θaN, θb1, . . . , θbN ∈ R är Lagrangemultiplikatorer. Denna

funk-tion visar sig vara en konvex funkfunk-tion eftersom man får att Hessianmatrisen har att alla element på diagonalen är negativa och resterande element noll [2], vilket innebär att den är positivt definit och medför att L(λ, λn) är en konvex funktion. Vi vet alltså att det optimum vi hittar kommer att vara ett globalt maximum. För att hitta denna punkt tar vi gradienten av para-metrar i Lagrangefunktionen för att finna den optimala lösningen. Vi får då följande ekvationssystem                                                        ∂L(λ, λn) ∂πi = 0, i = 1, . . . , N ∂L(λ, λn) ∂θπ = 0 ∂L(λ, λn) ∂aij = 0, i = 1, . . . , N j = 1, . . . , N ∂L(λ, λn) ∂θai = 0, i = 1, . . . , N ∂L(λ, λn) ∂bij = 0, i = 1, . . . , N j = 1, . . . , M ∂L(λ, λn) ∂θbj = 0, i = 1, . . . , N (4.18)

och som vi vill lösa för alla πi, aij och bij. Att ta fram dessa kräver en del beräkningar och kommer därför inte att visa alla steg. Hänvisar till [3] s.79-89 för en mer detaljerad beskrivning. När alla beräkningarna är gjorda kommer vi att intuitivt tolka resultaten.

4.3.1 Optimum för πi För att lösa ut πi behöver vi

(38)

∂L(λ, λn) ∂πi = ∂ ∂πi K X k=1 wk N X i=1 ln πiP (x1= si, Y(k)| λn) ! − θπi = 0 ⇐⇒ K X k=1 wk P (x1= si, Y(k)| λn) πi − θ_π_i = 0 (4.19) och ∂L(λ, λn) ∂θπ = 1 − N X i=1 πi= 0. (4.20)

Med hjälp av dessa kan man få ut att

∂L(λ, λn) ∂aij = ∂ ∂aij   K X k=1 wk N X i=1 N X j=1 Tk X t=2 ln aijP (xt−1 = si, xt= sj, Y(k)| λn)  − θai = 0 ⇐⇒ K X k=1 wk Tk X t=2 P (xt−1= si, xt= sj, Y(k)| λn) aij − θai = 0 (4.22)

(39)

och ∂L(λ, λn) ∂θai = 1 − N X i=1 aij = 0 (4.23)

vilket tillsammans ger att

a∗_ij = K X k=1 wk Tk X t=2 P (xt−1= si, xt= sj, Y(k)| λn) K X k=1 wk Tk X t=2 N X j=1 P (xt−1= si, xt= sj, Y(k)| λn) = K X k=1 wk Tk X t=2 P (xt−1= si, xt= sj, Y(k) | λn) K X k=1 wk Tk X t=2 P (xt−1= si, Y(k) | λn) = K X k=1 wkP (Y(k)| λn) Tk X t=2 P (xt−1= si, xt= sj | Y(k), λn) K X k=1 wkP (Y(k)| λn) Tk X t=2 P (xt−1= si| Y(k), λn) = K X k=1 wkP (Y(k)| λn) Tk−1 X t=1 P (xt= si, xt+1= sj | Y(k), λn) K X k=1 wkP (Y(k)| λn) Tk−1 X t=1 P (xt= si| Y(k), λn) (4.24) där i = 1, . . . , N och j = 1, . . . , N. 4.3.3 Optimum för bij För att ta fram bij behöver vi

∂L(λ, λn) ∂bij = ∂ ∂bij K X k=1 wk N X i=1 Tk X t=1 ln b iyt(k) P (xt= si, Y(k)| λn) ! − θ_b_i = 0 ⇐⇒ K X k=1 wk Tk X t=1 P (xt= si, Y(k)| λn)I(yt(k)= vj) bij − θ_b_i = 0 (4.25) där I(y_t(k)= vj) = ( 1 om y_t(k)= vj, 0 annars.

(40)

Vi behöver också ∂L(λ, λn) ∂θbi = 1 − N X i=1 bij = 0 (4.26)

som tillsammans ger att

b∗_ij = K X k=1 wk Tk X t=1 P (xt= si, Y(k)| λn)I(y(k)t = vj) K X k=1 wk Tk X t=1 P (xt= si, Y(k)| λn) = K X k=1 wkP (Y(k) | λn) Tk X t=1 P (xt= si | Y(k), λn)I(yt(k)= vj) K X k=1 wkP (Y(k)| λn) Tk X t=1 P (xt= si | Y(k), λn) (4.27) där i = 1, . . . , N och j = 1, . . . , M. 4.3.4 Att tolka resultaten

Nu när vi har funnit optimum till varje parameter är det lämpligt att försöka tolka vad resultaten egentligen säger. För att göra det enkelt börjar vi med fallet då vi enbart har en observationssekvens, alltså K = 1. Då reduceras de optimala formlerna till

π_i∗= P (x1 = si| Y, λn), i = 1, . . . , N a∗_ij = T −1 X t=1 P (xt= si, xt+1 = sj | Y, λn) T −1 X t=1 P (xt= si | Y, λn) , i = 1, . . . , N j = 1, . . . , N b∗_ij = T X t=1 P (xt= si | Y, λn)I(yt= vj) T X t=1 P (xt= si | Y, λn) i = 1, . . . , N j = 1, . . . , N

(41)

Notera att P (xt= si| Y, λn) förekommer i varje lösning, alltså sannolik-heten att befinna sig i tillstånd si givet observationssekvensen och modellen. Vi tolkar

T X

t=1

P (xt= si | Y, λn) =förväntat antal gånger man är i tillstånd si och om vi exkluderar sista tidpunkten T tolkar vi

T −1 X t=1

P (xt= si | Y, λn) =förväntat antal övergångar från tillstånd si. I formeln för aij har man även P (xt = si, xt+1 = sj | Y, λ) som betyder sannolikheten att vid tidpunkt t vara i tillstånd si och vid tidpunkt t + 1 vara i tillstånd sj givet observationssekvensen samt modellen. Vi kan tolka T −1

X t=1

P (xt= si, xt+1= sj | Y, λn) =förväntat antal övergångar från si till sj. Med dessa kan vi nu konstatera att i specialfallet med enbart en observa-tionssekvens blir

π_i∗=sannolikheten att vara i tillstånd si vid t=1

a∗_ij = förväntat antal övergångar från tillstånd si till tillstånd sj förväntat antal övergångar från tillstånd si

b∗_ij = förväntat antal gånger man är i tillstånd si då man observerat vj förväntat antal gånger man är i tillstånd si

Detta resultat kan vi enkelt att generalisera till fallet då vi har K stycken oberoendeobservationssekvenser. Vi kan se det som att vi beräknar de förvän-tade för varje observationssekvens och viktar dessa med 1/K och summerar alla resultat som kommer att ge det genomsnittliga. I fallet med beroende observationssekvenser viktas de förväntade för varje observationssekvens en-ligt (4.10) istället för 1/K. Man behöver också förstå hur dessa kan beräknas i praktiken. Som vi sett tidigare kan vi uttrycka

P (xt= si| Y, λn) = ˆ

αt(i) ˆβt(i) PN

i=1αˆt(i) ˆβt(i) .

(42)

ˆ

αt(i)och ˆβt(i). Men detta kan vi lösa eftersom P (xt= si, xt+1= sj | Y, λn) = P (xt= si, xt+1= sj, Y | λn) P (Y | λn) = P (xt= si, xt+1= sj, Y | λn) PN i=1 PN j=1P (xt= si, xt+1= sj, Y | λn) = _P_N αt(i)aijbjyt+1βt+1(j) i=1 PN j=1αt(i)aijbjyt+1βt+1(j) . Anledningen till det sista steget är att

P (xt= si, xt+1= sj, Y | λn) =P (y1, . . . , yt, xt= si, xt+1= sj | λn) P (yt+1, . . . , yT | xt= si, xt+1= sj, λn)

=P (y1, . . . , yt, xt= si | λn)P (xt+1= sj | xt= si) P (yt+2, . . . , yT | xt= si, λn)P (yt+1| xt+1= sj) =αt(i)aijβt+1(i)bjyt+1.

Till sist kan vi också visa att (se [3] s.45) αt(i)aijbjyt+1βt+1(j) PN i=1 PN j=1αt(i)aijbjyt+1βt+1(j) = αˆt(i)aijbjyt+1βˆt+1(j) PN i=1 PN j=1αˆt(i)aijbjyt+1βˆt+1(j) . (4.28) Eftersom vi har beräkningseffketiva metoder för att ta fram ˆαt(i) och ˆβt(i) variabler leder det till att man enkelt kan träna en modell. Vi beräknar alltså ˆ

αt(i) och ˆβt(i) för alla tidpunkter och alla tillstånd för en modell λn och använda de optimala uppdateringsformlerna som då ger den nya modellen λn+1 och upprepar denna process tills konvergens.

4.4 Baum-Welch alogritm

Låt oss nu sammanfatta resultaten från detta kapitel för hur man tränar en modell med multipla beroende observationssekvenser i termer av ˆαt(i) och ˆβt(i)variabler. Definitionen för multipla observationssekvenser och dess vikter wk kommer från (4.7), (4.8) och (4.9).

1. Välj den initiala modellen λ1 och beräkna wk för k = 1, . . . , K. Sätt n := 1.

2. Använd den aktuella modellen λn för att beräkna ˆα(k)t (i) och ˆβ (k) t (i) för k = 1, . . . , K och i = 1, . . . , N samt t = 1, . . . , T . Beräkna även P (Y(k)| λ_n) för k = 1, . . . , K med den skalade framåtalgoritmen.

(43)

3. Den nya modellen λn+1 beräknas enligt π∗_i = K X k=1 wkP (Y(k)| λn) Tk−1 X t=1 ˆ α(k)_t (i) ˆβ_t(k)(j) K X k=1 wkP (Y(k)| λn) a∗_ij = K X k=1 wkP (Y(k)| λn) Tk−1 X t=1 ˆ α(k)_t (i)aijb_jy(k) t+1 ˆ β(k)_t+1(j) K X k=1 wkP (Y(k)| λn) Tk−1 X t=1 ˆ α(k)_t (i) ˆβ_t(k)(j) b∗_ij = K X k=1 wkP (Y(k)| λn) Tk−1 X t=1 ˆ α(k)_t (i) ˆβ_t(k)(j)I(y_t(k)= vj) K X k=1 wkP (Y(k)| λn) Tk−1 X t=1 ˆ α(k)_t (i) ˆβ_t(k)(j) och sätt n := n + 1.

4. Om avslutningskriteriet inte är uppfyllt gå till steg 2, annars avsluta algoritmen. Ett avslutningskriterium kan exempelvis vara

P (Y | λn+1) − P (Y | λn) < där > 0 eller att man väljer att avsluta algoritmen efter n iterationer.

Exempel 4.4.1. Med Matlab’s inbyggda funktion hmmgenerate genererades 10 stycken oberoende observationssekvenser med längd 100 utifrån modellen från exempel 2.1.1. I detta fall är alltså wk = 1/K för alla k. Med dessa observationssekvenser vill vi se hur nära man kan komma den ursprungliga modellen från exempel 2.1.1 genom att använda Baum-Welch algoritm. Den initiala modellen λ1 sätts till aij ≈ 1/N och bij ≈ 1/M, närmare bestämt 1

A1 = 0.51 0.49 0.49 0.51 B1 = 0.35 0.33 0.32 0.32 0.33 0.35 . Med funktionen hmmtrain blir den optimala modellen λ∗

A∗= 0.9247 0.0753 0.1335 0.8665 B∗ = 0.4365 0.3692 0.1942 0.1605 0.3657 0.4737

där avslutningskriteriet är satt till = 10−6_{. I den optimala} emissionsmatri-sen B∗ _{ser vi en viss likhet med den “sanna” emissionsmatrisen. Detsamma} kan inte sägas om övergånssannolikhetsmatrisen A∗ _{som inte ser ut som} för-väntat.

1

Matlab’s inbyggda HMM-funktioner använder ej den initiala tillståndsfördelningsvek-torn π.

(44)

Litteraturförteckning

[1] Sean Borman, The Expectation Maximization Algorithm, A Short Tu-torial, 2004.

[2] Xiaolin Li, Marc Parizeau, Training Hidden Markov Models with Mul-tiple Observations - A Combinatorial Method, 2000.

[3] Mikael Nilsson, First Order Hidden Markov Model, Theory and Im-plementation Issues, 2005.

[4] Stephen Tu, Derivation of Baum-Welch Algorithm for Hidden Markov Models.