Preliminares: mezclas y cadenas de Markov

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]# Set your own path.
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.