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


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"), 
                           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?

  • $\begingroup$ I made some changes that I hope will clarify what I'm asking. Any clearer? $\endgroup$ Mar 9, 2017 at 19:17
  • 1
    $\begingroup$ So, are you looking for the joint probability of the entire observation sequence? (Why?) $\endgroup$ Mar 9, 2017 at 19:37
  • $\begingroup$ 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. $\endgroup$ Mar 9, 2017 at 20:04
  • $\begingroup$ 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$ $\endgroup$
    – Taylor
    Mar 9, 2017 at 20:39
  • $\begingroup$ @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! $\endgroup$ Mar 9, 2017 at 20:56

1 Answer 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.


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


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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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