feat: complete statistics trait coverage (Min/Max/Median/Mode)

[agene0001]

Closes #276. Independent of my other open PRs — branches directly off main.

Addresses both actionable bullets from your comment: implement the traits whose
value is unambiguous, and document the ambiguity where it isn't.

⚠️ One behavioural change, flagged deliberately

Dirichlet::pdf no longer panics for input off the simplex — it returns
0.0 (and ln_pdf returns NEG_INFINITY). This is the one thing here that
isn't purely additive, so please look at it first.

Why it's in scope: Multinomial::pmf already returns 0 for coordinates that
don't sum to n, and every univariate distribution returns 0 outside its
support. Dirichlet was the sole holdout, and I found it by writing
assert_eq!(d.pdf(&d.min()), 0.0) for the Min/Max tests below and watching
it panic.

What exactly changed:

• The set of inputs classified as out-of-support is unchanged — same
  (f64::MIN_POSITIVE..1.0) element predicate, same 1e-4 sum tolerance. Only
  the response changed, from panic! to 0.0.
• The dimension-mismatch panic is kept. That's a usage error with no
  meaningful density to return, and Multinomial panics on it too.
• Four #[should_panic] tests became behavioural tests
  (test_pdf_bad_input_range, test_pdf_bad_input_sum, and the two ln_pdf
  equivalents → test_pdf_out_of_support_is_zero,
  test_ln_pdf_out_of_support_is_neg_infinity). Those tests encoded the old
  contract, so their rewrite is the honest record of the break.

If you'd rather not take this, say so and I'll drop it to a separate PR — the
rest of the branch doesn't depend on it, beyond two assertions in the
Min/Max test.

Coverage

	Min	Max	Median	Mode
Multinomial	❌→✅	❌→✅	⊘	⊘
Dirichlet	❌→✅	❌→✅	⊘	❌→✅
MultivariateNormal	✅	✅	❌→✅	✅
MultivariateStudent	✅	✅	❌→✅	✅
Empirical	✅	✅	❌→✅	⊘
Categorical	✅	✅	✅	❌→✅
Beta Chi Erlang FisherSnedecor Gamma InverseGamma Hypergeometric NegativeBinomial	✅	✅	❌→✅	✅

⊘ = deliberately not implemented; reasons at the bottom.

Min and Max are now implemented for all 33 distributions. Median and
Mode are complete except for the four ⊘ cells.

Min/Max for the multivariate types

Multinomial: min = 0, max = n, in every coordinate
Dirichlet:   min = 0, max = 1, in every coordinate

Typed to each distribution's support, following the univariate precedent:
Multinomial's Discrete impl is over OVector<u64, D>, and every univariate
discrete distribution uses its integer support type (Min<u64> for Binomial,
Min<i64> for DiscreteUniform). Dirichlet gets OVector<f64, D>, agreeing
with Beta.

On your "Max and Min for multivariate distributions are not clear". A
partial order admits no unique greatest element, so these can only be
componentwise bounds — and the resulting vector need not be in the support:
Multinomial::max sums to k*n where an outcome must sum to n;
Dirichlet::max sums to k where a sample lies on the unit simplex. Both are
the corner of the smallest axis-aligned box containing the support. Each
individual bound is attained one coordinate at a time. Documented at each
impl and pinned by tests: pmf is zero at both bound vectors, nonzero at the
per-coordinate extremes, and n == 0 is covered as the sole case where the two
bounds coincide and are a real outcome.

Median for the multivariate types

MultivariateNormal and MultivariateStudent are symmetric about their
location. The competing notions of multivariate median — marginal, geometric
(spatial), halfspace — disagree in general, but at a centre of symmetry they all
coincide, along with the mean and the mode. So median == μ is the only value
any definition could give, which is what makes these safe.

For MultivariateStudent this holds for every ν, including ν <= 1 where
mean() is None because the first moment doesn't converge. A median needs no
moments; there's a test for that case.

Median for Empirical

The sample median: middle order statistic for odd counts, mean of the two
middles for even ones. That second half is a convention, not a derivation — it
matches NumPy's and R's defaults and is the only value equidistant from both —
so it's documented as such, including that the result need not be an observed
value. Multiplicity is respected, so [1, 1, 1, 2] has median 1, not 1.5.

Multiplicity is the part of the BTreeMap walk most likely to be wrong, so it's
checked against a naive sort-and-index reference across 156 generated data sets
plus hand-written edge cases. Panics on empty data, matching Empirical's
existing Min/Max.

Mode for Dirichlet and Categorical

Dirichlet mirrors Beta, of which it is the generalization:
(α_i - 1) / (α_0 - K) reduces to (α - 1) / (α + β - 2) at K == 2, and it
returns None unless every α_i > 1 exactly as Beta returns None unless
both shapes exceed 1. Mode<Option<_>> is the crate's norm — 26 of 27
univariate impls — so no new signature was needed. Unlike min/max the mode
is a support point: its coordinates sum to (α_0 - K)/(α_0 - K) = 1. A test
confirms the formula genuinely maximizes the density, by shifting mass between
coordinate pairs (which keeps the point on the simplex) and checking the density
never increases.

Categorical had Median but not Mode, though the mode is the simpler of the
two. Ties are real — a fair die has no unique mode — so this returns the
lowest maximizing index, the only choice that doesn't depend on iteration
order.

The eight new Median impls for univariate distributions

All are inverse_cdf(0.5). None of these has a closed-form median, and the
trait already admits both approximations and cdf inversion: Poisson uses
floor(λ + 1/3 - 0.02/λ), Binomial uses floor(n*p), and Categorical was
already inverse_cdf(0.5). Note also that Median's trait docs make no
closed-form claim, whereas Mode's explicitly do — see below.

The cost is not uniform, and each impl says so. Beta and FisherSnedecor
reduce to inv_beta_reg; Gamma, Erlang, Chi and InverseGamma to a
numerical search; Hypergeometric and NegativeBinomial to a discrete
bisection returning an exact integer. So median() is O(many cdf evaluations)
for these where it is O(1) elsewhere, with nothing in the signature to say so.
That asymmetry is the main reason to review this section — happy to drop any
subset.

Validated three ways in tests/median_consistency.rs, because self-consistency
alone would hide an error shared with the cdf:

1. cdf(median()) == 0.5 across 45 parameter sets, including extremes like
   Beta(1e3, 1e3), Beta(0.1, 0.1), Gamma(0.1, 1) and
   Hypergeometric(1000, 500, 100).
2. Against mpmath at 50 digits, using its own incomplete beta and gamma
   rather than statrs's: worst |cdf(median) - 0.5| = 2.5e-13, typically
   ~1e-15.
3. Against closed forms where the distributions coincide — Gamma(1, r) and
   Erlang(1, r) against Exp(r)'s ln(2)/r, Chi(1) against the normal
   0.75 quantile, Beta(1,1) against 0.5, InverseGamma(1,1) against
   1/ln(2). Plus the known bracket [shape - 1/3, shape] for Gamma, checked
   independently of the cdf.

Agreement is ~1e-13 relative, not to the last bit — a root-find isn't
correctly rounded the way a closed form is. The tests assert that tolerance
explicitly rather than pretending otherwise.

Deliberately not implemented

• Median for Multinomial and Dirichlet. The vector of marginal medians
  is not the multivariate median and, worse, is not even in the support:
  Dirichlet(1,1,1) has Beta(1,2) marginals with median 1 - sqrt(1/2), and
  three of those sum to 0.879, not 1. Publishing that as "the median" would
  be wrong, not merely approximate. The geometric median has no closed form
  for either.
• Mode for Multinomial. There is no closed form. The exact mode is the
  Jefferson/D'Hondt divisor apportionment of n among the p_i, found by a
  threshold search — and the Mode trait's own docs say it "specifies that an
  object has a closed form solution for its mode(s)", which every existing impl
  honours. Implementing it would be the first violation of that stated contract,
  so it belongs behind your decision rather than in this PR. Happy to add it,
  either as an inherent method or with the trait doc amended, whichever you
  prefer.
• Mode for Empirical. Every value of continuous sample data occurs once,
  so everything ties and the answer is arbitrary.

Also noticed, not touched

multinomial.rs:371 has a commented-out impl Skewness<Vec<f64>> for Multinomial referencing a Skewness trait that doesn't exist in
statistics::traits. Dead code, worth deleting whenever convenient.

Verification

cargo fmt --check, cargo clippy --all-features --all-targets clean,
RUSTDOCFLAGS=-D warnings cargo doc --all-features clean, cargo hack check --feature-powerset --no-dev-deps clean, cargo hack check --rust-version --locked clean against Cargo.lock.MSRV, and the full suite across
--all-features, default, --no-default-features, and --no-default-features --features rand.

🤖 Generated with Claude Code