Struct hmmm::HMM

pub struct HMM {
    pub a: Array2<f64>,
    pub b: Array2<f64>,
    pub pi: Array1<f64>,
}

This struct represents a trained HMM, including values for each parameter.

Math

The HMM is used to predict a sequence of observations:

$$Y=(Y_0=y_0, Y_1=y_1, \ldots, Y_{T-1}=y_{T-1})$$

...where each $y_t \in [0, K)$.

It accomplishes this with latent variables for hidden state $X=(X_0, \ldots, X_{T-1})$ where each $x_t \in [0, N)$.

A trained HMM has three parameters:

• $A$, the $N × N$ state transition matrix: $a_{ij}=P(X_t=j|X_{t-1}=i)$
• $B$, the $N × K$ observation matrix: $b_{ik}=P(Y_t=y_k|X_t=i)$
• $π$, the $N$-length initial state distribution: $π_i=P(X_1=i)$

Fields

a: Array2<f64>

b: Array2<f64>

pi: Array1<f64>

Implementations

impl HMM

pub fn new(a: Array2<f64>, b: Array2<f64>, pi: Array1<f64>) -> Self

Create a new HMM with the given parameters.

This could be useful for loading a saved trained model.

Panics if any of:

• Dimensions are invalid
• Probability distributions are invalid

pub fn n(&self) -> usize

$N$, the number of states in this HMM

pub fn k(&self) -> usize

$K$, the number of possible observations that this model can emit

pub fn sampler<'a, R: Rng + ?Sized>(
    &'a self,
    rng: &'a mut R
) -> HMMSampleIter<'_, R>

pub fn filter<I>(&self, ys: I) -> HMMFilterIter<'_, I::IntoIter> where
    I: IntoIterator<Item = usize>, 

Given an iterator of observations, this returns a new iterator that yields the probability of being in each hidden state at each future time step. This method is relatively efficient from the standpoint of memory and computation time.

If you can store the whole sequence in memory, there are more accurate ways to compute the probability of being in each state at a particular point in time, because it is possible to use the observations from the future to better inform the probability of being in each hidden state at any particular time t.

This is not closely related to the meaning of "filter" as in std::iter::Iterator::filter.

Panics if an observation is out of bounds.

pub fn predict(&self, p_states: Array1<f64>, n_time_steps: usize) -> Array1<f64>

Given a distribution over states, calculate the probable distribution of states at a time in the future.

This is currently only efficient for small values of n_time_steps. In order to be more efficient, we want to be able to efficiently raise self.a to a power, which means getting eigenvalues and eigenvectors.

Panics if:

• The length of p_states is invalid
• p_states is not probability distribution

pub fn smooth(&self, ys: &Array1<usize>) -> Array2<f64>

Given a sequence of observations, compute the probability of being in any given state at each point in time.

Return a $T × N$ matrix where element (t, k) is the probability that we are in state k at time t.

This is the forward-backward algorithm.

pub fn most_likely_sequence(&self, ys: &Array1<usize>) -> Array1<usize>

This is the Viterbi algorithm. Given a sequence of observations, return the most likely sequence of states.

It's possible to do this in log space but I normalized instead to make it feel more like the forwards-backwards algorithm.

pub fn train<R: Rng>(
    ys: &Array1<usize>,
    n: usize,
    k: usize,
    rng: &mut R
) -> Self

Find the maximum likelihood estimate for the parameters. Caveats:

• If there is not enough data, the MLE is undefined. This implementation will use a uniform prior for any parameters for which there isn't enough data.
• This is not guaranteed to find a global minimal, only a local minimum.
• Due to a lack of identifiability, an HMM with $N$ states has $N!$ equivalent solutions.
• Taking the most likely state at each point in time doesn't necessarily result in the most likely sequence of states, or even a possible sequence of states. If you want that, use HMM::most_likely_states.

Baum-Welch (Baum et. al. 1970) is a variant of the Expectation-Maximization algorithm for HMMs.

Let $α_i(t) = P(Y_0=y_0, \ldots, Y_t=y_t, X_t=i | θ)$

Let $β_i(t) = P(Y_{t+1}=y_{t+1}, \ldots, Y_T=y_T | X_t=i, θ)$

$$ γ_i(t) = P(X_t=i|Y,θ) = \frac{P(X_t=i,Y|θ)}{P(Y|θ)} = \frac{α_i(t)β_i(t)}{\sum_{j=1}^N α_j(t)β_j(t)} $$

$γ_i(t)$ is the same thing that the smooth method computes: the probability of being in state $i$ at time $t$, conditioned on both past and future observations.

$ξ_{ij}(t)$ is the probability of being in state $i$ at time $t$ and in state $j$ at time $t + 1$:

$$ ξ_{ij}(t) = P(X_t=i,X_{t+1}=j|Y,θ) = \frac{P(X_t=i,X_{t+1}=j,Y|θ)}{P(Y|θ)} = \frac{α_i(t) a_{ij} β_j(t+1) b_j(y_{t+1})} {\sum_{i=1}^N \sum_{j=1}^N α_i(t) a_{ij} β_j(t+1) b_j(y_{t+1}) } $$

Updating the counts:

$$ π_i^* = γ_i(0) $$

For $α$:

$$ α_{ij}^*=\frac{\sum^{T-1}_{t=1}ξ_{ij}(t)}{\sum^{T-1}_{t=1}γ_i(t)} $$

For $β$:

$$ β_i^*(v_k)=\frac{\sum^T_{t=1} 1_{y_t=v_k} γ_i(t)}{\sum^T_{t=1} γ_i(t)} $$

pub fn ll_given_states(&self, xs: &[usize], ys: &[usize]) -> f64

Return the log likelihood of a sequence of states and observations. This is not a typical use case, because often the vector of hidden states is not available.

Panics if:

• The number of states and observations is not equal
• A state or observation is out of bounds

Trait Implementations

impl Debug for HMM

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.