# \[ANN\] HMMBase.jl - A lightweight and efficient Hidden Markov Model abstraction

**URL:** <https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604>\
**Category:** Package Announcements\
**Tags:** statistics\
**Created:** [March 7, 2019, 9:30pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604 "2019-03-07T21:30:33Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![maxmouchet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxmouchet/32/7834_2.png) [@maxmouchet](https://discourse.julialang.org/u/maxmouchet)\
**Post date:** [March 7, 2019, 9:30pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/1 "2019-03-07T21:30:33Z")

</div>

Hi,

I’d like to introduce [HMMBase.jl](https://github.com/maxmouchet/HMMBase.jl). A package providing basic building blocks for hidden Markov models in Julia.

### Motivation

While there are many implementations of HMMs for Julia, and other languages, I found out that:

- Julia packages, such as [HiddenMarkovModels.jl](https://github.com/BenConnault/HiddenMarkovModels.jl) were not maintained anymore and thus not compatible with recent Julia version.
- Most HMMs libraries tend to directly implement the algorithms for a given probability distribution, most commonly discrete and normal distributions, and hence cannot easily support new distributions.

Instead, I chose to rely on the `Distribution` interface provided by [Distributions.jl](https://github.com/JuliaStats/Distributions.jl) to handle arbitrary distributions. As long as the distribution implements the `pdf` and `fit!` method, it is supported by HMMBase. Whether it is univariate or multivariate, a single distribution or a mixture model.

For example, we can sample a two states HMM with two different distributions as follows:

```julia-auto
hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Gamma(1,1)])
z, y = rand(hmm, 1000)

```

The only constraint being that each observation distribution must have the same dimension (e.g., it is not possible to mix 2d and 3d normal distributions).

Similarly, arbitrary _containers_ that conforms to the `AbstractArray` interface are supported to store the model parameters and the states/observations. This means that, for example, [ArrayFire.jl](https://github.com/JuliaComputing/ArrayFire.jl) could be used to perform some computations on the GPU. That said, I only have tested standard Julia arrays and [StaticArrays](https://github.com/JuliaArrays/StaticArrays.jl) for now.

### Goals

My goal is to provide well-tested and efficient implementations of the following algorithms:

- Forward-backward
- Viterbi
- Baum-Welch (MLE estimator)

The MLE estimator implements only the E part and relies on the `fit_mle` method of each distribution for the M part.

For now my (early) benchmark is against the [hmmlearn](https://github.com/hmmlearn/hmmlearn) Python library which implements the same algorithm in Cython (for normal distributions only). [This issue](https://github.com/maxmouchet/HMMBase.jl/issues/2) shows some early results.

 ![](https://global.discourse-cdn.com/julialang/original/3X/c/c/cceb0682dcebebc1db9d1046b2141ed8c05ca816.png)

### Availability

The package is available for Julia 1.1+ in my own registry:

```julia-auto
pkg> registry add https://github.com/maxmouchet/JuliaRegistry.git
pkg> add HMMBase

```

Feel free to consult the [documentation](https://maxmouchet.github.io/HMMBase.jl/stable/) for examples.

Right now the only package that depends on HMMBase, is my implementation of the Gibbs sampler for the hierarchical Dirichlet process hidden Markov model, [HDPHMM.jl](https://github.com/maxmouchet/HDPHMM.jl). Both of these packages will be maintained during my PhD as some of my work depends on them. Hopefully, HMMBase will be stable enough by the time I graduate that it does not need much work anymore.

I’m open to feedback, ideas, and contributions 🙂

---

<div class="post-metadata">

**Author:** ![dellison](https://avatars.discourse-cdn.com/v4/letter/d/ce73a5/32.png) [@dellison](https://discourse.julialang.org/u/dellison)\
**Post date:** [March 13, 2019, 3:16am UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/2 "2019-03-13T03:16:45Z")

</div>

This is cool, thanks for sharing! 🙂

To me, it looks like HMMBase.jl implements only first-order models at the moment (i.e. the transition probabilities are conditioned on only the previous one state of history). Is this right?

The reason that I’m asking is that in natural language processing, HMMs have been used for sequence labeling tasks like part-of-speech tagging. In this task, for example, a second-order HMM that conditions on the previous two tags generally seems to do better than a first-order model that just conditions on one previous tag.

(If this would need to be a separate subtype of `AbstractHMM`, perhaps this is an area where I could contribute!)

---

<div class="post-metadata">

**Author:** ![maxmouchet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxmouchet/32/7834_2.png) [@maxmouchet](https://discourse.julialang.org/u/maxmouchet)\
**Post date:** [March 21, 2019, 5:27pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/3 "2019-03-21T17:27:07Z")

</div>

Sorry for my late reply! Indeed it only implements first-order models.

I’m not very familiar with higher-order models, but right now I see two ways of implementing them:

- Writing specialized versions of the functions (e.g. `messages_forwards_2`);
- Implementing more generic algorithms (such as belief propagation) to handle models of arbitrary orders in a single function;

The second option seems cleaner but I’m worried about the performance implications. On the other hand, macros could be used to generate specialized function for an arbitrary order instead.

As for the types, I can add a new parameter like `AbstractHMM{F<:VariateForm, N}` where `N` represents the order (1 by default).

That said, feel free to draft something for second-order models and to open an issue/PR if you want. I’ll do some experiments also on my side.

---

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [May 10, 2022, 2:13pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/4 "2022-05-10T14:13:19Z")

</div>

Hi! Did you do end up doing any work on second-order HMMs? I have an interest in that.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [February 22, 2024, 8:37pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/5 "2024-02-22T20:37:55Z")

</div>

Hi! What do you mean by second-order HMMs? I’ve taken over from Max as mister HMM, so my package HiddenMarkovModels.jl is probably your best bet

---

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [February 22, 2024, 8:55pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/6 "2024-02-22T20:55:52Z")

</div>

Hah! I had just asked, and then deleted, a similar question in response to your HiddenMarkovModels discussion in another thread, because I thought it was an old conversation and thus not the right place, but then see you replied to my old question in this thread. I’m all turned around!

In any case, I meant/mean HMMs where the observed state depends on two hidden states.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [February 23, 2024, 5:43am UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/7 "2024-02-23T05:43:07Z")

</div>

I think it is possible to do it from a standard HMM by encoding a tuple of states `(s1, s2)` as a single state `s` and using a sparse transition matrix. I suggested something similar as an answer to [this issue](https://github.com/gdalle/HiddenMarkovModels.jl/issues/12), do you need help figuring out the details or is it sufficiently clear?

---

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [February 23, 2024, 2:16pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/8 "2024-02-23T14:16:58Z")

</div>

Interesting, I’ll look at that, thanks. The implementation I’m currently using calculates a separate transition matrix for each observed state and then selects which one to use for each observed state in the forward calculation.

Why do you say the transition matrix would be sparse?

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [February 23, 2024, 3:15pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/9 "2024-02-23T15:15:45Z")

</div>

> [@BioTurboNick](#):
>
> Why do you say the transition matrix would be sparse?

Because a transition would go from (the encoding of) `(s1, s2)` to (the encoding of) `(s1', s2')`. In particular, the only valid transitions are those for which `s2 = s1'`.

---

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [February 23, 2024, 6:32pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/10 "2024-02-23T18:32:15Z")

</div>

Thanks, I understand. I may have a more complicated situation then, or at least I’m not sure how it maps - the original code was written by someone else and I’m getting my head around HMMs.

The matrices being used are computed from physical parameters, expressing “given we observed state {s1, s2, …}, the probability that the hidden states transitioned from hidden states hx to hy are given by {M1, M2, …}”. Actually they might be called “emission matrices”? But I don’t see that word appear in your documentation. Very likely I’m not understanding something fundamental.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [February 23, 2024, 6:37pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/11 "2024-02-23T18:37:25Z")

</div>

Emissions are usually a synonym for observations in HMM parlance, so I picked one but the other term is equivalent. What you’re referring to as the “hidden state” is also just called the state, and what you’re calling “state” might be the observations?  
In any case, I do think there’s an issue with your mathematical reasoning (cause I don’t understand it), and it would be easier to clarify on a Zoom call. Wanna chat in the coming days?

---

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [February 23, 2024, 7:00pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/12 "2024-02-23T19:00:01Z")

</div>

Thanks for the offer, but it’s not critical. Let me just try this to make it clearer:

`n` Observations: 0, 1

`m` States: 0, 1, 2

Given some physical parameters of the process, compute `n`, `m x m` matrices, which are normalized by row across all `n`.

The forward log likelihood function takes a vector of observations, a vector of initial state probabilities, and the matrices. For each observation, the matching matrix is selected and combined with the previous likelihoods.

If still not understandable, no worries. Though if you know of a text or other resource that could help, I’d appreciate that!

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [February 24, 2024, 7:50am UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/13 "2024-02-24T07:50:59Z")

</div>

The framework used by HiddenMarkovModels.jl is that of controlled HMMs, it is described in section 4.2 of my PhD thesis: [Machine learning and combinatorial optimization algorithms, with applications to railway planning - PASTEL - Thèses en ligne de ParisTech](https://pastel.hal.science/tel-04053322).  
In the graphical model below, the state X\_t is hidden, while the controls u\_t and observations Y\_t are known. An arrow represents an influencing relation, roughly speaking: for instance the distribution of X\_t can be written \mathbb{P}(X\_t \mid X\_{t-1}, u\_t), and the distribution of Y\_t can be written \mathbb{P}(Y\_t \mid X\_t, u\_t).

 ![Screenshot 2024-02-24 at 08.44.55](https://global.discourse-cdn.com/julialang/original/3X/7/5/757ad55da96f06d512f71c5deb497b19869c29fa.png)

Let us denote by \theta the parameters of the model, which include the initial state distribution, state transition matrices and the distributions relating each state to the observations it can generate. Usually, when we compute (and maximize) the “likelihood” of an HMM, we are talking about the following quantity

\mathbb{P}\_{\theta}(Y\_{1:T} \mid u\_{1:T}) = \sum\_{X\_{1:T}} \mathbb{P}\_{\theta}(Y\_{1:T} \mid u\_{1:T}, X\_{1:T}) \mathbb{P}\_{\theta}(X\_{1:T} \mid u\_{1:T})

This sum has an exponential number of terms, but the special structure of HMMs means we can compute it efficiently with dynamic programming. That is one thing my package does.

From what you tell me, it seems that those binary “observations” you speak of are actually “controls”, in the sense that they influence the evolution of the state instead of being influenced by it. But if my interpretation is correct, it means that in your case, there is no Y\_t? How do you get any information about the states X\_t? What are those “physical parameters of the process” you mentioned?

---

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [February 26, 2024, 5:04pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/14 "2024-02-26T17:04:59Z")

</div>

We have videos of single-molecule blinks that are converted into a list of points in each frame, and the points are clustered and then assembled into a temporal trace of ON/OFF observations for each cluster. So each frame is a discretization of a continuous underlying process we want to simulate/learn.

Whether a spot is detected in a discrete frame depends the on rate and off rate of the underlying continuous physical process, but also on the minimum on time required for an ON observation to be recorded, and there’s also an adjustment for false positive detections.

These parameters are input into an inverse Laplace transformation to populate the matrices I was describing, and accounts for the possibility that one or many hidden state changes could occur during the discrete frame. One matrix describes the probability of an ON observation given that a given frame started in one state and ended in another state. The other does the same for an OFF observation.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [February 26, 2024, 5:12pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/15 "2024-02-26T17:12:49Z")

</div>

Do clusters evolve independently? What set does the “state” of a frame belong to? Do you observe this state at each frame, and they are just “hidden” in the continuous time evolution between two frames?

---

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [February 26, 2024, 5:27pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/16 "2024-02-26T17:27:31Z")

</div>

Yes, each cluster is independent - the input being analyzed or simulated is a vector of 0s and 1s from one cluster, for “spot not observed” and “spot observed”, or OFF/ON.

If I understand your last question properly, yes. It’s a series of state changes of variable duration that are integrated and discretized into the observation sequence. (though also there is a possibility of a permanent-off state that would be indistinguishable from the normal off state by observation, other than no longer seeing any ON observations after that point).

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [February 26, 2024, 7:12pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/17 "2024-02-26T19:12:35Z")

</div>

Okay, let’s see if I got this right.

Our random variables are:

- a process of hidden cluster states X(t) \in \mathcal{X} = \{1, \dots, m\} evolving in continuous time (t \in \mathbb{R})
- a sequence of frames of that state X\_k = X(t\_k) taken at discrete instants (k \in \mathbb{N}) but still hidden
- a sequence of on-off observations Y\_k \in \mathcal{Y} = \{0, 1\} associated with each frame

Our parameters are:

- a vector of initial state probabilities \pi \in \{0, 1\}^m such that

\mathbb{P}(X\_0 = i) = \pi\_i

- two emission matrices M\_0 and M\_1 in \mathbb{R}^{m \times m} such that

\begin{align\*} \mathbb{P}(Y\_k = 0 \mid X\_{k-1} = i, X\_{k} = j) &= M\_0[i,j] \\ \mathbb{P}(Y\_k = 1 \mid X\_{k-1} = i, X\_{k} = j) &= M\_1[i,j] \end{align\*}

The parameters are known, and for some reason we want to compute the loglikelihood

\log \mathbb{P}(Y\_0, \dots, Y\_K)

Does that sound right?  
If that is correct, we should be able to massage this model into an HMM, but we are still missing one specification: what is the distribution of the state evolution? In other words I do need

\mathbb{P}(X\_k = j \mid X\_{k-1} = i)

---

<div class="post-metadata">

**Author:** ![BioTurboNick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bioturbonick/32/6380_2.png) [@BioTurboNick](https://discourse.julialang.org/u/BioTurboNick)\
**Post date:** [February 26, 2024, 7:35pm UTC](https://discourse.julialang.org/t/ann-hmmbase-jl-a-lightweight-and-efficient-hidden-markov-model-abstraction/21604/18 "2024-02-26T19:35:23Z")

</div>

I will get back to you on this :-). Thanks for your help in the meantime! Translating to mathematical formalization is useful.