Documento Final

Este Documento final consta de la recopilación de los capitulos: \[ * Preliminaries: \quad Mixtures\quad and \quad Markov \quad Chains;\\ * Hidden \quad Markov \quad Models:\quad definition \quad and \quad properties;\\ * Estimation \quad by \quad direct \quad maximization \quad of \quad the \quad likelihood; \]

Capitulo 1: Preliminaries: Mixtures and Markov Chains

1.1Introducción

Los modelos ocultos de Markov (HMM) son modelos en los cuales se tiene la atribución de que el sistema que se está trabajando es un proceso de Markov, cuyos parametros no se conocen. Esto permite que el proceso se vuelva de proposito general y sirva para series de tiempo: univariadas y multivariadas.

Para ilustrar al lector sobre las aplicaciones que tienen las HMM, se usará un archivo llamado “EarthQuakeDataCH1.csv”; este archivo contiene informacion de terremotos de magnitud mayor a 7 desde 1900 al 2007.

Para ello se procede a leer el CSV con la librería “readr”:

library(readr)
library(htmltools)
data_terremotos <- data.frame(read_csv('EarthQuakeDataCH1.csv'))

Parsed with column specification:
cols(
  Anio = col_integer(),
  cantidadTerremotos = col_integer()
)

#graficacion de datos
plot(x = data_terremotos[,1],y = data_terremotos[,2], type = "b", col = "#33d6ff", lwd = 4, xlab = "Años", ylab = "Magnitud")

Figura 1.1

Para darnos una idea sobre que distrubucion puede seguir la data, utilizamos un histograma que se muestra a continuación:

#se utiliza la columna 2 debido a que alli se encuentran las magnitudes
magnitudes <- data_terremotos[,2]
#creacion del histograma
h <- hist(magnitudes, breaks = 10, density = 10,
          col = "#0099cc", xlab = "Accuracy", main = "Terremotos - Distribucion") 
#definicion de limites X y Y
xfit <- seq(min(magnitudes), max(magnitudes), length = 40) 
yfit <- dnorm(xfit, mean = mean(magnitudes), sd = sd(magnitudes)) 
yfit <- yfit * diff(h$mids[1:2]) * length(magnitudes) 
lines(xfit, yfit, col = "#ff4775", lwd = 2)

Tal como se puede observar la dispersión corresponde a un podelo de Poisson. Las cadenas ocultas de Markov(HMM) se han usado desde hace mas de tres decadas, el mayor campo de operación ha sido el analisis de señales, tales como:

• Reconocimiento: voz, facial, gestos, firmas.
• Ambiente: terremotos, direccion del viento, lluvia.
• Finanzas
• Biofisica: modelacion de canales de iones.

Una de las grandes caracteristicas del HMM es que es simple se realizar y debido a su trasfondo matematico tambien es adecuado para tener una respuesta computacional rapida.

1.2 Modelos de mezcla independientes

1.2.1 Definicion y propiedades

Considere nuevamente la serie de cuentas de terremotos ilustrada en la Figura 1.1. Un modelo estándar para cuentas ilimitadas es la distribución de Poisson, con su función de probabilidad \[p(x)=e^{-λ}*λ^x/x!\] y la propiedad que indica que la varianza equivale a la media. Sin embargo, para la serie de terremotos, la muestra de varianza, \(s^2≈52\), es mucho mas grande que la muestra de la media, \[\bar{x} ≈19\], lo cual indica una fuerte sobre dispersión relativa a la distribución de Poisson. La falta de semejanza es confirmada por la Figura 1.2, la cual ilustra la distribución de Poisson ajustada y un diagrama de barras de las frecuencias relativas de las cuentas. Un método para resolver observaciones sobre dispersas con una distribución binodo o (más generalmente) distribución multimodo es usar un modelo de mezcla. Los modelos de mezcla son diseñados para acomodar heterogeneidad inobservada en la población; eso es, la población puede consistir de grupos inobservados, cada uno teniendo una distribución distinta para la variable observada. Considérese el ejemplo, la distribución del número, X, de paquetes de cigarrillos comprados por los consumidores de un supermercado. Los consumidores pueden ser divididos en grupos, por ejemplo, no fumadores, fumadores ocasionales, y fumadores regulares. Ahora aun si el numero de paquetes comprado por los consumidores entre cada grupo fuera distribuible por Poisson, la distribución X no seria Poisson; seria sobre dispersa relativa a Poisson, y talvez aun multimodo. Análogamente, supóngase que cada cuenta en la serie de terremotos es generada por una de dos distribuciones de Poisson, con medias \(λ_1\) y \(λ_2\), donde la elección de media es determinada por otros mecanismos aleatorios que llamaremos el proceso de parámetro. Supóngase también que \(λ_1\) es seleccionado con probabilidad \(δ_1\) y \(λ_2\) con probabilidad \(δ_2=1-δ_1\). Veremos mas adelante en este capitulo que la varianza de la distribución resultante excede la media por \(δ_1 δ_2 (λ_1-λ_2 )^2\). Si el proceso de parámetro es una serie de variables aleatorias independientes, las cuentas también son independientes, de ahí el término “mezcla independiente”. En general, una distribución de mezcla independiente consiste de un numero finito, dígase \(m\), de distribuciones de componente y una “distribución de mezcla” la cual selecciona de estas componentes. Las distribuciones de componente pueden ser tanto discretas como continuas. En el caso de dos componentes, la mezcla de distribución depende de dos probabilidades o funciones de densidad:

\[ Componente \qquad \qquad\qquad\qquad\qquad1 \qquad\qquad 2 \\ Probabilidad o función de densidad \qquad p_1(x) \qquad p_2(x) \]

Para especificar la componente, una necesita una variable aleatoria discreta C la cual realiza la mezcla:

\[ C= \begin{cases} 1 \quad con \quad probabilidad \quad\delta_1\\ 2 \quad con \quad probabilidad \quad\delta_2 = 1-\delta_1 \end{cases} \]

La estructura de ese proceso para el caso de dos distribuciones de componentes continuas es ilustrada en la Figura 1.3. En ese ejemplo uno puede pensar de C como el resultado de lanzar una moneda con probabilidad 0.75 de “cara”: si el resultado es “cara”, entonces C = 1 y una observación es tomada de p_1: si es “cruz”, entonces C = 2 y una observación es tomada de \(p_2\). Suponemos que no sabemos el valor de C, eso es, cual de \(p_1\) o \(p_2\) estaba activo cuando la observación fue generada.

Figura 1.3

La extensión de m componentes es sencilla. Que \(δ_1,…,δ_m\) denote las probabilidades asignadas a las diferentes componentes, y que \(p_1,…,p_m\) denote sus probabilidades o funciones de densidad. Que X denote la variable aleatoria la cual tiene la distribución de mezcla. En el caso discreto, la función de probabilidad de X es dada por:

\[ p(x) = \sum_{i=1}^{m} Pr(X=x|C=i)Pr(C=i)\\ =\sum_{i=1}^{m}\delta_ip_i(x) \] El caso continuo es análogo. La expectativa de la mezcla puede ser dada en términos de las expectaciones de las distribuciones de componentes. Dejando que Y, denote la variable aleatoria con función de probabilidad \(p_i\), tenemos:

\[ E(X) = \sum_{i=1}^{m} Pr(C=i) E(X|C=i) = \sum_{i=1}^{m} \delta_iE(Y_i) \]

El mismo resultado aguarda para una mezcla de distribuciones continuas. Mas generalmente, para una mezcla del k-esimo momento del origen es simplemente una combinación lineal de los k-esimos momentos de sus componentes \(Y_i\):

\[ E(X^k) = \sum_{i=1}^{m} \delta_iE(Y_i^k), k=1,2,... \]

Nótese que el resultado análogo no cuenta para momentos centrales. En particular, la varianza X no es una combinación linear de las varianzas de sus componentes \(Yi\). El Ejercicio 1 le pide al lector probar eso, en el caso de dos componentes, la varianza de la mezcla está dada por:

\[ Var(X) = \delta_1Var(Y_1)+\delta_2Var(Y_2)+\delta_1\delta_2(E(Y_1)-E(Y_2))^2 \]

1.2.2 Estimación de parámetros

La estimación de parámetros de una distribución de mezclas es a menudo realizada por una probabilidad máxima (ML). La probabilidad de un modelo de mezcla con m componentes esta dada, para tanto casos discretos como continuos, por:

\[ L(\theta_1,...,\theta_m,\delta_1,...,\delta_m|x_1,...,x_n) = \prod_{j=1}^{n}\sum_{i=1}^{m}\delta_ip_i(x_j,\theta_i) \]

Aquí \(θ_1,…,θ_m\) son los vectores de parametros de las distribuciones de componentes, \(δ_1,…,δ_m\) son los parametros de mezcla, totalizando 1, y \(x_1,…,x_n\) son las n observaciones. Entonces, en el caso de distribuciones de componentes cada una especificada por un parámetro, \(2m-1\) parametros independientes tienen que ser estimados. Excepto talvez en casos especiales, maximización analítica de tal probabilidad no es posible, pero es generalmente sencillo de evaluarlo rápidamente; ver Ejercicio 3. Maximización numérica será ilustrada aquí al considerar el caso de una mezcla de distribuciones de Poisson. Supongase que m=2 y las dos componentes son distribuibles por Poisson con medias \(λ_1\) y \(λ_2\). Que \(δ_1\) y \(δ_2\) sean los parámetros de mezcla (con \(δ_1 + δ_2 = 1\). La distribución de mezcla p esta dada por:

\[ p(x) = \delta_1 \frac{\lambda_1^xe^{-\lambda_1}}{x!}+\delta_2 \frac{\lambda_2^xe^{-\lambda_2}}{x!} \] Desde que \(δ_2=1-δ_1\), hay solamente tres parámetros a estimar: \(λ_1\), \(λ_2\) y \(δ_1\). La probabilidad es:

\[ L(\lambda_1,\lambda_2,\delta_1|x_1,...,x_n) = \prod_{i=1}^n(\delta_1\frac{\lambda_1^{x_i}e^{-\lambda_1}}{x_i!}+(1-\delta_1)\frac{\lambda_2^{x_i}e^{-\lambda_2}}{x_i!}) \] La maximización analítica de L con respecto a \(λ_1\), \(λ_2\) y \(δ_1\) seria rara, ya que L es el producto de n factores, del cual cada uno es una suma. Primero tomar el logaritmo y luego la diferenciación no simplifica mucho las cosas tampoco. Por lo tanto, la estimación de parámetros es mas convenientemente llevada a cabo por maximización numérica directa de la probabilidad (o su logaritmo), aunque el logaritmo EM es comúnmente usado como alternativa; véase, por ejemplo, McLachlan y Peel (2000) o Fruhwirth-Schnatter (2006). Un paquete útil de R para estimación en modelos de mezcla es \(flexmix\) (Leisch, 2004). Sin embargo, es sencillo escribir nuestro propio código de R para evaluar, y luego maximizar, probabilidades de mezcla en casos sencillos. Esta probabilidad puede ser luego maximizada al usar (por ejemplo) la función de R nlm. Sin embargo, los parámetros \(\delta\) y \(\lambda\) son sujetos por (ecuación rara) y (para \(i = 1, …, m\)) \(δ_i > 0\) y \(λ_i>0\). Es entonces necesario Re parametrizar si uno así lo desea para utilizar un optimizador sin restricciones tal como nlm. Una posibilidad para maximizar la probabilidad con respecto a los 2m-1 “parámetros funcionales” sin restricción.

\[ \eta_i = log\lambda_i (i=1,...m) \] y

\[ \tau_i = log(\frac{\delta_i}{1-\sum_{j=2}^m \delta_j}) (i=2,...,m) \]

Se recuperan los “parámetros naturales” originales por medio de:

\[ \lambda_i = e^{\eta_i} \quad (i=1,...,m) \\ \delta_i = \frac{e^{\tau_i}}{1+\sum_{j=2}^me^{\eta_j}} \quad (i=2,...,m) \]

Y \(δ_1=1-\sum_{j=2}^{m}δ_i\) . El siguiente código implementa las ideas anteriores en orden para que encajen un modelo de cuatro distribuciones de Poisson a las cuentas de terremotos. Los resultados están dados por m = 1,2,3,4 en la Tabla 1.2

mllk <- function(wpar,x){ zzz <- w2n(wpar)
        -sum(log(outer(x,zzz$lambda,dpois)%*%zzz$delta)) }
n2w  <- function(lambda,delta)log(c(lambda,delta[-1]/(1-sum(delta[-1]))))
w2n  <- function(wpar){m <- (length(wpar)+1)/2
        lambda <- exp(wpar[1:m])
        delta  <- exp(c(0,wpar[(m+1):(2*m-1)]))
return(list(lambda=lambda,delta=delta/sum(delta))) }
x        <- read.table("/Users/bracruz/Documents/8vo. Trimestre IO/Sist. Expertos y Redes Prob./earthquakes.txt")[,2]
wpar     <- n2w(c(10,20,25,30),c(1,1,1,1)/4)
w2n(nlm(mllk,wpar,x)$estimate)

NA/Inf replaced by maximum positive value

$lambda
[1] 15.52750 10.58416 20.96911 32.07922

$delta
[1] 0.35429030 0.09302729 0.43662984 0.11605257

Notese como, en este código, el uso de la función outer hace posible el evaluar la probabilidad de mezcla de Poisson en una simple y compacta expresión en vez de en un ciclo. Pero si la distribución siendo mezclada fueran distribuciones con mas de un parámetro (dígase, normal), un enfoque un poco distinto seria necesario.

Los resultados anteriores pueden ser verificados en la siguiente tabla:

Tabla 1.2

1.2.3 Probabilidad ilimitada en mezclas

Hay un aspecto en mezclas de distribuciones continuas que difieren del caso discreto y vale la pena resaltarlo. Es el siguiente: puede pasar que: en la vecinidad de ciertas combinaciones de parámetros, la probabilidad es ilimitada. Por ejemplo, en el caso de una mezcla de distribuciones normales, la probabilidad se vuelve arbitrariamente grande si uno asigna la media componente igual a una de las observaciones y permite que la varianza correspondiente tienda a cero. El problema ha sido ampliamente discutido en la literatura de modelos de mezcla, y hay aquellos que sugieren que si la probabilidad es entonces ilimitada, la ML estima simplemente que “no existe”.

La fuente del problema, sin embargo, es solamente el uso de densidades en vez de probabilidades en la posibilidad; no surgiría si se remplazara cada valor de densidad en una posibilidad por la probabilidad del intervalo correspondiente al valor registrado. (Por ejemplo, una observación registrada como “12.4” esta asociada con el intervalo [12.35, 12.45]). En el contexto de mezclas independientes se remplaza la expresión:

\[ \prod_{j=1}^n\sum_{i=1}^m\delta_i p_i(x_j,\theta_i) \]

de la posibilidad (véase ecuación anterior) por la posibilidad dicreta:

\[ L = \prod_{j=1}^n\sum_{i=1}^m\delta_i \int_{a_j}^{b_j} p_i(x,\theta_i)dx \]

Donde el intervalo \((a_j, b_j)\) consiste de aquellos valores los cuales, si son observados, serian registrados como \(x_j\). Esto simplemente equivale a reconocer explícitamente la naturaleza del intervalo de todas las supuestas observaciones continuas. Mas generalmente, la posibilidad discreta de observaciones en un conjunto de variables aleatorias \(X_1, X_2,…, X_n\) es la probabilidad de la forma \(Pr(a_i < X_i < b_i \quad para\quad toda \quad i)\). Utilizaremos el termino posibilidad continua para la densidad conjunta evaluada en la observación. Otra manera de evitar este problema es el imponer un limite inferior en las varianzas y buscar el mejor sujeto máximo local para ese límite. Puede suceder, sin embargo, que se puede ser tan afortunado para evitar los “picos” de posibilidad al buscar por un máximo local; en este aspecto, valores buenos de inicio pueden ayudar. El fenómeno de posibilidades ilimitadas no surge para observaciones de valores discretos porque la posibilidad en ese caso es una probabilidad y por ende, limitada a valores entre 0 y 1.

1.2.4. Ejemplos de modelos de mezcla ajustados

Mezcla de distribuciones de Poisson

Si uno utiliza nlm para ajustar una mezcla de m distribuciones de Poisson (m=1,2,3,4) a los datos del terremoto, uno obtiene los resultados desplegados en la Tabla 1.2. Notese que hay una clara mejora en la posibilidad resultante de la adicion de una segunda componente, y muy poca mejora de la adicion de una cuarta – aparentemente insuficiente para justificar los dos parámetros adicionales. La Figura 1.4 presenta un histograma de las cuentas observadas y los cuatro modelos ajustados. Es claro que las mezclas encajan las observaciones mucho mejor a como lo hace una sola distribución de Poisson, y visualmente los modelos de tres y cuatro estados parecen adecuados. El mejor ajuste para las mezclas es también evidente de las varianzas de los cuatro modelos como esta presentada en la Tabla 1.2. Al computar las medias y varianzas de los modelos hemos usado \(E(X)=\sum_iδ_i λ_i\) y \(Var(X)=E(X^2 )-(E(X))^2\), con \(E(X^2 )=\sum_iδ_i (λ_i+λ_i^2)\). Para comparación también hemos usado el paquete de R flexmix para ajustar los cuatro modelos. Los resultados correspondieron mucho a excepción de laso del modelo de cuatro componentes, donde el valor de posibilidad mas grande que encontramos por flexmix fue 356.7759 y las medias de los componentes difirieron un poco. Notese también, que la discusión superior ignora la posibilidad de dependencia serial en los datos del terremoto, un punto el cual se retomara en el Capítulo 2.

Figura 1.4

1.3 Cadenas de Markov

Ahora introduciremos las cadenas de Markov. Nuestro análisis se restringe a aquellos aspectos de cadenas de Markov discretas los cuales necesitaremos. Asi, sin embargo haremos referencia a propiedades tales como la irreductibilidad y aperiodicidad, no lidiaremos con cuestiones técnicas. 1.3.1 Definiciones y ejemplos Una secuencia de variables aleatorias discretas \({C_t:t ∈ N} \)se dice que es una (tiempo discreto) Cadena de Markov (CM) si, para todo t ∈N, satisface la propiedad de Markov:

\[ Pr(C_{t+1} | C_t,...,C_1) = Pr((C_{t+1}|C_t) \] Eso es, condicionando a la historia del proceso hasta el tiempo t es equivalente a condicionar solo el valor mas reciente \(C_t\). Para hacerlo más compacto, definimos \(C^t\) como la historia \((C_1,C_2,…,C_t )\), en cual caso, la propiedad de Markov se puede escribir como:

\[ Pr(C_{t+1} | C^{(t)} = Pr(C_{t+1} | C_t) \] La propiedad de Markov puede ser interpretada como la primera relajación de la suposicion de independencia. Las variables aleatorias \({C_t}\) son dependientes de una manera que es matemáticamente conveniente, como es mostrado en el siguiente grafico en el cual el pasado y el futuro son dependientes solo a través del presente.

Cantidades importantes asociadas con una cadena de Markov son las probabilidades conditionales, también llamadas probabilidades de transición:

\[ Pr(C_{s+t}=j | C_s = i) \] Si estas probabilidades no dependen de s, la cadena de Markov es llamada homogénea, de otra manera seria no-homogenea. A no ser que haya una indicación explicita al contrario, asumiremos que la cadena de Markov en discusión es homogénea, en cual caso las probabilidades de transición serán denotadas por:

\[ \gamma_{ij}(t)= Pr(C_{s+t}|C_s=i) \] Notese que la notación \(γ_{ij} (t)\) no involucra a s. La matriz \(Γ(t)\) esta definida como la matriz con elementos \(γ_{ij} (t)\).

Una importante propiedad de todas las cadenas de Markov homogéneas y de estado y espacio finito es la que satisface las ecuaciones Chapman-Kolmogorov:

\[ \Gamma(t+u) = \Gamma(t) \Gamma(u) \] La solución requiere solamente la definición de probabilidad condicional y la aplicación de la propiedad de Markov. Las ecuaciónes Chapman-Kolmogorov implican que, para todo \(t ∈N\):

\[\Gamma(t)=\Gamma(1)^t\] Esto es, la matriz de probabilidades de transición con paso-t es la t-esima potencia de \(Γ(1)\), la matriz de probabilidades de transición de un paso. La matriz \(Γ(1)\), la cual será abreviada como \(Γ\), es una matriz cuadrada de probabilidades con la suma de sus filas igual a 1:

\[ \quad \] Donde m denota el numero de estados de la cadena de Markov. La declaración de que la suma de las filas es igual a 1 puede ser escrito como \(Γ1'=1'\); eso es, el vector columna \(1’\) es un eigenvector derecho de \(Γ\) y corresponde al eigenvalor 1. Nos referiremos a \(Γ\) como la matriz de transición de probabilidad (t.p.m). Muchos autores utilizan en vez del termino “matriz de transición”; evitaremos este termino debido a posibles confusiones con cuentas de matriz de transición, o intensidades de una matriz de transición.

Las probabilidades incondicionales \(Pra(C_t=j)\) de una cadena de Markov siendo en un estado dado en un tiempo dado t son a menudo de interés. Denotamos estos por el vector fila:

\[u(t) = (Pr(C_t=1),...,Pr(C_t=m)), t∈N\]

Nos referimos a u(1) como la distribución inicial de una cadena de Markov Para deducir la distribución en el tiempo t+1 desde el cual en t post multiplicamos por la matriz de transición de probabilidad \(Γ\):

\[u(t+1)=u(t)\Gamma \] Ejemplo. Imaginese que la declaración de días soleados y lluviosos es tal que el clima del dia depende solamente en el del dia anterior, y las probabilidades de transición están dadas por la siguiente tabla:

\[ \quad \]

Esto es, si hoy fue un dia lluvioso, la probabilidad que mañana sea lluvioso es de 0.9; si hoy estuvo soleado, la probabilidad de es 0.6. El clima es entonces una cadena de Markov homogénea de dos estados, con t.p.m. \(Γ\) dado por

\[ \Gamma = 0.9 \quad 0.1\\ \qquad0.6 \quad 0.4 \]

Ahora suponga que hoy (tiempo 1) es un dia soleado. Esto significa que la distribución del clima del dia de hoy es:

\[u(1) = (Pr(C_1 = 1),Pr(C_1 = 2)) = (0,1)\]

La distribución del clima de mañana, del pasado mañana, y asi, puede ser calculada post multiplicando repetidamente \(u(1)\) por
\(Γ\), el t.p.m:

\[ u(2) = (Pr(C_2 = 1),Pr(C_2 = 2)) = u(1)\Gamma = (0.6,0.4)\\ u(3) = (Pr(C_3 = 1),Pr(C_3 = 2)) = u(2)\Gamma = (0.78,0.22) \]

1.3.2. Distribuciones estacionarias

Una cadena de Markov con matriz de transición de probabilidad Γ se dice que tiene una distribución estacionaria δ (un vector fila con elementos no negativos) si δΓ=δ y δ1’ = 1. El primero de estos requerimientos expresa la estacionariedad, la segunda es el requerimiento que δ es una distribución de probabilidad. Por ejemplo, la cadena de Markov con t.p.m dado por

\[ \Gamma = 1/3 | \quad 1/3| \quad 1/3\\ \quad 2/3| \quad 0| \quad 1/3 \\ \quad 1/2| \quad 1/2| \quad 0 \] Tiene una distribución estacionaria \(δ=1/32 (15,9,8)\).

Desde que \(u(t+1) = u(t)Γ\), una cadena de Markov iniciada desde su distribución estacionaria continuara teniendo esa distribución en todos lo puntos de tiempo subsecuente, y nos referiremos a tal proceso como una cadena de Markov estacionaria. Es talvez valioso resaltar que esto asume mas que solamente homogeneidad. Solo la homogeneidad no seria suficiente para volver la cadena de Markov un proceso estacionario, y preferimos reservar el adjetivo “estacionario” para cadenas de Markov estacionarias que tienen la propiedad adicional que la distribución inicial \(u(1)\) es la distribución estacionaria y por ende procesos estacionarios.

Una cadena de Markov irreducible (homogénea, tiempo discreto, estado-espacio finito) tiene una distribución estacionaria, única y estrictamente positiva. Nótese sin embargo la suposición de irreductibilidad es necesaria para esta conclusión, la aperiodicidad no. Si, sin embargo, uno añade la suposición de aperiodicidad, una única distribución limitante existe, y es precisamente la distribución estacionaria. Desde que siempre debemos asumir aperiodicidad e irreductibilidad, los términos “distribución limitante” y “distribución estacionaria” son sinónimos para nuestros propósitos.

Un resultado general que puede ser convenientemente usado para computar una distribución estacionaria. El vector \(δ\) con elementos no negativos es una distribución estacionaria de la cadena de Markov con t.p.m Γ si y solo si:

\[\delta(I_m - \Gamma + U) = 1\]

Donde 1 es un vector fila de unos, \(I_m\) es la matriz identidad m x m, y U es la matriz m x m de unos. Alternativamente, una distribución estacionaria puede ser encontrada al borrar una de las ecuaciones en el sistema \(δΓ=δ\) y reemplazándola por \(\sum_iδ_i=1\).

1.3.3 Función de autocorrelación

Tendremos una oportunidad, por ejemplo en la Seccion 2.2.3 y en el Ejercicio 4(f) en el Capitulo 2, para comparar las funciones de autocorrelación (ACF) de uno modelo oculto de Markov con su cadena de Markov subyacente \({C_t}\), en los estados \(1, 2, …, m\). Asumiremos que estos estados son cuantitativos y no solamente categóricos. La ACF de \({C_t}\), asumida estacionaria e irreducible, puede ser obtenida de la siguiente manera. Primero, definiendo \(v = (1, 2, …, m)\) y \(V = diag(1, 2, …, m)\), tenemos, para todos los enteros k no negativos, \[ Cov(C_t,C_{t+k}) = \delta V\Gamma^kv' - (\delta v')^2 \] Segundo, si Γ es diagonizable, y sus eigenvalores (otro de 1) son denotados por w_2, w_3, …, w_m, entonces Γ puede ser escrito como:

\[\Gamma = U \Omega U^{-1} \] Donde \(Ω\) es \(diag(1, w_2, w_3, …, w_m)\) y las columnas de U son eigenvectores derechos correspondientes de \(Γ\). Entonces tenemos, para enteros no negativos k,

\[ Cov(C_t,C_{t+k}) = \delta V U \Omega^k U^{-1} v' - (\delta v')^2 \\ = a\Omega^k b' - a_1b_1 \\ =\sum_{i=2}^m a_ib_iw_i^k \] donde a = \(\delta VU\) y b’ = \(U^{-1}\). Por lo tanto \(Var(C_t )=\sum_{i=2}^m a_i b_i\) y, para enteros no negativos k:

\[ \rho(k) = Corr(C_t,C_{t+k}) = \sum_{i=2}^m a_ib_iw_i^k / \sum_{i=2}^m a_ib_i \]

Esto es un peso promedio de las k-esimas potencias de los eigenvalores \(w_2, w_3, …, w_m\), y de alguna manera similares a la ACF de un proceso auto regresivo Gaussiano de orden m-1. Nótese que la ecuación anterior implica en el caso m=2 que \(ρ(k)=ρ(1)^k\) para todos los enteros no negativos k, y que \(ρ(1)\) es el eigenvalor distinto a 1 de \(Γ\).

1.3.4. Estimando probabilidades de transición

Si nos dan una realización de una cadena de Markov, y deseamos estimar las probabilidades de transición, un enfoque – pero no solo el único – es el encontrar las cuentas de transición y estimar las probabilidades de transición como frecuencias relativas. Por ejemplo, si la CM tiene tres estados y la secuencia observada es:

( 2332111112 3132332122 3232332222 3132332212 3232132232 3132332223 3232331232 3232331222 3232132123 3132332121 )

Entonces la matriz de cuentas de transición es:

\[ (f_{ij}) = 4 \quad 7 \quad 6 \\ \qquad\qquad8 \quad 10 \quad 24 \\ \qquad\qquad 6 \quad 24 \quad 10 \]

Donde \(f_ij\) denota el numero de transiciones observadas del estado i al estado j. Desde que el numero de transiciones del estado 2 al estado 3 es 24, y el numero total de transiciones del estado 2 es 8+10+24, un estimado de frecuencia relativa de \(γ_23\) es 24/42. El t.p.m \(Γ\) es entonces estimado por:

\[ (f_{ij}) = 4/17 \quad 7/17 \quad 6/17 \\ \qquad\qquad8/42 \quad 10/42 \quad 24/42 \\ \qquad\qquad 6/40 \quad 24/40 \quad 10/40 \]

Ahora mostraremos que esto es de hecho el estimado condicional ML de \(Γ\), condicionada en la primera observación. Supongamos entonces que deseamos estimar los \(m^2-m\) parámetros \(γ_{ij} (i≠j)\) de una cadena e Markov de estado m \({C_t}\) de una realización c_1, c_2, …, c_T. La posibilidad condicionada en la primera observación es:

\[ L = \prod_{i=1}^m\prod_{j=1}^m \gamma_{ij}^{f_{ij}} \] El registro de probabilidad es entonces:

\[ l = \sum_{i=1}^m(\sum_{j=1}^m f_{ij} log \gamma_{ij}) = \sum_{i=1}^m l_i() \]

Y podemos maximizar l al maximizar cada \(l_i\) por separado. Sustituyendo \(1-\sum_{(k≠i)}γ_{ik}\) para \(γ_ii\) diferenciando \(l_i\) con respecto a una posibilidad de transición fuera de la diagonal \(γ_ij\) e igualando la derivada a cero rendimientos

\[ 0 = \frac{-f_{ii}}{1-\sum_{k≠i}\gamma_{ik} +\frac{f_{ij}}{\gamma_{ij}} = -\frac{f_{ii}}{\gamma_{ii}}+\frac{f_{ij}}{\gamma_{ij}}} \]

Entonces, a menos que un denominador sea cero en la ecuación superior,( f_{ij} γ_{ii}=f_{ii} γ_{ij} )y entonces \(γ_{ii} \sum_{(j=1)}^mf_{ij}=f_{ii}\). Esto implica que, en un máximo de la posibilidad,

\[ \gamma_{ii} = f_{ii}/\sum_{j=1}^m f_{ij} \quad y \quad \gamma_{ij} = f_{ij}/\sum_{j=1}^m f_{ij} \]

(Podríamos utilizar multiplicadores de Lagrange para expresar las restricciones \(\sum_{j=1}^mγ_{ij}=1\) sujeto al cual buscamos maximizar los términos \(l_i\) y por ende la posibilidad; es solo la probabilidad de transición empírica – es entonces visto a ser un estimador condicional ML de \(γ_ij\). El estimador \(Γ\) satisface el requerimiento que la suma de las filas es igual a 1.

La suposicion de estacionariedad de la cadena de Markov no fue usada en la derivación superior. Si deseamos asumir la estacionariedad, podemos usar la posibilidad incondicional. Esta es la posibilidad condicional mostrada anteriormente, multiplicada por la probabilidad estacionaria δ_c1. La posibilidad incondicional o su logaritmo pueden ser entonces maximizados numéricamente, sujeto a restricciones no negativas y suma de filas, en orden para estimar las probabilidades de transición \(γ_{ij}\). Bisgaard y Travis (1991) muestran en el caso de una cadena de Markov de dos estados que, excepto en algunos casos extremos, ecuaciones de posibilidad incondicional tienen una única solución. Para algunos casos especiales no triviales de la cadena de dos estados, también derivan expresiones explicitas para los estimados incondicionales de posibilidad máxima (MLEs) de las probabilidades de transición. Desde que usamos tal resultado luego, lo declaramos aquí. Suponga que la cadena de Markov \({C_t}\) toma los valores 0 y 1, y que deseamos estimar las probabilidades de transición \(γ_{ij}\) de una secuencia de observaciones en el cual hay \(f_{ij}\)transiciones del estado i al estado j (i,j=0,1), y \(f_{11}>0\) pero \(f_{00}=0\). Asi que en las observaciones un cero siempre es seguido por un uno. Definase \(c=f_{10}+(1-c_1)\) y \(d=f_{11}\). Entonces las MLEs incondicionales de las probabilidades de transición están dadas por:

\[ \gamma_{01} = 1 \quad y \quad \gamma_{10} = \frac{-(1+d)+((1+d))^2 + 4c(c+d-1))^{1/2}}{2(c+d+1)} \]

1.3.5. Cadenas de Markov de orden superior

En casos donde observaciones en un proceso con espacio y estado finito aparecen no satisfacer la propiedad de Markov, una posibilidad que sugiere a si misma es el utilizar una cadena de Markov de orden superior, esto es, un modelo \({C_t}\) satisfaciendo la siguiente generalización de la propiedad de Markov para algunos l≥2:

\[ Pr(C_t | C_{t-1} = j_1,...,C_{t-l}=j_l) = \sum_{i=1}^l \lambda_i q(j_i,j_o) \] Un registro de tal orden de cadenas de Markov puede ser encontrado, por ejemplo, el Lloyd (1980, Sección 19.9). Aunque tal modelo no es en un sentido usual una cadena de Markov, podemos redefinir el modelo de tal manera como para producir un proceso equivalente. Si permitimos que \(Y_t=(C_{(t-l+1)},C_{(t-l+2)},…,C_t)\), entonces {Y_t} es una cadena de Markov de primer orden en M^l, donde M es el espacio de estado de {C_t}. Aunque algunas propiedades pueden ser raras de establecer, ninguna nueva teoría es involucrada en analizar una cadena de Markov de orden superior mas que una de primer orden. Una cadena de Markov de segundo orden, si es estacionaria, esta caracterizada por las probabilidades de transición:

\[ \gamma(i,j,k) = Pr(C_t=k | C_{t-1} = j, C_{t-2} = i) \]

Y tiene una distribución estacionaria de variable doble \(u(j,k)=Pra(C_{t-1}=j,C_t=k)\) que satisface

\[ u(j,k) = \sum_{i=1}^m u(i,j)\gamma(i,j,k) \quad y \quad \sum_{j=1}^m\sum_{k=1}^m u(j,k) =1 \]

Por ejemplo, la cadena de Markov estacionaria de segundo orden mas general \({C_t}\), en los dos estados 1 y 2, esta caracterizada por las siguientes cuatro probabilidades de transición:

\[ a = Pr(C_t=2 | C_{t-1}=1, C_{t-2} = 1) \\ b = Pr(C_t=1 | C_{t-1}=2, C_{t-2} = 2) \\ c = Pr(C_t=1 | C_{t-1}=2, C_{t-2} = 1) \\ d = Pr(C_t=2 | C_{t-1}=1, C_{t-2} = 2) \\ \] El proceso \({Y_t }={(C_{t-1},C_t )}\) es entonces una cadena de Markov de primer orden, con los cuatro estados (1,1), (1,2), (2,1), (2,2), con matriz de transición de probabilidades

\[ \left(\begin{array}{cc} 1-a & a & 0 & 0\\ 0 & 0 & c & 1-c \\ 1-d & d & 0 & 0 \\ 0 & 0 & b & 1-b \end{array}\right) \]

Note los ceros estructurales que aparecen en la matriz. No es posible, por ejemplo, hacer una transición directa de (2,1) a (2,2); por lo tanto, el cero en la fila 3 y columna 4 en el t.p.m. Los parámetros a, b, c y d están limitados a 0 y 1 pero son sin restricciones en otro caso. La distribución estacionaria de \({Y_t}\) es proporcional al vector

\[(b(1-d),ab,ab,a(1-c))\]

Del cual se dice que la matriz \(u(j,k)\) de probabilidades bi-variable estacionaria para \({C_t }\) es:

\[ \frac{1}{b(1-d)+2ab+a(1-c)} \left(\begin{array}{cc} b(1-d) & ab\\ ab & a(1-c) \end{array}\right) \]

El uso de una cadena de Markov general de orden superior incrementa el número de parámetros del modelo; una cadena de Markov general de orden l en m estados tiene \(m^l (m-1)\) probabilidades de transicion independientes. Pegram (1980) y Raftery (1985) han por lo tanto propuesto ciertas clases de modelos para cadenas de orden superior. El modelo de Pegram tiene m+l-1 parametros, y aquellos de Raftery m(m-1)+l-1. Para m=2 los modelos son equivalentes, pero para m>2 aquellos de Raftery son mas generales y pueden representar un rango mas amplio de patrones de dependencia y estructuras de autocorrelación. En ambos casos un incremento de uno en el orden de la cadena de Markov requiere solamente de un parámetro adicional. Modelos de Raftery, a los cuales el se refiere como modelos de “distribución de transición de la mezcla” (MTD), son definidos como sigue. El proceso \({C_t}\) toma valores \(M={1,2,…,m}\) y satisface

\[Pr(C_t = j_0 | C_{t-1} = j_1,..., C_{t-l} = j_l) = \sum_{i-1}^l \lambda_i q(j_i,j_0)\] Donde \(\sum_{i=1}^l λ_i=1\), y \(Q=(q(j,k))\) es una matriz m×m con entradas no negativas y la suma de sus filas equivale a uno, tal que el lado derecho de la ecuación anterior esta limitada por cero y uno para todo \(j_o,j_1,…,j_l∈M\). Este último requerimiento, el cual general m^(l+1) pares de restricciones no lineares en los parámetros, asegura que las probabilidades condicionales en la ecuación anterior son en efecto probabilidades, y la condición en la suma de filas de Q asegura que la suma sobre \(j_0\) de estas probabilidades condicionales es uno. Note que Raftery no asume que los parámetros \(λ_i\), son no negativos.

---

Capítulo 2: Modelos de Markov ocultos: definición y propiedades

2.1. Un modelo de Markov oculto simple

Considérese nuevamente la serie de terremotos mostrada en el capitulo anterior. Las observaciones son cuentas ilimitadas, haciendo la distribución de Poisson una elección natural para describirlas, pero su distribución es claramente relativa sobre dispersa a la de Poisson. Observamos en el capítulo anterior que este hecho puede ser acomodado al usar una mezcla de distribuciones de Poisson con medias \(λ_1,λ_2,…,λ_m\). La elección de la media esta hecha por un segundo proceso aleatorio, el proceso de parámetro. La media \(λ_i\) es seleccionada con probabilidad \(δ_i\), donde \(i=1,2,…,m\) y \(\sum_{i=1}^mδ_i=1\).

Un modelo de mezcla independiente no resultara con la serie de terremotos porque, por definición, no permite la dependencia serial en las observaciones. La muestra de función de autocorrelación (ACF), desplegada en la Figura 2.1, claramente indica que las observaciones son dependientes en serie. Una manera de permitir la dependencia serial en las observaciones es relajar la suposición que el proceso de parámetro es independiente en serie. Una manera simple y matemáticamente conveniente para hacer esto es asumir que es una cadena de Markov. El modelo resultante para las observaciones es llamado un modelo de Markov oculto de Poisson.

library(forecast)
library(readr)
#leer la data 
data_terremotos_ch2 <- read.table("/Users/bracruz/Documents/8vo. Trimestre IO/Sist. Expertos y Redes Prob./earthquakes.txt")
#Obtener la ACF de la data anterior
ACF_ch2 <- acf(data_terremotos_ch2, plot=TRUE)

Figura 2.1

2.2 Lo básico

2.2.1 Definición y notación

Figura 2.2

Un modelo de Markov oculto \({X_t:t ∈N}\) es un tipo particular de mezcla dependiente. Con \(X^t\) y \(C^t\) representando las historias del tiempo 1 al tiempo t, uno puede resumir el modelo más simple de este tipo por:

\[ Pr(C_t | C^{t-1}) = Pr(C_t | C_{t-1}), t=2,3,... \\ Pr(X_t | X^{t-1},C^t) = Pr(X_t | C_t), t ∈ N \] El modelo consiste de dos partes: primero, un “proceso de parámetro” no observado \({C_t:t=1,2,…}\) satisfaciendo la propiedad de Markov; y como segundo, el “proceso estado-dependiente” \({X_t:t=1,2,…}\), en el cual la distribución de \(X_t\) depende solamente en el estado actual \(C_t\) y no de estados u observaciones previas. La estructura se representa en la Figura 2.2

Si la cadena de Markov \({C_t}\) tiene m estados, llamamos \({X_t}\) una HMM de m-estados. Aunque es la terminología usual en aplicaciones de procesamiento del habla, el nombre “modelo de Markov oculto” es por ningún medio el único usado para tales o similares.

Figura 2.3 - Proceso gerenando las observaciones en un estado doble de HMM. Siguiendo la cadena 2,1,1,1,2,1

El proceso generando las observaciones es nuevamente demostrado en la Figura anterior, para distribuciones estado-dependientes \(p_1\) y \(p_2\), distribución estacionaria \(δ=(0.75,0.25)\), y t.p.m. \(\Gamma= \left(\begin{array}{cc} 0.9 & 0.1\\ 0.3 & 0.7 \end{array}\right)\). En contraste al caso de una mezcla independiente, aquí la distribución de \(C_t\), el estado en tiempo t, sí depende de \(C_(t-1)\). También es cierto de las mezclas independientes, hay para cada estado una distribución diferente, discreta o continua. Ahora introduciremos una notación la cual cubrirá tanto observaciones valuadas discretas como continuas. En el caso de observaciones discretas, definimos, para \(i=1,2,…,m\)

\[p_i(x) = Pr(X_t = x | C_t = i)\] Esto es, \(p_i\), es la función de probabilidad en masa de \(X_t\) si la cadena de Markov está en el estado i en el tiempo t. El caso continuo es tratado de manera similar: aquí definimos \(p_i\) como la función de densidad de probabilidad de \(X_t\) asociada con el estado i. Nos referimos a las m distribuciones \(p_i\) como la distribución estado-dependiente del modelo.

2.2.2 Distribuciones marginales

A menudo necesitamos las distribuciones marginales de \(X_t\) y también distribuciones marginales de orden superior, tales como \((X_t,X_{t+k})\). Derivaremos los resultados para el caso en el cual la cadena de Markov es homogénea pero no necesariamente estacionaria, y luego veremos para el caso especial en donde la cadena de Markov es estacionaria. Por conveniencia la derivación es dada solo para las distribuciones estado-dependientes discretas; el caso continuo puede ser derivado análogamente.

Distribuciones de una variable

Para observaciones de valores discretos \(X_t\), definimos \(u_i(t)=Pra(C_t=i)\quad para\quad t=1,…,T\), tenemos

\[ Pr(X_t = x) = \sum_{i=1}^m Pr(C_t = i) Pr(X_t = x |C_t = i) \\ = \sum_{i=1}^m u_i(t)p_i(x) \] La expresión puede ser convenientemente rescrita en notacion de matriz

\[ Pr(X_t = x) = (u_1(t),...,u_m(t)) \left(\begin{array}{cc} p_1(x) & 0\\ 0 & p_m(x) \end{array}\right) \left(\begin{array}{cc} 1\\ 1 \end{array}\right) \\ =u(t) P(x)1' \] Donde \(P(x)\) es definida como la matriz diagonal con su i-esimo elemento diagonal \(p_i (x)\). Tomando en cuenta la ecuación del capítulo anterior que dicta que \(u(t)=u(1) Γ^{t-1}\), y por lo tanto

\[Pr(X_t = x) = u(1)\Gamma^{t-1} P(x)1'\] Esta ecuación funciona solamente si la cadena de Markov es puramente homogénea, y no necesariamente estacionaria. Si, como lo asumiremos usualmente, la cadena de Markov es estacionaria, con distribución estacionaria \(δ\), el resultado es aún más simple: en ese caso \(δΓ^{(t-1)}=δ\) para todo \(t∈N\), y entonces

\[Pr(X_t = x) = \delta P(x)1'\]

Distribuciones de dos variables

El cálculo de muchas de las distribuciones relacionadas a un HMM es realizado más fácilmente al notar de primero que, para un modelo grafico dirigido, la distribución conjunta de un set de variables aleatorias \(V_i\) esta dada por

\[Pr(V_1,V_2,...,V_n) = \prod_{i=1}^n Pr(V_i |pa(V_i))\]

Donde \(pa(V_i )\) denota el set de todos los “parientes” de \(V_i\) en el set \(V_1,V_2,…,V_n\). En el grafico dirigido de las cuatro variables aleatorias \(X_t,X_{(t+k)},C_t,C_{(t+k)}\) (para enteros positivos k), \(C_t\) no tiene parientes, \(pa(X_t )={C_t },pa(C_{(t+k)} )={C_t}\) y \(pa(X_{(t+k)} )={C_{(t+k)}}\). Y por tanto sigue que

\[ Pr(X_t,X_{t+k},C_t,C_{t+k}) = Pr(C_t)Pr(X_t|C_t) Pr(C_{t+k}|C_t) Pr(X_{t+k}|C_{t+k}) \] Y por lo tanto dicta que:

\[ Pr(X_t = v, X_{t+k}=w) = \sum_{i=1}^m \sum_{j=1}^m u_i(t)p_i(v)\gamma_{ij}(k)p_j(w) \]

(Aquí y en cualquier otro lado, \(γ_{ij} (k)\) denota los (i,j) elementos de \(Γ^k\).) Escribiendo la suma doble superior como el productor de matrices obtenemos

\[ Pr(X_t = v, X_{t+k}=w) = u(t) P(v)\Gamma^k P(w)1' \] Si la cadena de Markov es estacionaria, esto se reduce a

\[ Pr(X_t = v, X_{t+k}=w) = \delta P(v) \Gamma^k P(w)1' \]

Similarmente, uno puede obtener expresiones para las distribuciones marginales de orden superior; en el caso estacionario, la formula para una distribución de tres variables es, para enteros positivos k y l

\[ Pr(X_t=v, X_{t+k}=w, X_{t+k+l} =z) = \delta P(v) \Gamma^k P(w) \Gamma^l P(z)1' \]

2.2.3 Momentos

Primero notemos que

\[ E(X_t) = \sum_{i=1}^m E(X_t | C_t=i ) Pr(C_t = i) = \sum_{i=1}^m u_i(t) E(X_t|C_t = i) \] En el cual, el caso estacionario se reduce a

\[ E(X_t) = \sum_{i=1}^m \delta_i E(X_t| C_t = i) \] Mas generalmente, resultados analogos aguardan para \(E(g(X_t ))\) y \(E(g(X_t,X_{t+k} ))\), para cualquier función \(g\) para las cuales expectativas estado-dependientes relevantes existen. En el caso estacionario

\[ E(g(X_t)) = \sum_{i=1}^m \delta_i E(g(X_t) | C_t =i) \] y además

\[ E(g(X_t,X_{t+k})) = \sum_{i,j=1}^m E(g(X_t,X_{t+k})|C_t = i, C_{t+k} = j) \delta_i \gamma_{ij}(k) \] A menudo también estaremos interesados en una función g la cual factoriza a

\[g(X_t, X_{t+k}) = g_1(X_t) g_2(X_{t+k})\] En cual caso la ecuación resultante termina siendo

\[ E(g(X_t,X_{t+k})) = \sum_{i,j=1}^m E(g_1(X_t)| C_t = i) E(g_2(X_{t+k})|C_{t+k}=j) \delta_i \gamma_{ij}(k) \] Estas expresiones nos permiten encontrar covarianzas y correlaciones; expresiones explicitas convenientes existen para muchos casos. Para el caso de un HMM de Poisson estacionario de dos estados:

\[ E(X_t) = \delta_1 \lambda_1 + \delta_2 \lambda_2 \\ Var(X_t) = E(X_t) + \delta_1 \delta_2(\lambda_2 - \lambda_1)^2 \geq E(X_t) \\ Cov(X_t, X_{t+k}) = \delta_1 \delta_2(\lambda_2 - \lambda_1)^2 (1 - \gamma_{12} - \gamma_{21})^k, para \quad k \in N \]

Notese que la formula resultante para la correlación de \(X_t\) y \(X_{t+k}\) es de la forma \(ρ(k)=A(1-γ_{12}-γ_{21} )^k\) con \(A∈[0,1)\), y que A=0 si \(λ_1=λ_2\).

2.3. La Posibilidad

El objetivo de esta sección es desarrollar una formula conveniente para la posibilidad de \(L_T\) de T observaciones continuas \(x_1,x_2,…,x_T\) asumida a ser generada por un HMM de m-estados. Veremos que la computación de la posibilidad, consiste de la suma de \(m^T\) términos, cada cual es un producto de 2T factores, parece requerir O(Tm^T) operaciones. Nuestro propósito aquí es demostrar que \(L^T\) puede ser en general computado simplemente relacionado en \(O(Tm^2 )\) operaciones. La manera luego será ampliada para estimar parámetros por maximización numéricas de la posibilidad. Primero, la posibilidad de modelos de dos estados será explorada, y entonces la formula general será presentada.

2.3.1. La posibilidad de un HMM de Bernoulli de dos estados

Considere el HMM estacionario de dos estados con t.p.m.

\[ \Gamma = \left(\begin{array}{cc} 1/2 & 1/2\\ 1/4 & 3/4 \end{array}\right) \quad y \quad \delta = \left(\begin{array}{cc} 1/3, 2/3 \end{array}\right) \]

Y distribuciones estado-dependientes dadas por

\[ Pr(X-t =x | C_t = 1) = 1/2 \quad(para\quad x =0,1) \\ Pr(X-t =1 | C_t = 2) = 1 \]

Llamaremos un modelo de este tipo un HMM de Bernoulli. . Considere la probabilidad \(Pra(X_1=X_2=X_3=1\). Primero, note que, por ecuaciones anteriores,

\[ Pra(X_1,X_2,X_3,C_1,C_2,C_3 )=Pra(C_1 ) Pra(X_1│C_1 ) Pra(C_2│C_1 ) Pra(X_2│C_2 ) Pra(C_3│C_2 ) Pra(X_3 |C_3); \] luego sumarmos los valores asumidos por \(C_1,C_2,C_3\). El resultado es:

\[ Pr(X_1=1, X_2 = 1, X3=1) \\ = \sum_{i=1}^2 \sum_{j=1}^2 \sum_{k=1}^2 Pr(X_1 = 1,X_2=1,X_3=1, C_1 =i, C_2=j, C_3=k) \\ = \sum_{i=1}^2 \sum_{j=1}^2 \sum_{k=1}^2 \delta_i p_i(1) \gamma_{ij} p_j(1) \gamma_{jk} p_k(1) \] Note que la triple suma de la ecuación anterior tiene \(m^T=2^3\) términos, cada del cual es producto de \(2T=2 x 3 \)factores. Para evaluar la probabilidad requerida, las diferentes posibilidades para los valores de i,j y k pueden ser listados y la suma calculada como en la siguiente tabla

Probabilidad de la probabilidad de una HMM Bernoulli

La suma de la última columna de la tabla anterior nos dice que \(Pra(X_1=1,X_2=1,X_3=1)=29/48\) . Al observar notamos que el elemento más grande en esa columna es 3/8; la secuencia de estados que maximiza la probabilidad conjunta

\[Pr(X_1 = 1,X_2=1,X_3=1, C_1 =i, C_2=j, C_3=k)\] Es entonces la secuencia \(i=2,j=2,k=2\). Equivalentemente, maximiza la probabilidad condicional \(Pra(C_1=i,C_2=j,C_3=k| X_1=1,X_2=1,X_3=1)\). Pero una manera mas conveniente de presentar la sumar es utilizar notación de matriz. Definamos P(u) como \(diag(p_1 (u),p_2 (u))\). Entonces

\[ P(0) = \left(\begin{array}{cc} 1/2 & 0\\ 0 & 0 \end{array}\right) \quad y \quad P(1)= \left(\begin{array}{cc} 1/2 & 0\\ 0 & 1 \end{array}\right) \] Y la triple suma anterior puede ser escrita como el producto de una matriz

\[ δP(1)ΓP(1)ΓP(1)1' \]

2.3.2. La posibilidad en general

Aquí consideramos la posibilidad de un HMM en general. Suponemos hay una secuencia de observación \(x_1,x_2,…,x_T \)generada por tal modelo. Buscamos la probabilidad \(L_T\) de observar esa secuencia, según lo calculado bajo un HMM de m-estados el cual tiene una distribución inicial \(δ\) y t.p.m. \(Γ\) para la cadena de Markov, funciones de probabilidad estado-dependientes \(p_i\).

Proposicion 1: La probabilidad es dada por

\[ L_T = \delta P(x_1)\Gamma P(x_2) \Gamma P(x_3) \quad ... \quad \Gamma P(x_T)1' \] si \(\delta\), la distribucioón de \(C_1\), si es una distribucion estacionaria de la cadena de Markov, entonces:

\[ L_T = \delta \Gamma P(x_1)\Gamma P(x_2) \Gamma P(x_3) \quad ... \quad \Gamma P(x_T)1' \]

Una simple pero crucial consecuencia de la expresión matricial para la posibilidad es el “algoritmo de avance” para computación recursiva de la posibilidad. Tal computación recursiva juega un rol clave, no solo en evaluación de la posibilidad y por tanto estimación de parámetros, pero también en predecir, decodificar y revisar modelos. La naturaleza recursiva de tal evaluación de posibilidad vía la Proposición 1 es computacionalmente mucho mas efectivo que suma de fuerza bruta sobre todas las secuencias de estado posible. El hecho que tales esquemas computacionales recursivos no costosos pueden ser utilizados para abordar varias preguntas de interés es una característica clave de los HMMs. Para definir el algoritmo de avance, definimos el vector \(α_t\) para \(t=1,2,…,T\), por

\[ \alpha_t = \delta P(x_1)\Gamma P(x_2) \Gamma P(x_3) \quad ... \quad \Gamma P(x_t) = \delta P(x_1) \prod_{s=2}^t \Gamma P(x_s) \] Con la convención que un producto vacío es la matriz identidad. Sigue inmediatamente de esta definición que

\[ L_T = \alpha_T 1' \quad y \quad \alpha_t = \alpha_{t-1} \Gamma P(x_t) \quad para \quad t \ge 2 \] Acordemente, podemos convenientemente fijar como sigue las computaciones involucradas en la fórmula de posibilidad

\[ \alpha_1 = \delta P(x_1) \\ \alpha_t = \alpha_{t-1} \Gamma P(x_t) \quad para \quad t =2,3,...,T \\ L_T = \alpha_T 1' \]

El número de operaciones involucradas es de orden \(Tm^2\) puede ser deducido así. Para cada uno de los valores de t en el ciclo, hay m elementos de \(α_t\) por ser computados, y cada uno de esos elementos es una suma de m productos de tres cantidades: un elemento de \(α_{t-1}\), una probabilidad de transición \(γ_{ij}\), y una probabilidad (o densidad) estado-dependiente \(p_j (x_t)\).

El esquema de correspondiente para la computación de la Propuesta 1 es

\[ \alpha_0 = \delta \\ \alpha_t = \alpha_{t-1} \Gamma P(x_t) \quad para \quad t =1,2,...,T \\ L_T = \alpha_T 1' \] Los elementos del vector \(α_t\) son usualmente referidos como probabilidades de avance; donde mostramos que el j-esimo elemento de \(α_t\) es \(Pra(X^{(t)}=x^{(t)},C_t=j\). A continuación, mostramos código de R que utiliza el algoritmo de avance para evaluar la posibilidad de observaciones \(x_1,…x_T\) bajo un HMM de Poisson con por lo menos dos estados, t.p.m. \(Γ\), vector de medias estado-dependientes \(λ\), y distribución inicial \(δ\).

alpha <- delta*dpois(x[1], lambda)
for(i in 2:T) 
  alpha <- %*% Gamma*dpois(x[i],lambda)
sum(alpha)

En la discusión superior hemos utilizado la expresión de múltiples sumas para la posibilidad para así llegar a la expresión matricial, y luego utilizado la expresión matricial para llegar a la recursión de avance. Una ruta alternativa es definir el vector de avance de probabilidades \(α_t\) por

\[ \alpha_t(j) = Pr(X^{(t)} = x^{(t)}, C_t = j), j=1,2,...,m \] Y luego deducir la recursión de avance

\[ \alpha_t = \alpha_{t-1} \Gamma P(x_t) \] La expresión matricial es entonces una simple consecuencia de la recursión de avance.

2.3.3 HMMs no son procesos de Markov

HMMs no satisfacen en general la propiedad de Markov. Esto lo podemos demostrar mediante un simple contrajemplo. Basado en la seccion 2.3.1 sabemos que:

\[ Pr(X_1=1,X_2=1,X_3=1) = \frac{29}{48} \\ Pr(X_2=1) = \frac{5}{6} \\ entonces \\ Pr(X_1,X_2=1) = Pr(X_2=1,X_3=1) = \frac{17}{24} \\ por \quad tanto \quad se\quad deduce\quad que \\ Pr(X_3=1 | X_1=1,X_2=1) = \frac{Pr(X_1=1,X_2=1,X_3=1)}{Pr(X_1=1,X_2=1)} = \frac{29/48}{17/24} \\ y \quad eso \\ Pr(X_3=1|X_2=1) = \frac{Pr(X_2=1,X_3=1)}{Pr(X_2=1)} = \frac{17/24}{5/6} = 17/20 \]

Por lo tanto \(Pra(X_3=1|X_2=1)≠Pra(X_3=1|X_1=1,X_2=1\); este HMM no satisface la propiedad de Markov. Para las condiciones bajo las cuales un HMM si satisface la propiedad de Markov.

2.3.4. La posibilidad cuando hay datos faltantes

En el contexto de una serie de tiempo es potencialmente raro si algunos de los datos faltan. En el caso de series de tiempo de modelos ocultos de Markov, sin embargo, el ajuste que se necesita ser realizado para la computación de la posibilidad si hay datos faltantes resulta ser uno fácil. Suponga, por ejemplo, que uno tiene disponible las observaciones \(x_1,x_2,x_4,x_7,x_8,…,x_T\) de un HMM pero \(x_3,x_5 \quad y \quad x_6 \)hacen falta. Entonces la posibilidad de las observaciones esta dada por

\[ Pr(X_1=x1,X_2=x2,X_4=x_4,X_7=x_7,X_8=x_8,...,,X_T=x_T) = \\ \sum \delta_{c_1}\gamma_{c1,c2}\gamma_{c2,c4}(2) \gamma_{c4,c7} (3) \gamma_{c7,c8} \quad ... \quad \gamma_{c_{T-1},c_T} \\ \times p_{c_1} (x_1) p_{c2}(x_2) p_{c4}(x_4) p_{c7}(x_7)\quad ... \quad p_{c_T}(x_T) \]

Donde (como antes) \(γ_{ij} (k)\) denota una probabilidad de transición de paso k, y la suma se toma sobre todos los índices \(c_t\) diferentes de \(c_3,c_5 y c_6\). Pero esto es solamente

\[ \sum \delta_{c_1} p_{c_1} (x_1) \gamma_{c_1},c_2p_{c_2} (x_2) \gamma_{c_2},c_4(2)p_{c_4} (x_4) \gamma_{c_4},c_7(3)p_{c_7} (x_7) \gamma_{c_7} \quad \times...\times \quad \gamma_{c_T-1}, c_T p_{c_T}(x_T)= \\ \delta P(x_1) \Gamma P(x_2)\Gamma^2 P(x_4) \Gamma^3 P(x_7)...\Gamma P(x_T)1' \]

Con \(L_T^{-(3,5,6)}\) denotando la posibilidad de las observaciones (diferentes de las faltantes), se declara que

\[ L_T^{-(3,5,6)} = P(x_1) \Gamma P(x_2)\Gamma^2 P(x_4) \Gamma^3 P(x_7)...\Gamma P(x_T)1' \]

En general, la expresión para la posibilidad de matrices diagonales \(P(x_t)\) correspondiente a las observaciones faltantes \(x_t\) son remplazadas por la matriz identidad; eso es, las probabilidades estado-dependientes correspondientes \(p_i (x_t )\) son remplazadas por 1 para todos los estados i.

2.3.5. La posibilidad cuando las observaciones son censoradas por intervalo

Suponga que deseamos ajustar un HMM de Poisson a una serie de cuentas, algunas de las cuales son censoradas por intervalo. Por ejemplo, el valor de \(x_t\) puede ser conocido para \(4≤t≤T\), con la información \(x_1≤5,2≤x_2≤3\) y \(x_3>10\) disponible acerca de las observaciones restantes. Para simplificar, asumamos de primero que la cadena de Markov solo tiene dos estados. En ese caso, uno remplaza las matrices diagonales \(P(x_i )(i=1,2,3)\) en la expresión de posibilidad por las matrices

\[ diag(Pr(X_1 \le |C_1 =1), Pr(X_1 \le 5| C_1=2)), \\ diag(Pr(2 \le X_2 \le 3 | C_2=1),Pr(2 \le X_2 \le 3 | C_2=2)) \\ diag(Pr(X_3 \gt 10 | C_3 = 1), Pr(X_3 \gt 10 | C_3 = 2)) \]

Mas generalmente, suponga que \(a≤x_t≤b\) donde a puede ser \(-∞\) (aunque eso no es relevante para el caso de Poisson), b puede ser \(∞\), y la cadena de Markov tiene m estados. Uno remplaza \(P(x_t)\) en la posibilidad por la matriz diagonal \(m×m\) de la cual el elemento diagonal i-esimo es \(Pra(a≤X_t≤b|C_t=i)\)

---

Capítulo 3: Estimación de la posibilidad por maximización directa

3.1. Introducción

Vimos anteriormente que la posibilidad de un HMM está dada por

\[ L_T = Pr(X^{(T)} = x^{(T)}) = \delta P(x_1) \Gamma P(x_2)... \Gamma P(x_T)1' \] Donde \(δ\) es la distribución inicial y P(x) la matriz diagonal m×m con su i-esimo elemento diagonal siendo la probabilidad estado-dependiente o densidad \(p_i (x)\). En principio, podemos entonces computar \(L_T=α_T 1'\) recursivamente medianteDonde \(δ\) es la distribución inicial y P(x) la matriz diagonal m×m con su i-esimo elemento diagonal siendo la probabilidad estado-dependiente o densidad \(p_i (x)\). En principio, podemos entonces computar \(L_T=α_T 1'\) recursivamente mediante

\[ \alpha_1 = \delta P(x_1) \\ y\\ \alpha_t = \alpha_{t-1} \Gamma P(x_t), para \quad t=2,3,...,T \]

Si se asume que la cadena de Markov es estacionaria (en cuyo caso \(δ=δΓ\), podemos en vez usar

\[ \alpha_0 = \delta \\ y\\ \alpha_t = \alpha_{t-1} \Gamma P(x_t). para\quad t=1,2,...,T \]

Primero consideraremos el caso estacionario.

El número de operaciones involucradas es de orden \(Tm^2\), haciendo la evaluación de la posibilidad bastante factible incluso para un gran T. La estimación de parámetros puede ser entonces realizada por maximización numérica de la posibilidad respecto a los parámetros. Pero existen algunos problemas que necesitan ser abordados cuando los parámetros son estimados de esta manera. Los problemas principales son el desbordamiento numérico, restricciones en los parámetros, y múltiples máximos locales en la función de posibilidad.

3.2. Escalando la computación de la posibilidad

En el caso de distribuciones estado-dependientes discretas, los elementos de \(α_t\), siendo productos de probabilidades, se vuelven progresivamente más pequeño conforme t incremente, y son eventualmente redondeados a cero. De hecho, con una probabilidad de 1, la posibilidad se acerca a 0 (o posiblemente \(∞\) en el caso continuo) exponencialmente rápido. Desde que la posibilidad es un producto de matrices, no escalares, no es posible evitar el desbordamiento numero simplemente al computar el registro de la posibilidad como la suma de los registros de sus factores. En este caso, la computación de la posibilidad de un modelo de mezcla independiente es más fácil que el de un HMM. Para resolver este problema, Durbin (1998, p. 78) sugiere un método de computación que cuenta con la siguiente aproximación. Suponga que deseamos computar \(log(p+q)\), donde\(p>q\). Escriba \(log(p+q)\) como \[ log(p)+log(1+q/p)=log(p)+log(1+exp(q̃-p ̃)) \] , donde \(p ̃=log(p) \quad y \quad q ̃=log(q)\). La función \(log(1+e^x)\) es entonces aproximada por interpolación de una tabla de sus valores; aparentemente tal pequeña tabla dará un grado razonable de exactitud. Preferimos computar el logaritmo de \(L_T\) al usar una estrategia de escalar el vector de probabilidades de avance \(α_t\). Efectivamente escalamos el vector \(α_t\) en cada tiempo t para que sus elementos sumen a 1, llevando récord de la suma de los registros de los factores de escala aplicados. Defina, para \(t=0,1,…,T\) el vector

\[ \phi_t =\alpha_t/w_t\]

Donde \(w_t=\sum_i α_t (i)=α_t 1'\). Primero, notamos ciertas consecuencias inmediatas de las definiciones de \(ϕ_t\) y \(w_t\):

\[ w_o = \alpha_01' = \delta1' = 1 \\ \phi_0 = \delta \\ w_t \phi_t = w_{t-1} \phi_{t-1} \Gamma P(x_t) \\ L_T = \alpha_T 1' = w_T(\phi_T 1') = w_T \]

Por lo tanto \(L_T=w_T=∏_{t=1}^T(w_t/w_{t-1})\). Y de la ecuación anterior tenemos que \(w_t=w_{t-1} (ϕ_{t-1} Γ P(x_t)1')\), y asi concluimos que

\[ logL_T = \sum_{t=1}^T log(w_t/w_{t-1}) = \sum_{t=1}^T log(\phi_{t-1} \Gamma P(x_t) 1') \]

La computación de la posibilidad de registro esta sumarizada a continuación en la forma de un algoritmo. Note que \(Γ\) y \(P(x_t)\) son matrices de m×m, v y \(ϕ_t\) son vectores de longitud m, u es un escalar, y l es el escalar en el cual la posibilidad de registro es acumulada.

\[ \phi_0 \leftarrow \delta; l \leftarrow0 \\ for \quad t=1,2,...,T \\ v \leftarrow \phi_{t-1} \Gamma P(x_t) \\ u \leftarrow v1' \\ l \leftarrow l + log(u) \\ \phi_t \leftarrow v/u \\ return \quad l \]

La posibilidad de registro requerida, \(loga(L_T)\), es entonces dada por el valor final de l. Este procedimiento prevendrá casi siempre desbordamiento. Claramente, variaciones menores de esta técnica son posibles: el factor de escala \(w_t\) podría ser elegido en lugar para ser el elemento más grande del vector siendo escalado, o la media de sus elementos (opuesto a la suma). El algoritmo es fácilmente modificado para computar la posibilidad de registro sin asumir estacionariedad de la cadena de Markov. Con \(δ\) denotando la distribución inicial, el algoritmo más general es

\[ w_1 \leftarrow \delta P(x_1) 1' ;\quad \phi_1 \leftarrow \delta P(x_1)/w_1 ; \quad l \leftarrow log(w_1) \\ for \quad t=2,3,...,T \\ v \leftarrow \phi_{t-1} \Gamma P(x_t) \\ u \leftarrow v1' \\ l \leftarrow l + log(u)\\ \phi_t \leftarrow v/u \\ return \quad l \]

Si la distribución inicial ocurre ser la distribución inicial, este algoritmo todavía aplica. El siguiente código implementa esta ultima versión del algoritmo para asi computar la posibilidad de registro de las observaciones \(x_1,…,x_T\) bajo un HMM de Poisson con al menos dos estados, matriz de probabilidad de transición \(Γ\), vector de medias estado-dependientes \(λ\) y distribución inicial \(δ\).

alpha <- delta*dpois(x[1], lambda)
lscale <- log(sum(alpha))
alpha <- alpha/sum(alpha)
for(i in 2:T){
  alpha <- alpha %*% Gamma*dpois(x[i], lambda)
  lscale <- lscale + log(sum(alpha))
  alpha <- alpha/sum(alpha)
}
lscale

3.3. Maximización de la posibilidad sujeta a restricciones

3.3.1 Re parametrización para evitar restricciones

Los elementos de \(Γ\) y aquellos de \(λ\), el vector de medias estado-dependientes en un HMM de Poisson, están sujetos a no negatividad y otras restricciones. En particular, la suma de las filas de \(Γ\) es igual a 1. Los estimados de los parámetros debería satisfacer también tales restricciones. Por ende, al maximizar la posibilidad, necesitamos resolver un problema de optimización restringida. En general, existen dos grupos de restricciones: aquellas que aplican a los parámetros de las distribuciones estado-dependientes y aquellas que aplican a los parámetros de la cadena de Markov. El primer grupo de restricciones depende en cual(es) distribución(es) estado-dependiente(s) son elegidas.

En el caso de un HMM de Poisson las restricciones relevantes son:

• Las medias \(λ_i\) de las distribuciones estado-dependientes deben, para \(i=1,…,m\) ser no negativos.
• Las filas de la matriz de probabilidad de transición \(Γ\) deben sumar a 1, y todos los parámetros \(γ_{ij}\) deben ser no negativos.

Aquí las restricciones pueden ser impuestas al hacer transformaciones. La transformación de los parámetros \(λ_i\) es fácil. Defina \(η_i=loga(λ_i)\), para \(i=1,…,m\). Entonces \(η_i∈R\). Luego que hemos maximizado la posibilidad con respecto a los parámetros no restringidos, los parámetros estimados restringidos pueden ser obtenidos al transformar de vuelta: \((λ_i )=expa(η_i))\).Note que \(Γ \) tiene \(m^2\) entradas pero solo \(m(m-1)\) parametros libres, asi como hay m restricciones suma de filas

\[ \sum_{j=1}^m \gamma_{ij} =1 \quad (i=1,...,m) \]

Mostraremos una posible transformación entre las \(m^2\) probabilidades restringidas \(γ_{ij}\) y \(m(m-1)\) números reales sin restricción \(τ_{ij},i≠j\). Por razones de visibilidad se muestra el caso para m=3. Empezaremos por definir la matriz

\[ \left(\begin{array}{cc} - , \tau_{12},\tau_{13}\\ \tau_{21} , - , \tau_{23}\\ \tau_{31} , \tau_{32} , - \end{array}\right) \]

Una matriz on \(m(m-1)\) entradas \(τ_{ij}∈R\). Ahora que \(g:R→R^+\) sea una función estrictamente incremental, por ejemplo

\[ g(x) = e^x \quad o \quad g(x) = \begin{cases} e^x \quad x \le 0 \\ x+1 \quad x \ge 0 \end{cases} \\ Define \quad que \\ v_{ij} = \begin{cases} g(\tau_{ij}) \quad para \quad i \neq j \\ 1 \quad para \quad i = j \end{cases} \]

Entonces fijamos \(γ_{ij}=v_{ij}/\sum_{k=1}^m v_{ik}\) (para i,j=1,2,…,m) y \(Γ=(γ_{ij})\). Nos referiremos a los parámetros \(η_i\) y \(τ_{ij}\) como parámetros funcionales, y a los parámetros \(λ_i\) y \(γ_{ij}\) como parámetros naturales.

Usando las transformaciones anteriores de \(Γ\) y \(λ\), podemos realizar el calculo de los parámetros maximizantes de la posibilidad en dos pasos

• Maximizar \(L_T\) con respecto a los parámetros funcionales \(T={τ_{ij}}\) y \(η=(η_1,…,η_m)\). Todos estos son sin restricciones.
• Transformar los estimados de los parámetros funcionales a estimados de parámetros naturales:

\[ T \rightarrow \Gamma , \eta \rightarrow \lambda \]

Considere \(Γ\) para el caso \(g(x)=e^x\) y m general. Tenemos entonces

\[ \gamma_{ij} = \frac{exp(\tau_{ij})}{1+ \sum_{k \neq i} exp(\tau_{ik})},\quad para \quad i\neq j \] Y los elementos diagonales de \(Γ\) siguen de la suma de filas de 1. La transformación en la dirección opuesta es

\[ \tau_{ij} = log(\frac{\gamma_{ij}}{1-\sum_{k \neq i}\gamma_{ik}}) = log(\gamma_{ij}/\gamma_{ii}),\quad para \quad i \neq j \] Ahora mostramos un código relativamente simple que transformara parámetros naturales a funcionales y viceversa. Este código refiere a un HMM de Poisson con m>2 estados, en el cual la cadena de Markov puede, si es apropiado, ser asumida estacionaria. En ese caso la distribución estacionaria \(δ\) no es dada, pero es computada cuando sea necesitada del t.p.m. \(Γ\) al resolver \(δ(I_m-Γ+U)=1\). De otra manera, \(δ\) es tratada como un parámetro (natural) y transformada para así remover las restricciones \(δ_i≥0\) y \(\sum_iδ_i=1\).

# Transformar los parámetros naturales de Poisson en parámetros de trabajo.
pois.HMM.pn2pw <- function(m,lambda,gamma,delta=NULL,stationary=TRUE)
{
 tlambda <- log(lambda)
 foo     <- log(gamma/diag(gamma))
 tgamma  <- as.vector(foo[!diag(m)])
 if(stationary) {tdelta <-NULL} else {tdelta<-log(delta[-1]/delta[1])}
 parvect <- c(tlambda,tgamma,tdelta)
 return(parvect)
}

# Transformar los parámetros de trabajo de Poisson en parámetros naturales.

pois.HMM.pw2pn <- function(m,parvect,stationary=TRUE)
{
 lambda        <- exp(parvect[1:m])
 gamma         <- diag(m)
 gamma[!gamma] <- exp(parvect[(m+1):(m*m)])
 gamma         <- gamma/apply(gamma,1,sum)
 if(stationary) {delta<-solve(t(diag(m)-gamma+1),rep(1,m))} else
                {foo<-c(1,exp(parvect[(m*m+1):(m*m+m-1)]))
                delta<-foo/sum(foo)}
 return(list(lambda=lambda,gamma=gamma,delta=delta))
}

3.4. Otros problemas

3.4.1. Múltiples máximos en la posibilidad

La posibilidad de un HMM es una función complicada de los parámetros y frecuentemente tiene varios máximos locales. La meta por supuesto es encontrar el máximo global, pero no hay un método simple de determinar en general si un algoritmo de maximización numérica ha alcanzado el máximo global. Dependiendo de los valores iniciales, puede suceder que el algoritmo fácilmente identifique un local, pero no el máximo, global. Esto también aplica al principal método alterno de estimación, el algoritmo EM. Una estrategia sensible entonces es usar un rango de valores iniciales para la maximización, y ver si el mismo máximo es identificado en cada caso.

3.4.2. Valores iniciales para las iteraciones

Es a menudo fácil el encontrar posibles valores iniciales para algunos de los parámetros de un HMM: por ejemplo, si uno busca el fijar un HMM de Poisson con dos estados, y la media de muestra es 10, uno podría probar 8 y 12, o 5 y 15, para los valores de las dos medias estado-dependientes. Estrategias mas sistemáticas basadas en los cuantiles de las observaciones son posibles, sin embargo. Por ejemplo, si el modelo tiene tres estados, use como los valores iniciales de las medias estado-dependientes el cuartil menor, mediano y cuartil superior de las cuentas observadas. Es menos fácil el adivinar valores de las probabilidades de transición \(γ_{ij}\). Una estrategia es el asignar un valor inicial común a todas las posibilidades de transición fuera de la diagonal. Una consecuencia de tal decisión, quizá conveniente, es que la distribución estacionaria correspondiente es uniforme sobre los estados; esto sigue por simetría. Elegir buenos valores iniciales para los parámetros tiende a alejarlo a uno de inestabilidad numérica.

3.4.3. Posibilidad sin limites

En el caso de HMMs con distribuciones continuas estado-dependientes, puede suceder que la posibilidad es sin límites en la proximidad de ciertas combinaciones de parámetros.

Serie de Terremotos - Seccion izquierda: Distribuciones de Poisson - HMM con dos y tres estados. Seccion Derecha: medios dependientes del estado en comparacion con las observaciones

Como antes, sugerimos que, si esto crea dificultado, se maximize la posibilidad discreta en vez de la densidad conjunta.

3.5. Ejemplo: Terremotos

La figura anterior muestra el resultado de ajustar HMMs de Poisson (estacionarios) con dos a tres estados a la serie de terremotos por medias del optimizador no restringido nlm. El modelo de dos estados es

\[ \Gamma = \left(\begin{array}{cc} 0.9340 , 0.00660\\ 0.1285 , 0.8715 \end{array}\right) \]

Con \(δ=(0.6608,0.3392)\),\(λ=(15.472,26.125)\), y la posibilidad de registro dada por \(l=-342.3183\). Es claro que la mezcla (dependiente de Markov) ajustada de dos distribuciones de Poisson provee un mejor ajuste a la distribución marginal de las observaciones que lo hace una única distribución de Poisson, pero el ajuste puede seguir siendo mejorado al usar una mezcla de 3 o 4 distribuciones de Poisson.

library(HiddenMarkov)
library(readr)
#-----  Poisson Distribution  -----
#leer la data 
x_data<- read.table("/Users/bracruz/Documents/8vo. Trimestre IO/Sist. Expertos y Redes Prob./earthquakes.txt")[,2]
Pi <- matrix(c(0.9340, 0.066,
               0.1285, 0.8715),
             byrow=TRUE, nrow=2)
delta <- c(0.6608, 0.3392)
x <- dthmm(x_data, Pi, delta, "pois", list(lambda=c(15.472, 26.125)),
           discrete = TRUE)
#    use above parameter values as initial values
y <- BaumWelch(x)

iter = 1 
LL = -342.3180692 
diff = Inf 

iter = 2 
LL = -341.8839295 
diff = 0.4341397 

iter = 3 
LL = -341.8808283 
diff = 0.003101278 

iter = 4 
LL = -341.8799367 
diff = 0.0008916081 

iter = 5 
LL = -341.8794228 
diff = 0.0005138536 

iter = 6 
LL = -341.8791236 
diff = 0.0002991765 

iter = 7 
LL = -341.8789489 
diff = 0.0001747594 

iter = 8 
LL = -341.8788465 
diff = 0.000102314 

iter = 9 
LL = -341.8787865 
diff = 5.999939e-05 

iter = 10 
LL = -341.8787513 
diff = 3.522865e-05 

iter = 11 
LL = -341.8787306 
diff = 2.070391e-05 

iter = 12 
LL = -341.8787184 
diff = 1.217641e-05 

iter = 13 
LL = -341.8787113 
diff = 7.165123e-06 

print(summary(y))

$delta
[1] 1.000000e+00 2.259134e-31

$Pi
          [,1]       [,2]
[1,] 0.9284124 0.07158762
[2,] 0.1191548 0.88084524

$nonstat
[1] TRUE

$distn
[1] "pois"

$pm
$pm$lambda
[1] 15.42342 26.02398


$discrete
[1] TRUE

$n
[1] 107

#log-likelihood ideal
print(logLik(y))

[1] -341.8787

{h<-hist(residuals(y))
xfit <- seq(min(residuals(y)), max(residuals(y)), length = 40) 
yfit <- dnorm(xfit, mean = mean(residuals(y)), sd = sd(residuals(y))) 
yfit <- yfit * diff(h$mids[1:2]) * length(residuals(y)) 
lines(xfit, yfit, col = "#0066ff", lwd = 2)}

El modelo de tres estados es:

\[ \Gamma = \left(\begin{array}{cc} 0.9555 , 0.024 , 0.021\\ 0.050 , 0.8899 , 0.051 \\ 0.000 , 0.197 , 0.803 \end{array}\right) \]

Con \(δ=(0.4436,0.4045,0.1519)\),\(λ=(13.146,19.721,29.7144)\) y \(l=-329.4603\).

library(HiddenMarkov)
library(readr)
#-----  Poisson Distribution  -----
#leer la data 
x_data<- read.table("/Users/bracruz/Documents/8vo. Trimestre IO/Sist. Expertos y Redes Prob./earthquakes.txt")[,2]
Pi <- matrix(c(0.9555 , 0.024 , 0.021 ,
0.050 , 0.8899 , 0.051 , 
0.000 , 0.197 ,  0.803),
             byrow=TRUE, nrow=3)
delta <- c(0.4436,0.4045,0.1519)
x <- dthmm(x_data, Pi, delta, "pois", list(lambda=c(13.146,19.721,29.714)),
           discrete = TRUE)
#    use above parameter values as initial values
y <- BaumWelch(x)

iter = 1 
LL = -329.9251807 
diff = Inf 

iter = 2 
LL = -328.5465037 
diff = 1.378677 

iter = 3 
LL = -328.5278471 
diff = 0.01865661 

iter = 4 
LL = -328.5274950 
diff = 0.0003520581 

iter = 5 
LL = -328.5274853 
diff = 9.758334e-06 

print(summary(y))

$delta
[1] 9.999999e-01 9.681354e-08 6.078985e-21

$Pi
           [,1]       [,2]       [,3]
[1,] 0.93930599 0.03207267 0.02862135
[2,] 0.04039751 0.90637871 0.05322378
[3,] 0.00000000 0.19036700 0.80963300

$nonstat
[1] TRUE

$distn
[1] "pois"

$pm
$pm$lambda
[1] 13.13389 19.71215 29.70890


$discrete
[1] TRUE

$n
[1] 107

#log-likelihood ideal
print(logLik(y))

[1] -328.5275

{h<-hist(residuals(y))
xfit <- seq(min(residuals(y)), max(residuals(y)), length = 40) 
yfit <- dnorm(xfit, mean = mean(residuals(y)), sd = sd(residuals(y))) 
yfit <- yfit * diff(h$mids[1:2]) * length(residuals(y)) 
lines(xfit, yfit, col = "#ff0000", lwd = 2)}

El modelo de cuatro estados es: \[ \Gamma = \left(\begin{array}{cc} 0.805, 0.102, 0.093, 0.000 \\ 0.000, 0.976, 0.000, 0.024 \\ 0.050, 0.000, 0.902, 0.048 \\ 0.000, 0.000, 0.188, 0.812 \end{array}\right) \]

Con \(δ=(0.0936,0.3983,0.3643,0.1439)\),\(λ=(11.283,13.853,19.695,29.700)\) y \(l=-327.8316\).

library(HiddenMarkov)
library(readr)
#-----  Poisson Distribution  -----
#leer la data 
x_data<- read.table("/Users/bracruz/Documents/8vo. Trimestre IO/Sist. Expertos y Redes Prob./earthquakes.txt")[,2]
Pi <- matrix(c(0.805, 0.102, 0.093, 0.000,
0.000, 0.976, 0.000, 0.024,
0.050, 0.000, 0.902, 0.048,
0.000, 0.000, 0.188, 0.812),
             byrow=TRUE, nrow=4)
delta <- c(0.0936,0.3983,0.3643,0.1439)
x <- dthmm(x_data, Pi, delta, "pois", list(lambda=c(11.283,13.853,19.695,29.700)),
           discrete = TRUE)
#    use above parameter values as initial values
y <- BaumWelch(x)

iter = 1 
LL = -327.8302182 
diff = Inf 

iter = 2 
LL = -326.9581836 
diff = 0.8720346 

iter = 3 
LL = -326.9289508 
diff = 0.02923274 

iter = 4 
LL = -326.9204573 
diff = 0.008493573 

iter = 5 
LL = -326.9134060 
diff = 0.007051262 

iter = 6 
LL = -326.9074158 
diff = 0.005990187 

iter = 7 
LL = -326.9024761 
diff = 0.004939741 

iter = 8 
LL = -326.8985207 
diff = 0.003955416 

iter = 9 
LL = -326.8954326 
diff = 0.003088091 

iter = 10 
LL = -326.8930709 
diff = 0.002361707 

iter = 11 
LL = -326.8912938 
diff = 0.001777022 

iter = 12 
LL = -326.8899734 
diff = 0.00132043 

iter = 13 
LL = -326.8890015 
diff = 0.0009718974 

iter = 14 
LL = -326.8882912 
diff = 0.0007103225 

iter = 15 
LL = -326.8877747 
diff = 0.0005164488 

iter = 16 
LL = -326.8874007 
diff = 0.0003740631 

iter = 17 
LL = -326.8871305 
diff = 0.0002701856 

iter = 18 
LL = -326.8869357 
diff = 0.0001947661 

iter = 19 
LL = -326.8867955 
diff = 0.0001401981 

iter = 20 
LL = -326.8866947 
diff = 0.000100815 

iter = 21 
LL = -326.8866223 
diff = 7.244193e-05 

iter = 22 
LL = -326.8865702 
diff = 5.202691e-05 

iter = 23 
LL = -326.8865329 
diff = 3.735121e-05 

iter = 24 
LL = -326.8865061 
diff = 2.680813e-05 

iter = 25 
LL = -326.8864868 
diff = 1.923742e-05 

iter = 26 
LL = -326.8864730 
diff = 1.380287e-05 

iter = 27 
LL = -326.8864631 
diff = 9.902635e-06 

print(summary(y))

$delta
[1]  8.887835e-05  9.999111e-01  1.172905e-46 1.444604e-128

$Pi
           [,1]       [,2]      [,3]       [,4]
[1,] 0.80911120 0.08812885 0.1027600 0.00000000
[2,] 0.00000000 0.95349019 0.0000000 0.04650981
[3,] 0.04138593 0.00000000 0.9118154 0.04679867
[4,] 0.00000000 0.00000000 0.1778585 0.82214151

$nonstat
[1] TRUE

$distn
[1] "pois"

$pm
$pm$lambda
[1] 11.25002 13.80970 19.69526 29.66059


$discrete
[1] TRUE

$n
[1] 107

#log-likelihood ideal
print(logLik(y))

[1] -326.8865

{h<-hist(residuals(y))
xfit <- seq(min(residuals(y)), max(residuals(y)), length = 40) 
yfit <- dnorm(xfit, mean = mean(residuals(y)), sd = sd(residuals(y))) 
yfit <- yfit * diff(h$mids[1:2]) * length(residuals(y)) 
lines(xfit, yfit, col = "#00ffff", lwd = 2)}

En todos estos casos el ACF es solo una combinación lineal de las k-esimas potencias de los eigenvalores distintos de 1 de la matriz de probabilidad de transición. Un fenómeno que es notable cuando uno ajusta modelos con tres o mas estados a series relativamente cortas es que los estimados de una o mas de las probabilidades de transición resulta ser muy cercana a cero.

Este fenómeno puede ser explicado como sigue. En una cadena de Markov estacionaria, el número esperado de transiciones del estado i al estado j en una serie de T observaciones es \((T-1) δ_i\) \(γ_{ij}\). Para \(δ_3=0.152\) y T=107 (como en nuestro modelo de tres estados), esta expectativa será menor que 1 si \(γ_{31}<0.062\). En tal serie, por ende, es probable que si \(γ_{31}\) es bastante pequeño no habrá transiciones del estado 3 al estado 1, y entonces cuando busquemos estimar \(γ_{31}\) en un HMM, el estimado es probable a ser efectivamente cero. Como m incremente, las probabilidades \(δ_i\) y \(γ_{ij}\) se vuelve más pequeña en promedio; esto lo hace mas probable que al menos una probabilidad de transición sea efectivamente cero.

3.6. Errores estándares e intervalos de confianza

Relativamente poco es conocido acerca de las propiedades de estimadores de máxima posibilidad de HMMs; solo resultados asintóticos están disponibles. Para explotar estos resultados se requiere estimar la matriz de varianza y covarianza de los estimadores de los parámetros. Se puede estimar los errores estándares del Hessiano de la posibilidad de registro en el máximo, pero este enfoque se encuentra con dificultades cuando algunos de los parámetros están en la frontera del espacio de su parámetro, el cual ocurre bastante a menudo cuando los HMMs son ajustados.

3.6.1. Errores estándares mediante el Hessiano

Aunque los estimados de punto \(\Theta = (\Gamma,\lambda)\) son fáciles de computar, estimados exacto de intervalo no están disponibles. Cappe (2005, Capitulo 12) muestra que, bajo ciertas condiciones de regularidad, los MLEs de parámetros de un HMM son consistentes, asintóticamente normales y eficientes. Por ende, si podemos estimar los errores estándares de los MLEs entonces, usando normalidad asintótica, podemos también computar intervalos aproximados de confianza. Sin embargo, como señala Fruhwirth-Schnatter (2006, p. 53) en el contexto de modelos independientes de mezcla, “Las condiciones de regularidad son a menudo violados, incluyendo casos de gran preocupación práctica, entre ellos sets de datos pequeños, mezclas con pesos de componentes pequeños, y mezclas de excesos con muchos componentes”. Además, McLachlan y Peel (2000, p. 68) advierten: “En particular, para modelos de mezcla, es bien conocido que el tamaño de muestra n tiene que ser muy grande antes que la teoría asintótica de la posibilidad máxima se aplique”. Con las advertencias previas en mente, podemos, en orden para estimar los errores estándares de los MLEs de un HMM, utilizar el Hessiano aproximado de menos la posibilidad de registro en el mínimo. Podemos invertirlo y así estimar la matriz de varianza y covarianza asintótica de los estimadores de los parámetros. Un problema con esta sugerencia es que, si los parámetros han sido transformados, el Hessiano disponible será el cual se refiere a los parámetros funcionales \(ϕ_i\), no los originales, para interpretarlo de mejor manera, parámetros naturales \(\Theta_i\) (\(Γ\) y \(λ\) en el caso de un HMM de Poisson). La situación es entonces que tenemos, en el mínimo de \(-l\), el Hessiano con respecto a los parámetros funcionales

\[ H = -(\frac{\partial^2l}{\partial \phi_i \partial \phi_j}) \]

Pero lo que realmente necesitamos es el Hessiano en respecto a los parámetros naturales

\[ G = -(\frac{\partial^2l}{\partial \theta_i \partial \theta_j}) \]

Existe, sin embargo, la siguiente relación entre los dos Hessianos en el mínimo

\[ H = MGM' \quad y \quad G^{-1} = M'H^{-1}M \]

Donde M es definido por \(m_{ij}=∂θ_j/∂ϕ_i\). En el caso de un HMM de Poisson, los elementos de M son bastante simples. Con M a nuestra disposición, podemos usar la relación anterior para deducir \(G^{-1}\) de \(H^{-1}\), y utilizar \(G^{-1}\) para encontrar errores estándares para los parámetros naturales, tales parámetros provistos no están en la frontera del espacio del parámetro. Una ruta alternativa para los errores estándares con respecto a los parámetros naturales el cual a menudo resulta bien, y es menos laboriosa, es la siguiente. Primero encuentre el MLE al resolver el problema de optimización restringida, luego vuelva a correr la optimización sin restricciones, empezando en o muy cercano al MLE. Si el estimado resultante es el mismo que el MLE ya encontrado, el Hessiano correspondiente entonces provee directamente los errores estándares con respecto a los parámetros naturales. Pero si uno hace la suposición de normalidad y basa un intervalo de confianza en él, tal suposición de normalidad es mas probable, pero no garantizada.

Computación recursiva del Hessiano

Un método alternativo de computar el Hessiano es el de Lystig y Hughes (2002). Ellos presentan el algoritmo de avance \(α_t=α_{t-1} ΓP(x_t)\) de manera en la cual incorporan una escala automática o “natural”, y luego extienden este enfoque para así computar su Hessiano y gradiente con respecto a los parámetros naturales, aquellos que hemos descrito arriba por \(θ_i\). Turner (2008) ha usado este enfoque para así buscar las derivadas analíticas requeridas para maximizar posibilidades de HMMs directamente por el algoritmo Levenberg-Marquardt. Aunque esto puede ser un método mas eficiente y preciso de computar el Hessiano que utilizar la relación descrita anteriormente, no resuelve el problema fundamental que el uso del Hessiano para computar errores estándares (y de allí intervalos de confianza) no es de confianza si algunos de los parámetros están en o cerca de la frontera del espacio del parámetro.

3.6.2. Errores estándares Bootstrap e intervalos de confianza

Como una alternativa a las técnicas descritas en la sección anterior, uno puede usar el Bootstrap paramétrico. Hablando burdamente, la idea del Bootstrap paramétrico es el evaluar las propiedades del modelo con parámetros \(\Theta\)usando aquellos del modelo con parámetros \(\widehat\Theta\). Los siguientes pasos son realizados para estimar la matriz de varianza-covarianza de \(\widehat\Theta\).

1. Ajuste el modelo, compute \(\widehat\Theta\).
2. 1. Genere una muestra, llamada muestra Bootstrap, de observaciones del modelo ajustado, el modelo con parámetros \(\widehat\Theta\). La longitud debería ser la misma que el número original de observaciones.

• 2. Estime los parámetros \(\Theta\) por \(\widehat\Theta^*\)para la muestra Bootstrap.
• 3. Repita los pasos (a) y (b) B veces (con B “grande”) y registre los valores \(\widehat\Theta^*\).

La matriz de varianza-covarianza de \(\widehat\Theta\) es entonces estimada por la matriz de varianza—covarianza muestra de los estimados Bootstrap \(\widehat\Theta^*\) (b), \(b=1,2,…,B\):

\[ Var-Cov(\widehat\Theta) = \frac{1}{B-1} \sum_{b=1}^B (\widehat\Theta^*(b) - \widehat\Theta^*(\cdotp))'(\widehat\Theta^*(b) - \widehat\Theta^*(\cdotp)) \\ donde \\ \widehat\Theta^*(\cdotp) = B^{-1} \sum_{b=1}^B \Theta^*(b) \]

El Bootstrap paramétrico requiere código para generar realizaciones de un modelo ajustado. El método de Bootstrap puede ser usado para estimar intervalos de confianza directamente.