BinomialDistribution CDF slower than scipy; PDF near-parity — PMF summation and scalar lgamma are limiting factors

[OldCrow]

Summary

BinomialDistribution CDF is consistently slower than scipy across all batch sizes and both AVX-512 (Zen4) and AVX2+FMA (Kaby Lake) test machines. PDF/log_PDF are near parity (1.0–1.4× on Zen4, 1.9–3.5× on Kaby Lake) rather than the 6–35× advantage seen for continuous distributions. This makes Binomial the weakest-performing discrete distribution in the benchmark.

Evidence

From benchmarks/scipy_comparison.py (v0.3.2, scipy 1.18.0):

Binomial PDF speedup:

Machine	N=1k	N=10k	N=100k	N=1M
Zen4 / AVX-512	1.4×	1.1×	1.1×	1.2×
Kaby Lake / AVX2+FMA	3.5×	2.2×	1.9×	2.1×

Binomial CDF speedup:

Machine	N=1k	N=10k	N=100k	N=1M
Zen4 / AVX-512	0.4×	0.6×	0.8×	0.9×
Kaby Lake / AVX2+FMA	0.9×	0.5×	0.4×	0.4×

Absolute CDF throughput: 5–16 M/s on both machines.

For comparison, Poisson CDF (also discrete) achieves 1.6–10.7×. The Binomial underperformance is specific, not a general discrete-distribution issue.

Suspected root causes

1. CDF algorithm: Binomial CDF is implemented as PMF summation (Σ C(n,k)·pᵏ·(1−p)ⁿ⁻ᵏ in log space). scipy.stats.binom.cdf uses scipy.special.bdtr (Cephes bdtr), which computes the regularised incomplete-beta function I_{1-p}(n-k, k+1). For large n, the incomplete-beta approach converges much faster than summing individual PMF terms.

2. PDF near-parity: Binomial PMF uses lgamma log-space evaluation per element in a scalar loop (no SIMD log-gamma primitive). Continuous distributions with SIMD vector_log paths achieve 10–50× speedup; Binomial is limited to ~1–3× because each PMF call reduces to lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1) + k·log(p) + (n-k)·log(1-p) — a mix of lgamma and log that doesn't currently vectorise.

Investigation steps

1. CDF: Evaluate switching from PMF summation to the regularised incomplete-beta approach (already available in the codebase via Beta/StudentT infrastructure). Benchmark the crossover point where beta is faster than summation.
2. PMF (PDF/log_PDF): Assess whether a vector_lgamma primitive is feasible. The AGENTS.md notes this as indefinitely deferred (too complex, low immediate impact) — validate against the benchmark data whether that assessment still holds.
3. Dispatch thresholds: Verify whether the Binomial scalar loop is even reaching the SIMD dispatch threshold at N=1k, or whether the poor N=1k ratio (1.4× Zen4, 1.0× on N=10k) reflects a threshold that puts small batches in an inefficient path.

Related

• deferred: vector_lgamma — AGENTS.md §Deferred Items
• Poisson CDF is also scalar but shows better ratios; understand whether it uses a different algorithm
• Issue bessel.h Tier 2 fallback (A&S polynomial) limits VonMises accuracy to ~10⁻⁷ on macOS/AppleClang #47 (Bessel Tier 2) shows the pattern of a known-limited scalar path; Binomial CDF is analogous