Skip to content

Fix and improve negative binomial distribution #960

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 8 commits into
base: master
Choose a base branch
from
Open

Fix and improve negative binomial distribution #960

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
8a89874
Fix description of the neg binomial distribution in the F# interface
Arlofin Sep 29, 2022
d1fb730
Fix sampling from the negative binomial distribution
Arlofin Sep 29, 2022
8b4cd7a
Avoid code duplication in sampling from the negative binomial distrib…
Arlofin Sep 29, 2022
6a16def
Disallow success probability p = 0 for negative binomial distribution
Arlofin Sep 29, 2022
fe863e6
Handle p = 1 edge case in sampling from the negative binomial distrib…
Arlofin Sep 29, 2022
c8cebe2
Fix the parameter descriptions of the negative binomial distribution
Arlofin Sep 30, 2022
3bd99d4
Handle p = 1 and r = 0 edge cases in PMF of the negative binomial dis…
Arlofin Sep 30, 2022
d05fab7
Adjust the unit tests of the negative binomial distribution
Arlofin Oct 1, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/FSharp/Distributions.fs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
25 changes: 16 additions & 9 deletions src/Numerics.Tests/DistributionTests/Discrete/NegativeBinomialTests.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ public class NegativeBinomialTests
/// </summary>
/// <param name="r">Number of trials.</param>
/// <param name="p">Probability of success.</param>
[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)
Expand All @@ -60,6 +60,7 @@ public void CanCreateNegativeBinomial(double r, double p)
/// </summary>
/// <param name="r">Number of trials.</param>
/// <param name="p">Probability of success.</param>
[TestCase(0.0, 0.0)]
[TestCase(0.0, double.NaN)]
[TestCase(0.0, -1.0)]
[TestCase(0.0, 2.0)]
Expand Down Expand Up @@ -99,7 +100,7 @@ public void ValidateEntropyThrowsNotSupportedException()
/// </summary>
/// <param name="r">Number of trials.</param>
/// <param name="p">Probability of success.</param>
[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)
Expand All @@ -113,8 +114,8 @@ public void ValidateSkewness(double r, double p)
/// </summary>
/// <param name="r">Number of trials.</param>
/// <param name="p">Probability of success.</param>
[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)
{
Expand Down Expand Up @@ -165,15 +166,18 @@ public void ValidateMaximum()
/// <param name="r">Number of trials.</param>
/// <param name="p">Probability of success.</param>
/// <param name="x">Input X value.</param>
[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)]
[TestCase(1.0, 1.0, 5)]
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));
}

/// <summary>
Expand All @@ -182,15 +186,18 @@ public void ValidateProbability(double r, double p, int x)
/// <param name="r">Number of trials.</param>
/// <param name="p">Probability of success.</param>
/// <param name="x">Input X value.</param>
[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)]
[TestCase(1.0, 1.0, 5)]
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));
}

/// <summary>
Expand All @@ -199,7 +206,7 @@ public void ValidateProbabilityLn(double r, double p, int x)
/// <param name="r">Number of trials.</param>
/// <param name="p">Probability of success.</param>
/// <param name="x">Input X value.</param>
[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)]
Expand Down
65 changes: 31 additions & 34 deletions src/Numerics/Distributions/NegativeBinomial.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ public class NegativeBinomial : IDiscreteDistribution
/// Initializes a new instance of the <see cref="NegativeBinomial"/> class.
/// </summary>
/// <param name="r">The number of successes (r) required to stop the experiment. Range: r ≥ 0.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 &lt; p ≤ 1.</param>
public NegativeBinomial(double r, double p)
{
if (!IsValidParameterSet(r, p))
Expand All @@ -68,7 +68,7 @@ public NegativeBinomial(double r, double p)
/// Initializes a new instance of the <see cref="NegativeBinomial"/> class.
/// </summary>
/// <param name="r">The number of successes (r) required to stop the experiment. Range: r ≥ 0.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 &lt; p ≤ 1.</param>
/// <param name="randomSource">The random number generator which is used to draw random samples.</param>
public NegativeBinomial(double r, double p, System.Random randomSource)
{
Expand Down Expand Up @@ -97,10 +97,10 @@ public override string ToString()
/// Tests whether the provided values are valid parameters for this distribution.
/// </summary>
/// <param name="r">The number of successes (r) required to stop the experiment. Range: r ≥ 0.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 &lt; p ≤ 1.</param>
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;
}

/// <summary>
Expand All @@ -109,7 +109,7 @@ public static bool IsValidParameterSet(double r, double p)
public double R => _r;

/// <summary>
/// Gets the probability of success. Range: 0 ≤ p ≤ 1.
/// Gets the probability of success. Range: 0 &lt; p ≤ 1.
/// </summary>
public double P => _p;

Expand Down Expand Up @@ -202,10 +202,12 @@ public double CumulativeDistribution(double x)
/// </summary>
/// <param name="k">The location in the domain where we want to evaluate the probability mass function.</param>
/// <param name="r">The number of successes (r) required to stop the experiment. Range: r ≥ 0.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 &lt; p ≤ 1.</param>
/// <returns>the probability mass at location <paramref name="k"/>.</returns>
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));
}

Expand All @@ -214,15 +216,17 @@ public static double PMF(double r, double p, int k)
/// </summary>
/// <param name="k">The location in the domain where we want to evaluate the log probability mass function.</param>
/// <param name="r">The number of successes (r) required to stop the experiment. Range: r ≥ 0.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 &lt; p ≤ 1.</param>
/// <returns>the log probability mass at location <paramref name="k"/>.</returns>
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.");
}

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)
Expand All @@ -235,12 +239,12 @@ public static double PMFLn(double r, double p, int k)
/// </summary>
/// <param name="x">The location at which to compute the cumulative distribution function.</param>
/// <param name="r">The number of successes (r) required to stop the experiment. Range: r ≥ 0.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 &lt; p ≤ 1.</param>
/// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
/// <seealso cref="CumulativeDistribution"/>
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.");
}
Expand All @@ -253,21 +257,14 @@ public static double CDF(double r, double p, double x)
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="r">The number of successes (r) required to stop the experiment. Range: r ≥ 0.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 &lt; p ≤ 1.</param>
/// <returns>a sample from the distribution.</returns>
static int SampleUnchecked(System.Random rnd, double r, double p)
{
var lambda = Gamma.SampleUnchecked(rnd, r, 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;
if (p == 1.0) return 0;

var lambda = Gamma.SampleUnchecked(rnd, r, p/(1-p));
return Poisson.SampleUnchecked(rnd, lambda);
}

static void SamplesUnchecked(System.Random rnd, int[] values, double r, double p)
Expand Down Expand Up @@ -317,10 +314,10 @@ public IEnumerable<int> Samples()
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="r">The number of successes (r) required to stop the experiment. Range: r ≥ 0.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 &lt; p ≤ 1.</param>
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.");
}
Expand All @@ -333,10 +330,10 @@ public static int Sample(System.Random rnd, double r, double p)
/// </summary>
/// <param name="rnd">The random number generator to use.</param>
/// <param name="r">The number of successes (r) required to stop the experiment. Range: r ≥ 0.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 &lt; p ≤ 1.</param>
public static IEnumerable<int> 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.");
}
Expand All @@ -350,10 +347,10 @@ public static IEnumerable<int> Samples(System.Random rnd, double r, double p)
/// <param name="rnd">The random number generator to use.</param>
/// <param name="values">The array to fill with the samples.</param>
/// <param name="r">The number of successes (r) required to stop the experiment. Range: r ≥ 0.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 &lt; p ≤ 1.</param>
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.");
}
Expand All @@ -365,10 +362,10 @@ public static void Samples(System.Random rnd, int[] values, double r, double p)
/// Samples a random variable.
/// </summary>
/// <param name="r">The number of successes (r) required to stop the experiment. Range: r ≥ 0.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 &lt; p ≤ 1.</param>
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.");
}
Expand All @@ -380,10 +377,10 @@ public static int Sample(double r, double p)
/// Samples a sequence of this random variable.
/// </summary>
/// <param name="r">The number of successes (r) required to stop the experiment. Range: r ≥ 0.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 &lt; p ≤ 1.</param>
public static IEnumerable<int> 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.");
}
Expand All @@ -396,10 +393,10 @@ public static IEnumerable<int> Samples(double r, double p)
/// </summary>
/// <param name="values">The array to fill with the samples.</param>
/// <param name="r">The number of successes (r) required to stop the experiment. Range: r ≥ 0.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
/// <param name="p">The probability (p) of a trial resulting in success. Range: 0 &lt; p ≤ 1.</param>
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.");
}
Expand Down
6 changes: 3 additions & 3 deletions src/Numerics/Distributions/Poisson.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -246,12 +246,12 @@ public static double CDF(double lambda, double x)
/// <param name="rnd">The random source to use.</param>
/// <param name="lambda">The lambda (λ) parameter of the Poisson distribution. Range: λ > 0.</param>
/// <returns>A random sample from the Poisson distribution.</returns>
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)
{
Expand Down Expand Up @@ -300,7 +300,7 @@ static void SamplesUnchecked(System.Random rnd, int[] values, double lambda)
}
}

static IEnumerable<int> SamplesUnchecked(System.Random rnd, double lambda)
internal static IEnumerable<int> SamplesUnchecked(System.Random rnd, double lambda)
{
if (lambda < 30.0)
{
Expand Down