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?

I made some changes that I hope will clarify what I'm asking. Any clearer? — lrthistlethwaite, Mar 09 '17 at 19:17
So, are you looking for the joint probability of the entire observation sequence? (Why?) — Juho Kokkala, Mar 09 '17 at 19:37
I believe that is the correct terminology, yes. P(O1=o1, O2=o2, ..., ON=oN) = ? I need to know how rare it is to see a given "time series", given a model. — lrthistlethwaite, Mar 09 '17 at 20:04
I don't know why this question was closed. Calculating a likelihood is a standard application of HMMs. @areyoujokingme you can run the forward algorithm all the way to the end, and then sum out the $X$ — Taylor, Mar 09 '17 at 20:39
@Taylor, I think you're right. I found a place in my textbook that validates your answer. I would love to post it in response to my question, but I guess it's "on hold". How do I reverse the status of this question? This is definitely an important use case of HMMs, for sure! — lrthistlethwaite, Mar 09 '17 at 20:56
@areyoujokingme I think your post could use a few edits first. Maybe mention that this probability only has to do with the observations and not the hidden states. Also maybe use the word likelihood. I do have to admit that I had a difficult time with the last couple paragraphs. — Taylor, Mar 09 '17 at 22:34
@Taylor, I'm not an amazing statistician at this point in my life, but I gave it another go, though without "likelihood", as I'm not comfortable with the term yet. Any help would be very appreciated. I have a page from Ewens & Grant's Statistical Models in Bioinformatics that mimicks your answer in the comments above. — lrthistlethwaite, Mar 09 '17 at 22:41

Taylor · Accepted Answer · 2017-03-10T19:20:41.620

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.

Interpreting output of HMM

1 Answers1

1.)

2.)

3.)