feat: add HypergeometricDistribution, BetaBinomialDistribution, ZipfDistribution — additional discrete distributions

[OldCrow]

Summary

Add three discrete distributions that extend libstats' coverage of common discrete models. None has a SIMD batch path (all use lgamma/special-function scalar loops), but all are commonly needed.

HypergeometricDistribution

Parametrisation: population N, success states K, draws n. PMF: P(X=k) = C(K,k)·C(N-K,n-k)/C(N,n). Scipy: hypergeom.

• Mean = nK/N, Variance = nK(N-K)(N-n) / (N²(N-1))
• MLE: typically K and N are known from the population; n is the sample size. MLE for K from sample not commonly needed.
• Use cases: sampling without replacement (quality control, genetics, card game probability, clinical trial stratification), Fisher's exact test
• PMF implementation: lgamma-based log-space computation (same pattern as Binomial); CDF by summation
• Edge case: HypergeometricDistribution reduces to Binomial as N → ∞ with K/N → p

BetaBinomialDistribution

Parametrisation: trials n, concentration α > 0, β > 0. The Binomial with a Beta-distributed success probability p ~ Beta(α, β). Scipy: betabinom.

PMF(k) = C(n,k) · B(k+α, n-k+β) / B(α, β)
log PMF = log C(n,k) + log B(k+α, n-k+β) - log B(α, β)

All via lgamma. The normalisation constant B(α,β) is precomputed at construction.

• Mean = nα/(α+β), Variance = nαβ(α+β+n) / ((α+β)²(α+β+1))
• MLE: method of moments for α, β given n; or Newton–Raphson on log-likelihood
• Use cases: overdispersed count data (when Binomial variance is insufficient), click-through rate modelling with Beta prior, ecological count models

ZipfDistribution (Zeta / Power-Law)

Parametrisation: exponent a > 1. PMF: P(X=k) = k^{-a} / ζ(a) where ζ is the Riemann zeta function. Scipy: zipf (also zipfian for truncated/shifted variant).

• Mean = ζ(a-1)/ζ(a) for a > 2; Variance requires a > 3
• MLE: no closed form; Newton–Raphson on the log-likelihood which involves the digamma function
• The Riemann zeta function ζ(a) is a scalar constant computed at construction time; no per-element special function needed for PMF evaluation given precomputed ζ(a).
• Use cases: word frequency in natural language (Zipf's law), city population sizes, internet traffic, network degree distributions
• Note: support k = 1, 2, 3, ... (no k = 0); document clearly to avoid confusion with scipy's 1-indexed convention

Implementation notes

• All three use lgamma for PMF computation (same infrastructure as Binomial and NegativeBinomial)
• CDF by PMF summation (no closed-form for any of the three)
• Dispatch thresholds: NEVER for all ops in all tables (no SIMD benefit; scalar summation CDF path)
• Register all three in distribution_meta.h with is_discrete = true
• Zeta function: implement detail::riemann_zeta(double a) using the Euler–Maclaurin formula or existing math libraries; verify accuracy against reference values at construction time