Add hsm_mode: Half-Sample Mode for continuous data

[anurag-mds]

Summary

This PR adds hsm_mode(), an implementation of the half-sample mode (HSM) which is a robust estimator of the mode for continuous distributions.

It is introduced as a separate function ( not an overload of mode() ) to preserve existing behaviour while providing a statistically meaningful alternative for continuous data

This addresses and closes issue #957.

---

Motivation

StatsBase.mode() is frequency-based and works well for discrete data. For continuous distributions, however, samples are usually unique, which makes frequency counts unstable and highly variable in practice

Issue #957 documents this behaviour, particularly for heavy-tailed distributions, where mode() can show extreme variance.

This PR provides an estimator designed specifically for continuous data

---

Approach

hsm_mode() implements the standard half-sample method described in the literature:

• Non-finite values (NaN, Inf) are filtered

• The data are sorted

• The algorithm repeatedly selects the contiguous half-sample with the smallest width

• Once ≤ 2 points remain, the midpoint of the final interval is returned

• The midpoint may not be a sample value, but provides a stable estimate of the location of highest density.

• Time complexity is dominated by sorting (O(n log n)); space complexity is O(n).

• After sorting, the contraction loop operates on SubArray views to avoid allocations.

---

API Design

The estimator is exposed as a new function:

hsm_mode(x::AbstractVector{T}) where T<:Real

It is NOT added as an overload of mode() in order to:

• avoid changing existing semantics

• clearly distinguish frequency-based and density-based estimation

• let users choose the appropriate method explicitly

The return type is an AbstractFloat, promoted from the input element type (e.g. integers → Float64, Float32 → Float32).

---

Testing and Documentation

Tests cover basic correctness, edge cases, robustness to outliers, handling of non-finite values, and type behavior. All tests pass.

The docstring explains intended use cases, compares with mode(), documents complexity, provides examples, and cites the relevant literature.

---

References

Robertson & Cryer (1974), JASA

Bickel & Fruehwirth (2006), CSDA

---

Notes

This PR is intentionally small and focused. Extensions such as weighted HSM or support for missing values are left for future work.

I have attached images showing that the tests are passing, and I would like to know whether I should address the existing warnings in the codebase.

Images:

Feedback on naming or API placement is welcome.

[devmotion]

Was this PR AI-generated?

[anurag-mds]

The implementation, test and commits are mine. I reviewed existing StatsBase prs to match project conventions, and I used AI assistance only for wording clarity in the pr description and for general performance review and bugs made by me. it was not used for algorithm or code . since in ci pipeline documentation issues which I totally forgot was missing I am actively fixing that like using ceil(len/2) as per the definition. I am eager to explain any design or choices in detail

[ForceBru]

introduced as a separate function ( not an overload of mode() ) to preserve existing behaviour

I feel like the existing behavior (for floating-point numbers) is a bug. If your data are not integers, we must assume that they come from a continuous distribution and use the HSM algorithm or another estimator for continuous data.

Or perhaps introduce a keyword argument:

mode(x::AbstractVector{<:AbstractFloat}; discrete::Bool=false) =
  discrete ? mode_discrete(x) : mode_hsm(x)

mode(x::AbstractVector) = mode_discrete(x)

Here, mode_discrete is what's currently called mode. The keyword argument lets one use the high-variance estimator for AbstractFloat.

I'm not a contributor to Distributions.jl, though (although I'd like to become a contributor; my complete PR with hyperbolic distributions seems to have joined its peers, unfortunately), that's just my opinion.

[anurag-mds]

I think using frequency-based mode for floats can be misleading.

But making HSM the default for AbstractFloat might break cases where floats are actually categories or rounded values.

Maybe a keyword like discrete=false or a dispatch mode(x, HSM()) works better. I can change the API like that if the maintainers think it’s right.

[anurag-mds]

@ForceBru
You have made an excellent point about API design consistency. You're right
that distinguishing frequency-based vs density-based estimation matters.

Here's the evidence for why hsm_mode() should remain separate:

As you can see
Insight: Frequency-based mode() counts collisions on continuous
data where samples are usually unique. The "most frequent" element becomes
arbitrary and order-dependent. HSM finds the actual density peak instead.

Design:

• mode(): Frequency-based, correct for discrete data [keep as-is]
• hsm_mode(): Density-based, essential for continuous data [separate]

Users can explicitly choose the right tool. No silent behavior changes.

So are you okay with this separate function approach or shall we explore the keyword argument variant further ?

[anurag-mds]

All tests pass, edge cases handled, documentation added. @devmotion are there any remaining technical concerns with hsm_mode() or its API before approval?

[nalimilan]

Thanks for the PR!

I tend to agree that a keyword argument would make sense here. The mode is a statistical concept which can be estimated in various ways. The one currently used by mode is the simplest, and of course it's bad for continuous samples unless you apply some binning. But mode(Normal()) returns zero already so the concept makes sense in general.

The API could look like mode(x, method=:halfsample).

[anurag-mds]

Thanks for your detailed review
This is very helpful. I'll address the points you raised API via keyword, iterator support, middle , non-finite handling, test adjustments and fix the CI failures.

[anurag-mds]

Apologies for the back-and-forth and the CI noise I'm currently away from my main development setup, but I’ve noted all changes precisely and will push a consolidated update shortly once I’m back. I’ll comment again once the updates are in.
Thanks again for the detailed guidance.

[anurag-mds]

@nalimilan Does everything seems good?

[nalimilan]

Thanks!

[anurag-mds]

Thanks @nalimilan for the thorough review. I agree with the remaining points (doc wording, method naming, references, and test adjustments). I’m addressing them now and will push a final cleanup commit so that
everything is consolidated.

[anurag-mds]

@nalimilan I have made necessary changes if anything to fix, modify or add do let me know!
Thanks a lot for the guidance

[anurag-mds]

@nalimilan What do you think?

[nalimilan]

Test with 2:4 to actually test that argument.

[anurag-mds]

Everything Noted.
I will make sure that all issues are properly addressed...

[anurag-mds]

Rebasing took more time Is everything good @nalimilan ?

[anurag-mds]

Is there anything left @nalimilan ?

[devmotion]

I can't see this restriction in the other mode implementation above. So maybe that is only a limitation of half-sample mode?

[anurag-mds]

Yes, it is

As you can see for Regular Frequency Mode:

Consider [Nan, Inf, -Inf]
Here,

1. It can still count frequencies
2. Nan appears once, Inf appears once
3. Returns first one

Now, for HSM Mode:

Consider the same array,
Here,

1. Filters out all non-finite values
2. Nothing left to calculate
3. MUST throw error
4. Cannot find a cluster with no data!

[nalimilan]

Then the error message should be more explicit:

Suggested change

	len == 0 && throw(ArgumentError("mode is not defined for collections with no finite values"))
	len == 0 && throw(ArgumentError("mode with `method=:halfsample` is not defined " *
	"for collections with no finite values"))

[anurag-mds]

I've addressed the final performance suggestions (sort! and LazyString). Ready for final review.
@nalimilan @devmotion

[anurag-mds]

Can you please review the updated changes now?

[nalimilan]

Sorry for the delay. Looks almost ready to me.

@devmotion Do you have better ideas of a name for method=:frequency? The definition of the mode implies frequency whatever the estimation method, but I can't find a better name...

[anurag-mds]

@nalimilan Done applied the suggested changes
Extremely sorry for such a noisy PR
But I got to know many standards and style julia follows

[anurag-mds]

Fixed the remaining points at most typo, exact test value pinned (1.0750000000000002), error message updated to 'for collections containing only NaN values', filteredv view used throughout for type stability, frequency tests updated to use longer vector, docstring signature added.

One change I made proactively based on your comment about:

• Inf: changed the filter from isfinite to !isnan, so Inf and -Inf are now allowed through.

• mode([NaN, Inf, -Inf]) returns NaN rather than throwing.

Happy to revert this if you'd prefer the stricter behaviour.

[nalimilan]

• mode([NaN, Inf, -Inf]) returns NaN rather than throwing.

Please throw an ArgumentError instead in the presence of NaN, this seems more explicit to me and it's easier to change later if we want.

[nalimilan]

Need to use filteredv here and below, right? I'm a bit worried that this wasn't caught by tests. Isn't this covered?

[anurag-mds]

The reason the tests didn't catch it is that the previous test cases were too simple they didn't 'squeeze' the data enough to show the mistake. I've now updated the code to use filteredv everywhere and added a more complex test case ([1.0, 2.0, ... 13.0]) that would have failed with the old bug.

I also switched to a strict ArgumentError for any non-finite values as you suggested. Ready for another look!

[anurag-mds]

I also wanted to say I'm genuinely sorry for the extra work this PR has caused. When this PR was created it was my very first PR, and I'm honestly a bit embarrassed that I ended up increasing your workload instead of making it easier. I’m learning a lot about the project's standards through this process, and I really appreciate your patience and guidance in helping me get this implementation right.