Fit Hidden Markov Model

Ajustar un modelo oculto de Markov (HMM) a datos de observación secuenciales usando el algoritmo de expectativa-maximización de Baum-Welch, decodificar la secuencia de estados ocultos más probable vía Viterbi, y seleccionar el número óptimo de estados ocultos mediante criterios de información.

Cuándo Usar

• Observas una secuencia de emisiones pero los estados generativos subyacentes no son directamente observables
• Sospechas que tus datos son generados por un sistema que alterna entre un número finito de regímenes
• Necesitas segmentar una serie temporal en fases latentes (ej., regímenes de mercado, fonemas del habla, anotación de secuencias biológicas)
• Quieres calcular la probabilidad de una secuencia observada bajo un modelo generativo
• Necesitas la secuencia de estados ocultos más probable dadas las observaciones (decodificación)
• Estás comparando modelos con diferentes números de estados ocultos para el mejor compromiso complejidad-ajuste

Entradas

Requerido

observations

n_hidden_states

emission_type

"gaussian"

"discrete"

"poisson"

"multinomial"

Opcional

initial_params

n_restarts

max_iterations

convergence_tol

state_range

[n_hidden_states]

covariance_type

"full"

"full"

"diagonal"

"spherical"

regularization

Procedimiento

Paso 1: Definir Estados Ocultos y Modelo de Observación

1.1. Especificar el número de estados ocultos

K

1.2. Elegir la familia de distribución de emisión basándose en el tipo de datos:

• Datos continuos: Gaussiana (univariante o multivariante)
• Datos de conteo: Poisson o binomial negativa
• Datos categóricos: discreta/multinomial

1.3. Definir los componentes del modelo:

• Matriz de transición

  A

  de tamaño

  K x K

  :

  A[i,j] = P(z_t = j | z_{t-1} = i)

• Parámetros de emisión

  theta_k

  para cada estado

  k

  : específicos de la distribución (ej., media y covarianza para gaussiana)
• Distribución de estado inicial

  pi

  :

  pi[k] = P(z_1 = k)

1.4. Verificar que los datos de observación estén correctamente formateados: sin valores faltantes en la secuencia, dimensionalidad consistente y longitud suficiente relativa al número de parámetros.

Esperado: Una arquitectura HMM claramente especificada con

K

T >> K^2

En caso de fallo: Si los datos contienen valores faltantes, imputar o eliminar los segmentos afectados. Si

T

K

K

Paso 2: Inicializar Parámetros

2.1. Generar parámetros iniciales para cada uno de los

n_restarts

• Matriz de transición: Matriz estocástica aleatoria (cada fila extraída de una distribución de Dirichlet) o una matriz uniforme ligeramente perturbada.
• Parámetros de emisión: Usar clustering K-means sobre las observaciones para inicializar medias; calcular varianzas de clusters para emisiones gaussianas.
• Distribución inicial: Uniforme o proporcional a tamaños de clusters de K-means.

2.2. Para el primer reinicio, usar la inicialización informada por K-means (generalmente el inicio más fuerte). Para reinicios posteriores, usar perturbaciones aleatorias.

2.3. Verificar que todos los parámetros iniciales sean válidos:

• Las filas de la matriz de transición suman 1 con todas las entradas positivas.
• Los parámetros de emisión están en el dominio válido (ej., las matrices de covarianza son definidas positivas).
• La distribución inicial suma 1.

Esperado:

n_restarts

En caso de fallo: Si K-means falla en converger, usar inicialización puramente aleatoria con más reinicios. Si las matrices de covarianza son singulares, agregar la constante de regularización a la diagonal.

Paso 3: Ejecutar Baum-Welch EM para Estimación de Parámetros

3.1. Paso-E (algoritmo Forward-Backward):

• Calcular probabilidades forward

  alpha[t,k]

  = P(o_1,...,o_t, z_t=k | modelo) usando la recursión:

  • alpha[1,k] = pi[k] * b_k(o_1)

  • alpha[t,k] = sum_j(alpha[t-1,j] * A[j,k]) * b_k(o_t)

• Calcular probabilidades backward

  beta[t,k]

  = P(o_{t+1},...,o_T | z_t=k, modelo):

  • beta[T,k] = 1

  • beta[t,k] = sum_j(A[k,j] * b_j(o_{t+1}) * beta[t+1,j])

• Calcular posterior de estado

  gamma[t,k]

  = P(z_t=k | O, modelo):

  • gamma[t,k] = alpha[t,k] * beta[t,k] / P(O | modelo)

• Calcular posterior de transición

  xi[t,i,j]

  = P(z_t=i, z_{t+1}=j | O, modelo).

3.2. Paso-M (Re-estimación de parámetros):

• Actualizar matriz de transición:

  A[i,j] = sum_t(xi[t,i,j]) / sum_t(gamma[t,i])

• Actualizar parámetros de emisión usando estadísticos suficientes ponderados:
  • Media gaussiana:

    mu_k = sum_t(gamma[t,k] * o_t) / sum_t(gamma[t,k])

  • Covarianza gaussiana: matriz de dispersión ponderada más regularización
  • Discreta:

    b_k(v) = sum_t(gamma[t,k] * I(o_t=v)) / sum_t(gamma[t,k])

• Actualizar distribución inicial:

  pi[k] = gamma[1,k]

3.3. Calcular log-verosimilitud:

log P(O | modelo) = log sum_k(alpha[T,k])

3.4. Escalado: Usar variables forward-backward escaladas para prevenir desbordamiento numérico inferior en secuencias largas. Normalizar

alpha

3.5. Repetir Paso-E y Paso-M hasta que el cambio en log-verosimilitud esté por debajo de

convergence_tol

max_iterations

3.6. A través de todos los reinicios, mantener el conjunto de parámetros con la log-verosimilitud final más alta.

Esperado: Log-verosimilitud monótonamente no decreciente a través de las iteraciones, convergiendo dentro de

max_iterations

En caso de fallo: Si la log-verosimilitud decrece, hay un error en el Paso-E o Paso-M — verificar las fórmulas. Si la convergencia es muy lenta, intentar mejor inicialización o aumentar

max_iterations

Paso 4: Aplicar Decodificación de Viterbi para la Secuencia de Estados Más Probable

4.1. Inicializar variables de Viterbi:

• delta[1,k] = log(pi[k]) + log(b_k(o_1))

• psi[1,k] = 0

  (sin predecesor)

4.2. Recurrir hacia adelante para

t = 2,...,T

• delta[t,k] = max_j(delta[t-1,j] + log(A[j,k])) + log(b_k(o_t))

• psi[t,k] = argmax_j(delta[t-1,j] + log(A[j,k]))

4.3. Terminar:

• z*_T = argmax_k(delta[T,k])

• Log-probabilidad de la mejor ruta:

  max_k(delta[T,k])

4.4. Retroceder para

t = T-1,...,1

• z*_t = psi[t+1, z*_{t+1}]

4.5. Generar la secuencia de estados decodificada

z* = (z*_1, ..., z*_T)

4.6. Comparar la probabilidad de la ruta de Viterbi con la probabilidad total de la secuencia del algoritmo forward para evaluar cuán dominante es la mejor ruta.

Esperado: Una única secuencia de estados más probable de longitud

T

{1,...,K}

En caso de fallo: Si la ruta de Viterbi tiene log-probabilidad de infinito negativo, alguna probabilidad de transición o emisión es cero donde no debería serlo. Agregar valores mínimos para prevenir log(0).

Paso 5: Realizar Selección de Modelo (BIC/AIC a Través de Órdenes de Modelo)

5.1. Para cada número candidato de estados ocultos

K

state_range

5.2. Calcular el número de parámetros libres

p

• Matriz de transición:

  K * (K - 1)

  (cada fila es un símplex)
• Parámetros de emisión: depende de la familia (ej., gaussiana con covarianza completa en

  d

  dimensiones:

  K * (d + d*(d+1)/2)

  )
• Distribución inicial:

  K - 1

5.3. Calcular criterios de información:

• BIC = -2 * log_verosimilitud + p * log(T)

• AIC = -2 * log_verosimilitud + 2 * p

• AICc = AIC + 2*p*(p+1) / (T - p - 1)

  (corrección para muestras pequeñas)

5.4. Seleccionar el modelo con el BIC más bajo (preferido por consistencia) o AIC (preferido para predicción). Reportar ambos.

5.5. Tabular resultados: para cada

K

5.6. Si el

K

state_range

Esperado: Un mínimo claro en BIC/AIC que identifique el número óptimo de estados ocultos. El modelo seleccionado debería haber convergido y tener significados de estado interpretables.

En caso de fallo: Si no existe un mínimo claro (BIC monótonamente decreciente), el modelo puede estar mal especificado — considerar una familia de emisión diferente. Si todos los modelos tienen log-verosimilitud pobre, los datos pueden no seguir una estructura HMM.

Paso 6: Validar con Datos Retenidos y Decodificación Posterior

6.1. Dividir datos en conjuntos de entrenamiento y validación (ej., 80/20 o usar múltiples secuencias si están disponibles).

6.2. Ajustar el modelo con datos de entrenamiento. Calcular la log-verosimilitud en datos retenidos usando el algoritmo forward (no re-ajustar parámetros).

6.3. Decodificación posterior (alternativa a Viterbi):

• Para cada paso temporal, asignar el estado con la mayor probabilidad posterior:

  z^_t = argmax_k(gamma[t,k])

• Esto maximiza el número esperado de estados correctamente decodificados (vs. Viterbi que maximiza la probabilidad conjunta de la ruta).

6.4. Comparar decodificación de Viterbi y posterior:

• Calcular la tasa de acuerdo entre las dos secuencias decodificadas.
• Las regiones de desacuerdo indican asignaciones de estado ambiguas.

6.5. Evaluar la interpretabilidad de los estados:

• Examinar los parámetros de emisión para cada estado (medias, varianzas, distribuciones discretas).
• Verificar que los estados correspondan a regímenes significativos en el contexto del dominio.
• Verificar que los tiempos de permanencia en estado (implicados por la diagonal de

  A

  ) sean razonables.

6.6. Calcular la log-verosimilitud de datos retenidos por observación y comparar a través de órdenes de modelo para confirmar la selección de modelo del conjunto de entrenamiento.

Esperado: La log-verosimilitud de datos retenidos es razonablemente cercana a la log-verosimilitud de entrenamiento (sin sobreajuste severo). La decodificación de Viterbi y posterior coinciden en 90%+ de los pasos temporales. Los estados tienen distribuciones de emisión distintas e interpretables.

En caso de fallo: Si la verosimilitud retenida es mucho peor que la de entrenamiento, el modelo está sobreajustando — reducir

K

Validación

• La log-verosimilitud es monótonamente no decreciente a través de las iteraciones de Baum-Welch para cada reinicio
• La matriz de transición es fila-estocástica (filas suman 1, todas las entradas no negativas)
• Los parámetros de emisión están en el dominio válido (covarianzas definidas positivas, distribuciones de probabilidad válidas)
• La log-probabilidad de la ruta de Viterbi no excede la log-probabilidad total de la secuencia
• Las curvas de BIC/AIC muestran un mínimo claro en el orden de modelo seleccionado
• La log-verosimilitud de datos retenidos confirma que el modelo generaliza más allá del conjunto de entrenamiento
• Los cálculos de probabilidad forward y backward concuerdan:

  P(O) = sum_k(alpha[T,k]) = sum_k(pi[k] * b_k(o_1) * beta[1,k])

Errores Comunes

• Óptimos locales en EM: El algoritmo de Baum-Welch converge a un máximo local, no necesariamente al global. Siempre usar múltiples reinicios aleatorios y seleccionar el mejor.
• Desbordamiento numérico inferior: Las probabilidades forward-backward se reducen exponencialmente con la longitud de la secuencia. Usar cómputo en espacio logarítmico o variables escaladas para prevenir desbordamiento inferior a cero.
• Sobreajuste con demasiados estados: Cada estado oculto adicional agrega

  O(K + d^2)

  parámetros. Usar BIC (no solo verosimilitud) para selección de modelo y validar con datos retenidos.
• Intercambio de etiquetas: Los estados ocultos son identificables solo hasta permutación. Al comparar modelos entre reinicios, emparejar estados por parámetros de emisión, no por índice.
• Estados degenerados: Un estado puede colapsar para explicar una única observación (gaussiana con varianza cercana a cero). La regularización en matrices de covarianza previene esto.
• Confundir Viterbi y decodificación posterior: Viterbi da la única mejor ruta conjunta; la decodificación posterior da el mejor estado marginal en cada paso temporal. Responden preguntas diferentes y pueden diferir significativamente.
• Ignorar tiempos de permanencia en estado: La distribución geométrica de tiempo de permanencia implícita en los HMM estándar puede ser un ajuste pobre para datos con duraciones de régimen largas. Considerar modelos ocultos semi-Markov si los tiempos de permanencia no son geométricos.

Habilidades Relacionadas

• Model Markov Chain -- prerrequisito para entender la estructura de transición que subyace a la capa oculta
• Simulate Stochastic Process -- puede usarse para generar datos HMM sintéticos para pruebas y para simular a partir de un modelo ajustado para verificaciones predictivas posteriores