From 8a89874d4543ceb041e641390ced1b768aad8e6f Mon Sep 17 00:00:00 2001 From: Florian Schwerteck Date: Thu, 29 Sep 2022 22:21:23 +0200 Subject: [PATCH 1/8] Fix description of the neg binomial distribution in the F# interface --- src/FSharp/Distributions.fs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/FSharp/Distributions.fs b/src/FSharp/Distributions.fs index d6f6a8a9b..266ceb10c 100644 --- a/src/FSharp/Distributions.fs +++ b/src/FSharp/Distributions.fs @@ -118,7 +118,7 @@ module Sample = let logNormal mu sigma (rng:System.Random) = LogNormal.Sample(rng, mu, sigma) let logNormalSeq mu sigma (rng:System.Random) = LogNormal.Samples(rng, mu, sigma) - /// Negative-Binomial with number of failures (r) until the experiment stopped and probability (p) of a trial resulting in success. + /// Negative-Binomial with number of successes (r) until the experiment stopped and probability (p) of a trial resulting in success. let negativeBinomial r p (rng:System.Random) = NegativeBinomial.Sample(rng, r, p) let negativeBinomialSeq r p (rng:System.Random) = NegativeBinomial.Samples(rng, r, p) From d1fb730b4da0ea1ad19b602d4d7e2172d0e4e62f Mon Sep 17 00:00:00 2001 From: Florian Schwerteck Date: Thu, 29 Sep 2022 22:51:38 +0200 Subject: [PATCH 2/8] Fix sampling from the negative binomial distribution Fix the rate parameter of the gamma distribution used to represent the negative binomial distribution as Poisson mixture. --- src/Numerics/Distributions/NegativeBinomial.cs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Numerics/Distributions/NegativeBinomial.cs b/src/Numerics/Distributions/NegativeBinomial.cs index f6d177693..62aa424bd 100644 --- a/src/Numerics/Distributions/NegativeBinomial.cs +++ b/src/Numerics/Distributions/NegativeBinomial.cs @@ -257,7 +257,7 @@ public static double CDF(double r, double p, double x) /// a sample from the distribution. static int SampleUnchecked(System.Random rnd, double r, double p) { - var lambda = Gamma.SampleUnchecked(rnd, r, p); + var lambda = Gamma.SampleUnchecked(rnd, r, p/(1-p)); var c = Math.Exp(-lambda); var p1 = 1.0; var k = 0; From 8b4cd7a36236d7829ce6625e32971cd837b2272a Mon Sep 17 00:00:00 2001 From: Florian Schwerteck Date: Thu, 29 Sep 2022 23:14:17 +0200 Subject: [PATCH 3/8] Avoid code duplication in sampling from the negative binomial distribution Make the unchecked sample methods of the Poisson distribution internal instead of private. Within the negative binomial sampling, call the Poisson sample method instead of duplicating the implementation. --- src/Numerics/Distributions/NegativeBinomial.cs | 11 +---------- src/Numerics/Distributions/Poisson.cs | 6 +++--- 2 files changed, 4 insertions(+), 13 deletions(-) diff --git a/src/Numerics/Distributions/NegativeBinomial.cs b/src/Numerics/Distributions/NegativeBinomial.cs index 62aa424bd..60dabb36b 100644 --- a/src/Numerics/Distributions/NegativeBinomial.cs +++ b/src/Numerics/Distributions/NegativeBinomial.cs @@ -258,16 +258,7 @@ public static double CDF(double r, double p, double x) static int SampleUnchecked(System.Random rnd, double r, double p) { var lambda = Gamma.SampleUnchecked(rnd, r, p/(1-p)); - var c = Math.Exp(-lambda); - var p1 = 1.0; - var k = 0; - do - { - k = k + 1; - p1 = p1*rnd.NextDouble(); - } - while (p1 >= c); - return k - 1; + return Poisson.SampleUnchecked(rnd, lambda); } static void SamplesUnchecked(System.Random rnd, int[] values, double r, double p) diff --git a/src/Numerics/Distributions/Poisson.cs b/src/Numerics/Distributions/Poisson.cs index 16e0def9e..0c8291244 100644 --- a/src/Numerics/Distributions/Poisson.cs +++ b/src/Numerics/Distributions/Poisson.cs @@ -246,12 +246,12 @@ public static double CDF(double lambda, double x) /// The random source to use. /// The lambda (λ) parameter of the Poisson distribution. Range: λ > 0. /// A random sample from the Poisson distribution. - static int SampleUnchecked(System.Random rnd, double lambda) + internal static int SampleUnchecked(System.Random rnd, double lambda) { return (lambda < 30.0) ? DoSampleShort(rnd, lambda) : DoSampleLarge(rnd, lambda); } - static void SamplesUnchecked(System.Random rnd, int[] values, double lambda) + internal static void SamplesUnchecked(System.Random rnd, int[] values, double lambda) { if (lambda < 30.0) { @@ -300,7 +300,7 @@ static void SamplesUnchecked(System.Random rnd, int[] values, double lambda) } } - static IEnumerable SamplesUnchecked(System.Random rnd, double lambda) + internal static IEnumerable SamplesUnchecked(System.Random rnd, double lambda) { if (lambda < 30.0) { From 6a16def84e48f02c09b2020f2f23f39e9555fe63 Mon Sep 17 00:00:00 2001 From: Florian Schwerteck Date: Fri, 30 Sep 2022 01:17:10 +0200 Subject: [PATCH 4/8] Disallow success probability p = 0 for negative binomial distribution The negative binomial distribution is undefined for p = 0 (asymptotically, it reaches infinity, which cannot be represented as integer). This also avoids unhandled divisions by 0 or other undefined mathematical operations in the properties and an endless loop in the sampling method. --- .../Distributions/NegativeBinomial.cs | 46 +++++++++---------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/src/Numerics/Distributions/NegativeBinomial.cs b/src/Numerics/Distributions/NegativeBinomial.cs index 60dabb36b..6cf08ce9d 100644 --- a/src/Numerics/Distributions/NegativeBinomial.cs +++ b/src/Numerics/Distributions/NegativeBinomial.cs @@ -51,7 +51,7 @@ public class NegativeBinomial : IDiscreteDistribution /// Initializes a new instance of the class. /// /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public NegativeBinomial(double r, double p) { if (!IsValidParameterSet(r, p)) @@ -68,7 +68,7 @@ public NegativeBinomial(double r, double p) /// Initializes a new instance of the class. /// /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. /// The random number generator which is used to draw random samples. public NegativeBinomial(double r, double p, System.Random randomSource) { @@ -97,10 +97,10 @@ public override string ToString() /// Tests whether the provided values are valid parameters for this distribution. /// /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public static bool IsValidParameterSet(double r, double p) { - return r >= 0.0 && p >= 0.0 && p <= 1.0; + return r >= 0.0 && p > 0.0 && p <= 1.0; } /// @@ -109,7 +109,7 @@ public static bool IsValidParameterSet(double r, double p) public double R => _r; /// - /// Gets the probability of success. Range: 0 ≤ p ≤ 1. + /// Gets the probability of success. Range: 0 < p ≤ 1. /// public double P => _p; @@ -202,7 +202,7 @@ public double CumulativeDistribution(double x) /// /// The location in the domain where we want to evaluate the probability mass function. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. /// the probability mass at location . public static double PMF(double r, double p, int k) { @@ -214,11 +214,11 @@ public static double PMF(double r, double p, int k) /// /// The location in the domain where we want to evaluate the log probability mass function. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. /// the log probability mass at location . public static double PMFLn(double r, double p, int k) { - if (!(r >= 0.0 && p >= 0.0 && p <= 1.0)) + if (!(r >= 0.0 && p > 0.0 && p <= 1.0)) { throw new ArgumentException("Invalid parametrization for the distribution."); } @@ -235,12 +235,12 @@ public static double PMFLn(double r, double p, int k) /// /// The location at which to compute the cumulative distribution function. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. /// the cumulative distribution at location . /// public static double CDF(double r, double p, double x) { - if (!(r >= 0.0 && p >= 0.0 && p <= 1.0)) + if (!(r >= 0.0 && p > 0.0 && p <= 1.0)) { throw new ArgumentException("Invalid parametrization for the distribution."); } @@ -253,7 +253,7 @@ public static double CDF(double r, double p, double x) /// /// The random number generator to use. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. /// a sample from the distribution. static int SampleUnchecked(System.Random rnd, double r, double p) { @@ -308,10 +308,10 @@ public IEnumerable Samples() /// /// The random number generator to use. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public static int Sample(System.Random rnd, double r, double p) { - if (!(r >= 0.0 && p >= 0.0 && p <= 1.0)) + if (!(r >= 0.0 && p > 0.0 && p <= 1.0)) { throw new ArgumentException("Invalid parametrization for the distribution."); } @@ -324,10 +324,10 @@ public static int Sample(System.Random rnd, double r, double p) /// /// The random number generator to use. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public static IEnumerable Samples(System.Random rnd, double r, double p) { - if (!(r >= 0.0 && p >= 0.0 && p <= 1.0)) + if (!(r >= 0.0 && p > 0.0 && p <= 1.0)) { throw new ArgumentException("Invalid parametrization for the distribution."); } @@ -341,10 +341,10 @@ public static IEnumerable Samples(System.Random rnd, double r, double p) /// The random number generator to use. /// The array to fill with the samples. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public static void Samples(System.Random rnd, int[] values, double r, double p) { - if (!(r >= 0.0 && p >= 0.0 && p <= 1.0)) + if (!(r >= 0.0 && p > 0.0 && p <= 1.0)) { throw new ArgumentException("Invalid parametrization for the distribution."); } @@ -356,10 +356,10 @@ public static void Samples(System.Random rnd, int[] values, double r, double p) /// Samples a random variable. /// /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public static int Sample(double r, double p) { - if (!(r >= 0.0 && p >= 0.0 && p <= 1.0)) + if (!(r >= 0.0 && p > 0.0 && p <= 1.0)) { throw new ArgumentException("Invalid parametrization for the distribution."); } @@ -371,10 +371,10 @@ public static int Sample(double r, double p) /// Samples a sequence of this random variable. /// /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public static IEnumerable Samples(double r, double p) { - if (!(r >= 0.0 && p >= 0.0 && p <= 1.0)) + if (!(r >= 0.0 && p > 0.0 && p <= 1.0)) { throw new ArgumentException("Invalid parametrization for the distribution."); } @@ -387,10 +387,10 @@ public static IEnumerable Samples(double r, double p) /// /// The array to fill with the samples. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public static void Samples(int[] values, double r, double p) { - if (!(r >= 0.0 && p >= 0.0 && p <= 1.0)) + if (!(r >= 0.0 && p > 0.0 && p <= 1.0)) { throw new ArgumentException("Invalid parametrization for the distribution."); } From fe863e6a23fd4ef882aee1b8b8f7491970608a5f Mon Sep 17 00:00:00 2001 From: Florian Schwerteck Date: Fri, 30 Sep 2022 01:23:22 +0200 Subject: [PATCH 5/8] Handle p = 1 edge case in sampling from the negative binomial distribution --- src/Numerics/Distributions/NegativeBinomial.cs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Numerics/Distributions/NegativeBinomial.cs b/src/Numerics/Distributions/NegativeBinomial.cs index 6cf08ce9d..e872999fb 100644 --- a/src/Numerics/Distributions/NegativeBinomial.cs +++ b/src/Numerics/Distributions/NegativeBinomial.cs @@ -257,6 +257,8 @@ public static double CDF(double r, double p, double x) /// a sample from the distribution. static int SampleUnchecked(System.Random rnd, double r, double p) { + if (p == 1.0) return 0; + var lambda = Gamma.SampleUnchecked(rnd, r, p/(1-p)); return Poisson.SampleUnchecked(rnd, lambda); } From c8cebe2097d6c757bf71f293d978493e8dc6f192 Mon Sep 17 00:00:00 2001 From: Florian Schwerteck Date: Sat, 1 Oct 2022 01:12:08 +0200 Subject: [PATCH 6/8] Fix the parameter descriptions of the negative binomial distribution Use correct XML encoding of the less than operator for the description of the range of the p parameter. --- .../Distributions/NegativeBinomial.cs | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/Numerics/Distributions/NegativeBinomial.cs b/src/Numerics/Distributions/NegativeBinomial.cs index e872999fb..ba74d833c 100644 --- a/src/Numerics/Distributions/NegativeBinomial.cs +++ b/src/Numerics/Distributions/NegativeBinomial.cs @@ -51,7 +51,7 @@ public class NegativeBinomial : IDiscreteDistribution /// Initializes a new instance of the class. /// /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public NegativeBinomial(double r, double p) { if (!IsValidParameterSet(r, p)) @@ -68,7 +68,7 @@ public NegativeBinomial(double r, double p) /// Initializes a new instance of the class. /// /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. /// The random number generator which is used to draw random samples. public NegativeBinomial(double r, double p, System.Random randomSource) { @@ -97,7 +97,7 @@ public override string ToString() /// Tests whether the provided values are valid parameters for this distribution. /// /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public static bool IsValidParameterSet(double r, double p) { return r >= 0.0 && p > 0.0 && p <= 1.0; @@ -109,7 +109,7 @@ public static bool IsValidParameterSet(double r, double p) public double R => _r; /// - /// Gets the probability of success. Range: 0 < p ≤ 1. + /// Gets the probability of success. Range: 0 < p ≤ 1. /// public double P => _p; @@ -202,7 +202,7 @@ public double CumulativeDistribution(double x) /// /// The location in the domain where we want to evaluate the probability mass function. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. /// the probability mass at location . public static double PMF(double r, double p, int k) { @@ -214,7 +214,7 @@ public static double PMF(double r, double p, int k) /// /// The location in the domain where we want to evaluate the log probability mass function. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. /// the log probability mass at location . public static double PMFLn(double r, double p, int k) { @@ -235,7 +235,7 @@ public static double PMFLn(double r, double p, int k) /// /// The location at which to compute the cumulative distribution function. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. /// the cumulative distribution at location . /// public static double CDF(double r, double p, double x) @@ -253,7 +253,7 @@ public static double CDF(double r, double p, double x) /// /// The random number generator to use. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. /// a sample from the distribution. static int SampleUnchecked(System.Random rnd, double r, double p) { @@ -310,7 +310,7 @@ public IEnumerable Samples() /// /// The random number generator to use. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public static int Sample(System.Random rnd, double r, double p) { if (!(r >= 0.0 && p > 0.0 && p <= 1.0)) @@ -326,7 +326,7 @@ public static int Sample(System.Random rnd, double r, double p) /// /// The random number generator to use. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public static IEnumerable Samples(System.Random rnd, double r, double p) { if (!(r >= 0.0 && p > 0.0 && p <= 1.0)) @@ -343,7 +343,7 @@ public static IEnumerable Samples(System.Random rnd, double r, double p) /// The random number generator to use. /// The array to fill with the samples. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public static void Samples(System.Random rnd, int[] values, double r, double p) { if (!(r >= 0.0 && p > 0.0 && p <= 1.0)) @@ -358,7 +358,7 @@ public static void Samples(System.Random rnd, int[] values, double r, double p) /// Samples a random variable. /// /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public static int Sample(double r, double p) { if (!(r >= 0.0 && p > 0.0 && p <= 1.0)) @@ -373,7 +373,7 @@ public static int Sample(double r, double p) /// Samples a sequence of this random variable. /// /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public static IEnumerable Samples(double r, double p) { if (!(r >= 0.0 && p > 0.0 && p <= 1.0)) @@ -389,7 +389,7 @@ public static IEnumerable Samples(double r, double p) /// /// The array to fill with the samples. /// The number of successes (r) required to stop the experiment. Range: r ≥ 0. - /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. + /// The probability (p) of a trial resulting in success. Range: 0 < p ≤ 1. public static void Samples(int[] values, double r, double p) { if (!(r >= 0.0 && p > 0.0 && p <= 1.0)) From 3bd99d4e4550ca5f36ec90cc9d3ee3ddad664443 Mon Sep 17 00:00:00 2001 From: Florian Schwerteck Date: Sat, 1 Oct 2022 00:31:15 +0200 Subject: [PATCH 7/8] Handle p = 1 and r = 0 edge cases in PMF of the negative binomial distribution --- src/Numerics/Distributions/NegativeBinomial.cs | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/Numerics/Distributions/NegativeBinomial.cs b/src/Numerics/Distributions/NegativeBinomial.cs index ba74d833c..3f36de68f 100644 --- a/src/Numerics/Distributions/NegativeBinomial.cs +++ b/src/Numerics/Distributions/NegativeBinomial.cs @@ -206,6 +206,8 @@ public double CumulativeDistribution(double x) /// the probability mass at location . public static double PMF(double r, double p, int k) { + if (p == 1.0 || r == 0.0) return k == 0 ? 1.0 : 0.0; + return Math.Exp(PMFLn(r, p, k)); } @@ -223,6 +225,8 @@ public static double PMFLn(double r, double p, int k) throw new ArgumentException("Invalid parametrization for the distribution."); } + if (p == 1.0 || r == 0.0) return k == 0 ? 0.0 : Double.NegativeInfinity; + return SpecialFunctions.GammaLn(r + k) - SpecialFunctions.GammaLn(r) - SpecialFunctions.GammaLn(k + 1.0) From d05fab707572ab1b888492b7b5a530f31c8997b5 Mon Sep 17 00:00:00 2001 From: Florian Schwerteck Date: Sat, 1 Oct 2022 02:02:19 +0200 Subject: [PATCH 8/8] Adjust the unit tests of the negative binomial distribution Adjust the unit tests after parameter restriction and proper handling of edge cases: - Use p = 0.1 instead of p = 0.0 - Add a test case with p = 0.0 in test for invalid parameters - Add special case in probability validation for p = 1.0 and r = 0.0 --- .../Discrete/NegativeBinomialTests.cs | 25 ++++++++++++------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/src/Numerics.Tests/DistributionTests/Discrete/NegativeBinomialTests.cs b/src/Numerics.Tests/DistributionTests/Discrete/NegativeBinomialTests.cs index 547aaee65..6eca7de38 100644 --- a/src/Numerics.Tests/DistributionTests/Discrete/NegativeBinomialTests.cs +++ b/src/Numerics.Tests/DistributionTests/Discrete/NegativeBinomialTests.cs @@ -45,7 +45,7 @@ public class NegativeBinomialTests /// /// Number of trials. /// Probability of success. - [TestCase(0.0, 0.0)] + [TestCase(0.0, 0.1)] [TestCase(0.1, 0.3)] [TestCase(1.0, 1.0)] public void CanCreateNegativeBinomial(double r, double p) @@ -60,6 +60,7 @@ public void CanCreateNegativeBinomial(double r, double p) /// /// Number of trials. /// Probability of success. + [TestCase(0.0, 0.0)] [TestCase(0.0, double.NaN)] [TestCase(0.0, -1.0)] [TestCase(0.0, 2.0)] @@ -99,7 +100,7 @@ public void ValidateEntropyThrowsNotSupportedException() /// /// Number of trials. /// Probability of success. - [TestCase(0.0, 0.0)] + [TestCase(0.0, 0.1)] [TestCase(0.1, 0.3)] [TestCase(1.0, 1.0)] public void ValidateSkewness(double r, double p) @@ -113,8 +114,8 @@ public void ValidateSkewness(double r, double p) /// /// Number of trials. /// Probability of success. - [TestCase(0.0, 0.0)] - [TestCase(0.3, 0.0)] + [TestCase(0.0, 0.1)] + [TestCase(0.3, 0.1)] [TestCase(1.0, 1.0)] public void ValidateMode(double r, double p) { @@ -165,7 +166,7 @@ public void ValidateMaximum() /// Number of trials. /// Probability of success. /// Input X value. - [TestCase(0.0, 0.0, 0)] + [TestCase(0.0, 0.1, 0)] [TestCase(0.1, 0.3, 1)] [TestCase(1.0, 1.0, 2)] [TestCase(1.0, 1.0, 3)] @@ -173,7 +174,10 @@ public void ValidateMaximum() public void ValidateProbability(double r, double p, int x) { var d = new NegativeBinomial(r, p); - Assert.AreEqual(Math.Exp(SpecialFunctions.GammaLn(r + x) - SpecialFunctions.GammaLn(r) - SpecialFunctions.GammaLn(x + 1.0) + (r * Math.Log(p)) + (x * Math.Log(1.0 - p))), d.Probability(x)); + if (p == 1.0 || r == 0.0) + Assert.AreEqual(x == 0 ? 1.0 : 0.0, d.Probability(x)); + else + Assert.AreEqual(Math.Exp(SpecialFunctions.GammaLn(r + x) - SpecialFunctions.GammaLn(r) - SpecialFunctions.GammaLn(x + 1.0) + (r * Math.Log(p)) + (x * Math.Log(1.0 - p))), d.Probability(x)); } /// @@ -182,7 +186,7 @@ public void ValidateProbability(double r, double p, int x) /// Number of trials. /// Probability of success. /// Input X value. - [TestCase(0.0, 0.0, 0)] + [TestCase(0.0, 0.1, 0)] [TestCase(0.1, 0.3, 1)] [TestCase(1.0, 1.0, 2)] [TestCase(1.0, 1.0, 3)] @@ -190,7 +194,10 @@ public void ValidateProbability(double r, double p, int x) public void ValidateProbabilityLn(double r, double p, int x) { var d = new NegativeBinomial(r, p); - Assert.AreEqual(SpecialFunctions.GammaLn(r + x) - SpecialFunctions.GammaLn(r) - SpecialFunctions.GammaLn(x + 1.0) + (r * Math.Log(p)) + (x * Math.Log(1.0 - p)), d.ProbabilityLn(x)); + if (p == 1.0 || r == 0.0) + Assert.AreEqual(x == 0 ? 0.0 : Double.NegativeInfinity, d.ProbabilityLn(x)); + else + Assert.AreEqual(SpecialFunctions.GammaLn(r + x) - SpecialFunctions.GammaLn(r) - SpecialFunctions.GammaLn(x + 1.0) + (r * Math.Log(p)) + (x * Math.Log(1.0 - p)), d.ProbabilityLn(x)); } /// @@ -199,7 +206,7 @@ public void ValidateProbabilityLn(double r, double p, int x) /// Number of trials. /// Probability of success. /// Input X value. - [TestCase(0.0, 0.0, 0)] + [TestCase(0.0, 0.1, 0)] [TestCase(0.1, 0.3, 1)] [TestCase(1.0, 1.0, 2)] [TestCase(1.0, 1.0, 3)]