Estimación por maximización directa de la probabilidad.

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

3.1. Introducción

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

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

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

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

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

Primero consideraremos el caso estacionario.

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

3.2. Escalando la computación de la posibilidad

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

\[ \phi_t =\alpha_t/w_t\]

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

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

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

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

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

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

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

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

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

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

3.3. Maximización de la posibilidad sujeta a restricciones

3.3.1 Re parametrización para evitar restricciones

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.4. Otros problemas

3.4.1. Múltiples máximos en la posibilidad

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

3.4.2. Valores iniciales para las iteraciones

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

3.4.3. Posibilidad sin limites

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

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

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

3.5. Ejemplo: Terremotos

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

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

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

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

iter = 1 
LL = -342.3180692 
diff = Inf 

iter = 2 
LL = -341.8839295 
diff = 0.4341397 

iter = 3 
LL = -341.8808283 
diff = 0.003101278 

iter = 4 
LL = -341.8799367 
diff = 0.0008916081 

iter = 5 
LL = -341.8794228 
diff = 0.0005138536 

iter = 6 
LL = -341.8791236 
diff = 0.0002991765 

iter = 7 
LL = -341.8789489 
diff = 0.0001747594 

iter = 8 
LL = -341.8788465 
diff = 0.000102314 

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

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

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

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

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

print(summary(y))

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

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

$nonstat
[1] TRUE

$distn
[1] "pois"

$pm
$pm$lambda
[1] 15.42342 26.02398


$discrete
[1] TRUE

$n
[1] 107

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

[1] -341.8787

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

El modelo de tres estados es:

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

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

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

iter = 1 
LL = -329.9251807 
diff = Inf 

iter = 2 
LL = -328.5465037 
diff = 1.378677 

iter = 3 
LL = -328.5278471 
diff = 0.01865661 

iter = 4 
LL = -328.5274950 
diff = 0.0003520581 

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

print(summary(y))

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

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

$nonstat
[1] TRUE

$distn
[1] "pois"

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


$discrete
[1] TRUE

$n
[1] 107

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

[1] -328.5275

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

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

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

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

iter = 1 
LL = -327.8302182 
diff = Inf 

iter = 2 
LL = -326.9581836 
diff = 0.8720346 

iter = 3 
LL = -326.9289508 
diff = 0.02923274 

iter = 4 
LL = -326.9204573 
diff = 0.008493573 

iter = 5 
LL = -326.9134060 
diff = 0.007051262 

iter = 6 
LL = -326.9074158 
diff = 0.005990187 

iter = 7 
LL = -326.9024761 
diff = 0.004939741 

iter = 8 
LL = -326.8985207 
diff = 0.003955416 

iter = 9 
LL = -326.8954326 
diff = 0.003088091 

iter = 10 
LL = -326.8930709 
diff = 0.002361707 

iter = 11 
LL = -326.8912938 
diff = 0.001777022 

iter = 12 
LL = -326.8899734 
diff = 0.00132043 

iter = 13 
LL = -326.8890015 
diff = 0.0009718974 

iter = 14 
LL = -326.8882912 
diff = 0.0007103225 

iter = 15 
LL = -326.8877747 
diff = 0.0005164488 

iter = 16 
LL = -326.8874007 
diff = 0.0003740631 

iter = 17 
LL = -326.8871305 
diff = 0.0002701856 

iter = 18 
LL = -326.8869357 
diff = 0.0001947661 

iter = 19 
LL = -326.8867955 
diff = 0.0001401981 

iter = 20 
LL = -326.8866947 
diff = 0.000100815 

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

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

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

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

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

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

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

print(summary(y))

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

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

$nonstat
[1] TRUE

$distn
[1] "pois"

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


$discrete
[1] TRUE

$n
[1] 107

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

[1] -326.8865

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

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

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

3.6. Errores estándares e intervalos de confianza

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

3.6.1. Errores estándares mediante el Hessiano

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

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

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

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

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

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

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

Computación recursiva del Hessiano

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

3.6.2. Errores estándares Bootstrap e intervalos de confianza

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

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

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

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

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

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