Interpreting output of HMM

Question

I have a 3 state HMM model ("state1", "state2", "state3"), with an alphabet of "1" ("hit") and "0" ("miss"). Here are the parameters (these are just for example's sake) that this HMM is defined by:

Start State Probabilities

startProbs=c(0.5,0.5,0)

State Transition Probability Matrix

transProbs=rbind(c(0.667,0.333,0), c(0.333,0.333,0.333), c(0,0.333,0.667))

Emission Probability Matrix

emissionProbs=cbind(c(0.95,0.05), c(0.75,0.25), c(0.05,0.95))

I'm using an R package called HMM. It makes it easy to initialize the model with the following call:

hmmModelApriori <- initHMM(States=c("S1", "S2", "S3"), Symbols=c("1", "0"), 
                           startProbs=c(0.5,0.5,0), 
                           transProbs=rbind(c(0.667, 0.333, 0),
                                            c(0.333, 0.333, 0.333),
                                            c(0,     0.333, 0.667)), 
                           emissionProbs=cbind(c(0.95,0.05), c(0.75,0.25), c(0.05,0.95)))

Using the package, I can calculate the posterior of a given observation sequence, such as:

obs = c("1","0","0","0","0","0","0","0","0","0","0","0","0","0","0","0","0","0","1")

The HMM model outputs probabilities for each "time point", for all 3 states, where the probabilities of all states for a given time point sum to 1.

p = posterior(hmmModelApriori, obs)

My question is, how do I get the joint probability for my test sequence (the variable "obs" shown above)?

P(obs) = ?

The output of the function "posterior" gives you a matrix of probabilities with dimension #states x #observations. However, I want a single probability for the entire observation sequence, obs. How do I calculate P(obs) from the output of HMM's posterior() function, stored as a variable "p" in the above code?

Answer 1

1.)

You want the likelihood: $p(O_1,\ldots,O_n|\lambda)$, where $\{O_i\}_i$ are the observations, and $\lambda = (A,B,\pi)$ are the model parameters (assumed known here). In your case, $A$ is your "emission probability matrix," $B$ is your "state transition probability matrix," and $\pi$ is(are) your "start state probabilities." In his tutorial, Rabiner calls this problem 1 of the three basic problems of HMMs.

2.)

I've never used the HMM package, but I think the only function they have to help you get this quantity is the function forward(), which gives you forward probabilities. Forward probabilities, denoted $\alpha_t(i)$, are defined as $$ \alpha_t(i) = p(O_1,\ldots,O_t,q_t=S_i|\lambda), $$ where $q_t$ is the time $t$ state variable.

There are recursive formulas for these, and they're pretty straightforward to verify (just use conditional independence assumed by your model). Let $a_{i,j}$ be the $(i,j)th$ element of the state transition matrix, and $b_j(O_{t+1}) = p(O_{t+1}|q_t=S_j)$. Then: \begin{align*} \alpha_{t+1}(j) &= p(O_1,\ldots, O_t,O_{t+1},q_{t+1} = S_j)\\ &= \sum_{i=1}^3 p(O_1,\ldots, O_t,O_{t+1},q_t = S_i, q_{t+1} = S_j)\\ &= \sum_{i=1}^3 p(O_{t+1}|q_{t+1}=S_j) p(q_{t+1}=S_j|q_t=S_i) p(O_1,\ldots, O_t,q_t = S_i)\\ &= \sum_{i=1}^3 \alpha_t(i) a_{i,j}b_j(O_{t+1}). \end{align*}

3.)

How do you use this? Well, if you get $\alpha_{T}(i)$, the $\alpha$ distribution for your final time point $T$, and you get this for $i=1,2,3$, then you can "sum out" the state.

$$ p(O_1,\ldots,O_T|\lambda) = \sum_{i=1}^3 \alpha_T(i). $$

For code, check out this thread.