12 - Datos secuenciales y modelos markovianos de estados ocultos
## Warning: package 'ggplot2' was built under R version 3.1.3Cuando los datos están ordenadas en sucesiones (temporales o espaciales, por ejemplo), usualmente encontramos estructuras de dependencia dadas por la proximidad. Por ejemplo, en series de tiempo observaciones de un tiempo dado están correlacionadas con observaciones en tiempos cercanos.
Una estrategia para modelar este tipo de datos es considerar modelos que explícitamente tratan con esta dependencia temporal entre las observaciones (por ejemplo, modelos ARMA). El enfoque que consideramos aquí es uno distinto que se basa en variables latentes.
En el caso de series de tiempo, estos modelos de variables latentes se llaman modelos de espacio de estado, cuyos dos ejemplos más populares son:
Modelos markovianos de estados escondidos (Hidden Markov Models): cuando las variables latentes son discretas. Las observadas pueden ser numéricas o discretas.
Modelos dinámicos lineales: cuando las variables latentes son gaussianas. Las observadas generalmente son numéricas.
Ejemplo: serie de terremotos. La siguiente serie cuenta el número de temblores de magnitud 7 o más en el mundo, desde 1900 a 2006 (Hidden Markov Models for Time Series, de Zucchini y MacDonald):
terremotos <- data.frame(num = c(13,14,8,10,16,26,32,27,18,32,36,24,22,23,22,18,
25,21,21,14,8,11,14,23,18,17,19,20,22,19,13,26,13,14,22,24,21,22,26,21,23,24,
27,41,31,27,35,26,28,36,39,21,17,22,17,19,15,34,10,15,22,18,15,20,15,22,19,16,
30,27,29,23,20,16,21,21,25,16,18,15,18,14,10,15,8,15,6,11,8,7,18,16,13,12,13,
20,15,16,12,18,15,16,13,15,16,11,11), anio = 1900:2006)
ggplot(terremotos, aes(anio, num)) + geom_line()
Es natural intentar modelar estos conteos con una distribución Poisson. Sin embargo, un modelo Poisson para observaciones independientes no ajusta a los datos:
library(ggplot2)
library(plyr)
library(dplyr)##
## Attaching package: 'dplyr'
##
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(nullabor)
set.seed(12875498)
lambda_1 <- mean(terremotos$num)
# pruebas de hipótesis visuales
terremotos_null <- lineup(null_dist('num', dist = 'poisson',
params = list(lambda = lambda_1)),
n = 20, terremotos)## decrypt("n20q ilOl TJ IH4TOTHJ 5m")ggplot(terremotos_null, aes(x = num)) +
geom_histogram(binwidth = 3) +
facet_wrap(~.sample)
Vemos que los datos presentan sobredispersión en relación a la distribución del modelo. Para lidiar con esta sobredispersión, podemos usar un modelo de clases latentes.
Consideramos un modelo \(C\to X\), donde \(C\) es la clase latente y \(X\) es la observación. La observación la suponemos Poison con una media dependiendo de la clase. Ajustamos por EM y seleccionamos un modelo:
library(flexmix)## Loading required package: latticemodel_1 <- FLXMRglm(family = "poisson")
# usa algoritmo EM para ajustar (5 repeticiones con distintos arranques):
# número de clases latentes: de 1 a 6.
modelos_fit <- initFlexmix(num ~ 1, data = terremotos, k = 1:6, model = model_1,
nrep = 5)## 1 : * * * * *
## 2 : * * * * *
## 3 : * * * * *
## 4 : * * * * *
## 5 : * * * * *
## 6 : * * * * *
## * * * * * *modelos_fit##
## Call:
## initFlexmix(num ~ 1, data = terremotos, model = model_1,
## k = 1:6, nrep = 5)
##
## iter converged k k0 logLik AIC BIC ICL
## 1 2 TRUE 1 1 -392 786 789 789
## 2 34 TRUE 2 2 -360 727 735 767
## 3 43 TRUE 3 3 -357 724 737 794
## 4 37 TRUE 4 4 -357 728 746 878
## 5 33 TRUE 5 5 -357 732 756 930
## 6 32 TRUE 6 6 -357 736 765 973Así que podríamos escoger, por ejemplo, el modelo de tres clases:
mod_1 <- getModel(modelos_fit, which='3')
mod_1##
## Call:
## initFlexmix(num ~ 1, data = terremotos, model = model_1,
## k = 3, nrep = 5)
##
## Cluster sizes:
## 1 2 3
## 26 11 70
##
## convergence after 43 iterationsprobs_clase <- mod_1@prior
medias_pois <- mod_1@components
probs_clase## [1] 0.30 0.13 0.58medias_pois## $Comp.1
## $Comp.1[[1]]
## $coef
## (Intercept)
## 2.6
##
##
##
## $Comp.2
## $Comp.2[[1]]
## $coef
## (Intercept)
## 3.5
##
##
##
## $Comp.3
## $Comp.3[[1]]
## $coef
## (Intercept)
## 3Para entender cómo quedó el ajuste, podemos graficar la serie original junto con la media poisson condicional de cada clase, por ejemplo:
terremotos$clase_1 <- clusters(mod_1)
terremotos$medias_est <- sapply(terremotos$clase_1, function(cl){
exp(mod_1@components[[cl]][[1]]@parameters$coef[1])
})
ggplot(terremotos, aes(anio, num)) +
geom_line(alpha = 0.8) +
geom_line(aes(anio, medias_est), color = "red", alpha= 0.8)
Simulamos de estos datos y vemos que ya no hay sobredispersión (es decir, este modelo ajusta mejor que el de una sola clase):
terremotos_null <- lineup(null_dist('num', dist = 'poisson',
params = list(lambda = terremotos$medias_est)),
n = 20, terremotos)## decrypt("n20q ilOl TJ IH4TOTHJ 1")ggplot(terremotos_null, aes(x = num)) +
geom_histogram(binwidth = 3) +
facet_wrap(~.sample)
Si una variable \(Y\) es Poisson, su media es igual a su varianza. Ahora supón que \(Y\) es una mezcla discreta Poisson: si la clase latente la denotamos por \(S\) y la observación por \(Y\), entonces \(Y|S=s\) es Poisson con media \(\lambda_s\), para \(s=1\ldots,M\). Muestra que en general la varianza de \(Y\) es más grande que su media (está sobredispersa). Explica en palabras de dónde viene esa sobredispersión.
Dependencia temporal
Existe un aspecto adicional que no hemos considerado: dependencia de observaciones contiguas en las serie. Podríamos gráficar para la serie de terremotos las pares \((X_t, X_{t+1})\):
ggplot(terremotos, aes(lag(num), num)) +
geom_point() +
geom_smooth()## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (geom_point).
Esta gráfica explica por qué si hacemos simulaciones de nuestro modelo ajustado (tomando en cuenta el orden temporal), estas simulaciones son claramente diferentes que los datos observados:
ggplot(terremotos_null, aes(anio, num)) +
geom_line(size=0.3) +
facet_wrap(~.sample)
Si nos interesa hacer predicciones, el modelo de mezcla no es un buen camino, pues no utiliza adecuadamente la información más reciente para hacer los pronósticos.
Otra manera de verificar que las observaciones tienen mayor grado de correlación temporal que los simulados del modelo de mezcla es a través de los siguientes cálculos.
terremotos_null %>%
group_by(.sample) %>%
mutate(num_lag = lag(num)) %>%
filter(!is.na(num_lag)) %>%
summarise(cor = cor(num, num_lag)) %>%
arrange(desc(cor))## Source: local data frame [20 x 2]
##
## .sample cor
## 1 4 0.576
## 2 8 0.407
## 3 14 0.378
## 4 18 0.312
## 5 7 0.312
## 6 10 0.304
## 7 15 0.251
## 8 12 0.246
## 9 16 0.237
## 10 9 0.227
## 11 19 0.221
## 12 2 0.211
## 13 13 0.204
## 14 11 0.175
## 15 5 0.146
## 16 17 0.117
## 17 20 0.099
## 18 1 0.084
## 19 3 0.078
## 20 6 0.032Y notamos que justamente los datos con mayor correlación corresponden a los datos observados (.sample = decrypt("n20q ilOl TJ IH4TOTHJ 1")).
Ahora evaluamos la tasa de error de predicción a un paso de este modelo usando validación cruzada para datos temporales. En este caso, la predicción es simplemente la media Poisson de la clase más popular:
preds_1 <- sapply(50:106 ,function(i){
mod_1 <- flexmix(num~1, data = terremotos[1:i, , drop = FALSE], k = 3,
model = model_1)
mod_1
clase_1 <- which.max(table(mod_1@cluster))
mod_1@components[[clase_1]][[1]]@parameters[[1]]
})Error medio absoluto:
mean(abs(exp(preds_1) - terremotos[51:107, "num"]))## [1] 5.7Autocorrelación muestral
Una manera popular de diagnosticar observaciones temporalmente correlacionadas es usando la función de autocorrelación muestral. Usualmente la correlación la calculamos para pares de observaciones \((X_i,Y_i)\). En el caso de la autocorrelación, calculamos la correlación entre una serie y versiones rezagadas de la misma variable, es decir, consideramos los pares \((X_i,X_{i-k})\) para distintos valores de \(k=1,2,\ldots\). La autcorrelación muestral se define:
\[r_k=\frac{\sum_{t=k+1}^{T}(x_t-\overline{x})(x_{t-k}-\overline{x})}{\sum_{t=1}^T(x_t-\overline{x})},\]
podemos graficar la serie \(r_k\) para descubrir estructura en la serie que estamos analizando. En nuestro ejemplo:
acf(terremotos$num)
El primer valor (\(k=1\)) siempre es igual a 1. Observamos sin embargo autocorrelacciones considerables de orden 1 a 5.
¿Cómo se ve esta gráfica cuando tenemos observaciones independientes? Esperamos observar coeficientes de correlación relativamente chicos:
set.seed(280572)
par(mfrow=c(1,2))
acf(rnorm(150))
acf(rpois(150, 20))
Modelando la dependencia temporal
Si queremos capturar la estrucutra temporal de los datos, o hacer predicciones para datos futuros, es necesario modelar la estructura temporal explícita o implícitamente. Por ejemplo,podríamos intentar construir el modelo:
O de manera más general, tomando en cuanto dependencias más largas:
Sin embargo, utilizaremos un enfoque de variables latentes. La idea es introducir un estado latente, a partir del cual generamos las observaciones. Por ejemplo:
En un modelo como este, pueden existir dependencias más largas entre las \(X_t\), aún cuando la estructura es relativamente simple.
¿cuáles son las independencias condicionales mostradas en la gráfica? ¿Las \(X_i\) son independientes entre ellas?}
Modelos markovianos de estados ocultos
Comenzamos con un modelo simple, que consiste de dos variables, una observada y una latente:
- \(X_t\) cantidad de lluvia en mm.
- \(S_t\) variable discreta con dos valores posibles: soleado o nublado.
El modelo está definido en dos grandes partes
Modelo de observaciones: \[P(X_t|S_t).\] Este modelo está dado con dos densidades: una para la cantidad de lluvia cuando \(S_t=nublado\) y otra para \(S_t=soleado\).
Modelo de transición de estados: \[P(S_{t+1}|S_t).\] Estos modelos están dados por una matriz de transición, pues supondremos que estas probabildades no cambian con \(t\) (supuesto de homogeneidad). Por ejemplo, podríamos tener:
p_trans <- matrix(c(0.8, 0.2, 0.3, 0.7), byrow = T, ncol = 2)
p_trans## [,1] [,2]
## [1,] 0.8 0.2
## [2,] 0.3 0.7Y el modelo completo se puede representar como la última gráfica de la sección anterior.
En este ejmplo, el modelo de observaciones es \(X_t\) normal con media 20 y desviación estándar 5 cuando se trata de un día soleado, y \(X_t\) es normal con media 100 y desviación estándar 20 cuando se trata de un día nublado.
set.seed(2805)
estadosLluvia <- function(n){
estado <- character(n)
estado[1] <- sample(c('soleado', 'nublado'), 1)
transicion <- matrix(c(0.8, 0.2, 0.2, 0.8), 2, byrow = T)
rownames(transicion) <- c('soleado', 'nublado')
colnames(transicion) <- c('soleado', 'nublado')
for(j in 2:n){
estado[j] <- sample(c('soleado', 'nublado'), 1,
prob = transicion[estado[j - 1], ])
}
estado
}
estadosLluvia(10)## [1] "soleado" "soleado" "soleado" "soleado" "soleado" "soleado" "soleado"
## [8] "nublado" "nublado" "nublado"prob_obs <- data.frame(edo = c("soleado", "nublado"), medias = c(20, 100),
desv = c(10, 50))
obsPrecip <- function(long_serie){
data.frame(dia = 1:long_serie, edo = estadosLluvia(long_serie)) %>%
left_join(prob_obs) %>%
group_by(dia) %>%
mutate(precipit = max(rnorm(1, medias, desv), 0))
}
ej_1 <- obsPrecip(10)## Warning: joining factors with different levels, coercing to character
## vectorSeries simuladas de este modelo se ven como sigue:
sims_1 <- rdply(10, obsPrecip(50))
ggplot(sims_1, aes(x = dia, ymax = precipit, ymin = 0, y = precipit, color = edo)) +
geom_linerange(alpha = 0.6) +
geom_point(size=0.8) +
facet_wrap(~.n, nrow = 2)
Nótese que hay dependencia secuencial en estos datos. Si sólo tuviéramos los datos observados, podríamos construir un modelo más complejo intentando entender las correlaciones secuenciales:
ej_1 <- obsPrecip(1000)## Joining by: "edo"acf(ej_1$precipit)
Ajusta con EM un modelo a estos datos usando el paquete depmixS4. Prueba con distintos tamaños de muestra. Por ejemplo:
library(depmixS4)
datos <- sims_1
# ntimes: son 10 series de longitud 50 cada una:
mod_hmm <- depmix(precipit~1, data=datos, nstates=2,
family=gaussian(), ntimes=rep(50,10))
fit_1 <- fit(mod_hmm, emcontrol=em.control(maxit=200))## iteration 0 logLik: -2683
## iteration 5 logLik: -2539
## iteration 10 logLik: -2441
## iteration 15 logLik: -2441
## converged at iteration 19 with logLik: -2441summary(fit_1)## Initial state probabilties model
## pr1 pr2
## 0.56 0.44
##
##
## Transition matrix
## toS1 toS2
## fromS1 0.79 0.21
## fromS2 0.23 0.77
##
## Response parameters
## Resp 1 : gaussian
## Re1.(Intercept) Re1.sd
## St1 20 10
## St2 101 46Observa que como tenemos varias series, podemos estimar las probabilidades iniciales. Cuando tenemos una sola serie, típicamente suponemos que la serie comienza en un estado, por ejemplo el primero.
Consideramos datos secuenciales \(X_1,\ldots,X_T\). Un modelo markoviano de estado oculto se factoriza como (es decir, cumple las relaciones de independencia condicional indicadas en la gráfica):
Las variables \(S_t\) (estado oculto) son latentes.
Las variables \(X_t\) no son independientes, pero son condicionalmente independientes dados los estados ocultos.
El modelo está dado por una factorización que solo incluye términos de la forma \(P(S_t|S_{t-1})\), \(P(X_t|S_t)\) y \(P(S_1)\).
En los modelos de transición de estados suponemos homogeneidad de la cadena de Markov subyacente, es decir,\[P(S_t=j|S_{t-1}=i)=p_{ij}\]no depende de \(t\).
- En los modelos de respuesta o de observación, las variables observadas \(X_t\) pueden ser discretas o continuas. Si \(X_t\) son discretas, cada \(P(X_t|S_t)\) está dada por una tabla, que también suponemos homogénea (no depende de \(t\)). Si \(X_t\) son continuas, entonces estas probabilidades están dadas por densidades condicionales.
Observación: una generalización usual de homogeneidad es incluir covariables en los modelos de respuesta o de observación, de forma que especificamos \(P(X_t|S_t,Z_t )\) con modelos de regresión, por ejemplo. El análisis es condicional a los valores de \(Z_t\) (no modelamos su dinámica, por ejemplo). También las probabilidades de transición pueden depender de covariables.}
Ejemplo: serie de terremotos. Ahora podemos regresar al ejemplo de los terremotos. Ajustamos los modelos de estado oculto a la serie (de 1 a 4 estados):
library(depmixS4)
terremotos <- data.frame(num = c(13,14,8,10,16,26,32,27,18,32,36,24,22,23,22,18,
25,21,21,14,8,11,14,23,18,17,19,20,22,19,13,26,13,14,22,24,21,22,26,21,23,24,
27,41,31,27,35,26,28,36,39,21,17,22,17,19,15,34,10,15,22,18,15,20,15,22,19,16,
30,27,29,23,20,16,21,21,25,16,18,15,18,14,10,15,8,15,6,11,8,7,18,16,13,12,13,
20,15,16,12,18,15,16,13,15,16,11,11), anio = 1900:2006)
head(terremotos)## num anio
## 1 13 1900
## 2 14 1901
## 3 8 1902
## 4 10 1903
## 5 16 1904
## 6 26 1905modelos_hmm <- lapply(1:4, function(i){
mod_hmm <- depmix(num~1, data = terremotos, nstates = i, family = poisson(),
ntimes = 107)
fit_hmm <- fit(mod_hmm, emcontrol = em.control(maxit = 200))
fit_hmm
})## iteration 0 logLik: -392
## converged at iteration 1 with logLik: -392
## iteration 0 logLik: -389
## iteration 5 logLik: -343
## iteration 10 logLik: -342
## iteration 15 logLik: -342
## iteration 20 logLik: -342
## iteration 25 logLik: -342
## converged at iteration 27 with logLik: -342
## iteration 0 logLik: -366
## iteration 5 logLik: -342
## iteration 10 logLik: -336
## iteration 15 logLik: -329
## iteration 20 logLik: -329
## iteration 25 logLik: -329
## converged at iteration 29 with logLik: -329
## iteration 0 logLik: -390
## iteration 5 logLik: -336
## iteration 10 logLik: -329
## iteration 15 logLik: -328
## iteration 20 logLik: -328
## iteration 25 logLik: -328
## iteration 30 logLik: -328
## iteration 35 logLik: -328
## iteration 40 logLik: -328
## iteration 45 logLik: -328
## iteration 50 logLik: -328
## iteration 55 logLik: -328
## iteration 60 logLik: -328
## iteration 65 logLik: -328
## iteration 70 logLik: -328
## iteration 75 logLik: -328
## iteration 80 logLik: -328
## iteration 85 logLik: -328
## iteration 90 logLik: -328
## iteration 95 logLik: -328
## iteration 100 logLik: -328
## iteration 105 logLik: -328
## iteration 110 logLik: -328
## iteration 115 logLik: -328
## iteration 120 logLik: -328
## iteration 125 logLik: -328
## iteration 130 logLik: -328
## iteration 135 logLik: -328
## iteration 140 logLik: -328
## iteration 145 logLik: -328
## iteration 150 logLik: -328
## iteration 155 logLik: -328
## iteration 160 logLik: -328
## iteration 165 logLik: -328
## iteration 170 logLik: -328
## iteration 175 logLik: -328
## iteration 180 logLik: -328
## iteration 185 logLik: -328
## iteration 190 logLik: -328
## iteration 195 logLik: -328
## iteration 200 logLik: -328sapply(modelos_hmm, BIC)## [1] 789 707 708 744sapply(modelos_hmm, AIC)## [1] 786 694 679 693Podemos seleccionar el modelo con 3 estados.
mod_1 <- modelos_hmm[[3]]
summary(mod_1)## Initial state probabilties model
## pr1 pr2 pr3
## 0 1 0
##
##
## Transition matrix
## toS1 toS2 toS3
## fromS1 0.906 4.0e-02 0.053
## fromS2 0.032 9.4e-01 0.029
## fromS3 0.190 1.1e-11 0.810
##
## Response parameters
## Resp 1 : poisson
## Re1.(Intercept)
## St1 3.0
## St2 2.6
## St3 3.4Simulaciones de este modelo dan los siguiente resultados, en donde es difícil distinguir los datos observados:
class(mod_1) <- 'depmix'
set.seed(2805)
sims_1 <- simulate(mod_1, nsim = 19)
terremotos_null <- data.frame(num = sims_1@response[[1]][[1]]@y,
anio = rep(terremotos$anio, 19)) %>%
rbind(terremotos)
codigo <- sample(1:20)
terremotos_null$tipo <- rep(codigo, each = 107)
ggplot(terremotos_null, aes(x = anio, y = num)) +
geom_line(size=0.3) +
facet_wrap(~tipo)
# verdaderas están en: codigo[20]Vemos que las autocorrelaciones de grado 1 son similares:
terremotos_null %>%
group_by(tipo) %>%
mutate(num_lag = lag(num)) %>%
filter(!is.na(num_lag)) %>%
summarise(cor = cor(num, num_lag))## Source: local data frame [20 x 2]
##
## tipo cor
## 1 1 0.38
## 2 2 0.56
## 3 3 0.58
## 4 4 0.52
## 5 5 0.47
## 6 6 0.51
## 7 7 0.35
## 8 8 0.50
## 9 9 0.35
## 10 10 0.63
## 11 11 0.44
## 12 12 0.38
## 13 13 0.49
## 14 14 0.43
## 15 15 0.57
## 16 16 0.64
## 17 17 0.57
## 18 18 0.47
## 19 19 0.57
## 20 20 0.50E incluso podemos comparar las gráficas de autocorrelación, que no muestran desajuste importante:
set.seed(227342)
# solo 10 gráficas
sims_1 <- simulate(mod_1, nsim = 9)
terremotos_null <- data.frame(num = sims_1@response[[1]][[1]]@y,
anio = rep(terremotos$anio, 9)) %>%
rbind(terremotos)
codigo <- sample(1:10)
terremotos_null$tipo <- rep(codigo, each = 107)
par(mfrow=c(2, 5))
sims_x <- ddply(terremotos_null, 'tipo', function(df){
x <- df$num
acf(x, lag.max=12, main="")
x
})
Podemos ver el desempeño de estos modelos en predicción:
set.seed(72)
preds_2 <- sapply(50:106 ,function(i){
dat.1 <- terremotos[1:i,,drop=FALSE]
mod.1 <- depmix(num~1, data=dat.1, nstates=3,
family=poisson(), ntimes=i)
fit.mod.1 <- fit(mod.1, verbose = FALSE, emcontrol=em.control(maxit=600))
probs.1 <- fit.mod.1@posterior[i,2:4]
estado.i <- sample(1:3, 1, prob=probs.1)
estado.pred <- which.max(predict(fit.mod.1@transition[[estado.i]]))
fit.mod.1@response[[estado.pred]][[1]]@parameters[1]$coefficients
})## converged at iteration 36 with logLik: -155
## converged at iteration 152 with logLik: -159
## converged at iteration 65 with logLik: -163
## converged at iteration 69 with logLik: -167
## converged at iteration 81 with logLik: -169
## converged at iteration 110 with logLik: -172
## converged at iteration 65 with logLik: -175
## converged at iteration 45 with logLik: -178
## converged at iteration 119 with logLik: -183
## converged at iteration 67 with logLik: -189
## converged at iteration 66 with logLik: -191
## converged at iteration 42 with logLik: -195
## converged at iteration 42 with logLik: -197
## converged at iteration 59 with logLik: -200
## converged at iteration 57 with logLik: -203
## converged at iteration 71 with logLik: -206
## converged at iteration 71 with logLik: -209
## converged at iteration 92 with logLik: -211
## converged at iteration 113 with logLik: -214
## converged at iteration 64 with logLik: -218
## converged at iteration 73 with logLik: -222
## converged at iteration 69 with logLik: -225
## converged at iteration 95 with logLik: -228
## converged at iteration 87 with logLik: -231
## converged at iteration 94 with logLik: -234
## converged at iteration 100 with logLik: -236
## converged at iteration 130 with logLik: -239
## converged at iteration 113 with logLik: -242
## converged at iteration 89 with logLik: -245
## converged at iteration 68 with logLik: -248
## converged at iteration 76 with logLik: -251
## converged at iteration 60 with logLik: -253
## converged at iteration 49 with logLik: -256
## converged at iteration 45 with logLik: -260
## converged at iteration 44 with logLik: -263
## converged at iteration 37 with logLik: -267
## converged at iteration 35 with logLik: -269
## converged at iteration 34 with logLik: -273
## converged at iteration 29 with logLik: -276
## converged at iteration 31 with logLik: -278
## converged at iteration 34 with logLik: -281
## converged at iteration 46 with logLik: -285
## converged at iteration 50 with logLik: -288
## converged at iteration 35 with logLik: -291
## converged at iteration 26 with logLik: -293
## converged at iteration 27 with logLik: -296
## converged at iteration 34 with logLik: -300
## converged at iteration 34 with logLik: -302
## converged at iteration 33 with logLik: -305
## converged at iteration 26 with logLik: -308
## converged at iteration 28 with logLik: -311
## converged at iteration 27 with logLik: -314
## converged at iteration 26 with logLik: -316
## converged at iteration 29 with logLik: -319
## converged at iteration 26 with logLik: -321
## converged at iteration 25 with logLik: -324
## converged at iteration 26 with logLik: -326## error medio absoluto
mean(abs(exp(preds_2) - terremotos[51:107,"num"]))## [1] 4.4Finalmente, comparamos con la predicción de tomar el valor anterior:
mean(abs(terremotos[51:107, "num"] - terremotos[50:106, "num"]))## [1] 5Ejemplo: prototipo reconocimiento de trazo de dígitos. Los modelos de estados escondidos se usan en varias partes de procesamiento de señales. En este ejemplo, veremos un enfoque para reconocer dígitos con base en los trazos que se hacen para escribirlos, este es un enfoque distinto al usual en donde se trabaja con una matriz de pixeles. Utilizaremos los datos pendigits.
load('./pendigits/digitos_lista.RData')
load('./pendigits/digitos.RData')
length(results.list)## [1] 7494digito <- results.list[[1]]
digito## [,1] [,2]
## init 267 333
## 267 336
## 267 339
## 266 342
## 264 344
## 262 346
## 259 346
## 254 347
## 249 346
## 244 345
## 239 343
## 235 340
## 233 333
## 233 327
## 235 318
## 239 309
## 245 297
## 251 285
## 258 272
## 265 259
## 272 247
## 277 235
## 280 223
## 284 208
## 284 199
## 281 192
## 278 184
## 272 178
## 265 171
## 258 166
## 249 162
## 241 158
## 233 155
## 225 153
## 219 152
## 214 152
## 210 153
## 207 155
## 204 158
## 202 163
## 201 169
## 201 174
## 202 181
## 205 188
## 208 195
## 213 202
## 220 209
## 227 216
## 236 223
## 246 229
## 258 237
## 268 245
## 279 252
## 290 260
## 301 268
## 310 277
## 317 286
## 323 293
## 327 300
## 329 306
## 330 312
## 330 317
## 330 320
## 329 323
## 327 325
## 325 328
## 322 330
## 317 333
## 313 336
## 306 337
## 300 338
## 292 339
## 285 340
## 277 339
## 272 338
## 265 336
## 262 333
## 259 330floor(nrow(digito)/9)## [1] 8par(mfrow=c(3,3))
for(i in 1:9){
plot(digito[1:(8*(i-1)+2),])
}
Nos interesa construir un modelo que describa la escritura de estos dígitos. Proponemos un modelo como sigue:
El espacio de estados ocultos tiene 16 estados distintos.
La matriz de transición entre estados es de izquierda a derecha. Es decir, la serie comienza en el estado 1, y solamente transiciona de 1 a 2, de 2 a 3 y así sucesivamente. Puede permanecer un número arbitrario de tiempos en cada estado.
Las observaciones son bivariadas: tenemos un ángulo de trazo (8 posibles direcciones), y una velocidad de trazo (que supondremos gaussiano). Supondremos estas dos variables son condicionalmente independientes dado el estado.
La idea general es:
Transformar los patrones de escritura a observaciones de ángulo y velocidad.
Ajustar un modelo distinto para cada uno de los dígitos.
Cuando observamos un dígito nuevo, calculamos la verosimilitud de la observación bajo los distintos modelos.
Clasificamos al dígito con verosimilitud más alta.
En estas notas veremos los primeros dos pasos.
La ventaja de usar HMM en este contexto es, en primer lugar, que no es necesario normalizar los dígitos a una longitud fija (aunque esto puede ayudar en la predicción), y mejor aún, que el modelo es flexible en cuanto a escalamientos de distintas partes de la escritura de los dígitos, por ejemplo, dibujar más lentamente unas partes, hacer un ocho con la parte de arriba más grande que la de abajo, etc.
Primero construimos funciones auxiliares:
# convertirSerie: recibe las coordenadas y devuelve el ángulo y la velocidad
convertirSerie <- function(patron){
d_1 <- rbind(patron[-1, ], c(NA,NA)) - patron
d_2 <- d_1[-nrow(d_1), ]
angulo <- acos(d_2[, 1] / (apply(d_2, 1, function(x) sqrt(sum(x ^ 2))))) *
sign(d_2[, 2])
velocidad <- apply(d_2, 1, function(x) sqrt(sum(x ^ 2)))
dat_out<- data.frame(angulo = angulo[], velocidad = velocidad[])
dat_out[!is.na(dat_out$angulo), ]
}
# trazar: recibe ángulo y velocidad y devuelve coordenadas
trazar <- function(datos){
angulos <- datos[, 1]
velocidad <- datos[, 2]
longitud <- length(angulos) + 1
puntos <- data.frame(x = rep(NA, longitud), y = rep(NA, longitud))
puntos[1, ] <- c(0, 0)
for(i in 2:length(angulos)){
puntos[i,] <- puntos[i - 1, ] + velocidad[i - 1] *
c(cos(angulos[i - 1]), sin(angulos[i - 1]))
}
puntos
}
con_ang <- convertirSerie(results.list[[6]])
con_ang## angulo velocidad
## 1 0.695 7.8
## 2 0.675 6.4
## 3 1.571 4.0
## 4 0.000 2.0
## 5 -2.521 8.6
## 6 -2.467 12.8
## 7 -2.356 17.0
## 8 -2.260 22.0
## 9 -2.196 22.2
## 10 -2.185 20.8
## 11 -2.013 21.0
## 12 -1.930 17.1
## 13 -1.642 14.0
## 14 -1.305 11.4
## 15 -0.559 9.4
## 16 -0.091 11.0
## 17 0.000 11.0
## 18 0.180 11.2
## 19 0.464 8.9
## 20 0.785 8.5
## 21 1.373 5.1
## 22 2.246 6.4
## 23 2.622 8.1
## 24 2.923 9.2
## 25 3.051 11.0
## 26 -2.863 7.3
## 27 -2.976 6.1plot(trazar(con_ang))
plot(results.list[[6]])
Comenzamos seleccionando los dígitos 8, convertimos los datos a una tabla, y categorizamos los ángulos:
library(Hmisc)
filtro <- sapply(digitos, function(x) x==8)
resultados <- results.list[filtro]
dat_8 <- ldply(1:700, function(i){
conv_1 <- convertirSerie(resultados[[i]])
data.frame(angulo = conv_1[,1], velocidad = conv_1[,2], serie=i,
tiempo=1:nrow(conv_1))
})
head(dat_8)## angulo velocidad serie tiempo
## 1 1.6 3.0 1 1
## 2 1.6 3.0 1 2
## 3 1.9 3.2 1 3
## 4 2.4 2.8 1 4
## 5 2.4 2.8 1 5
## 6 0.0 3.0 1 6# Longitud de los dígitos
longitudes <- ddply(dat_8, 'serie', summarise, long=max(tiempo))
# Creamos 8 categorías usando cuantiles
dat_8$angulo_cat <- factor(cut2(dat_8$angulo, g = 8, levels.mean = TRUE))
levs <- as.numeric(levels(dat_8$angulo_cat))Ahora construimos nuestro modelo. Una diferencia importante en relación a los ejempos anteriores es que aquí tenemos varias series de tiempo observadas (cuyas longitudes están en el argumento ntimes abajo):
mod <- depmix(list(angulo_cat ~ 1, velocidad ~ 1), data = dat_8,
nstates = 16,
family = list(multinomial("identity"), gaussian()),
ntimes = longitudes$long)El modelo está incializado de la siguiente forma
summary(mod)## Initial state probabilties model
## pr1 pr2 pr3 pr4 pr5 pr6 pr7 pr8 pr9 pr10 pr11 pr12
## 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## pr13 pr14 pr15 pr16
## 0.062 0.062 0.062 0.062
##
##
## Transition matrix
## toS1 toS2 toS3 toS4 toS5 toS6 toS7 toS8 toS9 toS10 toS11
## fromS1 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS2 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS3 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS4 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS5 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS6 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS7 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS8 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS9 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS10 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS11 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS12 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS13 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS14 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS15 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## fromS16 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062 0.062
## toS12 toS13 toS14 toS15 toS16
## fromS1 0.062 0.062 0.062 0.062 0.062
## fromS2 0.062 0.062 0.062 0.062 0.062
## fromS3 0.062 0.062 0.062 0.062 0.062
## fromS4 0.062 0.062 0.062 0.062 0.062
## fromS5 0.062 0.062 0.062 0.062 0.062
## fromS6 0.062 0.062 0.062 0.062 0.062
## fromS7 0.062 0.062 0.062 0.062 0.062
## fromS8 0.062 0.062 0.062 0.062 0.062
## fromS9 0.062 0.062 0.062 0.062 0.062
## fromS10 0.062 0.062 0.062 0.062 0.062
## fromS11 0.062 0.062 0.062 0.062 0.062
## fromS12 0.062 0.062 0.062 0.062 0.062
## fromS13 0.062 0.062 0.062 0.062 0.062
## fromS14 0.062 0.062 0.062 0.062 0.062
## fromS15 0.062 0.062 0.062 0.062 0.062
## fromS16 0.062 0.062 0.062 0.062 0.062
##
## Response parameters
## Resp 1 : multinomial
## Resp 2 : gaussian
## Re1.pr1 Re1.pr2 Re1.pr3 Re1.pr4 Re1.pr5 Re1.pr6 Re1.pr7 Re1.pr8
## St1 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St2 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St3 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St4 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St5 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St6 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St7 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St8 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St9 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St10 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St11 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St13 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St14 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St15 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St16 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## Re2.(Intercept) Re2.sd
## St1 0 1
## St2 0 1
## St3 0 1
## St4 0 1
## St5 0 1
## St6 0 1
## St7 0 1
## St8 0 1
## St9 0 1
## St10 0 1
## St11 0 1
## St12 0 1
## St13 0 1
## St14 0 1
## St15 0 1
## St16 0 1pero antes de ajustarlo, es necesario restringir los parámetros como describimos arriba (para tener un modelo de izquierda a derecha)
# usamos setpars para ver el orden de los parámetros y poder escribir las
# restricciones
setpars(mod, value=1:npar(mod))## Initial state probabilties model
## pr1 pr2 pr3 pr4 pr5 pr6 pr7 pr8 pr9 pr10 pr11 pr12 pr13 pr14 pr15
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## pr16
## 16
##
##
## Transition matrix
## toS1 toS2 toS3 toS4 toS5 toS6 toS7 toS8 toS9 toS10 toS11 toS12
## fromS1 17 18 19 20 21 22 23 24 25 26 27 28
## fromS2 33 34 35 36 37 38 39 40 41 42 43 44
## fromS3 49 50 51 52 53 54 55 56 57 58 59 60
## fromS4 65 66 67 68 69 70 71 72 73 74 75 76
## fromS5 81 82 83 84 85 86 87 88 89 90 91 92
## fromS6 97 98 99 100 101 102 103 104 105 106 107 108
## fromS7 113 114 115 116 117 118 119 120 121 122 123 124
## fromS8 129 130 131 132 133 134 135 136 137 138 139 140
## fromS9 145 146 147 148 149 150 151 152 153 154 155 156
## fromS10 161 162 163 164 165 166 167 168 169 170 171 172
## fromS11 177 178 179 180 181 182 183 184 185 186 187 188
## fromS12 193 194 195 196 197 198 199 200 201 202 203 204
## fromS13 209 210 211 212 213 214 215 216 217 218 219 220
## fromS14 225 226 227 228 229 230 231 232 233 234 235 236
## fromS15 241 242 243 244 245 246 247 248 249 250 251 252
## fromS16 257 258 259 260 261 262 263 264 265 266 267 268
## toS13 toS14 toS15 toS16
## fromS1 29 30 31 32
## fromS2 45 46 47 48
## fromS3 61 62 63 64
## fromS4 77 78 79 80
## fromS5 93 94 95 96
## fromS6 109 110 111 112
## fromS7 125 126 127 128
## fromS8 141 142 143 144
## fromS9 157 158 159 160
## fromS10 173 174 175 176
## fromS11 189 190 191 192
## fromS12 205 206 207 208
## fromS13 221 222 223 224
## fromS14 237 238 239 240
## fromS15 253 254 255 256
## fromS16 269 270 271 272
##
## Response parameters
## Resp 1 : multinomial
## Resp 2 : gaussian
## Re1.pr1 Re1.pr2 Re1.pr3 Re1.pr4 Re1.pr5 Re1.pr6 Re1.pr7 Re1.pr8
## St1 273 274 275 276 277 278 279 280
## St2 283 284 285 286 287 288 289 290
## St3 293 294 295 296 297 298 299 300
## St4 303 304 305 306 307 308 309 310
## St5 313 314 315 316 317 318 319 320
## St6 323 324 325 326 327 328 329 330
## St7 333 334 335 336 337 338 339 340
## St8 343 344 345 346 347 348 349 350
## St9 353 354 355 356 357 358 359 360
## St10 363 364 365 366 367 368 369 370
## St11 373 374 375 376 377 378 379 380
## St12 383 384 385 386 387 388 389 390
## St13 393 394 395 396 397 398 399 400
## St14 403 404 405 406 407 408 409 410
## St15 413 414 415 416 417 418 419 420
## St16 423 424 425 426 427 428 429 430
## Re2.(Intercept) Re2.sd
## St1 281 282
## St2 291 292
## St3 301 302
## St4 311 312
## St5 321 322
## St6 331 332
## St7 341 342
## St8 351 352
## St9 361 362
## St10 371 372
## St11 381 382
## St12 391 392
## St13 401 402
## St14 411 412
## St15 421 422
## St16 431 432pars <- c(unlist(getpars(mod)))
# Restringimos a solo poder ir del estado uno al dos, dos al tres, ...
# esto se hace modificando la inical de la matriz de transición
pars[17:272] <- 0
pars[sapply(1:16, function(i){0+17*i})] <- 1/2
pars[sapply(1:15, function(i){1+17*i})] <- 1/2
pars[255:256] <- 0.5
pars[272] <- 1
# restringimos estado inicial al 1
pars[2:16] <- 0
pars[1] <- 1
mod <- setpars(mod, pars)
summary(mod)## Initial state probabilties model
## pr1 pr2 pr3 pr4 pr5 pr6 pr7 pr8 pr9 pr10 pr11 pr12 pr13 pr14 pr15
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## pr16
## 0
##
##
## Transition matrix
## toS1 toS2 toS3 toS4 toS5 toS6 toS7 toS8 toS9 toS10 toS11 toS12
## fromS1 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## fromS2 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## fromS3 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## fromS4 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## fromS5 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0
## fromS6 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0
## fromS7 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0
## fromS8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0
## fromS9 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0
## fromS10 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.0
## fromS11 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.5
## fromS12 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5
## fromS13 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## fromS14 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## fromS15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## fromS16 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## toS13 toS14 toS15 toS16
## fromS1 0.0 0.0 0.0 0.0
## fromS2 0.0 0.0 0.0 0.0
## fromS3 0.0 0.0 0.0 0.0
## fromS4 0.0 0.0 0.0 0.0
## fromS5 0.0 0.0 0.0 0.0
## fromS6 0.0 0.0 0.0 0.0
## fromS7 0.0 0.0 0.0 0.0
## fromS8 0.0 0.0 0.0 0.0
## fromS9 0.0 0.0 0.0 0.0
## fromS10 0.0 0.0 0.0 0.0
## fromS11 0.0 0.0 0.0 0.0
## fromS12 0.5 0.0 0.0 0.0
## fromS13 0.5 0.5 0.0 0.0
## fromS14 0.0 0.5 0.5 0.0
## fromS15 0.0 0.0 0.5 0.5
## fromS16 0.0 0.0 0.0 1.0
##
## Response parameters
## Resp 1 : multinomial
## Resp 2 : gaussian
## Re1.pr1 Re1.pr2 Re1.pr3 Re1.pr4 Re1.pr5 Re1.pr6 Re1.pr7 Re1.pr8
## St1 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St2 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St3 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St4 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St5 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St6 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St7 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St8 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St9 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St10 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St11 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St13 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St14 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St15 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## St16 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
## Re2.(Intercept) Re2.sd
## St1 0 1
## St2 0 1
## St3 0 1
## St4 0 1
## St5 0 1
## St6 0 1
## St7 0 1
## St8 0 1
## St9 0 1
## St10 0 1
## St11 0 1
## St12 0 1
## St13 0 1
## St14 0 1
## St15 0 1
## St16 0 1Ahora ajustamos usando EM:
set.seed(280572)
fm <- fit(mod, emcontrol=em.control(maxit=100))## iteration 0 logLik: -147614
## iteration 5 logLik: -125759
## iteration 10 logLik: -124984
## iteration 15 logLik: -124584
## iteration 20 logLik: -124487
## iteration 25 logLik: -124415
## iteration 30 logLik: -123993
## iteration 35 logLik: -123867
## iteration 40 logLik: -123859
## iteration 45 logLik: -123858
## iteration 50 logLik: -123850
## iteration 55 logLik: -123828
## iteration 60 logLik: -123821
## iteration 65 logLik: -123821
## iteration 70 logLik: -123821
## iteration 75 logLik: -123821
## iteration 80 logLik: -123821
## converged at iteration 82 with logLik: -123821summary(fm)## Initial state probabilties model
## pr1 pr2 pr3 pr4 pr5 pr6 pr7 pr8 pr9 pr10 pr11 pr12 pr13 pr14 pr15
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## pr16
## 0
##
##
## Transition matrix
## toS1 toS2 toS3 toS4 toS5 toS6 toS7 toS8 toS9 toS10 toS11 toS12
## fromS1 0.78 0.22 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## fromS2 0.00 0.55 0.45 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## fromS3 0.00 0.00 0.64 0.36 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## fromS4 0.00 0.00 0.00 0.57 0.43 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## fromS5 0.00 0.00 0.00 0.00 0.46 0.54 0.00 0.00 0.00 0.00 0.00 0.00
## fromS6 0.00 0.00 0.00 0.00 0.00 0.68 0.32 0.00 0.00 0.00 0.00 0.00
## fromS7 0.00 0.00 0.00 0.00 0.00 0.00 0.63 0.37 0.00 0.00 0.00 0.00
## fromS8 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.41 0.59 0.00 0.00 0.00
## fromS9 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.55 0.45 0.00 0.00
## fromS10 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.58 0.42 0.00
## fromS11 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.58 0.42
## fromS12 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.58
## fromS13 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## fromS14 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## fromS15 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## fromS16 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## toS13 toS14 toS15 toS16
## fromS1 0.00 0.00 0.00 0.0
## fromS2 0.00 0.00 0.00 0.0
## fromS3 0.00 0.00 0.00 0.0
## fromS4 0.00 0.00 0.00 0.0
## fromS5 0.00 0.00 0.00 0.0
## fromS6 0.00 0.00 0.00 0.0
## fromS7 0.00 0.00 0.00 0.0
## fromS8 0.00 0.00 0.00 0.0
## fromS9 0.00 0.00 0.00 0.0
## fromS10 0.00 0.00 0.00 0.0
## fromS11 0.00 0.00 0.00 0.0
## fromS12 0.42 0.00 0.00 0.0
## fromS13 0.68 0.32 0.00 0.0
## fromS14 0.00 0.84 0.16 0.0
## fromS15 0.00 0.00 0.70 0.3
## fromS16 0.00 0.00 0.00 1.0
##
## Response parameters
## Resp 1 : multinomial
## Resp 2 : gaussian
## Re1.pr1 Re1.pr2 Re1.pr3 Re1.pr4 Re1.pr5 Re1.pr6 Re1.pr7 Re1.pr8
## St1 0.00939 0.0201 0.0084 0.06457 0.1847 0.32372 0.3057 0.0834
## St2 0.05431 0.0324 0.0145 0.31046 0.0304 0.00294 0.0081 0.5470
## St3 0.63192 0.0549 0.0274 0.09398 0.0540 0.06518 0.0257 0.0468
## St4 0.00177 0.8782 0.0723 0.04708 0.0001 0.00000 0.0000 0.0006
## St5 0.00073 0.0804 0.8506 0.06823 0.0000 0.00000 0.0000 0.0000
## St6 0.00000 0.0347 0.2024 0.76285 0.0000 0.00000 0.0000 0.0000
## St7 0.08032 0.1913 0.3309 0.15949 0.1096 0.05282 0.0345 0.0411
## St8 0.00000 0.2245 0.7476 0.02787 0.0000 0.00000 0.0000 0.0000
## St9 0.00478 0.9409 0.0258 0.02857 0.0000 0.00000 0.0000 0.0000
## St10 0.83349 0.0049 0.0064 0.15140 0.0000 0.00000 0.0000 0.0038
## St11 0.00000 0.0000 0.0000 0.13273 0.0109 0.00000 0.0027 0.8537
## St12 0.00000 0.0000 0.0000 0.00738 0.0668 0.03626 0.8668 0.0227
## St13 0.00000 0.0000 0.0000 0.00018 0.0501 0.94972 0.0000 0.0000
## St14 0.00000 0.0000 0.0000 0.00000 0.8797 0.12025 0.0000 0.0000
## St15 0.00000 0.0000 0.0000 0.02549 0.0355 0.47455 0.4645 0.0000
## St16 0.38437 0.1515 0.0068 0.09699 0.0000 0.00039 0.0053 0.3546
## Re2.(Intercept) Re2.sd
## St1 5.7 3.0
## St2 7.3 3.2
## St3 10.9 4.9
## St4 10.1 3.9
## St5 12.8 4.4
## St6 16.6 4.7
## St7 18.1 8.9
## St8 12.2 2.9
## St9 9.6 2.7
## St10 10.3 3.5
## St11 10.4 3.5
## St12 10.6 4.1
## St13 14.8 5.5
## St14 15.4 5.4
## St15 7.9 2.9
## St16 10.4 4.2Y podemos simular algunas trayectorias:
set.seed(2805)
fm_1 <- fm
class(fm_1) <- 'depmix'
sim_1 <- simulate(fm_1)
resp_1 <- sim_1@response
estados <- sim_1@states
angulosSim <- lapply(resp_1[estados], function(x){
params <- x[[1]]@parameters$coefficients
c(sample(levs, 1, prob=params ), rnorm(1, x[[2]]@parameters[1][[1]], sd=x[[2]]@parameters[2]$sd))
# rnorm(1, params[[1]], sd=params$sd )}
}
)
dat_x <- data.frame(Reduce(rbind, angulosSim))
names(dat_x) <- c('angulo','velocidad')
# estados
head(dat_x, 6)## angulo velocidad
## 1 -0.48 3.7
## 2 2.77 7.5
## 3 -0.48 9.4
## 4 -0.48 12.4
## 5 2.77 4.4
## 6 -2.78 9.1plot(trazar(dat_x[1:78,]), type='l')
plot(trazar(dat_x[79:(38+78),]), type='l')
La ruta de máxima probabilidad (estados) se puede ver en el objeto posterior. Esta secuencia indica los estados más probables en el sentido de máxima probabilidad posterior (MAP).
plot(fm@posterior[1:60,1])
respuesta <- fm@response[fm@posterior[1:60,1]]
post <- ldply(respuesta, function(comp){
angulo.probs <- comp[[1]]@parameters$coefficients
velocidad <- comp[[2]]@parameters$coefficients
c(levs[which.max(angulo.probs)], velocidad)
})
names(post) <- c('angulo.cut','velocidad')
plot(trazar(post))
En el caso de arriba, se inició en el estado 1 (por construcción), permaneció ahí por 5 tiempos, pasó al estado 2,…
Veamos el caso del dígito 3:
load('./pendigits/modelo_digitos_3.RData')
filtro <- sapply(digitos, function(x) x==3)
resultados <- results.list[filtro]
dat.6 <- ldply(1:500, function(i){
conv.1 <- convertirSerie(resultados[[i]])
data.frame(angulo=conv.1[,1], velocidad=conv.1[,2], serie=i, tiempo=1:nrow(conv.1))
})
longitudes <- ddply(dat.6, 'serie', summarise, long=max(tiempo))
dat.6$angulo.cut <- factor(cut2(dat.6$angulo, g=8, levels.mean=TRUE))
levs <- as.numeric(levels(dat.6$angulo.cut))
plot(fm@posterior[1:44,1])
respuesta <- fm@response[fm@posterior[1:44,1]]
post <- ldply(respuesta, function(comp){
angulo.probs <- comp[[1]]@parameters$coefficients
velocidad <- comp[[2]]@parameters$coefficients
c(levs[which.max(angulo.probs)], velocidad)
})
names(post) <- c('angulo.cut','velocidad')
plot(trazar(post))
En HMM, la variable \(x_{n+1}\) es independiente de:
\(x_1,...,x_{n}\)
\(x_1,...,x_{n-1}\)
\(x_{n-1}\)
Ninguna de las anteriores.
Sea \(A\) la matriz de transición en un HMM. Un modelo de mezclas para datos iid corresponde a un caso especial de HMM donde:
\(A_{jk}\) toman el mismo valor para toda \(j\).
\(A_{jk}\) toman el mismo valor para toda \(j\) y para toda \(k\).
\(A_{jk}\) toman el mismo valor para toda \(k\).
- Ninguna de las anteriores.
Algoritmo EM para modelos markovianos de estados ocultos
Ahora veremos cómo se construye el paso E y M del algoritmo esperanza- maximización para HMM.
Similar al caso de varaibles latentes, comenzamos escribiendo la verosimilitud con datos completos. La serie de datos observados la denotamos como \[X=X_1,\ldots, X_T\] y la de datos latentes \[S=S_1,\ldots, S_T.\]
Supongamos que tenemos una observación \(x=x_1,\ldots, x_T\) con \(s=s_1,\ldots, s_T\). Si suponemos que \(P(S_1=s_1)=1\) la verosimilitud (usando nuestro modelo HMM) es
\[P(x,S=s)=P(x_1|S_1=s_1)P(S_2=s_2|S_1=s_1)P(x_2|S_2=s_2)\cdots P(S_{T}=s_T|S_{T-1}=s_{t-1})P(x_T|S_T=s_T).\]
Si denotamos \(p_{ij}=P(S_t=j|S_{t-1}=i)\) y \(p_j=P(x_t|S_t=j)\). Tomando logaritmo y reacomodando:
\[\log P(x,S=s) =\sum_{t=2}^t \log p_{s_{t-1}, s_{t}} + \sum_{t=1}^T \log p_{s_t}(x_t). \]
Nos conviene reescribir esta ecuación de otra manera. Introducimos entonces las variables indicadoras \(u_j(t)\), tales que \(u_j(t)=1\) si y sólo si \(s_t=j\), y \(v_{ij}(t)\), tales que \(v_{ij}(t)=1\) si y sólo si \(s_{t-1}=i, s_{t}=j\). Entonces podemos reescribr
\[\log P(x,S=s) = \sum_{i=1}^M\sum_{j=1}^M \left( \sum_{t=2}^T v_{ij} (t) \right) \log p_{ij} + \sum_{i=1}^M \sum_{t=1}^T u_i(t)\log p_i(x_t). \]
A partir de esta ecuación es fácil ver la forma del algoritmo.
- Paso E: Necesitamos calcular el valor esperado condicional (dados los datos observados \(X\)) de la expresión de arriba. Basta entonces calcular \[\delta_j (t)=E(u_j(t)|X)=P(S_t=j|X)\] y \[\gamma_{ij}(t) =E(v_{ij}(t)|X)=P(S_{t-1}=i, S_{t}=j|X).\]
Una vez que calculamos estas cantidades (usando el algoritmo hacia adelante-hacia atrás, o de Baum-Welch, como veremos más adelante), procedemos al paso M.
- Paso M: Buscamos maximizar \[\sum_{i=1}^M\sum_{j=1}^M \left( \sum_{t=2}^T \hat{\gamma}_{ij} (t) \right) \log p_{ij} + \sum_{i=1}^M \sum_{t=1}^T \hat{\delta}_i(t)\log p_i(x_t). \]
Notemos que la optimización de los parámetros de \(p_i\) se puede separar de la correspondiente a \(\hat{p}_{ij}\). Más aun, la optimización de \(\hat{p}_{ij}\) se puede separar para cada \(i\).
Es fácil ver que \[\hat{p}_{ij}=\frac{\hat{\gamma}_{ij}}{\sum_l \hat{\gamma}_{il}}\] donde \[\hat{\gamma}_{ij}=\sum_{t=2}^T \hat{\gamma}_{ij} (t),\] que son conteos ponderados por las probabilidades que resultan del paso E.
Para el segundo término, tenemos que usar la forma particular de las densidades condicionales \(p_j(x)\). Si suponemos \(p_i\) densidad Normal \(\hat{\mu_i}\) y \(\hat{\Sigma_i}\) resultan: \[\hat{\mu_i}=\frac{\sum_{t=1}^T \hat{\delta_i}(t) x_t}{\sum_{t=1}^T \hat{\delta_i}(t)}\] \[\hat{\Sigma_i}=\frac{\sum_{t=1}^T \hat{\delta_i}(t) (x_t-\hat{\mu}_i)(x_t-\hat{\mu}_i)^T}{\sum_{t=1}^T \hat{\delta}(t)}\] es el mismo caso que en mezcla de normales. Como ejercicio, calcula para el modelo Poisson.
Algoritmo hacia adelante-hacia atrás (forward-backward)
Utilizamos el algoritmo forward-backward para calcular \(\hat{\gamma}_{ij}\) y \(\hat{\delta}_j\). Recordemos que
\[ \begin{aligned} \hat{\delta}_{j}&=P(S_{t}=j|x)=\frac{P(x, S_t=j)}{P(x)}\\ \hat{\gamma}_{ij}&=P(S_{t-1}=i,S_t=j|x)=\frac{P(x,S_{t-1}=i,S_t=j)}{P(x)} \end{aligned} \]
El cálculo de cada una de las probabilidades de arriba es computacionalmente intensivo, por ejemplo para calcular \(P(x)\):
\[P(x)=\sum_{\mathcal{S}} p_{s_1}(x_1)p_{s_1,s_2}p_{s_2}(x_2)\cdots p_{s_{T-1},s_T}p_{s_T}(x_T)\]
donde \(\mathcal{S}\) son las combinaciones de posibles estados (\(M^T\) posibilidades) por tanto esta aproximación no es factible. Es por esto que surge la necesidad de un algoritmo más eficiente.
El algoritmo hacia adelante-hacia atrás usa el principio de programación dinámica(recursión inteligente) para calcular \(\hat{\gamma}_{ij}\) y \(\hat{\delta}_j\) en tiempo lineal (\(M^2T\)), consta de dos pasos y explota las independencias condicionales del modelo.
Probabilidad hacia adelante
Definimos la probabilidad hacia adelante \(\alpha_i(t)\) como la probabilidad conjunta de observar las primeras \(t\) observaciones \(x^j\) (\(j=1,...,t\)) y siendo \(i\) el estado al tiempo \(t\):
\[\alpha_i(t)=P(X_1=x_1,...,X_T=x_t,S_t=i)\]
La probabilidad se puede evaluar de manera recursiva siguiendo la fórmula:
\(\alpha_i(1) = \pi_k p_i(x_1)\) para \(i=1,...,M\)
\(\alpha_i(t) = p_i(x_t)\sum_{j=1}^M \alpha_j(t-1)p_{j,i}\) para \(t=2,...,T\) e \(i=1,...,M\).
Prueba:
La idea clave es usar \((S_t,X_t) \perp (X_1,...,X_{t-1})|S_{t-1}\) \[ \begin{aligned} \alpha_i(t)&=P(x_1,...,x_t,S_t=i)\\ &=\sum_{j=1}^M P(x_1,...,x_t,S_t=i, S_{t-1}=j)\\ &=\sum_{j=1}^M P(x_1,...,x_{t-1},S_{t-1}=j)P(x_t|S_t=i)P(S_t=1|S_{t-1}=j)\\ &=\sum_{j=1}^M \alpha_j(t-1)p_i(x_t)p_{j,i}\\ \end{aligned} \]
Probabilidad hacia atrás
Definimos la probabilidad hacia atrás \(\beta_i(t)\) como la probabilidad condicional de las observaciones posteriores al tiempo \(t\) (\(x_{t+1},...,x_T\)) dado que el estado al tiempo \(t\) es \(i\).
\[\beta_i(t)=P(x_{t+1},...,x_T|S_t=i)\] para \(t=1,...T-1\).
La recursión de la probabilidad hacia atrás se evalúa como:
\(\beta_i(T)=1\), para \(i=1,...,M\).
\(\beta_i(t)=\sum_{i=1}^M p_{i,j}p_j(x_{t+1})\beta_i(t+1)\) para \(t = 1,...,T-1\).
Prueba:
La idea clave es usar \(X_{t+1} \perp (X_{t+2},...,X_T)|S_{t+1}\)
\[ \begin{aligned} \beta_i(t)&=P(x_{t+1},...,x_T|S_t=i) \\ &=\sum_{j=1}^M P(x_{t+1},...,x_T,S_{t+1}=j |S_t=i)\\ &=\sum_{j=1}^M P(S_{t+1}=j|S_t=i) P(x_{t+1},...,x_T|S_{t+1}=j)\\ &=\sum_{j=1}^M p_{i,j}P(x_{t+1},...,x_T|S_{t+1}=j)\\ &=\sum_{j=1}^M p_{i,j}P(x_{t+1}|S_{t+1}=j)P(x_{t+2},...,x_T|s_{t+1}=j)\\ &=\sum_{j=1}^M p_{i,j}p_j(x_{t+1})\beta_j(t+1) \end{aligned} \]
Escribimos \(\delta\) y \(\gamma\)
Ahora vemos como escrbir \(\delta_j\) y \(\gamma_{i,j}\) usando las probabilidades hacia adelante y hacia atrás:
\[ \begin{aligned} \hat{\delta}_{j}(t)&=P(S_{t}=j|x)\\ &=\frac{P(x,S_{t}=j)}{P(x)}\\ &=\frac{\alpha_j(t)\beta_j(t)}{\sum_{i=1}^M \alpha_i(T)} \end{aligned} \]
Prueba: \[ \begin{aligned} P(x_1,...,x_T,S_t=j)&=P(x_1,...,x_t,S_t=j)P(x_{t+1},...,x_T|S_t=j)\\ &=\alpha_j(t)\beta_j(t) \end{aligned} \]
Para el denominador notemos que:
\[ \begin{aligned} P(x) &= \sum_i^M P(x,S_{t}=i)\\ &= \sum_{i=1}^M \alpha_i(t)\beta_i(t) \end{aligned} \] esto se cumple para cualquier \(t\), así que si tomamos \(t=T\):
\[P(x) = \sum_{i=1}^M \alpha_i(T)\]
En el caso de \(\gamma_{i,j}\) tenemos:
\[ \begin{aligned} \hat{\gamma}_{ij}&=P(S_{t-1}=i,S_t=j|x)\\ &=\frac{P(x,S_{t-1}=i,S_t=j)}{P(x)}\\ &=\frac{\alpha_i(t-1)\beta_j(t)p_{i,j}p_j(x_{t})}{P(x)} \end{aligned} \]
Prueba: \[ \begin{aligned} P(x_1,...,x_T,S_{t-1}=i,S_t=j)&=P(x_1,...,x_{t-1},S_{t-1}=i)P(x_{t+1},...,X_T|S_t=j)P(S_t=j|S_{t-1}=i)P(x_t|S_t=j)\\ &=\alpha_i(t-1)\beta_j(t)p_{i,j}p_j(x_{t}) \end{aligned} \]
Resumen de algoritmo de estimación
Entonces, el algoritmo de estimación itera de la siguiente manera:
Comenzamos con valores inciales \(\hat{\delta}_j\) y \(\hat{\gamma}_{ij}\).
Actualizamos los parámetros \[\hat{p}_{ij}=\frac{\hat{\gamma}_{ij}}{\sum_l \hat{\gamma}_{il}}\] y los correspondientes a las densidades \(p_j(x)\).
Utilizando el conjunto de parámetros actuales (\(\hat{p}_{ij}\) y los correspondientes a \(p_j(x)\)) calculamos \(\hat{\delta}_j\) y \(\hat{\gamma}_{ij}\) a través del algoritmo hacia adelante-hacia atrás.
Iteramos entre 2 y 3.
Algoritmo de Viterbi
En muchas de las aplicaciones de HMM nos interesa hacer inferencia de la secuencia de estados \(\{S_1,...,S_T\}\), en este caso el criterio de optimización es:
\[MAP(S|x)=argmax_s P(S|x) = argmax_s P(S,x)\]
Aqui estamos buscando el camino más probable. Si consideramos un algoritmo de fuerza bruta, esto es, realizamos búsqueda exahustiva sobre todas las posibles secuencias tendríamos que considerar \(M^T\) casos. Es por ello que nuevamente recurrimos a un algoritmo de programación dinámica.