Pasos de EM para un modelo ocultado de Markov cuando observables son agregados de los Estados

Question

Se dan varias cadenas de Markov que se ejecutan en paralelo, supongamos que uno no sabe que la cadena se mueve de un estado a otro en cada momento, pero en lugar de observar el monto total de las cadenas en cada estado en cada momento.

¿Cómo se puede encontrar el Modelo Oculto de Markov en este entorno?

Ejemplo:

Tomar dos cadenas y dos estados $s_1,s_2$:

• En el momento $t=1$, $|S_1|=1, |S_2|=1$ lo que significa que hay una persona de la cadena en cada estado
• En el momento $t=2$, $|S_1|=2, |S_2|=0$ lo que significa que ambas cadenas se trasladó a en el primer estado y ninguno en el otro
• En el momento $t=3$, $|S_2|=0, |S_1|=2$ por lo tanto ambas cadenas se trasladó a la segundo estado y no quedaba nada en el otro

Más preguntas relacionadas con:

1. Lo que si algunos estados no ir a trabajar?

2. Lo que si la probabilidad de transición está cambiando con el tiempo?

Answer 1

Inicialización: Si $\mathbf{M}=(m_{ij})_{i,j=1,\ldots,k}$ $k\times k$ de la matriz de transición y si $(X_t^c)_{t=1,\ldots,T}^{i=1,\ldots,C}$ indica el $C$ cadenas de Markov, la distribución (y por lo tanto la probabilidad) de la muestra es $$\prod_{c=1}^C\prod_{t=2}^T m_{x_{t-1}^cx_{t-1}^c}$$ con un log-probabilidad igual a $$\sum_{c=1}^C\sum_{t=2}^T \log(m_{x_{t-1}^cx_{t-1}^c})=\sum_{i=1}^k\sum_{j=1}^k\log(m_{ij})\sum_{c=1}^C\sum_{t=2}^T\mathbb{I}_i(x_{t-1}^c)\mathbb{I}_j(x_{t}^c)=\sum_{i=1}^k\sum_{j=1}^k\log(m_{ij})\sum_{t=2}^T n_{ij}^t$$

E-step: Given that $(X_t^c)_{t=1,\ldots,T}^{i=c,\ldots,C}$ is not observed but only $S_t^j=\sum_{c=1}^C \mathbb{I}_j(x_{t}^c)$ for $j=1,\ldots,k$ and $t=1,\ldots,T$, el E-paso consiste en computación $$\varrho_{ij}^t=\mathbb{E}_{\mathbf{M}}\left[\sum_{c=1}^C\mathbb{I}_i(x_{t-1}^c)\mathbb{I}_j(x_{t}^c)\Big|\mathbf{S}\right]=\sum_{c=1}^C \mathbb{E}_{\mathbf{M}}\left[\mathbb{I}_i(x_{t-1}^c)\mathbb{I}_j(x_{t}^c)\Big|\mathbf{S}\right]=C\mathbb{E}_{\mathbf{M}}\left[\mathbb{I}_i(x_{t-1}^c)\mathbb{I}_j(x_{t}^c)\Big|\mathbf{S}\right]$$Now, the distribution of $(n_{ij}^t)_{i,j=1,\ldots,k}$ given $\mathbf{S}$ es restringido Multinomial $$\prod_{i=1}^k{S_t^i\choose n_{i1}\cdots n_{ik}}\prod_{j=1}^k m_{ij}^{n_{ij}^t}\times\prod_{j=1}^k\mathbb{I}_{\sum_{i=1}^k n^t_{ij}=S_{t+1}^j}$$ y los vectores $(n_{ij}^t)_{i,j=1,\ldots,k}$ son independientes para todos los $t$'s. Yo no creo que exista una forma cerrada de la fórmula en este caso para $\varrho_{ij}^t=\mathbb{E}_{\mathbf{M}}\left[n^{t}_{ij}\big|\mathbf{S}\right]$ y usted puede calcular las expectativas, al sumar sólo al $k$ es muy pequeña. Por lo tanto, la solución razonable parece ser un Monte Carlo aproximación a las cantidades.

M-paso: la Optimización de la$$\sum_{i=1}^k\sum_{j=1}^k\log(m_{ij})\sum_{t=2}^T\varrho_{ij}^t$$ in $\mathbf{M}$ es sencillo: $$\hat{m}_{ij}=\sum_{t=1}^T\varrho_{ij}^t\Big/\sum_{\ell=1}^k \sum_{t=1}^T\varrho_{i\ell}^t$$

Nota: En el ejemplo proporcionado en la pregunta, el oculto de la cadena de Markov no se oculta como es posible reconstruir la individual cadenas a través de los tiempos $t=2,3$, hasta una permutación de los índices de $c$.

Sin ninguna pretensión de eficiencia, aquí es un R de aplicación del código de Monte Carlo EM (o más precisamente, de una MCMC-EM) concepto en esta configuración:

k=3  #number of states
T=10 #number of time steps
C=20 #number of parallel Markov chains

#true matrix M
M=matrix(0,k,k)
for (i in 1:k){
  M[i,]=rgamma(k,i)
  M[i,]=M[i,]/sum(M[i,])}
#production of C parallel chains
cha=matrix(0,C,T)
cha[,1]=sample(1:k,C,rep=TRUE)
for (t in 2:T)
for (c in 1:C)
  cha[c,t]=sample(1:k,1,prob=M[cha[c,t-1],])
#observable summaries
S=matrix(0,k,T)
for (i in 1:k)
  S[i,]=apply(cha==i,2,sum)

#complete likelihood
complik <-function(mov,M){
  (min(mov)>-1)*sum(mov*log(M)-lfactorial(mov))}

#MCMCEM function (L stands for the number of simulations)
MCMCEM <- function(S,L=1e3){

M0=matrix(0,k,k) #arbitrary starting Markov matrix
for (i in 1:k){
  M0[i,]=rgamma(k,i)
  M0[i,]=M0[i,]/sum(M0[i,])}

run=0 #stopping rule for EM
while (run==0){

#E step
en=matrix(0,k,k) #matrix of n_{ij}
for (t in 2:T){

  #starting realisation complying with constraints
  ave=cur=matrix(0,k,k)
  check=0 #are constraints verified?
  while (check==0){

     for (u in 1:(k-1)) 
       cur[u,]=rmultinom(1,S[u,t-1],prob=S[,t]*M0[u,])
     #filling the last row
     cur[k,]=S[,t]-apply(cur[-k,],2,sum)
     check=(min(cur[k,])>=0)}
     ave=ave+cur

  for (l in 1:L){

     #random change under the constraint
     prop=cur
     #changing a row or a column of current (n_ij) matrix
     if (runif(1)<.5){
         u=sample(1:k,2)
         check=0
         while (check==0){
           prop[u[1],]=rmultinom(1,S[u[1],t-1],prob=S[,t]*M0[u[1],])
           prop[u[2],]=S[,t]-apply(prop[-u[2],],2,sum)
           check=(min(prop[u[2],])>=0)}
     }else{
         u=sample(1:k,2)
         check=0
         while (check==0){
           prop[,u[1]]=rmultinom(1,S[u[1],t],prob=S[,t]*M0[,u[1]])
           prop[,u[2]]=S[,t-1]-apply(prop[,-u[2]],1,sum)
           check=(min(prop[,u[2]])>=0)}

     # Metropolis acceptance step
     if (log(runif(1))<complik(prop,M)-complik(cur,M)) 
         cur=prop
     ave=ave+cur}

  #Monte Carlo & temporal average    
  en=en+ave/L
} #t

#M step
for (u in 1:k) en[u,]=en[u,]/sum(en[u,])
diff=max(abs(M0-en))
#new value of Markov transition matrix
M0=en
#stopping rule
run=(diff<1e-2)}

return(M0)}

Mientras que la Metrópoli paso no es totalmente correcto en el que yo no se dividen por la propuesta de distribución, aquí está la comparación de la solución producida por esta función con el real M:

> MCMCEM(S)
          [,1]      [,2]      [,3]
[1,] 0.6719400 0.2228070 0.1052530
[2,] 0.2558264 0.5534466 0.1907271
[3,] 0.1326174 0.5417582 0.3256244
> M
          [,1]      [,2]      [,3]
[1,] 0.6363601 0.2500658 0.1135741
[2,] 0.2488576 0.5786064 0.1725360
[3,] 0.1229905 0.5530896 0.3239198