From 1314de85a01639a5a44dfb6c50f78f3b66babcfe Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 7 Oct 2026 06:27:58 +0000 Subject: [PATCH] Keep the binomial upper tail precise instead of rounding it to zero [patch] BinomialDistributionProvider fell back to IDistribution's 1 - Cdf for its survival function, so any upper tail below about 1e-16 came back as zero. Evaluate it directly as I(p; k + 1, n - k), as Poisson already does with its own complementary incomplete function. Fixes #50 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01Vpup29Ycn8SGi7LvBjAq1c --- .../BinomialDistributionProvider.cs | 17 ++++++++++++++++- Essentials.Tests/DistributionProviderTests.cs | 19 ++++++++++++++++++- 2 files changed, 34 insertions(+), 2 deletions(-) diff --git a/Essentials.DistributionProviders.Binomial/BinomialDistributionProvider.cs b/Essentials.DistributionProviders.Binomial/BinomialDistributionProvider.cs index d856106..3377d9f 100644 --- a/Essentials.DistributionProviders.Binomial/BinomialDistributionProvider.cs +++ b/Essentials.DistributionProviders.Binomial/BinomialDistributionProvider.cs @@ -1,4 +1,4 @@ -// Copyright (c) 2023-2026 ktsu-dev contributors +// Copyright (c) 2023-2026 ktsu-dev contributors namespace ktsu.Essentials.DistributionProviders.Binomial; @@ -118,6 +118,21 @@ public double Cdf(int value) : SpecialFunctions.RegularizedIncompleteBeta(Trials - value, value + 1.0, failureProbability); } + /// + public double SurvivalFunction(int value) + { + if (value < 0) + { + return 1.0; + } + + // The upper tail is P(X > k) = I(p; k + 1, n - k). Taking it as one minus the CDF would round + // every tail below about 1e-16 to zero, and those small tails are usually the reason to ask. + return value >= Trials + ? 0.0 + : SpecialFunctions.RegularizedIncompleteBeta(value + 1.0, Trials - value, Probability); + } + /// public int Sample(IRandomProvider random) { diff --git a/Essentials.Tests/DistributionProviderTests.cs b/Essentials.Tests/DistributionProviderTests.cs index 8e00dfc..1545c4a 100644 --- a/Essentials.Tests/DistributionProviderTests.cs +++ b/Essentials.Tests/DistributionProviderTests.cs @@ -1,4 +1,4 @@ -// Copyright (c) 2023-2026 ktsu-dev contributors +// Copyright (c) 2023-2026 ktsu-dev contributors namespace ktsu.Essentials.Tests; @@ -622,6 +622,23 @@ public void Binomial_Stays_Finite_Where_The_Coefficient_Would_Overflow() Assert.AreEqual(1.0, (2.0 * distribution.Cdf(999)) + distribution.Pmf(1000), 1e-11); } + [TestMethod] + public void Binomial_Keeps_Its_Precision_Far_Out_In_The_Upper_Tail() + { + // Both values were computed exactly, as rationals, outside this codebase. One minus the CDF rounds + // the first to zero and keeps only eight digits of the second. + IDiscreteDistribution distribution = new BinomialDistributionProvider(100, 0.01); + Assert.AreEqual(9.576655593219767e-22, distribution.SurvivalFunction(20), 9.576655593219767e-22 * 1e-12); + Assert.AreEqual(6.255518382834089e-09, distribution.SurvivalFunction(10), 6.255518382834089e-09 * 1e-12); + + Assert.AreEqual(1.0, distribution.SurvivalFunction(-1), Tight); + Assert.AreEqual(0.0, distribution.SurvivalFunction(100), Tight); + IDiscreteDistribution never = new BinomialDistributionProvider(5, 0.0); + Assert.AreEqual(0.0, never.SurvivalFunction(0), Tight); + IDiscreteDistribution always = new BinomialDistributionProvider(5, 1.0); + Assert.AreEqual(1.0, always.SurvivalFunction(4), Tight); + } + [TestMethod] public void Poisson_Matches_Values_Summed_Independently() {