### Start Julia REPL for Comparison Environment

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/libs/HMMComparison/README.md

Launch a single-threaded Julia REPL within the specific project environment for running comparisons.

```bash
julia -t 1 --project=libs/HMMComparison
```

--------------------------------

### Install HiddenMarkovModels.jl

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/README.md

Install the HiddenMarkovModels package using Julia's package manager.

```julia
pkg> add HiddenMarkovModels
```

--------------------------------

### Implement Custom HMM with AbstractHMM Interface

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Create custom HMM types by implementing the `AbstractHMM` interface. This example shows a controlled HMM where observation distributions depend on a control input sequence.

```julia
using Distributions, HiddenMarkovModels, LinearAlgebra, StatsAPI
import HiddenMarkovModels as HMMs

# Define a controlled HMM where transitions depend on control input
struct ControlledGaussianHMM{T} <: AbstractHMM
    init::Vector{T}
    trans::Matrix{T}
    dist_coeffs::Vector{Vector{T}}  # Linear regression coefficients per state
end

# Required interface methods
HMMs.initialization(hmm::ControlledGaussianHMM) = hmm.init

function HMMs.transition_matrix(hmm::ControlledGaussianHMM, control::AbstractVector)
    return hmm.trans  # Could also depend on control
end

function HMMs.obs_distributions(hmm::ControlledGaussianHMM, control::AbstractVector)
    # Observation mean is a linear function of control
    return [Normal(dot(hmm.dist_coeffs[i], control), 1.0) for i in 1:length(hmm)]
end

# Create instance
d = 3  # Control dimension
init = [0.6, 0.4]
trans = [0.7 0.3; 0.2 0.8]
dist_coeffs = [-ones(d), ones(d)]  # State 1: negative response, State 2: positive
hmm = ControlledGaussianHMM(init, trans, dist_coeffs)

# Simulate with control sequence
control_seq = [randn(d) for _ in 1:100]
result = rand(hmm, control_seq)
state_seq, obs_seq = result.state_seq, result.obs_seq

# Inference with controls
best_states, logL = viterbi(hmm, obs_seq, control_seq)
smoothed, _ = forward_backward(hmm, obs_seq, control_seq)

println("Decoded states: ", best_states[1:10])
println("Joint logL: ", only(logL))
```

--------------------------------

### Get Sequence Limits from Concatenated Data

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

The `seq_limits` function helps identify the start and end indices of individual sequences within a single concatenated observation vector, which is useful for batch processing.

```julia
using Distributions, HiddenMarkovModels

# Setup
hmm = HMM([0.6, 0.4], [0.7 0.3; 0.2 0.8], [Normal(-1.0), Normal(1.0)])

# Generate multiple sequences of different lengths
obs_seqs = [last(rand(hmm, len)) for len in [50, 75, 100, 60]]

# Concatenate for batch processing
obs_concat = reduce(vcat, obs_seqs)
seq_ends = cumsum(length.(obs_seqs))
println("Sequence ends: ", seq_ends)  # [50, 125, 225, 285]

# Get limits for sequence k
k = 2
t1, t2 = seq_limits(seq_ends, k)
println("Sequence $k spans indices $t1 to $t2")  # 51 to 125

# Extract individual sequence results after batch inference
all_states, all_logL = viterbi(hmm, obs_concat; seq_ends)

# Get results for sequence k
states_k = all_states[t1:t2]
logL_k = all_logL[k]
println("Sequence $k has $(length(states_k)) states, logL = $(round(logL_k, digits=2))")

# Verify against individual inference
states_individual, logL_individual = viterbi(hmm, obs_seqs[k])
println("Individual inference matches: ", states_k == states_individual)

# Iterate over all sequences
for k in eachindex(seq_ends)
    t1, t2 = seq_limits(seq_ends, k)
    println("Sequence $k: indices $t1:$t2, length $(t2-t1+1), logL $(round(all_logL[k], digits=2))")
end
```

--------------------------------

### Create a Basic HMM in Julia

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/README.md

Demonstrates how to initialize a Hidden Markov Model (HMM) with initial state probabilities, transition matrix, and emission distributions.

```julia
using Distributions, HiddenMarkovModels
init = [0.6, 0.4]
trans = [0.7 0.3; 0.2 0.8]
dists = [Normal(-1.0), Normal(1.0)]
hmm = HMM(init, trans, dists)
```

--------------------------------

### In-Place HMM Operations

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Illustrates setting up an HMM and generating data, preparing for performance-critical applications that require allocation-free in-place operations.

```julia
using Distributions, HiddenMarkovModels
import HiddenMarkovModels as HMMs

# Setup
init = [0.6, 0.4]
trans = [0.7 0.3; 0.2 0.8]
dists = [Normal(-1.0), Normal(1.0)]
hmm = HMM(init, trans, dists)

# Generate data
_, obs_seq = rand(hmm, 1000)
control_seq = fill(nothing, length(obs_seq))
seq_ends = (length(obs_seq),)

```

--------------------------------

### Initialize and Run HMM Algorithms In-Place

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Pre-allocate storage for HMM algorithms once, then run them in-place to avoid further allocations. This is useful for streaming or online inference.

```julia
forward_storage = HMMs.initialize_forward(hmm, obs_seq, control_seq; seq_ends)
viterbi_storage = HMMs.initialize_viterbi(hmm, obs_seq, control_seq; seq_ends)
fb_storage = HMMs.initialize_forward_backward(hmm, obs_seq, control_seq; seq_ends)

HMMs.forward!(forward_storage, hmm, obs_seq, control_seq; seq_ends)
println("Forward logL: ", sum(forward_storage.logL))

HMMs.viterbi!(viterbi_storage, hmm, obs_seq, control_seq; seq_ends)
println("Viterbi logL: ", sum(viterbi_storage.logL))

HMMs.forward_backward!(fb_storage, hmm, obs_seq, control_seq; seq_ends)
println("FB logL: ", sum(fb_storage.logL))
```

--------------------------------

### Create a Hidden Markov Model (HMM)

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Defines the core HMM data structure with initial probabilities, transition matrix, and observation distributions. Supports various observation types like Normal and Multivariate Normal, and sparse transition matrices.

```julia
using Distributions, HiddenMarkovModels

# Create a simple 2-state HMM with Normal observation distributions
init = [0.6, 0.4]                           # Initial state probabilities
trans = [0.7 0.3; 0.2 0.8]                  # State transition matrix
dists = [Normal(-1.0), Normal(1.0)]         # Observation distributions per state
hmm = HMM(init, trans, dists)

# Output:
# Hidden Markov Model with:
#  - initialization: [0.6, 0.4]
#  - transition matrix: [0.7 0.3; 0.2 0.8]
#  - observation distributions: [Normal{Float64}(μ=-1.0, σ=1.0), Normal{Float64}(μ=1.0, σ=1.0)]
```

```julia
# HMM with multivariate Normal observations
using LinearAlgebra
init = [0.6, 0.4]
trans = [0.7 0.3; 0.2 0.8]
dists = [MvNormal([-0.5, -0.8], I), MvNormal([0.5, 0.8], I)]
hmm_mv = HMM(init, trans, dists)
```

```julia
# HMM with sparse transition matrix for large state spaces
using SparseArrays
trans_sparse = sparse([
    0.7 0.3 0
    0 0.7 0.3
    0.3 0 0.7
])
init_sparse = [0.2, 0.6, 0.2]
dists_sparse = [Normal(1.0), Normal(2.0), Normal(3.0)]
hmm_sparse = HMM(init_sparse, trans_sparse, dists_sparse)
```

--------------------------------

### Define Custom Observation Distribution

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Demonstrates how to create a custom observation distribution `StuffDist` by implementing `rand`, `DensityInterface.DensityKind`, `DensityInterface.logdensityof`, and `StatsAPI.fit!`. This allows HMMs to work with user-defined data types.

```julia
using DensityInterface, HiddenMarkovModels, Random, StatsAPI

# Custom observation type
struct Stuff{T}
    quantity::T
end

# Custom distribution
mutable struct StuffDist{T}
    quantity_mean::T
end

# Required: sampling
function Random.rand(rng::Random.AbstractRNG, dist::StuffDist)
    quantity = dist.quantity_mean + randn(rng)
    return Stuff(quantity)
end

# Required: declare density exists
DensityInterface.DensityKind(::StuffDist) = HasDensity()

# Required: log-density computation
function DensityInterface.logdensityof(dist::StuffDist, obs::Stuff)
    return -abs2(obs.quantity - dist.quantity_mean) / 2
end

# Optional but needed for Baum-Welch: in-place fitting
function StatsAPI.fit!(
    dist::StuffDist,
    obs_seq::AbstractVector{<:Stuff},
    weight_seq::AbstractVector{<:Real}
)
    dist.quantity_mean = sum(w * o.quantity for (w, o) in zip(weight_seq, obs_seq)) / sum(weight_seq)
    return nothing
end
```

```julia
# Use custom distribution in HMM
init = [0.6, 0.4]
trans = [0.7 0.3; 0.2 0.8]
dists = [StuffDist(-1.0), StuffDist(+1.0)]
hmm = HMM(init, trans, dists)

# All algorithms work with custom distributions
state_seq, obs_seq = rand(hmm, 100)
println("Observation type: ", typeof(obs_seq[1]))  # Stuff{Float64}

best_states, _ = viterbi(hmm, obs_seq)
smoothed, _ = forward_backward(hmm, obs_seq)

# Learning also works
hmm_guess = HMM([0.5, 0.5], [0.6 0.4; 0.3 0.7], [StuffDist(-0.5), StuffDist(0.5)])
hmm_est, logL_evo = baum_welch(hmm_guess, obs_seq)
println("Estimated means: ", [d.quantity_mean for d in obs_distributions(hmm_est)])
```

--------------------------------

### Include Experiment Scripts in Julia

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/libs/HMMComparison/README.md

Execute the necessary Julia scripts to run the measurements and generate plots for the experiments.

```julia
include("libs/HMMComparison/experiments/measurements.jl")
include("libs/HMMComparison/experiments/plots.jl")
```

--------------------------------

### Simulate State and Observation Sequences with `rand`

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Generates state and observation sequences from a given HMM for a specified number of time steps. Supports custom random number generators for reproducibility and handles multivariate observations.

```julia
using Distributions, HiddenMarkovModels, Random

# Create HMM
init = [0.6, 0.4]
trans = [0.7 0.3; 0.2 0.8]
dists = [Normal(-1.0), Normal(1.0)]
hmm = HMM(init, trans, dists)

# Simulate for T time steps
T = 100
state_seq, obs_seq = rand(hmm, T)

# state_seq is a Vector{Int} of hidden states
println("First 5 states: ", state_seq[1:5])    # e.g., [1, 1, 2, 2, 2]

# obs_seq is a Vector matching the output type of the distributions
println("First 5 observations: ", round.(obs_seq[1:5], digits=2))  # e.g., [-0.87, -1.23, 0.95, 1.12, 0.78]
```

```julia
# Simulate with a specific RNG for reproducibility
using StableRNGs
rng = StableRNG(42)
state_seq, obs_seq = rand(rng, hmm, T)
```

```julia
# For multivariate observations, each element is a vector
hmm_mv = HMM([0.5, 0.5], [0.8 0.2; 0.3 0.7], [MvNormal([0,0], I), MvNormal([3,3], I)])
state_seq_mv, obs_seq_mv = rand(hmm_mv, 50)
println("First observation vector: ", obs_seq_mv[1])  # e.g., [0.12, -0.45]
```

--------------------------------

### Estimate HMM Parameters with Baum-Welch

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Use the Baum-Welch algorithm to estimate HMM parameters from observation sequences. It can handle single or multiple sequences and allows for custom convergence criteria.

```julia
using Distributions, HiddenMarkovModels

# True model
init_true = [0.6, 0.4]
trans_true = [0.7 0.3; 0.2 0.8]
dists_true = [Normal(-1.0), Normal(1.0)]
hmm_true = HMM(init_true, trans_true, dists_true)

# Generate training data
_, obs_seq = rand(hmm_true, 1000)

# Initial guess (must be close enough for EM to converge)
init_guess = [0.5, 0.5]
trans_guess = [0.6 0.4; 0.3 0.7]
dists_guess = [Normal(-0.5), Normal(0.5)]
hmm_guess = HMM(init_guess, trans_guess, dists_guess)

# Run Baum-Welch
hmm_est, logL_evolution = baum_welch(hmm_guess, obs_seq)

# Check convergence
println("Initial logL: ", round(first(logL_evolution), digits=2))
println("Final logL: ", round(last(logL_evolution), digits=2))
println("Iterations: ", length(logL_evolution))

# Compare estimated vs true parameters
println("\nTransition matrix:")
println("True: ", trans_true)
println("Estimated: ", round.(transition_matrix(hmm_est), digits=3))

println("\nObservation means:")
println("True: ", [mean(d) for d in dists_true])
println("Estimated: ", [round(mean(d), digits=3) for d in obs_distributions(hmm_est)])

# Baum-Welch with multiple sequences (more robust estimation)
obs_seqs = [last(rand(hmm_true, rand(100:200))) for _ in 1:100]
obs_concat = reduce(vcat, obs_seqs)
seq_ends = cumsum(length.(obs_seqs))

hmm_est_multi, _ = baum_welch(hmm_guess, obs_concat; seq_ends)
println("\nMulti-sequence estimation:")
println("Transition: ", round.(transition_matrix(hmm_est_multi), digits=3))

# Custom convergence criteria
hmm_est_strict, _ = baum_welch(
    hmm_guess, obs_seq;
    atol=1e-8,              # Stricter convergence threshold
    max_iterations=500,     # Allow more iterations
    loglikelihood_increasing=true  # Error if logL decreases
)
```

--------------------------------

### Learning Algorithms

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/docs/src/api.md

Functions for training and learning parameters of Hidden Markov Models.

```APIDOC
## Learning

```@docs
baum_welch
fit!
```
```

--------------------------------

### Sequence Formatting

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/docs/src/api.md

Details on how to format observation and control sequences for HMM algorithms. Supports single and multiple sequences.

```APIDOC
## Sequence Formatting

Most algorithms ingest data with `obs_seq` (mandatory) and `control_seq` (optional), and a keyword argument `seq_ends` (optional).

- **Single Sequence**: If data is a single sequence, `obs_seq` and `control_seq` are the corresponding vectors. `seq_ends` is not needed.
- **Multiple Sequences**: If data consists of multiple sequences, `obs_seq` and `control_seq` are concatenations of several vectors. The end indices of these sequences are provided by `seq_ends`.

To prepare multiple sequences:
```julia
obs_seq = reduce(vcat, obs_seqs)
control_seq = reduce(vcat, control_seqs)
seq_ends = cumsum(length.(obs_seqs))
```
```

--------------------------------

### Differentiable HMM with ForwardDiff and Zygote

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Sets up a `DiffusionHMM` struct and demonstrates automatic differentiation of its log-likelihood using `ForwardDiff.jl` for forward-mode and `Zygote.jl` for reverse-mode gradients.

```julia
using ComponentArrays, Distributions, HiddenMarkovModels, ForwardDiff, Zygote
import HiddenMarkovModels as HMMs

# Define a differentiable HMM
struct DiffusionHMM{V1,M2,V3} <: AbstractHMM
    init::V1
    trans::M2
    means::V3
end

HMMs.initialization(hmm::DiffusionHMM) = hmm.init

function HMMs.transition_matrix(hmm::DiffusionHMM, λ::Number)
    N = length(hmm)
    return (1 - λ) * hmm.trans + λ * ones(N, N) / N
end

function HMMs.obs_distributions(hmm::DiffusionHMM, λ::Number)
    return [Normal((1 - λ) * hmm.means[i]) for i in 1:length(hmm)]
end

# Setup
init = [0.6, 0.4]
trans = [0.7 0.3; 0.3 0.7]
means = [-1.0, 1.0]
hmm = DiffusionHMM(init, trans, means)

# Generate data with controls between 0 and 1
control_seq = rand(50)
obs_seq = rand(hmm, control_seq).obs_seq
seq_ends = (length(obs_seq),)

# Pack parameters for differentiation
parameters = ComponentVector(; init, trans, means)

function f(params, obs, ctrl; seq_ends)
    new_hmm = DiffusionHMM(params.init, params.trans, params.means)
    return logdensityof(new_hmm, obs, ctrl; seq_ends)
end

# Forward-mode differentiation (ForwardDiff.jl)
grad_params = ForwardDiff.gradient(p -> f(p, obs_seq, control_seq; seq_ends), parameters)
println("Gradient w.r.t. parameters: ", round.(grad_params, digits=4))

grad_obs = ForwardDiff.gradient(o -> f(parameters, o, control_seq; seq_ends), obs_seq)
println("Gradient w.r.t. observations (first 5): ", round.(grad_obs[1:5], digits=4))

# Reverse-mode differentiation (Zygote.jl) - more efficient for many parameters
grads = Zygote.gradient((p, o, c) -> f(p, o, c; seq_ends), parameters, obs_seq, control_seq)
grad_params_zygote, grad_obs_zygote, grad_ctrl_zygote = grads

println("\nZygote gradients match ForwardDiff: ", grad_params_zygote ≈ grad_params)
```

--------------------------------

### Clone HiddenMarkovModels.jl Repository

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/libs/HMMComparison/README.md

Use this command to clone the project repository from GitHub.

```bash
git clone https://github.com/gdalle/HiddenMarkovModels.jl
cd HiddenMarkovModels.jl
```

--------------------------------

### In-place Versions of Algorithms

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/docs/src/api.md

In-place versions of common HMM algorithms for memory efficiency.

```APIDOC
## In-place versions

### Forward

```@docs
HiddenMarkovModels.ForwardStorage
HiddenMarkovModels.initialize_forward
HiddenMarkovModels.forward!
```

### Viterbi

```@docs
HiddenMarkovModels.ViterbiStorage
HiddenMarkovModels.initialize_viterbi
HiddenMarkovModels.viterbi!
```

### Forward-backward

```@docs
HiddenMarkovModels.ForwardBackwardStorage
HiddenMarkovModels.initialize_forward_backward
HiddenMarkovModels.forward_backward!
```

### Baum-Welch

```@docs
HiddenMarkovModels.baum_welch!
```
```

--------------------------------

### Access HMM Algorithm Results from Storage

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Retrieve results such as filtered marginals, smoothed marginals, and the best state sequence from the storage objects after running the HMM algorithms.

```julia
filtered_marginals = forward_storage.α
smoothed_marginals = fb_storage.γ
best_states = viterbi_storage.q

println("Last filtered state: ", filtered_marginals[:, end])
println("Last smoothed state: ", smoothed_marginals[:, end])
println("Viterbi final state: ", best_states[end])
```

--------------------------------

### Compute Observation Sequence Log-Likelihood with `logdensityof`

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Calculates the log-likelihood of an observation sequence by integrating over all possible state sequences. Useful for comparing different HMM models. Requires `Distributions` and `HiddenMarkovModels`.

```julia
using Distributions, HiddenMarkovModels

# Setup
init = [0.6, 0.4]
trans = [0.7 0.3; 0.2 0.8]
dists = [Normal(-1.0), Normal(1.0)]
hmm = HMM(init, trans, dists)
_, obs_seq = rand(hmm, 100)

# Compute log-likelihood P(Y_1:T)
logL = logdensityof(hmm, obs_seq)
println("Log-likelihood: ", logL)  # e.g., -145.23

# Compare models using log-likelihood
hmm_wrong = HMM([0.5, 0.5], [0.5 0.5; 0.5 0.5], [Normal(0.0), Normal(0.0)])
logL_correct = logdensityof(hmm, obs_seq)
logL_wrong = logdensityof(hmm_wrong, obs_seq)
println("Correct model logL: ", round(logL_correct, digits=2))
println("Wrong model logL: ", round(logL_wrong, digits=2))
# Correct model should have higher log-likelihood

# Multiple sequences
obs_seqs = [last(rand(hmm, rand(50:100))) for _ in 1:10]
obs_concat = reduce(vcat, obs_seqs)
seq_ends = cumsum(length.(obs_seqs))
total_logL = logdensityof(hmm, obs_concat; seq_ends)
println("Total log-likelihood for 10 sequences: ", round(total_logL, digits=2))
```

--------------------------------

### Fit Controlled Gaussian HMM

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Implements the `fit!` function for ControlledGaussianHMM to update transition and observation parameters using provided observation and control sequences. Requires `StatsAPI`.

```julia
function StatsAPI.fit!(
    hmm::ControlledGaussianHMM{T},
    fb_storage::HMMs.ForwardBackwardStorage,
    obs_seq::AbstractVector,
    control_seq::AbstractVector;
    seq_ends
) where {T}
    (; γ, ξ) = fb_storage
    N = length(hmm)

    # Fit transition parameters
    hmm.init .= 0
    hmm.trans .= 0
    for k in eachindex(seq_ends)
        t1, t2 = HMMs.seq_limits(seq_ends, k)
        hmm.init .+= γ[:, t1]
        hmm.trans .+= sum(ξ[t1:t2])
    end
    hmm.init ./= sum(hmm.init)
    for row in eachrow(hmm.trans)
        row ./= sum(row)
    end

    # Fit observation coefficients via weighted least squares
    U = reduce(hcat, control_seq)'
    y = obs_seq
    for i in 1:N
        W = sqrt.(Diagonal(γ[i, :]))
        hmm.dist_coeffs[i] = (W * U) \ (W * y)
    end
end
```

--------------------------------

### Core Types

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/docs/src/api.md

Documentation for the core abstract types and concrete types used in HiddenMarkovModels.jl.

```APIDOC
## Types

```@docs
AbstractHMM
HMM
```
```

--------------------------------

### Compute Smoothed State Marginals with `forward_backward`

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Calculates the probability of each state given ALL observations, providing more accurate state inference. Compares filtered vs. smoothed marginals. Requires `Distributions` and `HiddenMarkovModels`.

```julia
using Distributions, HiddenMarkovModels

# Setup
init = [0.6, 0.4]
trans = [0.7 0.3; 0.2 0.8]
dists = [Normal(-1.0), Normal(1.0)]
hmm = HMM(init, trans, dists)
_, obs_seq = rand(hmm, 100)

# Run forward-backward algorithm
smoothed_marginals, logL = forward_backward(hmm, obs_seq)

# smoothed_marginals[i, t] = P(X_t = i | Y_1:T)
# Uses ALL observations, not just up to time t
println("Size of smoothed marginals: ", size(smoothed_marginals))  # (2, 100)
println("Log-likelihood: ", only(logL))

# Compare filtered vs smoothed at an early time point
filtered, _ = forward(hmm, obs_seq)
println("Filtered P(X_10 | Y_1:10): ", round.(filtered[:, 10], digits=3))
println("Smoothed P(X_10 | Y_1:T): ", round.(smoothed_marginals[:, 10], digits=3))

# Note: At the final time step, filtered and smoothed are identical
println("Filtered at T: ", round.(filtered[:, end], digits=3))
println("Smoothed at T: ", round.(smoothed_marginals[:, end], digits=3))

# Most likely state at each time (soft assignment)
most_likely_states = [argmax(smoothed_marginals[:, t]) for t in 1:size(smoothed_marginals, 2)]
```

--------------------------------

### Interface Functions

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/docs/src/api.md

Functions that define the interface for interacting with HMM objects.

```APIDOC
## Interface

```@docs
initialization
transition_matrix
obs_distributions
```
```

--------------------------------

### Inference Algorithms

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/docs/src/api.md

Functions for performing inference on Hidden Markov Models.

```APIDOC
## Inference

```@docs
logdensityof
joint_logdensityof
forward
viterbi
forward_backward
```
```

--------------------------------

### Utility Functions

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/docs/src/api.md

General utility functions for HMMs, including length, random generation, and element type.

```APIDOC
## Utils

```@docs
length
rand
eltype
seq_limits
```
```

--------------------------------

### Compute Filtered State Marginals with `forward`

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Calculates the probability of each state given observations up to the current time. Useful for online inference. Requires `Distributions` and `HiddenMarkovModels` packages.

```julia
using Distributions, HiddenMarkovModels

# Setup
init = [0.6, 0.4]
trans = [0.7 0.3; 0.2 0.8]
dists = [Normal(-1.0), Normal(1.0)]
hmm = HMM(init, trans, dists)
_, obs_seq = rand(hmm, 100)

# Run forward algorithm
filtered_marginals, logL = forward(hmm, obs_seq)

# filtered_marginals[i, t] = P(X_t = i | Y_1:t)
# Matrix of size (num_states, T)
println("Size of marginals: ", size(filtered_marginals))  # (2, 100)

# logL is a Vector containing one log-likelihood per sequence
println("Observation log-likelihood: ", only(logL))  # e.g., -143.27

# Final state distribution given all observations
println("P(X_T | Y_1:T): ", filtered_marginals[:, end])  # e.g., [0.23, 0.77]

# Filtered marginals at time t only use observations up to t
# Useful for online/streaming inference
for t in [1, 10, 50, 100]
    println("P(X_$t | Y_1:$t): ", round.(filtered_marginals[:, t], digits=3))
end
```

--------------------------------

### Miscellaneous Functions

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/docs/src/api.md

Various helper and miscellaneous functions for HMMs.

```APIDOC
## Miscellaneous

```@docs
HiddenMarkovModels.valid_hmm
HiddenMarkovModels.rand_prob_vec
HiddenMarkovModels.rand_trans_mat
HiddenMarkovModels.fit_in_sequence!
```
```

--------------------------------

### Compute Joint Log-Likelihood with `joint_logdensityof`

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Calculates the joint log-likelihood of observations and a specific state sequence. Useful for comparing true vs. inferred state sequences (e.g., Viterbi). Requires `Distributions` and `HiddenMarkovModels`.

```julia
using Distributions, HiddenMarkovModels

# Setup
init = [0.6, 0.4]
trans = [0.7 0.3; 0.2 0.8]
dists = [Normal(-1.0), Normal(1.0)]
hmm = HMM(init, trans, dists)
true_state_seq, obs_seq = rand(hmm, 100)

# Compute joint log-likelihood P(X_1:T, Y_1:T)
joint_logL = joint_logdensityof(hmm, obs_seq, true_state_seq)
println("Joint log-likelihood: ", joint_logL)  # e.g., -156.34

# Viterbi finds the state sequence with highest joint log-likelihood
viterbi_states, viterbi_logL = viterbi(hmm, obs_seq)
joint_logL_viterbi = joint_logdensityof(hmm, obs_seq, viterbi_states)
joint_logL_true = joint_logdensityof(hmm, obs_seq, true_state_seq)

println("Viterbi sequence joint logL: ", round(joint_logL_viterbi, digits=2))
println("True sequence joint logL: ", round(joint_logL_true, digits=2))
# Viterbi sequence should have >= joint log-likelihood (it's the argmax)

# Relationship: P(Y) = sum over X of P(X, Y)
# So: logdensityof(hmm, obs_seq) >= joint_logdensityof(hmm, obs_seq, any_state_seq)
marginal_logL = logdensityof(hmm, obs_seq)
println("Marginal logL: ", round(marginal_logL, digits=2))
println("Joint logL: ", round(joint_logL_viterbi, digits=2))
```

--------------------------------

### Find Most Likely State Sequence with `viterbi`

Source: https://context7.com/gdalle/hiddenmarkovmodels.jl/llms.txt

Implements the Viterbi algorithm to find the most likely sequence of hidden states given a sequence of observations. Returns the decoded state sequence and its joint log-likelihood. Can process multiple concatenated observation sequences.

```julia
using Distributions, HiddenMarkovModels

# Setup HMM and generate observations
init = [0.6, 0.4]
trans = [0.7 0.3; 0.2 0.8]
dists = [Normal(-1.0), Normal(1.0)]
hmm = HMM(init, trans, dists)

# Generate test data
true_states, obs_seq = rand(hmm, 50)

# Run Viterbi decoding
best_state_seq, best_joint_logL = viterbi(hmm, obs_seq)

# best_state_seq: Vector{Int} of most likely states
# best_joint_logL: Vector with one joint log-likelihood value
println("Decoded states (first 10): ", best_state_seq[1:10])
println("Joint log-likelihood: ", only(best_joint_logL))  # e.g., -72.45

# Check decoding accuracy
accuracy = mean(true_states .== best_state_seq)
println("Decoding accuracy: ", round(accuracy * 100, digits=1), "%")  # Often 80-95%
```

```julia
# Viterbi with multiple observation sequences
obs_seqs = [last(rand(hmm, rand(50:100))) for _ in 1:5]
obs_seq_concat = reduce(vcat, obs_seqs)
seq_ends = cumsum(length.(obs_seqs))

best_states_all, logL_all = viterbi(hmm, obs_seq_concat; seq_ends)
println("Number of sequences processed: ", length(logL_all))  # 5
println("Total states decoded: ", length(best_states_all))
```

--------------------------------

### Internal Functions

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/docs/src/api.md

Internal functions and types used within the HiddenMarkovModels.jl package.

```APIDOC
## Internals

```@docs
HiddenMarkovModels.LightDiagNormal
HiddenMarkovModels.LightCategorical
HiddenMarkovModels.log_initialization
HiddenMarkovModels.log_transition_matrix
HiddenMarkovModels.mul_rows_cols!
HiddenMarkovModels.argmaxplus_transmul!
```
```

--------------------------------

### Concatenate Sequence Data

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/docs/src/api.md

Use this snippet when dealing with multiple sequences to concatenate observation and control sequences and calculate their end indices.

```julia
obs_seq = reduce(vcat, obs_seqs)
control_seq = reduce(vcat, control_seqs)
seq_ends = cumsum(length.(obs_seqs))
```

--------------------------------

### HMM BibTeX Citation

Source: https://github.com/gdalle/hiddenmarkovmodels.jl/blob/main/README.md

Provides the BibTeX entry for citing the HiddenMarkovModels.jl package in research.

```bibtex
@article{
    Dalle2024,
    doi = {10.21105/joss.06436},
    url = {https://doi.org/10.21105/joss.06436},
    year = {2024},
    publisher = {The Open Journal},
    volume = {9},
    number = {96},
    pages = {6436},
    author = {Guillaume Dalle},
    title = {HiddenMarkovModels.jl: generic, fast and reliable state space modeling},
    journal = {Journal of Open Source Software}
}
```

=== COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.