API reference

HiddenMarkovModels

A Julia package for HMM modeling, simulation, inference and learning.

Exports

• AbstractHMM
• AbstractHSMM
• ControlBoundEmission
• ControlBoundEmissionVector
• ControlledEmission
• ControlledEmissionHMM
• HMM
• HSMM
• baum_welch
• duration_distributions
• fit!
• forward
• forward_backward
• initialization
• joint_logdensityof
• logdensityof
• obs_distributions
• seq_limits
• transition_matrix
• viterbi

AbstractHMM

Abstract supertype for an HMM amenable to simulation, inference and learning.

Interface

To create your own subtype of AbstractHMM, you need to implement the following methods:

• initialization
• transition_matrix
• obs_distributions
• fit! (for learning)

Applicable functions

Any AbstractHMM which satisfies the interface can be given to the following functions:

• rand
• logdensityof
• forward
• viterbi
• forward_backward
• baum_welch (if [fit!](@ref) is implemented)

struct HMM{V<:(AbstractVector), M<:(AbstractMatrix), VD<:(AbstractVector), Vl<:(AbstractVector), Ml<:(AbstractMatrix)} <: AbstractHMM

Basic implementation of an HMM.

Fields

• init::AbstractVector: initial state probabilities

• trans::AbstractMatrix: state transition probabilities

• dists::AbstractVector: observation distributions

• loginit::AbstractVector: logarithms of initial state probabilities

• logtrans::AbstractMatrix: logarithms of state transition probabilities

struct ControlledEmissionHMM{V<:(AbstractVector), M<:(AbstractMatrix), VD<:(AbstractVector{<:ControlledEmission}), Vl<:(AbstractVector), Ml<:(AbstractMatrix)} <: AbstractHMM

An AbstractHMM where control variables affect only the emission distributions, not the transition dynamics. The init and trans fields are control-independent; each element of dists must be a ControlledEmission subtype, supporting the non-standard signatures logdensityof(dist, obs, control) and rand(rng, dist, control).

At inference time, obs_distributions(hmm, control) returns a lazy ControlBoundEmissionVector that binds each emission to control, exposing the standard logdensityof(dist, obs) / rand(rng, dist) interface expected by the inference algorithms.

A nothing control is forwarded to the emissions like any other control value; for genuinely uncontrolled models, prefer HMM.

• init::AbstractVector: initial state probabilities

• trans::AbstractMatrix: state transition probabilities (control-independent)

• dists::AbstractVector{<:ControlledEmission}: control-aware emission distributions; each is a ControlledEmission subtype supporting logdensityof(dist, obs, control) and rand(rng, dist, control)

• loginit::AbstractVector: logarithms of initial state probabilities

• logtrans::AbstractMatrix: logarithms of state transition probabilities

initialization(model)

Return the vector of initial state probabilities for model (any AbstractLatentStateModel).

transition_matrix(model)
transition_matrix(model, control)

Return the matrix of state transition probabilities for model (possibly when control is applied).

obs_distributions(model)
obs_distributions(model, control)

Return a vector of observation distributions, one for each state of model (possibly when control is applied).

These distribution objects should implement

• Random.rand(rng, dist) for sampling
• DensityInterface.logdensityof(dist, obs) for inference
• StatsAPI.fit!(dist, obs_seq, weight_seq) for learning

length(model)

Return the number of states of model (any AbstractLatentStateModel).

rand([rng,] hsmm::AbstractHSMM, T)
rand([rng,] hsmm::AbstractHSMM, control_seq)

Simulate hsmm for T time steps, or when the sequence control_seq is applied.

Sojourn lengths are drawn from the model's duration distributions (via rand_duration), so consecutive timesteps share the same state until the sojourn elapses, after which a new state is drawn under the no-self-transition constraint.

Return a named tuple (; state_seq, obs_seq, duration_seq), where duration_seq[t] is the sojourn length originally sampled for the segment containing t (it may exceed the remaining horizon if the sequence was truncated at T).

rand([rng,] hmm, T)
rand([rng,] hmm, control_seq)

Simulate hmm for T time steps, or when the sequence control_seq is applied.

Return a named tuple (; state_seq, obs_seq).

eltype(model, obs, control)

Return a type that can accommodate forward-backward computations for model on observations similar to obs.

It is typically a promotion between the element type of the initialization, the element type of the transition matrix, and the type of an observation logdensity evaluated at obs. Subtypes that carry additional numeric ingredients extend this promotion (e.g. AbstractHSMM folds in the duration logdensity type).

seq_limits(seq_ends, k)

Return a tuple (t1, t2) giving the begin and end indices of subsequence k within a set of sequences ending at seq_ends.

logdensityof(model)

Return the prior loglikelihood associated with the parameters of model.

logdensityof(hmm, obs_seq; ...)
logdensityof(hmm, obs_seq, control_seq; seq_ends)

Run the forward algorithm to compute the loglikelihood of obs_seq for hmm, integrating over all possible state sequences.

logdensityof(hsmm, obs_seq; ...)
logdensityof(
    hsmm,
    obs_seq,
    control_seq;
    seq_ends,
    max_duration
)

Run the forward algorithm to compute the loglikelihood of obs_seq for hsmm, integrating over all possible state sequences.

Unlike forward, this never throws when the loglikelihood is not finite, and returns -Inf instead.

max_duration limits the sojourn lengths considered, see forward. The default is exact but costs O(N T²) per sequence of length T, so consider a smaller value for long sequences.

joint_logdensityof(hmm, obs_seq, state_seq; ...)
joint_logdensityof(
    hmm,
    obs_seq,
    state_seq,
    control_seq;
    seq_ends
)

Run the forward algorithm to compute the the joint loglikelihood of obs_seq and state_seq for hmm.

joint_logdensityof(hsmm, obs_seq, state_seq; ...)
joint_logdensityof(
    hsmm,
    obs_seq,
    state_seq,
    control_seq;
    seq_ends
)

Compute the joint loglikelihood of obs_seq and state_seq for hsmm.

Each constant-state segment contributes its sojourn length to its duration distribution. The control at the segment's first timestep determines its duration distribution, and the control at the first timestep of the next segment determines the transition into it.

The final segment is right-censored and therefore contributes duration_logsurvival rather than duration_logdensityof.

forward(hmm, obs_seq; ...)
forward(
    hmm,
    obs_seq,
    control_seq;
    seq_ends,
    error_if_not_finite
)

Apply the forward algorithm to infer the current state after sequence obs_seq for hmm.

Return a tuple (storage.α, storage.logL) where storage is of type ForwardStorage.

forward(hsmm, obs_seq; ...)
forward(
    hsmm,
    obs_seq,
    control_seq;
    seq_ends,
    max_duration,
    error_if_not_finite
)

Apply the forward algorithm to infer the current state after sequence obs_seq for hsmm. Uses the explicit-duration formulation of (Yu, 2010), Section 3.1.

Return a tuple (storage.α, storage.logL) where storage is of type HSMMForwardStorage.

storage.α[:, t] contains the filtered state marginals. The segment covering t is treated as right-censored.

max_duration limits the sojourn lengths considered. It defaults to the longest sequence, which gives the exact result but costs O(N T²) per sequence of length T. For long sequences, a smaller value brings this down to O(N T max_duration), at the price of underestimating the loglikelihood.

viterbi(hmm, obs_seq; ...)
viterbi(hmm, obs_seq, control_seq; seq_ends)

Apply the Viterbi algorithm to infer the most likely state sequence corresponding to obs_seq for hmm.

Return a tuple (storage.q, storage.logL) where storage is of type ViterbiStorage.

forward_backward(hmm, obs_seq; ...)
forward_backward(hmm, obs_seq, control_seq; seq_ends)

Apply the forward-backward algorithm to infer the posterior state and transition marginals during sequence obs_seq for hmm.

Return a tuple (storage.γ, storage.logL) where storage is of type ForwardBackwardStorage.

baum_welch(hmm_guess, obs_seq; ...)
baum_welch(
    hmm_guess,
    obs_seq,
    control_seq;
    seq_ends,
    atol,
    max_iterations,
    loglikelihood_increasing,
    progress
)

Apply the Baum-Welch algorithm to estimate the parameters of an HMM on obs_seq, starting from hmm_guess.

Return a tuple (hmm_est, loglikelihood_evolution) where hmm_est is the estimated HMM and loglikelihood_evolution is a vector of loglikelihood values, one per iteration of the algorithm.

Keyword arguments

• atol: minimum loglikelihood increase at an iteration of the algorithm (otherwise the algorithm is deemed to have converged)
• max_iterations: maximum number of iterations of the algorithm
• loglikelihood_increasing: whether to throw an error if the loglikelihood decreases
• progress: whether to emit progress log records (rendered as a bar in REPL with TerminalLoggers, natively in Pluto, Jupyter and VS Code); off by default to keep the loop allocation-free

StatsAPI.fit!(
    hmm, fb_storage::ForwardBackwardStorage,
    obs_seq, [control_seq]; seq_ends,
)

Update hmm in-place based on information generated during forward-backward.

This function is allowed to reuse fb_storage as a scratch space, so its contents should not be trusted afterwards.

struct ForwardStorage{R}

Fields

Only the fields with a description are part of the public API.

• α::Matrix: posterior last state marginals α[i] = ℙ(X[T]=i | Y[1:T])

• logL::Vector: one loglikelihood per observation sequence

• B::Matrix

• c::Vector

initialize_forward(hmm, obs_seq, control_seq; seq_ends)

initialize_forward(
    hsmm,
    obs_seq,
    control_seq;
    seq_ends,
    max_duration
)

forward!(
    storage,
    hmm,
    obs_seq,
    control_seq;
    seq_ends,
    error_if_not_finite
)

forward!(
    storage,
    hsmm,
    obs_seq,
    control_seq;
    seq_ends,
    error_if_not_finite
)

struct HSMMForwardStorage{R}

The per-sequence buffers are sized from the seq_ends given to initialize_forward, so forward! must be called with the same seq_ends.

Fields

Only the fields with a description are part of the public API.

• α::Matrix: filtered state marginals α[i,t] = ℙ(X[t]=i | Y[1:t]), segment at t right-censored

• logL::Vector: one loglikelihood per observation sequence

• max_duration::Int64: longest sojourn duration considered (longer sojourns are truncated); each sequence uses min(max_duration, sequence length)

• log_ends::Matrix

• log_ongoing::Matrix

• log_prefix::Vector

• cum_log_obs::Matrix

• obs_zeros::Matrix{Int64}

• log_dur::Array{Matrix{R}, 1} where R

• log_surv::Array{Matrix{R}, 1} where R

• incoming::Array{Vector{R}, 1} where R

struct ViterbiStorage{R}

Fields

Only the fields with a description are part of the public API.

• q::Vector{Int64}: most likely state sequence q[t] = argmaxᵢ ℙ(X[t]=i | Y[1:T])

• logL::Vector: one joint loglikelihood per pair of observation sequence and most likely state sequence

• logB::Matrix

• ϕ::Matrix

• ψ::Matrix{Int64}

initialize_viterbi(hmm, obs_seq, control_seq; seq_ends)

viterbi!(storage, hmm, obs_seq, control_seq; seq_ends)

struct ForwardBackwardStorage{R, M<:AbstractArray{R, 2}}

Fields

Only the fields with a description are part of the public API.

• γ::Matrix: posterior state marginals γ[i,t] = ℙ(X[t]=i | Y[1:T])

• ξ::Vector{M} where {R, M<:AbstractMatrix{R}}: posterior transition marginals ξ[t][i,j] = ℙ(X[t]=i, X[t+1]=j | Y[1:T])

• logL::Vector: one loglikelihood per observation sequence

• B::Matrix

• α::Matrix

• c::Vector

• β::Matrix

• Bβ::Matrix

initialize_forward_backward(
    hmm,
    obs_seq,
    control_seq;
    seq_ends,
    transition_marginals
)

forward_backward!(
    storage,
    hmm,
    obs_seq,
    control_seq;
    seq_ends,
    transition_marginals
)

baum_welch!(
    fb_storage,
    logL_evolution,
    hmm,
    obs_seq,
    control_seq;
    seq_ends,
    atol,
    max_iterations,
    loglikelihood_increasing,
    progress
)

valid_hmm(hmm)

Perform some checks to rule out obvious inconsistencies with an AbstractHMM object.

Dispatches on AbstractLatentStateModel so the same structural checks (probability initialization, stochastic transition matrix, matching dimensions, density-bearing observation distributions) back both valid_hmm and valid_hsmm.

rand_prob_vec([rng, ::Type{R},] N)

Generate a random probability distribution of size N with normalized uniform entries.

rand_trans_mat([rng, ::Type{R},] N)

Generate a random transition matrix of size (N, N) with normalized uniform entries.

fit_in_sequence!(dists, i, x, w)

Modify the i-th element of dists by fitting it to an observation sequence x with associated weight sequence w.

Default behavior:

fit!(dists[i], x, w)

Override for Distributions.jl (in the package extension)

dists[i] = fit(eltype(dists), turn_into_vector(x), w)

struct LightDiagNormal{T1, T2, T3, V1<:AbstractArray{T1, 1}, V2<:AbstractArray{T2, 1}, V3<:AbstractArray{T3, 1}}

An HMMs-compatible implementation of a multivariate normal distribution with diagonal covariance, enabling allocation-free in-place estimation.

This is not part of the public API and is expected to change.

Fields

• μ::AbstractVector: means

• σ::AbstractVector: standard deviations

• logσ::AbstractVector: log standard deviations

struct LightCategorical{T1, T2, V1<:AbstractArray{T1, 1}, V2<:AbstractArray{T2, 1}}

An HMMs-compatible implementation of a discrete categorical distribution, enabling allocation-free in-place estimation.

This is not part of the public API and is expected to change.

Fields

• p::AbstractVector: class probabilities

• logp::AbstractVector: log class probabilities

log_initialization(model)

Return the vector of initial state log-probabilities for model.

Falls back on initialization.

log_transition_matrix(model)
log_transition_matrix(model, control)

Return the matrix of state transition log-probabilities for model (possibly when control is applied).

Falls back on transition_matrix, allocating a new matrix on every call. Inference routines such as Viterbi call it at every time step, even for models without control, so override it when this matters.

mul_rows_cols!(B, l, A, r)

Perform the in-place operation B .= l .* A .* transpose(r).

argmaxplus_transmul!(y, ind, A, x)

Perform the in-place multiplication transpose(A) * x in the sense of max-plus algebra, store the result in y, and store the index of the maximum for each component of y in ind.

abstract type ControlledEmission

Abstract supertype for control-aware emission distributions used by ControlledEmissionHMM.

A concrete subtype D <: ControlledEmission describes an emission whose density, sampling and fitting all depend on an external control value. It must implement the three non-standard, control-aware methods:

• DensityInterface.logdensityof(d::D, obs, control) for inference
• Random.rand(rng::AbstractRNG, d::D, control) for sampling
• StatsAPI.fit!(d::D, obs_seq, control_seq, weights) for learning

Subtyping ControlledEmission makes these requirements explicit and provides DensityInterface.DensityKind(::D) = HasDensity() automatically.

struct ControlBoundEmission{D<:ControlledEmission, C}

A single ControlledEmission bound to a specific control value. Wraps a control-aware emission together with one control, exposing the standard logdensityof(_, obs) and rand(rng, _) interface expected by the inference algorithms.

• dist::ControlledEmission: Control-aware emission distribution (a ControlledEmission subtype).

• control::Any: The control value bound to this emission; passed as the third argument to logdensityof and rand.

struct ControlBoundEmissionVector{D<:ControlledEmission, VD<:AbstractArray{D<:ControlledEmission, 1}, C} <: AbstractArray{ControlBoundEmission{D<:ControlledEmission, C}, 1}

A lazy vector pairing each ControlledEmission in dists with the same control value, yielding ControlBoundEmission elements. No allocation occurs until individual elements are accessed.

This is returned by obs_distributions(hmm::ControlledEmissionHMM, control) so that the standard contract is satisfied: the returned distributions are already bound to control and respond to logdensityof(dist, obs) and rand(rng, dist).

• dists::AbstractVector{D} where D<:ControlledEmission: Vector of control-aware emission distributions (ControlledEmission subtypes).

• control::Any: The control value shared across all emissions; passed to each wrapped distribution.

duration_logdensityof(dist, k)

Log-probability that a state with sojourn-time distribution dist is held for exactly k consecutive timesteps.

dist is supplied by the user with support on {0, 1, 2, ...} and is interpreted as the law of (sojourn time - 1). The sojourn time itself lives on {1, 2, 3, ...}, so this returns logdensityof(dist, k - 1). Durations k < 1 fall outside the user's support and yield -Inf through dist (e.g. a Distributions.Distribution returns -Inf for a negative argument).

See also rand_duration and duration_logsurvival.

rand_duration(rng, dist)

Sample a sojourn time on {1, 2, 3, ...} from the duration distribution dist.

dist is interpreted as the law of (sojourn time - 1) (see duration_logdensityof), so this returns rand(rng, dist) + 1.

duration_logsurvival(dist, k)

Return log ℙ(sojourn time >= k) for the duration distribution dist, following the convention of duration_logdensityof (the sojourn time lives on {1, 2, 3, ...}).

This is the right-censoring term for the final segment of a sequence, whose sojourn is only known to last at least k timesteps. Since every sojourn lasts at least one timestep, the result must be 0 = log(1) for k <= 1.

A closed-form method based on logccdf is provided for every Distributions.DiscreteUnivariateDistribution when Distributions.jl is loaded. Custom duration distributions must implement this function themselves, for instance as

function HiddenMarkovModels.duration_logsurvival(dist::MyDuration, k::Integer)
    k <= 1 && return 0.0
    return log1p(-sum(exp(duration_logdensityof(dist, j)) for j in 1:(k - 1)))
end

(but taking more precautions with numerical instabilities).

StateSegments(state_seq, t1, t2)

Non-allocating iterator over the maximal constant-state runs of state_seq[t1:t2].

Each element is a tuple (t_start, t_end) such that state_seq[t_start:t_end] is constant and differs from its neighbors within t1:t2. Consecutive segments tile t1:t2 in order.

AbstractLatentStateModel

Abstract supertype for a discrete-latent-state sequence model amenable to simulation, inference and learning.

This type factors out the methods shared by AbstractHMM and AbstractHSMM (state initialization, transition matrix, observation distributions, and the derived length, eltype, log-accessors and no-control fallbacks).

Interface

Every subtype must implement:

• initialization
• transition_matrix
• obs_distributions
• fit! (for learning)

Subtypes may add further requirements (e.g. AbstractHSMM additionally requires duration_distributions).

AbstractHSMM

Abstract supertype for a Hidden Semi-Markov Model amenable to simulation, inference and learning.

An AbstractHSMM explicitly models state durations through a per-state duration distribution. Sojourn lengths are drawn from those distributions, so self-transitions in the transition matrix are forbidden (the diagonal must be zero).

Interface

To create your own subtype of AbstractHSMM, you need to implement the following methods:

• initialization
• transition_matrix
• obs_distributions
• duration_distributions
• fit! (for learning)

Applicable functions

Any AbstractHSMM which satisfies the interface can be given to the following functions:

• rand
• logdensityof
• joint_logdensityof
• forward

struct HSMM{V<:(AbstractVector), M<:(AbstractMatrix), VD<:(AbstractVector), VDur<:(AbstractVector), Vl<:(AbstractVector), Ml<:(AbstractMatrix)} <: AbstractHSMM

Basic implementation of an HSMM with explicit state durations.

Fields

• init::AbstractVector: initial state probabilities

• trans::AbstractMatrix: state transition probabilities (zero diagonal, no self-transitions)

• dists::AbstractVector: observation distributions

• dur_dists::AbstractVector: state duration distributions, interpreted as the law of (sojourn - 1)

• loginit::AbstractVector: logarithms of initial state probabilities

• logtrans::AbstractMatrix: logarithms of state transition probabilities

duration_distributions(model)
duration_distributions(model, control)

Return a vector of sojourn-time distributions, one for each state of model (possibly when control is applied).

Following the convention of duration_logdensityof and rand_duration, each distribution has support on {0, 1, 2, ...} and is interpreted as the law of (sojourn time - 1). Any density-bearing object works (e.g. a Distributions.Geometric or Distributions.Poisson); it should implement

• Random.rand(rng, dist) for sampling
• DensityInterface.logdensityof(dist, k) for inference
• duration_logsurvival(dist, k) for right-censoring the final segment of a sequence

The last method is already provided for every Distributions.DiscreteUnivariateDistribution.

valid_hsmm(hsmm, control=nothing)

Perform some checks to rule out obvious inconsistencies with an AbstractHSMM object.

On top of the valid_hmm-style checks (probability initialization, stochastic transition matrix, matching dimensions, density-bearing observation distributions), this verifies:

• the transition matrix has a zero diagonal (no self-transitions);
• there is one density-bearing duration distribution per state.