Add HSMM forward algorithm by rsenne · Pull Request #168 · JuliaStats/HiddenMarkovModels.jl

rsenne · 2026-08-27T15:40:57Z

Adds the segment-based forward algorithm for the HSMM mirrored as closely as possible to the HMM counterpart

forward follows the HMM interface and returns (α, logL). As with ForwardStorage.α, α[:, t] contains the filtered state marginals on the probability scale.
To support the right-censoring convention introduced in #167, the recursion tracks two quantities:

log_ends[i, t]: probability mass for segments that end at t, used by the transition step.
log_ongoing[i, t]: probability mass for segments that cover t but may continue beyond it.

cc: @simonsteiger

codecov · 2026-08-27T15:44:57Z

Codecov Report

❌ Patch coverage is 99.41176% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 97.06%. Comparing base (055ef57) to head (a0949e9).

Files with missing lines	Patch %	Lines
src/inference/forward_hsmm.jl	99.40%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #168      +/-   ##
==========================================
+ Coverage   96.50%   97.06%   +0.55%     
==========================================
  Files          24       25       +1     
  Lines         716      886     +170     
==========================================
+ Hits          691      860     +169     
- Misses         25       26       +1

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

rsenne · 2026-08-28T01:39:05Z

+    if seq_ends isa NTuple{1}
+        for k in eachindex(seq_ends)
+            _forward!(storage, hsmm, obs_seq, control_seq, seq_ends, k; error_if_not_finite)
+        end
+    else
+        @threads for k in eachindex(seq_ends)
+            _forward!(storage, hsmm, obs_seq, control_seq, seq_ends, k; error_if_not_finite)
+        end
+    end


This could be shared via something like this:

Suggested change

if seq_ends isa NTuple{1}

for k in eachindex(seq_ends)

_forward!(storage, hsmm, obs_seq, control_seq, seq_ends, k; error_if_not_finite)

end

else

@threads for k in eachindex(seq_ends)

_forward!(storage, hsmm, obs_seq, control_seq, seq_ends, k; error_if_not_finite)

end

end

function foreach_sequence(kernel::F, seq_ends::AbstractVectorOrNTuple{Int}) where {F}

if seq_ends isa NTuple{1}

for k in eachindex(seq_ends)

kernel(k)

end

else

@threads for k in eachindex(seq_ends)

kernel(k)

end

end

return nothing

end

unclear if its worth the hassle or not

The reason I put this workaround in place is because @threads was type-unstable, and I wanted a type-stable path at least in the single-sequence case. Not sure this is still necessary, and perhaps OhMyThreads.jl also has a better way to unify these cases

Good point. I have switched @threads to use OhMyThreads.jl on StateSpaceDynamics.jl and it works very well. I'll redo this with that in mind then. Mayeb I'll do a separte PR for the HMM case and you can review before that

rsenne · 2026-08-28T01:47:00Z

+$(TYPEDFIELDS)
+"""
+struct HSMMForwardStorage{R}
+    "posterior last state marginals `α[i] = ℙ(X[T]=i | Y[1:T])`"


i copied this form the equivalent struct for the HMM case. Would this be more accurate to say: "filtered state marginals α[i,t] = ℙ(X[t]=i | Y[1:t])". This somewhat to me implies only the last column

Segment-based forward recursion for `AbstractHSMM`, with the same public interface as the HMM version: `forward` returns `(α, logL)` where `α[:, t]` are the filtered state marginals `ℙ(X[t]=i | Y[1:t])` on the probability scale, exactly as for an `AbstractHMM`. Two quantities are propagated per timestep. `log_ends[i, t]` is the mass of segments that end at `t` in state `i` and drives the recursion; `log_ongoing[i, t]` is the mass of segments that merely cover `t`, i.e. are right-censored there. Every prefix of a sequence ends mid-segment under the censoring convention of `rand`, so the likelihood and the state marginals are read from `log_ongoing` while the transition sweep consumes `log_ends`. Both are accumulated in the same duration sweep, so censoring costs no extra duration evaluations. `max_duration` is a keyword defaulting to the longest sequence length, which makes the recursion exact. Lowering it discards every segment longer than the cutoff, so the result is a strict lower bound that tightens as the cutoff grows, rather than an error of unpredictable sign. The observation logdensity prefix sums live in an `N × T` matrix matching the layout of `ForwardStorage.B`, and each sequence writes only inside its own columns, so threaded sequences never share an entry. The transition matrix is hoisted out of the segment loop when the controls are uniform, since for a model whose log-transitions are not cached, looking it up per timestep allocates a matrix per timestep. Tested by brute-force marginalization of `joint_logdensityof` over all state sequences, by normalization to one over all short observation sequences, and against the equivalent HMM under geometric sojourns, where both the loglikelihoods and the filtered marginals agree elementwise.

gdalle · 2026-09-08T07:26:48Z

+    cum_log_obs = Matrix{R}(undef, N, T)
+    obs_zeros = Matrix{Int}(undef, N, T)
+
+    # Per-sequence scratch space keeps parallel calls independent.


Back when I coded this package I handled this manually, but I think there are better ways now with OhMyThreads.jl

https://juliafolds2.github.io/OhMyThreads.jl/stable/literate/tls/tls/

gdalle · 2026-09-08T07:28:19Z

+    return reachable
+end
+
+function extend_segments!(


This has lots of positional arguments, the risk of a mixup is non-negligible

Switch to kwargs?

gdalle · 2026-09-08T07:29:23Z

+    if seq_ends isa NTuple{1}
+        for k in eachindex(seq_ends)
+            _forward!(storage, hsmm, obs_seq, control_seq, seq_ends, k; error_if_not_finite)
+        end
+    else
+        @threads for k in eachindex(seq_ends)
+            _forward!(storage, hsmm, obs_seq, control_seq, seq_ends, k; error_if_not_finite)
+        end
+    end


The reason I put this workaround in place is because @threads was type-unstable, and I wanted a type-stable path at least in the single-sequence case. Not sure this is still necessary, and perhaps OhMyThreads.jl also has a better way to unify these cases

gdalle · 2026-09-08T07:29:43Z

+"""
+$(SIGNATURES)
+"""
+function forward!(


Do you have a paper reference for this algorithm, so that I can try to follow the math?

https://www.sciencedirect.com/science/article/pii/S0004370209001416
https://dspace.mit.edu/entities/publication/fb421d82-a37b-4656-8c00-88e846365927

We specifically are using the Explicit Duration HMM variant

(we can reuse less code, but space complexity is far better)

rsenne · 2026-09-13T16:43:12Z

I'm thinking of punting the OhMyThreads.jl addition to its own PR as it will impact both HSMMs and HMMs. That will keeo this one focused. Any objection?

There is also an issue with testing related to 1.13 where it seems the norm keyword is not being honored. I did not have this issue in 1.12.5. I did a quick infnorm thing as a workaround.

gdalle · 2026-09-14T07:34:01Z

Yeah, let's separate the threading updates from this new algo

Julia 1.13 ignores the `norm` keyword of `isapprox` for arrays, so the coherence tests silently compared with the Frobenius norm instead of the intended infinity norm, making the tolerance stricter than intended. Replace the six call sites with an explicit `infapprox` helper.

rsenne · 2026-09-15T16:59:47Z

just want to confirm this--I'm again not seeing a review/approval, but GH is letting me merge this, but I do not want to merge this prior to your go ahead. Can you confirm or deny an ongoing review on this?

gdalle · 2026-09-15T18:03:58Z

I haven't reviewed that one again and I would like to!

rsenne changed the title ~~Add HSMM forward algorithm and logdensityof~~ Add HSMM forward algorithm Aug 27, 2026

rsenne commented Aug 28, 2026

View reviewed changes

rsenne mentioned this pull request Aug 28, 2026

Add HSMM joint_logdensityof and duration_logsurvival #167

Merged

rsenne force-pushed the part3b-hsmm-forward branch from 56d2596 to 6d6064e Compare August 28, 2026 20:34

rsenne force-pushed the part3b-hsmm-forward branch 2 times, most recently from e78bcc9 to 76ea431 Compare September 7, 2026 23:05

Base automatically changed from part3a-hsmm-joint-logdensity to main September 8, 2026 07:20

rsenne added 2 commits September 8, 2026 09:20

Changes for self review

912cc26

gdalle force-pushed the part3b-hsmm-forward branch from 76ea431 to 912cc26 Compare September 8, 2026 07:20

gdalle reviewed Sep 8, 2026

View reviewed changes

rsenne added 2 commits September 10, 2026 10:17

small changes

8151273

add formulas

59badaf

rsenne requested a review from gdalle September 13, 2026 16:43

rsenne mentioned this pull request Sep 13, 2026

Fix isapprox norm kwarg no longer honored on Julia 1.13 #170

Merged

delete

e30fb0f

Uh oh!

Conversation

rsenne commented Aug 27, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

codecov Bot commented Aug 27, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rsenne commented Sep 13, 2026

Uh oh!

gdalle commented Sep 14, 2026

Uh oh!

rsenne commented Sep 15, 2026

Uh oh!

gdalle commented Sep 15, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

rsenne commented Aug 27, 2026 •

edited

Loading

codecov Bot commented Aug 27, 2026 •

edited

Loading