diff --git a/CLAUDE.md b/CLAUDE.md index f2fdc5c..dc50833 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -45,7 +45,7 @@ This is a .NET library (`ktsu.Essentials`) providing high-performance interfaces - `Essentials/IContinuousDistribution.cs` - `IDistribution` plus `Pdf`/`LogPdf` - `Essentials/IDiscreteDistribution.cs` - `IDistribution` plus `Pmf`/`LogPmf`, and a default `Quantile` that bisects the CDF. That default is an explicit reimplementation of `IDistribution.Quantile`, which CA1033 misreads as an explicit implementation on a class; the suppression on it says why - `Shared/DeflateBound.cs` - zlib's `deflateBound` plus container overhead, linked into the Deflate, Gzip and ZLib providers as their `GetMaxCompressedLength` -- `Shared/SpecialFunctions.cs` - Incomplete gamma, incomplete beta, error function, normal quantile and cancellation-free `ExpM1`/`Log1P`, linked into the Normal, LogNormal, Exponential, Binomial and Poisson providers rather than placed in the interfaces-only package +- `Shared/SpecialFunctions.cs` - Incomplete gamma, incomplete beta, error function, normal quantile and cancellation-free `ExpM1`/`Log1P`, linked into the Normal, LogNormal, Exponential, Binomial, Poisson and Geometric providers rather than placed in the interfaces-only package - `Shared/DistributionArguments.cs` - Parameter validation shared by every distribution provider, linked in the same way - `Shared/PersistenceDefaults.cs` - The `IPersistenceProvider` default `RetrieveOrCreateAsync` body, linked into the FileSystem and Temp providers, which declare that member themselves and so hide the default - `Essentials/ISerializationProvider.cs` - Object serialization/deserialization interface diff --git a/Essentials.DistributionProviders.Geometric/Essentials.DistributionProviders.Geometric.csproj b/Essentials.DistributionProviders.Geometric/Essentials.DistributionProviders.Geometric.csproj index 1be4173..cdeaed3 100644 --- a/Essentials.DistributionProviders.Geometric/Essentials.DistributionProviders.Geometric.csproj +++ b/Essentials.DistributionProviders.Geometric/Essentials.DistributionProviders.Geometric.csproj @@ -13,6 +13,7 @@ + diff --git a/Essentials.DistributionProviders.Geometric/GeometricDistributionProvider.cs b/Essentials.DistributionProviders.Geometric/GeometricDistributionProvider.cs index 80197ec..d45eba5 100644 --- a/Essentials.DistributionProviders.Geometric/GeometricDistributionProvider.cs +++ b/Essentials.DistributionProviders.Geometric/GeometricDistributionProvider.cs @@ -1,4 +1,4 @@ -// Copyright (c) 2023-2026 ktsu-dev contributors +// Copyright (c) 2023-2026 ktsu-dev contributors namespace ktsu.Essentials.DistributionProviders.Geometric; @@ -23,6 +23,11 @@ namespace ktsu.Essentials.DistributionProviders.Geometric; /// public sealed class GeometricDistributionProvider : IDiscreteDistribution { + /// + /// How far either way the closed-form quantile may be walked before giving up on it and bisecting. + /// + private const int MaxQuantileCorrection = 2; + private readonly double failureProbability; private readonly double logFailureProbability; @@ -52,7 +57,11 @@ public GeometricDistributionProvider(double probability) Probability = probability; failureProbability = 1.0 - probability; - logFailureProbability = Math.Log(failureProbability); + + // Taken as log1p(-p) rather than the log of 1 - p, which rounds to exactly one below p of about + // 1e-16 and leaves a logarithm of zero: every CDF value would then be zero and every quantile a + // walk to int.MaxValue. Above that it still costs digits in proportion to how small p is. + logFailureProbability = SpecialFunctions.Log1P(-probability); } /// @@ -73,10 +82,29 @@ public GeometricDistributionProvider(double probability) public double Variance => failureProbability / (Probability * Probability); /// - public double Pmf(int value) => value < 0 ? 0.0 : Math.Pow(failureProbability, value) * Probability; + public double Pmf(int value) + { + if (value < 0) + { + return 0.0; + } + + // Zero failures is the success probability itself. Taking it through the logarithm would multiply a + // zero count by the infinite logarithm of a certain success, which is a NaN rather than one. + if (value == 0) + { + return Probability; + } + + return Math.Exp(value * logFailureProbability) * Probability; + } /// - public double Cdf(int value) => value < 0 ? 0.0 : 1.0 - Math.Pow(failureProbability, value + 1.0); + /// + /// One minus (1 - p)^(k + 1), written as -expm1((k + 1) log1p(-p)) so that neither the + /// failure probability nor the subtraction from one throws away the digits a small CDF lives in. + /// + public double Cdf(int value) => value < 0 ? 0.0 : -SpecialFunctions.ExpM1((value + 1.0) * logFailureProbability); /// /// Evaluates the survival function, the probability of more than failures. @@ -87,7 +115,7 @@ public GeometricDistributionProvider(double probability) /// /// The failure count to evaluate at. /// A probability in the range [0, 1]. - public double SurvivalFunction(int value) => value < 0 ? 1.0 : Math.Pow(failureProbability, value + 1.0); + public double SurvivalFunction(int value) => value < 0 ? 1.0 : Math.Exp((value + 1.0) * logFailureProbability); /// public int Quantile(double probability) @@ -106,23 +134,51 @@ public int Quantile(double probability) // The smallest count k whose CDF reaches the target satisfies (k + 1) >= log(1 - p) / log(1 - q), // with the inequality flipping because the logarithm of the failure probability is negative. - double raw = Math.Ceiling(Math.Log(1.0 - probability) / logFailureProbability) - 1.0; + double raw = Math.Ceiling(SpecialFunctions.Log1P(-probability) / logFailureProbability) - 1.0; int count = raw <= 0.0 ? 0 : raw >= int.MaxValue ? int.MaxValue : (int)raw; // That form is exact in real arithmetic, but where the target is itself a CDF value the ratio of // the two logarithms lands a rounding step either side of a whole number, and the ceiling turns - // a step into a whole outcome. One comparison against the CDF in each direction puts it back, and - // neither loop can run more than once because the estimate is never off by more than one. - while (count > 0 && Cdf(count - 1) >= probability) + // a step into a whole outcome. A comparison against the CDF in each direction puts it back. The + // estimate should never be off by more than one, but the walks are bounded all the same: an + // unbounded one is what turned a bad estimate into two billion steps, and bisecting the whole + // support costs only 31 CDF evaluations. + for (int step = 0; step < MaxQuantileCorrection && count > 0 && Cdf(count - 1) >= probability; step++) { count--; } - while (count < int.MaxValue && Cdf(count) < probability) + for (int step = 0; step < MaxQuantileCorrection && count < int.MaxValue && Cdf(count) < probability; step++) { count++; } - return count; + bool settled = (count == int.MaxValue || Cdf(count) >= probability) && (count == 0 || Cdf(count - 1) < probability); + return settled ? count : Bisect(probability); + } + + /// + /// Finds the smallest count whose CDF reaches by bisecting the support. + /// + /// The target probability. + /// The quantile. + internal int Bisect(double probability) + { + int low = 0; + int high = int.MaxValue; + while (low < high) + { + int mid = low + ((high - low) / 2); + if (Cdf(mid) >= probability) + { + high = mid; + } + else + { + low = mid + 1; + } + } + + return low; } } diff --git a/Essentials.Tests/DistributionProviderTests.cs b/Essentials.Tests/DistributionProviderTests.cs index 8e00dfc..91f99e4 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; @@ -672,6 +672,44 @@ public void Geometric_Matches_Its_Closed_Form() Assert.AreEqual(0, certain.Quantile(1.0)); } + [TestMethod] + public void Geometric_Keeps_Its_Precision_For_A_Rare_Success() + { + // Below p of about 1e-16, 1 - p rounds to exactly one, so a CDF built on it is zero everywhere and + // the quantile walks one outcome at a time to int.MaxValue. The reference values were computed to + // 60 digits outside this codebase. + GeometricDistributionProvider rare = new(1e-17); + Assert.AreEqual(1.0000000095e-9, rare.Cdf(100_000_000), 1.0000000095e-9 * 1e-12); + Assert.AreEqual(1.0 - 1.0000000095e-9, rare.SurvivalFunction(100_000_000), 1e-15); + Assert.AreEqual(100_000_000, rare.Quantile(1e-9)); + + // Above that threshold 1 - p survives but carries a relative error of about ε/p into the CDF. + GeometricDistributionProvider uncommon = new(1e-12); + Assert.AreEqual(9.995001676245083e-4, uncommon.Cdf(1_000_000_000), 9.995001676245083e-4 * 1e-12); + + // At the bottom of the double range the CDF is (k + 1) p to every digit that can be stored. + GeometricDistributionProvider vanishing = new(1e-300); + Assert.AreEqual(1001e-300, vanishing.Cdf(1000), 1001e-300 * 1e-15); + Assert.AreEqual(1e-300, vanishing.Pmf(1000), 1e-300 * 1e-15); + Assert.AreEqual(0.0, vanishing.Pmf(-1), Tight); + } + + [TestMethod] + public void Geometric_Bisection_Fallback_Agrees_With_The_Closed_Form_Quantile() + { + // The quantile falls back to bisecting the support only when its closed-form estimate is off by more + // than the bounded correction allows, which no ordinary parameter provokes. Checked directly, so the + // fallback is known to land on the same outcome whenever it is needed. + foreach (double p in new[] { 0.25, 1e-6, 1e-17 }) + { + GeometricDistributionProvider distribution = new(p); + foreach (double q in new[] { 0.0, 1e-9, 0.25, 0.5, 0.999 }) + { + Assert.AreEqual(distribution.Quantile(q), distribution.Bisect(q), $"Geometric({p}) at {q}"); + } + } + } + [TestMethod] public void Categorical_Normalises_Its_Weights() {