3
$\begingroup$

Given several Markov chains running in parallel, suppose one doesn't know which chain is moving from one state to another at each time, but instead observe the aggregate amount of chains in each state at each time.

How can one find the Hidden Markov Model in this setting?

Example:

Take two individual chains and two states $s_1,s_2$:

  • At time $t=1$, $|S_1|=1, |S_2|=1$ which means there is one individual chain in each state
  • At time $t=2$, $|S_1|=2, |S_2|=0$ which means both chains moved to the first state and none to in the other
  • At time $t=3$, $|S_2|=0, |S_1|=2$ hence both chains moved to the second state and none remained in the other

More related questions:

  1. What if some states do not commute?

  2. What if transition probability is changing over time?

$\endgroup$
  • $\begingroup$ Could you write down your model? I don't understand how 'individuals' relates to the state space. Also you haven't told us what $S_1$ and $S_2$ are. $\endgroup$ – Taylor Jan 19 '17 at 0:25
  • $\begingroup$ @Taylor see if you understand this $\endgroup$ – ZHU Jan 19 '17 at 0:31
  • $\begingroup$ @Xi'an can you give an E step and an M step? $\endgroup$ – ZHU Jan 19 '17 at 16:18
  • $\begingroup$ @Xi'an first of all how to specify a likelihood function of a Markov Chain $\endgroup$ – ZHU Jan 21 '17 at 4:33
4
+50
$\begingroup$

Initialisation: If $\mathbf{M}=(m_{ij})_{i,j=1,\ldots,k}$ is the $k\times k$ transition matrix and if $(X_t^c)_{t=1,\ldots,T}^{i=1,\ldots,C}$ denotes the $C$ Markov chains, the distribution (and hence the likelihood) of the sample is $$\prod_{c=1}^C\prod_{t=2}^T m_{x_{t-1}^cx_{t-1}^c}$$ with a log-likelihood equal to $$\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$, the E-step consists in computing $$\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}$ is a constrained 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}$$ and the vectors $(n_{ij}^t)_{i,j=1,\ldots,k}$ are independent for all $t$'s. I do not think there is a closed-form formula in this case for $\varrho_{ij}^t=\mathbb{E}_{\mathbf{M}}\left[n^{t}_{ij}\big|\mathbf{S}\right]$ and you can compute those expectations by summing up only when $k$ is very small. Therefore, the reasonable solution seems to be a Monte Carlo approximation to those quantities.

M-step: Optimising$$\sum_{i=1}^k\sum_{j=1}^k\log(m_{ij})\sum_{t=2}^T\varrho_{ij}^t$$ in $\mathbf{M}$ is straightforward: $$\hat{m}_{ij}=\sum_{t=1}^T\varrho_{ij}^t\Big/\sum_{\ell=1}^k \sum_{t=1}^T\varrho_{i\ell}^t$$

Note: In the example provided in the question, the hidden Markov chain is not hidden as it is possible to reconstruct the individual chains across the times $t=2,3$, up to a permutation of the indices $c$.

Without any claim to efficiency, here is an R code implementing the Monte Carlo EM (or more precisely an MCMC-EM) concept in this setting:

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)}

While the Metropolis step is not totally right in that I do not divide by the proposal distribution, here is the comparison of the solution produced by this function with the actual 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
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Can you add in more detail in explaining the meaning of the variables such as $C$ and the decomposition of likelihood function into indicator function and the meaning of the indicator functions? And can you continue on finishing the $M$-step? $\endgroup$ – ZHU Jan 22 '17 at 6:35
  • $\begingroup$ Does this apply to higher dimensions? i.e. more individuals and more states? How to calculate $Q_{i,l}$ actually? $\endgroup$ – ZHU Jan 23 '17 at 7:19
  • $\begingroup$ @ZHU: As updated above, the EM step is unlikely to produce a closed-form expression. After a day or so, I gave up. The alternative is to use simulation to produce an approximate expectation, this is called MCEM. $\endgroup$ – Xi'an Jan 25 '17 at 7:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.