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)
 
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
         /// </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)
@@ -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)]
@@ -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)
@@ -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)
         {
@@ -165,7 +166,7 @@ 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)]
@@ -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));
         }
 
         /// <summary>
@@ -182,7 +186,7 @@ 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)]
@@ -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));
         }
 
         /// <summary>
@@ -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)]
diff --git a/src/Numerics/Distributions/NegativeBinomial.cs b/src/Numerics/Distributions/NegativeBinomial.cs
index f6d177693..3f36de68f 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 <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))
@@ -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)
         {
@@ -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>
@@ -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;
 
@@ -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));
         }
 
@@ -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)
@@ -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.");
             }
@@ -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)
@@ -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.");
             }
@@ -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.");
             }
@@ -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.");
             }
@@ -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.");
             }
@@ -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.");
             }
@@ -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.");
             }
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)
         /// <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)
             {
@@ -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)
             {