hmm.pdf

``Wst ep do obliczeniowej biologii molekularnej''

(J. Tiuryn, wyklad nr.11, 18 stycznia 2006)

Spis tresci

7 Ukryte modele Markowa 75

7.1 Algorytm Viterbiego . . . . . . . . . . . . . . . . . . . . . . . 77

7.2 Prawdopodobienstwo wyemitowania zadanego slowa . . . . . . 78

7.2.1 Algorytm preksowy . . . . . . . . . . . . . . . . . . . 78

7.2.2 Algorytm suksowy . . . . . . . . . . . . . . . . . . . . 78

7.3 Estymacje parametrow . . . . . . . . . . . . . . . . . . . . . . 79

7.3.1 Algorytm Bauma-Welcha . . . . . . . . . . . . . . . . . 80

7 Ukryte modele Markowa

Zaczniemy od denicji lancucha Markowa. B edziemy si e wyl acznie zaj-

mowac skonczonymi lancuchami Markowa dzialaj acymi w dyskretnym cza-

sie. Niech Q 6=;, b edzie skonczonym zbiorem (zbior stanow). Pewien

stan k 0 2Q jest wyrozniony jako stan pocz atkowy. Lancuch Markowa jest

zadany przez macierz przejsc M = (p k;l ) k;l2Q , ktora dla k; l2Q podaje

prawdopodobienstwo p k;l przejscia ze stanu k do stanu l. M musi spelniac

nast epuj acy warunek (tzn. macierz M musi byc stochastyczna): dla kazdego

k2Q mamy

p k;l = 1:

l2Q

Lancuch Markowa opisuje pewien uklad , ktory w kazdym momencie moze

si e znajdowac tylko w jednym ze stanow k2Q. Uklad ten obserwujemy w

dyskretnych chwilach czasowych t = 0; 1; : : :. Przyjmujemy, ze na pocz atku

uklad znajduje si e w stanie pocz atkowym k 0 . Jesli w danym momencie t,

uklad znajduje si e w stanie k, to w momencie t + 1 przechodzi on do stanu l

z prawdopodobienstwem p k;l . Podstawow a cech a lancucha Markowa jest to,

ze nast epny stan zalezy tylko od obecnego stanu i nie zalezy od wartosci t,

ani od historii dojscia do obecnego stanu.

Ukryte modele Markowa (HMM, od \hidden Markov models") stanowi a

pewne rozszerzenie denicji lancucha Markowa. Niech b edzie alfabetem.

B edziemy rozwazac lancuchy Markowa mog ace si e komunikowac z otoczeniem

poprzez emitowanie ci agow liter z alfabetu . Tak wi ec maj ac dany HMM

b ed acy w stanie k2Q, emituje on symbol x2 z prawdopodobienstwem

e k (x) oraz przechodzi do stanu l z prawdopodobienstwem p k;l . Ukryte mod-

ele Markowa znane s a w literaturze informatycznej rowniez pod nazw a prob-

abilistycznych automatow z wyjsciem. Jesli chcemy aby w kazdym stanie

k2Q emitowany byl jakis symbol, to musimy przyj ac

x2 e k (x)1. Przyjmujemy, ze to co daje si e obserwowac to symbole emi-

towane przez uklad, a nie stany wewn etrzne ukladu (st ad nazwa \ukryte").

Przyklad 7.0.1 Przykladem ukrytego lancucha Markowa jest nieuczciwe

kasyno gry, ktore ma dwa rodzaje kosci: uczciw a kostk e do gry (z praw-

dopodobienstwem 1/6 wyrzuca si e kazda z szesciu mozliwych wartosci) oraz

nieuczciw a kostk e (dla ktorej prawdopodobienstwo wyrzucenia szostki wynosi

1/2, a dla pozostalych liczb 1/10). Tak wi ec mamy dwa stany: F (uczciwa

kostka) i L (nieuczciwa). Uklad moze zmieniac swoj stan z pewnym praw-

dopodobienstwem, ale my stanu nie mozemy zaobserwowac. Jedynie widzimy

ci ag liczb b ed acych wynikiem rzutow kostk a.

Mozemy, na przyklad, miec do czynienia z nast epuj acym modelem: p F;F =

0:95; p L;L = 0:9; p F;L = 0:05; p L;F = 0:1; ponadto prawdopodobienstwo

emisji jest zdeniowane nast epuj aco: e F (x) = 1=6 dla x2f1; 2; 3; 4; 5; 6g

oraz e L (6) = 0:5 i e L (x) = 0:1 dla x2f1; 2; 3; 4; 5g

Przyklad 7.0.2 (Wyspy CpG) Wiadomo, ze w ludzkim genomie, jesli po-

jawi si e para nukleotydow CG (tzn G pojawia si e tuz za C w jednej nici, w

kierunku 5' do 3'), to nukleotyd C zwykle ulega zmianie (pod wplywem mety-

lacji) i z duzym prawdopodobienstwem zostaje zmutowany na T. Tak wi ec

pary CpG 1 zwykle pojawiaj a si e rzadziej w genomie czlowieka. Z pewnych

biologicznych powodow proces metylacji jest zatrzymywany na krotkich od-

cinkach genomu w poblizu miejsc rozpoczynaj acych rejony koduj ace bialko

(lub tez w poblizu tzw. odcinkow promotorowych). Takie odcinki s a nazy-

wane wyspami CpG. Rozpoznawanie wysp CpG moze byc wykorzystane jako

wskazowka przy poszukiwaniu biologicznie interesuj acych fragmentow genomu.

Mamy tutaj do czynienia z ukrytym modelem Markowa. Ma on osiem

stanow A + ; C + ; G + ; T + ; A ; C ; G ; T . Stan A + oznacza, ze znajdujemy si e

w rejonie wyspy CpG i czytamy nukleotyd A. Podobnie, stan G oznacza, ze

znajdujemy si e poza rejonem wyspy CpG i czytamy G. Emitowane symbole:

b ed ac w stanie A (gdzie 2f; +g), emitujemy A z prawdopodobienstwem

1 Litera 'p' pomi edzy C oraz G oznacza, ze chodzi o kolejne nukleotydy w tej samej

nici, a nie odpowiadaj ace sobie pary komplementarnych nukleotydow w DNA.

x2 e k (x) = 1.

W ogolnosci, jesli dopuszczamy sytuacje, ze w pewnych stanach (z pewnym

prawdopodobienstwem, nic nie jest emitowane, to przyjmujemy slabsze zalozenie

1; kazdy inny symbol jest emitowany w stanie A z prawdopodobienstwem

0. Podobnie dla innych stanow.

Niech S2 oraz 2Q b ed a niepustymi ci agami rownej dlugosci

n =jSj=jj> 0. Prawdopodobienstwo tego ze S zostanie wyemitowane

oraz uklad b edzie zmienial stany wedlug kolejnosci wynosi

P (S; ) =

e (t+1) (S(t + 1))p (t);(t+1) ;

t=0

gdzie w powyzszym wzorze przyjmujemy, ze (0) = k 0 , jest stanem pocz atkowym. 2

7.1 Algorytm Viterbiego

Dane slowo S2 zaobserwowane na wyjsciu HMM. Chcemy znalezc na-

jbardziej prawdopodobny ci ag stanow, ktory doprowadzil do wyemitowania

S, czyli znalezc takie, ze

P (S; ) = maxfP (S; )j2Q ;jj=jSjg:

Zastosujemy metod e dynamicznego programowania. Dla 0 < ijSjoraz

k2Q, niech v(i; k) b edzie prawdopodobienstwem najbardziej optymalnej

drogi koncz acej si e w stanie k, ktora doprowadzila do wyemitowania preksu

S[1::i] z maksymalnym prawdopodobienstwem. Tzn.

v(i; k) = maxfP (S[1::i]; )j2Q i ; (i) = kg:

Funkcj e v rozszerzamy dla przypadku i = 0, deniuj ac

(

1 gdy k = k 0 ;

0 gdy k 6= k 0 :

v(0; k) =

Twierdzenie 7.1.1 Dla 0 < ijSjoraz k2Q zachodzi

v(i; k) = e k (S(i))max

l2Q

[v(i1; l)p l;k ]:

Wowczas maksymalne prawdopodobienstwo P (S; ) znajduje si e ze wzoru

P (S; ) = max

k2Q [v(jSj; k)]:

2 Czasami przyjmuje si e, ze HMM ma wyrozniony stan koncowy, do ktorego uklad

przechodzi po wyemitowaniu slowa S. Wprowadzenie takiego stanu niewiele zmienia

rozwazania, wi ec dla prostoty b edziemy go opuszczac.

Optymaln a drog e (moze ich byc wi ecej niz jedna) znajduje si e przez za-

pami etywanie wskaznikow, dla ktorych maksimum jest realizowane w rownaniu

rekurencyjnym z Twierdzenia 7.1.1. Twierdzenie to w naturalny sposob

prowadzi do algorytmu wyznaczania najbardziej prawdopodobnej sciezki emi-

tuj acej zadane slowo. Jest to tzw. algorytm Viterbiego.

Twierdzenie 7.1.2 Algorytm Viterbiego wyznacza najbardziej prawdopodobn a

sciezk e emituj ac a przez HMM dane slowo S w czasie O(jSjjQj 2 ) i przestrzeni

O(jSjjQj), gdzie Q jest zbiorem stanow modelu HMM.

7.2 Prawdopodobienstwo wyemitowania zadanego slowa

Zajmiemy si e algorytmem obliczania prawdopodobienstwa P (S), wyemitowa-

nia slowa S przez dany HMM. Mamy

P (S) =

P (S; );

gdzie w powyzszej sumie przyjmujemy, ze P (S; ) = 0 dla drog , takich ze

jSj6=jj.

7.2.1 Algorytm preksowy

Dla 0ijSjoraz k2Q niech f (i; k) b edzie prawdopodobienstwem

wyemitowania preksu S[1::i] przy dodatkowym zalozeniu, ze droga prowadz aca

do tej emisji konczy si e w stanie k. Czyli

f (i; k) =

P (S[1::i]; ):

f j (i)=kg

Podobnie jak dla algorytmu Viterbiego mamy

(

1 gdy k = k 0 ;

0 gdy k 6= k 0 :

f (0; k) =

Ponizsze twierdzenie sugeruje pewien sposob obliczania wartosci f (i; k)

oraz prawdopodobienstwa P (S).

Twierdzenie 7.2.1 Dla 0 < ijSjoraz k2Q mamy

f (i; k) = e k (S(i))

f (i1; l)p l;k :

l2Q

Prawdopodobienstwo wyemitowania slowa S wynosi wowczas

P (S) =

f (jSj; k):

k2Q

7.2.2 Algorytm suksowy

Prawdopodobienstwo P (S) mozemy rowniez obliczyc uzywaj ac suksow. Niech

0ijSj= m i niech b(i; k) oznacza prawdopodobienstwo wyemitowania

suksu S[i+1::m], zakladaj ac, ze stan pocz atkowy modelu jest k. Mamy wi ec

b(jSj; k) = 1, dla k2Q (bo oczywiscie prawdopodobienstwo wyemitowania

pustego slowa przy dowolnym stanie pocz atkowym jest rowne 1). Ponadto

oczywiscie mamy

P (S) = b(0; k 0 ):

Ponizsze twierdzenie podaje sposob obliczania wartosci b(i; k).

Twierdzenie 7.2.2 Dla 0i <jSjoraz k2Q mamy

b(i; k) =

p k;l e l (S(i + 1))b(i + 1; l):

l2Q

Twierdzenie 7.2.3 Opieraj ac si e na rownaniach z Twierdzen 7.2.1 oraz

7.2.2, wartosci f (i; k) oraz b(i; k) mozna obliczyc w czasie O(jSjjQj 2 ) i

w pami eci O(jSjjQj).

Uzywaj ac funkcji f oraz b mozemy wyrazic warunkowe prawdopodobienstwo

tego, ze na i-tym kroku dany HMM byl w stanie k, pod warunkiem ze wyemi-

towane zostalo slowo S:

P ((i) = kjS) = P ((i) = k & S)

P (S)

= f (i; k)b(i; k)

P (S)

(7.1)

7.3 Estymacje parametrow

Problem estymacji parametrow polega na jak najlepszym doborze praw-

dopodobienstw przejsc i emisji w oparciu o zaobserwowany zbior slow S 1 ; : : : ; S n ,

wyemitowanych przez pewien HMM. Slowa te s a nazywane ucz acymi b adz

treningowymi sekwencjami. Zakladamy, ze znamy alfabet nad ktorym

tworzone s a slowa S 1 ; : : : ; S n oraz zbior Q stanow modelu HMM. Nie znamy

natomiast prawdopodobienstw p k;l oraz e k (x), dla k; l2Q oraz x2.

Tak wi ec, naszym zadaniem jest dobranie tych prawdopodobienstw tak, aby

zmaksymalizowac prawdopodobienstwo wyemitowania S 1 ; : : : ; S n . NiechM

oznacza pewien HMM nad alfabetem i zbiorem stanow Q. Przez P M (S 1 & : : : &S n )

b edziemy oznaczac prawdopodobienstwo wyemitowania slow S 1 ; : : : ; S n w

modeluM. Poniewaz slowa s a emitowane niezaleznie, to

P M (S 1 & : : : & S n ) =

P M (S j ):

j=1

Plik z chomika:

Inne pliki z tego folderu:

Inne foldery tego chomika: