1 Modelos ocultos de Markov (HMM)Matías Ison Patricia Muzulin Laura Kamenetzky Introducción a la Biología Computacional 5 de Octubre de 2004

2 Organización de la charlaFundamentos teóricos Algunos conceptos de probabilidades HMMS: Formalismo Preguntas a resolver Algoritmos utilizados Aplicaciones Búsqueda de genes Uso del Programa HMMER Búsqueda de familias de proteínas

3 Empecemos con un algunos conceptos de probabilidadesLa probabilidad de ocurrencia de un evento X es P(X) Si queremos calcular la probabilidad de un conjunto de eventos independientes: Las probabilidades se multiplican Probabilidad condicional P(X|Y) : La probabilidad de que ocurra X dado Y Probabilidad conjunta P(X,Y) = P(X|Y)P(Y) P(X,Y|Z) = P(X|Y,Z)P(Y|Z) Probabilidad marginal P(X) = S Y P(X|Y)P(Y)

4 Probabilidad a posterioriUsualmente queremos conocer la probabilidad de una observación O dada una suposición (modelo) M; P(O|M) Problema inverso: Dado O queremos conocer la probabilidad de que M sea correcto: La probabilidad a posteriori P(M|O) Teorema de Bayes: Para dos eventos aleatorios e independientes X, Y P(X|Y) = P(Y|X)P(X)/P(Y) P(M|O) = P(O|M)P(M)/P(O)

5 Definición de un modelo oculto de MarkovUn modelo oculto de Markov (HMM) es un conjunto finito de estados. Las transiciones entre estados están dadas por un conjunto de probabilidades de transición. En cualquier estado particular, la observación puede ser generada, de acuerdo a la distribución de probabilidades de emisión. Sólo el resultado observable, no el estado, es visible a un observador externo por lo que los estados están “ocultos”. Est. 2 Est.1

6 Ejemplos de modelos ocultos de MarkovTexto escrito por Shakespeare que en algunas partes ha sido editado por un mono Un casino tiene dos dados, uno cargado y el otro no. Alterna entre ambos Una secuencia de DNA con segmentos que codifican y otros que no

7 Ejemplos (continuación)Caso Observables Estado oculto Texto alfabeto Shakespeare/mono Dado dado justo (F) /dado cargado (L) DNA A,C,G,T codifica/no-codifica (bases)

8 Ejemplo: El Casino deshonestoUn casino tiene dos dados: Dado “justo” P(1) = P(2) = P(3) = P(5) = P(6) = 1/6 Dado cargado P(1) = P(2) = P(3) = P(5) = 1/10 P(6) = 1/2 El casino alterna entre el dado justo y el cargado una vez cada 20 turnos Juego: Apostamos $1 Tiramos (siempre con un dado justo) El Casino tira (tal vez con un dado justo, tal vez con uno cargado) EL número mas alto gana $2

9 Pregunta # 1 – EvaluaciónDada Una secuencia de tiradas del casino PREGUNTA ¿Cuan probable es esta secuencia, dado nuestro modelo de como opera el casino? Este es el problema de EVALUACIÓN en HMMs

10 Pregunta # 2 – DecodificaciónDada Una secuencia de tiradas del casino PREGUNTA ¿Que partes de la secuencia fueron generadas por el dado cargado y cuales por el dado justo? Este es el problema de DECODIFICACIÓN en HMMs

11 Pregunta # 3 – AprendizajeDada Una secuencia de tiradas del casino PREGUNTA ¿Cómo está cargado el dado cargado? ¿Cada cuanto cambia el casino entre el dado cargado y el justo? Este es el problema de APRENDIZAJE en HMMs

12 El HMM del casino deshonestoJUSTO CARGADO 0.05 0.95 P(1|F) = 1/6 P(2|F) = 1/6 P(3|F) = 1/6 P(4|F) = 1/6 P(5|F) = 1/6 P(6|F) = 1/6 P(1|L) = 1/10 P(2|L) = 1/10 P(3|L) = 1/10 P(4|L) = 1/10 P(5|L) = 1/10 P(6|L) = 1/2

13 Definición de un modelo oculto de MarkovDefinición: Un modelo oculto de Markov (HMM) Alfabeto  = { b1, b2, …, bM } Conjunto de estados  = { 1, ..., K } Probabilidades de transición entre dos estados cualesquiera aij = prob. de transición del estado i al estado j ai1 + … + aiK = 1, para todos los estados i = 1…K Probabilidades iniciales a0i a01 + … + a0K = 1 Probabilidades de emisión dentro de cada estado ei(b) = P( xi = b | i = k) ei(b1) + … + ei(bM) = 1, para todos los estados i = 1…K 1 2 K …

14 Un HMM no tiene memoria En cada paso de tiempo t,Lo único que afecta los futuros estados es el estado actual t P(t+1 = k | “cualquier cosa que pasó”) = P(t+1 = k | 1, 2, …, t, x1, x2, …, xt) = P(t+1 = k | t) 1 2 K …

15 Secuencias de observables y camino de estadosDada una secuencia de observaciones x = x1……xN, Y la secuencia de estados  = 1, ……, N 1 2 K … 1 1 2 K … 1 2 K … 1 2 K … … 2 2 K x1 x2 x3 xN

16 Verosimilitud de un camino de estadosDados una secuencia x = x1……xN Y un camino de estados  = 1, ……, N, Para averiguar cuan probable es el Camino: (dado nuestro HMM) P(x, ) = P(x1, …, xN, 1, ……, N) = P(xN, N | N-1) P(xN-1, N-1 | N-2)……P(x2, 2 | 1) P(x1, 1) = P(xN | N) P(N | N-1) ……P(x2 | 2) P(2 | 1) P(x1 | 1) P(1) = a01 a12……aN-1N e1(x1)……eN(xN) 1 2 K … x1 x2 x3 xN

17 Ejemplo: El casino deshonestoImaginemos que la secuencia de dados es: x = 1, 2, 1, 5, 6, 2, 1, 6, 2, 4 Entonces, cuan probable es = Fair, Fair, Fair, Fair, Fair, Fair, Fair, Fair, Fair, Fair? (imaginos inicialmente a0Fair = ½, aoLoaded = ½) ½  P(1 | Fair) P(Fair | Fair) P(2 | Fair) P(Fair | Fair) … P(4 | Fair) = ½  (1/6)10  (0.95)9 = = 0.5  10-9

18 Ejemplo: El casino deshonestoEntonces, la probabilidad de que el dado sea justo En toda la corrida es solamente  10-9 Pero…¿Cual es la probabilidad de = Loaded, Loaded, Loaded, Loaded, Loaded, Loaded, Loaded, Loaded, Loaded, Loaded? ½  P(1 | Loaded) P(Loaded, Loaded) … P(4 | Loaded) = ½  (1/10)8  (1/2)2 (0.95)9 = = 7.9  10-10 Entonces, despues de todo es 6.59 veces mas probable que el dado haya sido siempre justo, a que haya sido siempre trucho.

19 Ejemplo: El casino deshonestoImaginemos que la secuencia de dados es: x = 1, 6, 6, 5, 6, 2, 6, 6, 3, 6 Entonces, cuan probable es  = F, F, …, F? ½  (1/6)10  (0.95)9 = 0.5  10-9, igual que antes ¿Cuan probable es = L, L, …, L? ½  (1/10)4  (1/2)6 (0.95)9 = = 0.5  10-7 Entonces, es 100 veces mas probable que el dado esté cargado

20 Las tres grandes preguntas sobre HMMsEvaluación DADO un HMM M, y una secuencia x, ENCUENTRE Prob[ x | M ] Decodificación ENCUENTRE la secuencia de estados  que maximiza P[ x,  | M ] Aprendizaje DADOS un HMM M, con probs transición/emisión desconocidas, y una secuencia x, ENCUENTRE los parámetros  = (ei(.), aij) que maximizan P[ x |  ]

21 Que no nos confunda la notación!P[ x | M ]: La probabilidad de que la secuencia x haya sido generada por el modelo El modelo es: arquitectura (#estados, etc) + parámetros  = aij, ei(.) Entonces, P[ x |  ], y P[ x ] son lo mismo, en lo que se refiere a la arquitectura y al modelo entero De igual forma, P[ x,  | M ] y P[ x,  ] son lo mismo En el problema de APRENDIZAJE siempre escribimos P[ x |  ] para enfatizar que buscamos el  que maximiza P[ x |  ]

22 HMM: Línea temporal xT x1 xt xt-1 xt+1 1 t+1 T t t-1 tiempo Las flechas indican dependencias estocásticas. Los ’s son estados ocultos, acada uno depende solamente del estado previo. Proceso Markoviano. Los x’s son observaciones, dependen solamente en el estado oculto correspondiente.

23 Probabilidad de una observaciónxT x1 xt xt-1 xt+1 1 t+1 T t t-1 Dados una secuencia de observaciones y un modelo, vimos que era fácil calcular la probabilidad de obtener dicha secuencia...simplemente multiplicamos las probabilidades individuales...entonces basta con calcular las probabilidades para los distintos caminos y elegir la máxima

24 HMM ¿Por qué la charla no termina acá?Para una dada secuencia de longitud T tenemos alrededor de 2T cálculos. Si N es el número de estados en el gráfico. Hay NT secuencias de estados posibles. Complejidad : O(2TNT ). Pero esto es fuerza bruta...podemos intentar algo mejor.

25 Problema 2: DecodificandoEncontrando la mejor secuencia de estados

26 Ejemplo de decodificaciónDada una secuencia de observaciones X,encuentre la secuencia de est.  . (1) Texto Shakespeare (s) vs mono (m) x = ..aefjkuhrgnandshefoundhappinesssdmcamoe…  = ..mmmmmmssssssssssssssssssssssssssssmmmmmm… (2) Dado fair (F) vs loaded (L) x = … …  = …LLLLLLLLLLLLFFFFFFFFFFFFFFFFLLLLLLLLLLLLLLLLLL… (3) DNA coding (C) vs non-coding (N) x = …AACCTTCCGCGCAATATAGGTAACCCCGG…  = …NNCCCCCCCCCCCCCCCCCNNNNNNNN…

27 Decodificando DADO x = x1x2……xN Queremos encontrar  = 1, ……, N,tal que P[ x,  ] esté maximizado * = argmax P[ x,  ] Podemos usar programación dinámica Sea Vk(i) = max{1,…,i-1} P[x1…xi-1, 1, …, i-1, xi, i = k] = Probabilidad de la secuencia de estados mas verosimil que termina en el estado i = k 1 2 K … x1 x2 x3 xK

28 El algoritmo de ViterbiInput: x = x1……xN Inicialización: V0(0) = 1 (0 es la posición inicial) Vk(0) = 0, para todo k > 0 Iteración: Vj(i) = ej(xi)  maxk akj Vk(i-1) ptrj(i) = argmaxk akj Vk(i-1) Terminación: P(x, *) = maxk Vk(N) Rastreo: N* = argmaxk Vk(N) i-1* = ptri (i)

29 El algoritmo de Viterbix1 x2 x3 ……. xi ……………………………..xN Estado 1 Vi(j) j K Es similar a “alinear” un conjunto de estados de una secuencia Complejidad temporal: O(K2N) Complejidad espacial: O(KN)

30 Volvamos al ejemplo del casinoJUSTO CARGADO 0.05 0.1 0.95 0.9 P(1|F) = 1/6 P(2|F) = 1/6 P(3|F) = 1/6 P(4|F) = 1/6 P(5|F) = 1/6 P(6|F) = 1/6 P(1|L) = 1/10 P(2|L) = 1/10 P(3|L) = 1/10 P(4|L) = 1/10 P(5|L) = 1/10 P(6|L) = 1/2

31 Ejemplo: El casino deshonestoDada la secuencia de dados: x = 3, 1, 5, 1, 1, 6, 4, 6, 4, 4, … V0(0) = 1 (estado Begin) Vj(1) = ej(xi)  maxk akj Vk(i-1) VF(1) = 1/6 * maxk [ ½ * 1, 0.95*0, 0.1*0] = 1/12 VL(1) = 1/10 * maxk [ ½ * 1, 0.9*0, 0.05*0] = 1/20 ptrF(1) = Estado Begin ptrL(1) = Estado Begin VF(2) = 1/6 * maxk [ 0.95 * 1/12, 0.1*1/20] = 0.079*1/6=0.013 VL(2) = 1/10 * maxk [0.90 * 1/20, 0.05*1/12] = 0.045*1/6=0.0045 ptrf(2) = Estado F ptrf(2) = Estado L

32 Ejemplo: El casino deshonestoDada la secuencia de dados: x = 3, 1, 5, 1, 1, 6, 4, 6, 4, 4, … Estado/obs. ……… xN Begin 1 .08 .013 .05 .005 F L

33 Testeo del algoritmo de ViterbiUna secuencia de 300 tiradas de dados justos (F) o cargados (L) ﾺ

34 El algoritmo de Viterbi – detalle prácticoLos errores de underflow pueden ser un problema P[ x1,…., xi, 1, …, i ] = a01 a12……ai e1(x1)……ei(xi) Números muy chicos! Solución: Tomar el logaritmo a todos los valores Vl(i) = log ek(xi) + maxk [ Vk(i-1) + log akl ]

35 Encontrando la verosimilitud de una secuencia generada por el modeloProblema 1: Evaluación Encontrando la verosimilitud de una secuencia generada por el modelo

36 El algoritmo Forward Queremos calcularP(x) = probabilidad de la secuencia x, dado el HMM Sumamos sobre todos las posibles formas de generar x: P(x) =  P(x, ) =  P(x | ) P() Para evitar sumar sobre un número de caminos exponencial , definimos fk(i) = P(x1…xi, i = k) (La probabilidad forward)

37 Derivación del algoritmo ForwardProbabilidad forward: fl(i) = P(x1…xi, i = l) = 1…i-1 P(x1…xi-1, 1,…, i-1, i = l) el(xi) = k 1…i-2 P(x1…xi-1, 1,…, i-2, i-1 = k) akl el(xi) = el(xi) k fk(i-1) akl Relación de recurrencia entre fl(i) y fk(i-1)

38 El algoritmo Forward Podemos calcular fk(i) para todo k, i, usando programación dinámica Inicialización: f0(0) = 1 fk(0) = 0, para todo k > 0 Iteración: fl(i) = el(xi) k fk(i-1) akl Terminación: P(x) = k fk(N) ak0 Donde, ak0 es la probabilidad de que el estado que termina sea k

39 Relación entre Forward y ViterbiInicialización: V0(0) = 1 Vk(0) = 0, para todo k > 0 Iteración: Vj(i) = ej(xi) maxk Vk(i-1) akj Terminación: P(x, *) = maxk Vk(N) FORWARD Inicialización: f0(0) = 1 fk(0) = 0, para todo k > 0 Iteración: fl(i) = el(xi) k fk(i-1) akl Terminación P(x) = k fk(N) ak0

40 Motivación para el algoritmo BackwardEl algoritmo de Viterbi encuentra el camino más probable a través del modelo pero eso no es lo único que nos interesa. Estamos interesados en calcular, por ejemplo, el estado más probable para cualquier observación x, o en general queremos saber, dada la secuencia observada x, la probabilidad de que el sistema esté en el estado k P(i = k | x), Empezamos por calcular P(i = k, x) = P(x1…xi, i = k, xi+1…xN) = P(x1…xi, i = k) P(xi+1…xN | x1…xi, i = k) = P(x1…xi, i = k) P(xi+1…xN | i = k)

41 Derivación del algoritmo BackwardDefinamos la probabilidad Backward: bk(i) = P(xi+1…xN | i = k) = i+1…N P(xi+1,xi+2, …, xN, i+1, …, N | i = k) = l i+1…N P(xi+1,xi+2, …, xN, i+1 = l, i+2, …, N | i = k) = l el(xi+1) akl i+1…N P(xi+2, …, xN, i+2, …, N | i+1 = l) = l el(xi+1) akl bl(i+1)

42 El algoritmo Backward Podemos calcular bk(i) para todo k, i, otra vez...usando programación dinámica Inicialización: bk(N) = ak0, para todo k Iteración: bk(i) = l el(xi+1) akl bl(i+1) Terminación: P(x) = l a0l el(x1) bl(1)

43 Complejidad computacionalComplejidad de los algoritmos Forward y Backward Tiempo: O(K2N) Espacio: O(KN) Técnica práctica para evitar errores de underflow Viterbi: suma de logaritmos Forward/Backward: rescaleo en cada posición multiplicando por una constante

44 Decodificación a posterioriPodemos ahora calcular fk(i) bk(i) P(i = k | x) = ––––––– P(x) Entonces, podemos preguntar ¿Cúal es el estado más probable en la posición i de la secuencia x?: Definamos ^ por decodificación a posteriori: ^i = argmaxk P(i = k | x)

45 Ejemplo: Decodificación a posterioriSecuencia de 300 tiradas al azar Areas grises: Dado cargado Areas blancas: Dado justo Línea llena: Probabilidad a posteriori de que el dado sea justo

46 Ejemplo: Decodificación a posterioriSecuencia de 1000 tiradas al azar, pero usando una probabilidad de transición justo-cargado de 0.01 En este caso, el camino más probable encontrado por el algoritmo de Viterbi nunca visita el estado “dado cargado”

47 Decodificación a posterioriPara cada estado, La decodificación a posteriori nos da una curva de similitud del estado para cada posición Algunas veces esto resulta más útil que el camino de Viterbi * Pero la decodificación a posteriori puede darnos caminos no validos

48 Rastreo de pesos máximosOtro approach es encontrar la secuencia de estados bajo alguna restricción, y maximizar la exactitud esperada de la asignación de estados Aj(i) = maxk tal que la condicion(k, j) Ak(i-1) + P(i = j | x)

49 Problema 3: AprendizajeReestimando los parámetros del modelo basado en los datos de entrenamiento

50 Estimación de parámetros – caso ITenemos un conjunto de secuencias de ejemplo del tipo de las que queremos que el modelo ajuste (secuencias de entrenamiento), que suponemos independientes. Si conocieramos el camino de estados que recorrió el modelo, los estados no están ocultos (el modelo oculto de Markov se transforma en una cadena de Markov), en la cual los estimadores de máxima verosimilitud para las frecuencias de emisión y transición se obtienen a partir de las frecuencias de observaciones. Si tenemos información (biológica o física) que nos aporte información previa a la distribución de probabilidades podemos agregarsela al modelo como pseudocuentas.

51 Estimación de parámetros – caso IIObjetivo: Dada una secuencia de observaciones, encuentre el modelo más probable que genere esa secuencia. Problema: No conocemos las frecuencias relativas de los estados ocultos visitados. No se conocen soluciones analíticas Nos acercamos a la solución por sucesivas aproximaciones. El problema es ahora de optimización, por lo que se pueden usar muchas heurísticas (simulated annealing, algoritmos genéticos, etc)

52 El algoritmo de Baum-WelchEste es el algoritmo de Expectation-Maximization (EM) para la estimación de parámetros. Aplicable a cualquier proceso estocástico, en teoría. Encuentra las frecuencias esperadas de los posibles valores de las variables ocultas. Calcula las distribuciones de maxima verosimilitud de las variables ocultas en base a las probabilidades forward y backward. Repite estos pasos hasta satisfacer algún criterio de “convergencia.”

53 Máxima probabilidad a posterioriSupongamos que hay una distribución de probabilidades P() de los parámetros. Entonces, por el teorema de Bayes, dada una observación x, la probabilidad a posteriori P(|x) = P(x|) P()/P(x) Como P(x) es independiente de , el mejor estimador está dado por la máxima probabilidad a posteriori MAP= argmax P(x|) P()

54 Máxima verosimilitud Denotemos con  los parámetros (probabilidades de transición y emisión) del HMM. Para una observación x, determinamos  usando el criterio de máxima verosimilitud ML = argmax P(x|) Si se usa para generar un conjunto de observables{xi}, entonces xi P(xi|) log P(xi|) Se maximiza para = de esta forma podemos encontrar ML por iteración (algoritmo de Baum-Welch).

55 Ejemplo de Baum-Welch Modelo real Modelo estimado (300 tiradas)

56 Algoritmo de Baum-Welch – ComentariosComplejidad temporal: # iteraciones  O(N2T) Garantiza incrementar la (log) verosimilitud del modelo No garantiza encontrar los mejores parámetros GLOBALES Converge a un óptimo LOCAL, dependiente de las condiciones iniciales Si hay demasiados parámetros/ el modelo es muy largo - Sobre- entrenamiento Algunas veces se reemplaza la reestimación de Baum-Welch por un entrenamiento basado sólo en el camino más probable (Viterbi training)

57 Aplicaciones biológicas de los HMMSA continuación… Aplicaciones biológicas de los HMMS