HK07 – Les 5

HK07 – Les 5

Verborgen Markov modellen

Yves Moreau 3de jr. Burg. Ir. Elektrotechniek Dataverwerking & Automatisatie 2001-2002

Schatting van frequentiematrices

 Schatting van waarschijnlijkheden op basis van tellingen

 Zie bvb. Positie-Specifieke ScoringsMatrix in PSI-BLAST

 Voorbeeld: matrix model van een locaal motief

GACGTG CTCGAG CGCGTG AACGTG CACGTG





















. . . . . .

. . . . . .

. . . . . .

. . . . . .

T G C A

Tel het aantal instanties in de kolom

 Indien er veel (N>>) gealigneerde sites zijn, kunnen we de frequenties schatten als

 Dit is de maximum aannemelijkheidschatting voor  N

n N

n N

n N

n_{A C C G G T T}

A  / ,  / ,  / ,  /



N n n

P

n P n

N n

N n

N n

N P

ML

T G C

A T

T G

G C

C A

A















)

| ( max arg

)

| ( )

, , ,

| ,

, ,

(

















‘Pseudocounts’

 Als er maar een beperkt aantal tellingen is, is de

lmaximumaannemelijkheidschatting niet betrouwbaar (bvb., voor symbolen die niet geobserveerd zijn in de gegevens)

 In zo een geval willen we de observaties combineren met prior kennis

 Stel dat we voor  een Dirichlet prior gebruiken:

 Laten we de Bayesiaanse update berekenen )

( )

| ) (

|

( P n

n n P

P   

D(  |  )

) (

)

( θ |α

P  D 

)

|

)

(

( ) ( ) (

) (

) ( ) ( ) ( ) 1

| (

1

) 1

( 



 

 ^ 









 n

n M Z

n P

n Z n

M Z

n n P

P ^K

i

n

i ^{i i}

D



 ^K 

i

i ⁱ

Z ₁

) 1 (

) ( ) 1

|

(  ^

 

D



 ^K

i

n i ⁱ

n n M

P

) 1

( ) 1

( 

) ( )

| ) (

|

( P n

n n P

P   

D(  |  )

)

| (

)

|

( ^{n }



P

D

 





k

n k i

i PME

i d

n d Z

n   ^{k k} 

 







 ^{ 1}

) (

) 1

| D(

Bayesiaanse update

=1 omdat beide verdelingen genormalizeerd zijn

Berekening van de posterior gemiddelde schatting

A N

n n

Z n

Z _{i i i}

PME

i 

 





  





 

) (

) (

Normalizatie-integrale Z(.)

) 0 ,..., 0 , 1 , 0 ,..., 0

 (



‘Pseudocounts’

 Pseudocounts

 De prior levert een contributie in de schatting in de vorm van pseudogegevens

 Als weinig gegevens beschikbaar zijn, dan speelt de prior een belangrijke rol

 Als veel gegevens beschikbaar zijn, dan spelen de pseudogegevens een verwasloorbare rol



 

 ⁱ  ⁱ _{i i}

PME

i A

A N

n  

 with

Dirichlet mengeling

 Soms worden de gegevens gegenereerd door een heterogeen process (bvb., hydrophobische vs. hydrophilische domeinen in proteïnen, AT-rijke vs. GC-rijke gebieden in DNA)

 In dergelijke situaties zouden we verschillende priors willen gebruiken afhankelijk van de context

 Maar we kennen niet noodzakelijk de context op voorhand

 Een mogelijkheid is het gebruik van een Dirichlet mengeling (de parameters kunnen voortkomen uit m verschillende bronnen)

)

| ( )

,...,

|

( ^{1 m} q_k ^k

P    ^

 ^D

Dirichlet mengeling

 Posterior

 Via de regel van Bayes

nt) (pseudocou )

| ( )

| (

e) (disjuncti )

| (

) ,

| ( )

| (

k k

k

k k

k

n n

P

n P

n P

n P

























D





l

l l

k k k

n P q

n P n q

P ( | )

)

| ) (

|

( 

 





l

k l

l

k k

k k

Z n

Z q

Z n

Z n q

P ( ) / ( )

) (

/ ) ) (

|

(  



 

Dirichlet mengeling

 Integratie om de posterior gemiddelde schatting te berekenen

 De verschillende componenten van de Dirichlet mengeling worden eerst als aparte pseudocounts beschouwd

 Daarna worden die gecombineerd met een gewicht afhankelijk van de aannemelijkheid van de Dirichlet component



 _k ^{k i ik}

PME

i N A

n n

P  

 ( | ) ^



lk lk

l

k k

k k

Z n

Z q

Z n

Z n q

P ( )/ ( )

) (

/ ) ) (

|

(  



 

) /(

)

(n_i _i^k N  A

)

|

( n

P ^k

Schatting van de prior

 De prior ¹, ..., ^k kan geschat worden via maximum likelihood op een grote referentie dataset

 Waarom maximum likelihood en geen Bayesiaanse schatting?

 Onmogelijk om steeds verder priors op te bouwen

 Start met zelf-gekozen prior

 Kennis gebaseerd

 Niet informatief

 Bouw prior op met ML methoden

) ,...,

; ,...,

| ( max

arg )

,...,

; ,...,

( ¹ ₁

) 1 ,...,

; ,..., 1 (

1

1 1 m

M m l

l r

m r

m q q P n r r

m m  















Overzicht

 Verborgen Markov modellen

 Schatting van sequentie- en toestandwaarschijnlijkheden

 Viterbi schatting van het beste toestandpad

 Het voorwaartse algoritme voor de schatting van de waarschijnlijkheid van een sequentie

 Het achterwaartse algoritme voor de schatting van toestandwaarschijnlijkheden

 Parameterschatting voor verborgen Markov modellen

 Parameterschatting met gekende toestandpaden

 Parameterschatting met onbekende toestandpaden

 Viterbi training

 Baum-Welch training

Regular expressions

 Alignering

 Regular expression

 Probleem: regular expression kan niet het verschil maken tussen

 Uitzonderlijk TGCTAGG

 Consensus ACACATC

ACA---ATG TCAACTATC ACAC--AGC AGA---ATC ACCG--ATC

[AT][CG][AC][ACGT]*A[TG][GC]

Verborgen Markov model

A .8 C 0 G 0 T .2

A 0 C .8 G .2 T 0

A .8 C .2 G 0 T 0

A 1 C 0 G 0 T 0

A 0 C 0 G .2 T .8

A 0 C .8 G .2 T 0 A .2

C .4 G .2 T .2

1.0 1.0

.6

.4

.6 .4

1.0 1.0

Sequentiescore

0472 0

8 0 1 8 0 1 1 6 0 4 0

6 0 8 0 1 8 0 1 8 0 )

(

.

. .

. .

. .

. .

P

























 ACACATC 

Transitiewaarschijnlijkheden

Emissiewaarschijnlijkheden

Log odds

 In praktijk, logaritmes voor scalering en normalizatie t.o.v. random

 Log odds voor sequentie S:

A: 1.16 T:-0.22

C: 1.16 G:-0.22

A: 1.16

C:-0.22 A: 1.39 G:-0.22

T: 1.16

C: 1.16 G:-0.22 A:-0.22

C: 0.47 G:-0.22 T:-0.22

0 0

-0.51

-0.92

-0.51 -0.92

0 0

25 . 0 log )

( 25 log

. 0

)

log P(S P S L

L  

Log odds

65 . 6

16 . 1 0 16 . 1 39 . 1 51 . 0 47 . 0

51 . 0 16 . 1 0 16 . 1 0 16 . 1 )

( odds log





























ACACATC

Sequentie Log odds

ACAC--ATC (consensus) 6.7

ACA---ATG 4.9

TCAACTATC 3.0

ACAC--AGC 5.3

AGA---ATC 4.9

ACCG--ATC 4.6

TGCT--AGG (uitzonderlijk) -0.97

Markov keten

 Sequentie:

 Voorbeeld van Markov keten

 Probabilistisch model van een DNA sequentie

A

C G

T a_st  P(x_i  t | x_i1  s)

Transitiewaarschijnlijkheden





with ,...,

, ₂

(e.g.,

1 x x x A C G T

x

x  _{L i} 

A A 

Markov eigenschap

 Waarschijnlijkheid van een sequentie via de regel van Bayes

 Markov eigenschap

 “De toekomst is enkel functie van het heden en niet van het verleden”

) ( )

,...,

| (

) ,...,

| (

) ,..., ,

( )

(

1 1

2 1

1 1

1 1

x P x

x x

P x x

x P

x x

x P x

P

L L

L L

L L























 



L

i

x x

L L

L L

i

a i

x P

x P x x

P x

x P x

x P x

P

2 1

1 1

2 2

1 1

) 1

(

) ( )

| ( )

| (

)

| ( )

( 

Begin en einde van een sequentie

 Waarschijnlijkheidsberekening is niet homogeen

 Lengteverschillen worden niet gemodeleerd

 Oplossing

 Modelering van begin en einde van de sequentie

 De waarschijnlijkheid om een sequentie van een bepaalde lengte te observeren

daalt exponentieel met de lengte van de sequentie

)

| (

) ( ₁

t x

P a

s x

P a

L t

s









 



A

C G

T

 

Verborgen Markov model

A .8 C 0 G 0 T .2

A 0 C .8 G .2 T 0

A .8 C .2 G 0 T 0

A 1 C 0 G 0 T 0

A 0 C 0 G .2 T .8

A 0 C .8 G .2 T 0 A .2

C .4 G .2 T .2

1.0 1.0

.6

.4

.6 .4

1.0 1.0

Sequentiescore

0472 0

8 0 1 8 0 1 1 6 0 4 0

6 0 8 0 1 8 0 1 8 0 )

(

.

. .

. .

. .

. .

P

























 ACACATC 

Transitiewaarschijnlijkheden

Emissiewaarschijnlijkheden

Verborgen Markov model

 In een verborgen Markov model, observeren wij de symboolsequentie x maar we willen de verborgen toestandsequentie (het pad ) reconstrueren

 Transitiewaarschijnlijkheden

 Emissiewaarschijnlijkheden

 Gezamenlijke waarschijnlijkheid van sequentie en pad )

|

( l ₁ k

P

a_kl  _i  _i 

)

| (

)

(b P x b k

e_k  _i  _i 



 ^L

i

i _{i i}

i x a

e a

x P

1

0 ₁ ( ) ₁

) ,

(  _{   }

Casino (I) - probleemstelling

 De casino gebruikt meestal een faire dobbelsteen maar schakelt soms over naar een vervalste dobbelsteen

 We observeren de uitkomst x van de opeenvolgende worpen maar willen weten wanneer de dobbelsteen fair of vervalst was (pad )

1: 1/6 2: 1/6 3: 1/6 4: 1/6 5: 1/6 6: 1/6

1: 1/10 2: 1/10 3: 1/10 4: 1/10 5: 1/10 6: 1/2 0.05

0.1 0.95 0.9

Fair Vervalst

Schatting van sequentie- en

toestandwaarschijnlijkheden

Het Viterbi algoritme

 We zoeken het meest waarschijnlijke pad ^*

 Dit probleem kan aangepakt worden via dynamisch programmeren

 Laten we vk(i) definieren als de waarschijnlijkheid van het meest waarschijnlijke pad die in toestand k eindigt met symbool i

 Dan kunnen we deze waarschijnlijkheid recursief berekenen als

)

| ( max arg

) , ( max

* arg P x  P  x

    

) )

( ( max )

( )

1

( _{1 k kl}

i k l

l i e x v i a

v   _

) ,

,..., (

max )

( ₁

,..., ₁

1

k x

x P i

v_{k i i}

i







 



Het Viterbi algoritme

 Viterbi algoritme groeit dynamisch het beste pad van begin tot einde

 Initiele conditie : sequentie in begintoestand

 Traceback pointers volgen het beste pad (= decodering)

) ( ptr :

) 1 ,..., (

Traceback

) )

( ( max arg

) )

( ( max )

, ( :

Einde

) )

1 (

( max arg

) ( ptr

) )

1 (

( max )

( )

( :

) ,..., 1 (

Recursie

0 ,

0 )

0 ( , 1 ) 0 ( :

) 0 (

tie Initialisa

*

* 1

0

*

0

* 0

i i i

k k k

L

k k k

kl k k

i

kl k k

i l l

k

L i

a L v

a L v x

P

a i

v l

a i

v x

e i

v L

i

k v

v i







































Casino (II) - Viterbi

Het voorwaartse algoritme

 Het voorwaartse algoritme laat toe de waarschijnlijkheid P(x) van een sequentie te schatten t.o.v. van een verborgen Markov model

 Dit is belangrijk voor het berekenen van posterior

waarschijnlijkheden en voor het vergelijken van Markov modellen

 De som over alle paden (exponentieel veel) kan berekend worden via dynamish programmeren

 Laten we fk(i) definieren als de waarschijnlijkheid van de sequentie voor de paden die in toestand k eindigen bij symbool i

 Dan kunnen we deze waarschijnlijkheid recursief berekenen als







 ) , ( )

(x P x

P

) ,

,..., (

)

(i P x₁ x k

f_k  _i _i 

Het voorwaartse algoritme

 Het voorwaartse algoritme groeit de totale waarschijnlijkheid dynamisch voorwaarts (van begin tot einde)

 Initiele conditie : sequentie in begintoestand

 Einde : alle toestanden convergeren naar eindtoestand























k

k k

k

kl k

i l l

k

a L f

x P

a i

f x

e i

f L

i

k f

f i

0 0

) ( )

( :

Einde

) 1 (

) ( )

( :

) ,..., 1 (

Recursie

0 ,

0 )

0 ( ,

1 ) 0 ( :

) 0 (

tie Initialisa

Het achterwaartse algoritme

 Het achterwaartse algoritme laat toe de waarschijnlijkheid te

berekenen van de gehele sequentie gezamenlijk met de conditie dat symbool x_i in toestand k zit

 Dit is belangrijk om de waarschijnlijkheid van een bepaalde toestand op symbool x_i te berekenen

 P(x1,...,xi,_i=k) kan berekend worden via voorwaartse algoritme fk(i)

 Laten we b_k(i) definieren als de waarschijnlijkheid van de rest van de sequentie voor de paden die in toestand k door symbool xi passeren

)

| ,...,

( ) ,

,..., (

) ,

,...,

| ,...,

( ) ,

,..., (

) ,

(

1 1

1 1

1

k x

x P k x

x P

k x

x x

x P k x

x P k

x P

i L

i i

i

i i L

i i

i i





























)

| ,...,

( )

(i P x ₁ x k

b_k  _{i L} _i 

Het achterwaartse algoritme

 Het achterwaartse algoritme groeit de waarschijnlijkheid bk(i) dynamisch achterwaarts (van einde tot begin)

 Grensconditie : start in eindtoestand

 Eens beide voorwaartse en achterwaartse waarschijnlijkhedeb

beschikbaar zijn, kunnen we de posterior waarschijnlijkheid van een toestand berekenen























l

l l

l l

l i

l kl k

k k

b x e a x

P

i b x

e a i

b L

i

k a

L b L

i

) 1 ( ) ( )

( :

Einde

) 1 (

) (

) ( :

) 1 ,..., 1 (

Recursie

, )

( :

) (

tie Initialisa

1 0

1 0

) ( ) ) (

|

( f i b i

x k

P _i   ^{k k}

Posterior decodering

 In plaats van het meest waarschijnlijke pad te gebruiken voor de decodering (Viterbi), kunnen we het pad van de meest waarschijnlijke toestanden gebruiken

 Het pad ^{^} kan “illegaal” zijn (P(^{^}|x)=0)

 Deze aanpak kan gebruikt worden als men

geinteresseerd is in een functie g(k) van de toestand (e.g., labeling)

)

| (

max

ˆ arg P _i k x

i  k  







k

i k x g k

P x

i

G( | ) ( | ) ( )

Casino (III) – posterior decodering

 Posterior waarschijnlijkheid van de toestand “fair” t.o.v.

de dobbelsteenworpen

Casino (IV) – posterior decodering

 Nieuwe situatie : P(x_i+1 = FAIR | x_i = FAIR) = 0.99

 Viterbi decodering kan de vervalsing niet detecteren uit 1000 worpen, posterior decodering wel

1: 1/6 2: 1/6 3: 1/6 4: 1/6 5: 1/6 6: 1/6

1: 1/10 2: 1/10 3: 1/10 4: 1/10 5: 1/10 6: 1/2 0.01

0.1 0.99 0.9

Fair Vervalst

Parameterschatting

voor verborgen Markov modellen

Keuze van de topologie

 Voor de parameterschatting, veronderstellen we dat de topologie van het verborgen Markov model gekend is

 Keuze van topologie is een belangrijke designkeuze

 Duurmodellen

 “Silent states” voor gaps

Parameterschatting met gekende paden

 Verborgen Markov model met parameters  (transitie- en emissiewaarschijnlijkheden)

 Trainingsverzameling

D

 Score van het model is aannemelijkheid van de parameters gezien de gegevensverzameling





 ^N

j

j

N P x

x x

P

1

1,..., | ) log ( | ) (

log )

, (

Score

D

Parameterschatting met gekende paden

 Indien de toestandpaden gekend zijn, worden de

parameters geschat via tellingen (hoe vaak wordt een transitie gebruikt, hoe vaak wordt een symbool

geproduceerd in een toestand)

 Gebruik van ‘pseudocounts’ indien nodig

 Akl = aantal transities van k naar l in trainingsverzameling + pseudocount rkl

 Ek(b) = aantal emissies van b uit k in trainingsverzameling + pseudocount r (b)







b

k k k

l

l k kl kl

b E

b b E

A e a A

) (

) ) (

( and

Parameterschatting met onbekende paden : Viterbi training

 Strategie : iteratieve methode

 Veronderstel dat de parameters gekend zijn en vind het beste pad

 Gebruik de Viterbi decodering om de parameters te schatten

 Itereer tot convergentie

 Viterbi training maximizeert de aannemelijkheid van de parameters niet

 Viterbi training convergeert exact in een eindig aantal stappen

)) (

),..., (

,

| ,...,

( max

arg ^{1 * 1 *}

Vit P x xN _ x _ xN

   

 

Parameterschatting met onbekende paden : Baum-Welch training

 Strategie : gelijklopend met Viterbi maar de verwachte waarden voor het aantal transities en emissies wordt gebruikt (in plaats van enkel het beste pad te gebruiken)

 Voor de transities

 Voor de emissies

)

| (

) 1 (

) (

) ) (

,

| ,

( ₁ ¹

 



 P x

i b x

e a i x f

l k

P _i  _i   ^{k kl l i l} 

 



 _{ }

j i

j l j i l kl j

j k

j i

j i

i

kl f i a e x b i

x x P

l k

P

A ( ) ( ) ( 1)

)

| ( ) 1

,

| ,

( _{1 1}

 





 

 







j i x b

j k j

j k

j ix b

j i

k

j j

i b i x f

x P k P

b E

}

| { }

| {

) ( ) ) (

| ( ) 1

,

| (

)

(   

Parameterschatting met onbekende paden : Baum-Welch training

 Initialisatie : Kies willekeurige modelparameters

 Recursie :

 Zet alle transitie- en emissievariabelen op hun pseudocount

 Voor alle sequenties j = 1,...,n

 Bereken fk(i) voor sequentie j met het voorwaartse algoritme

 Bereken bk(i) voor sequentie j met het achterwaartse algoritme

 Voeg de bijdrage toe aan A en E

 Bereken de nieuwe modelparameters akl =A_kl/kl’ en ek(b)

 Bereken de log-aannemelijkheid van het model

 Einde : stop wanneer de log-aannemelijkheid niet meer verandert dan een gegeven drempel of het maximum aantal iteraties is overschreden

Casino (V) – Baum-Welch training

1: 0.19 2: 0.19 3: 0.23 4: 0.08 5: 0.23 6: 0.08

1: 0.07 2: 0.10 3: 0.10 4: 0.17 5: 0.05 6: 0.52 0.27

0.29 0.73 0.71

Fair Vervalst

1: 0.17 2: 0.17 3: 0.17 4: 0.17 5: 0.17 6: 0.15

1: 0.10 2: 0.11 3: 0.10 4: 0.11 5: 0.10 6: 0.48 0.07

0.12 0.93 0.88

Fair Vervalst

1: 1/6 2: 1/6 3: 1/6 4: 1/6 5: 1/6 6: 1/6

1: 1/10 2: 1/10 3: 1/10 4: 1/10 5: 1/10 6: 1/2 0.05

0.1

0.9 0.95

Fair Vervalst

Oorspronkelijk model

300 worpen 30000 worpen

Numerieke stabiliteit

 Veel uitdrukkingen bevatten producten van waarschijnlijkheden

 Dit veroorzaakt underflows als we deze uitdrukkingen berekenen

 Voor Viterbi kan dit opgelost worden door met logaritmes te werken

 Voor de voorwaartse en achterwaartse algoritmen kan er met logaritmes gewerkt worden via een benadering of kan er met een herschalering van de variabelen gewerkt worden

Samenvatting

 Verborgen Markov modellen

 Schatting van sequentie- en toestandwaarschijnlijkheden

 Viterbi schatting van het beste toestandpad

 Het voorwaartse algoritme voor de schatting van de waarschijnlijkheid van een sequentie

 Het achterwaartse algoritme voor de schatting van toestandwaarschijnlijkheden

 Parameterschatting voor verborgen Markov modellen

 Parameterschatting met gekende toestandpaden

 Parameterschatting met onbekende toestandpaden

 Viterbi training

 Baum-Welch training