diff --git a/src/Numerics.Tests/DistributionTests/Continuous/BetaTests.cs b/src/Numerics.Tests/DistributionTests/Continuous/BetaTests.cs
index 3351107c..fa28e6af 100644
--- a/src/Numerics.Tests/DistributionTests/Continuous/BetaTests.cs
+++ b/src/Numerics.Tests/DistributionTests/Continuous/BetaTests.cs
@@ -405,5 +405,13 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
             Assert.That(dist.InverseCumulativeDistribution(p), Is.EqualTo(x).Within(1e-6));
             Assert.That(Beta.InvCDF(a, b, p), Is.EqualTo(x).Within(1e-6));
         }
+
+        [TestCase(0.001)]
+        public void ProbableNaNWhenABBothSmall(double ab)
+        {
+            Beta dist = new Beta(ab, ab);
+            for (byte i = 0; i < 100; ++i)
+                Assert.That(!double.IsNaN(dist.Sample()),"Generate NaN");
+        }
     }
 }
diff --git a/src/Numerics/Distributions/Beta.cs b/src/Numerics/Distributions/Beta.cs
index 591317fe..67b07709 100644
--- a/src/Numerics/Distributions/Beta.cs
+++ b/src/Numerics/Distributions/Beta.cs
@@ -392,9 +392,22 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a random number from the Beta distribution.</returns>
         internal static double SampleUnchecked(System.Random rnd, double a, double b)
         {
-            var x = Gamma.SampleUnchecked(rnd, a, 1.0);
-            var y = Gamma.SampleUnchecked(rnd, b, 1.0);
-            return x/(x + y);
+            double x, y;
+            if (a == b)
+            {
+                x = Gamma.SampleUnchecked(rnd, a, 1.0);
+                y = Gamma.SampleUnchecked(rnd, b, 1.0);
+                //When a==b (and possibly a==b==0), return value is equally possible to be 0 or 1
+                if (x == 0 && y == 0)
+                    return Bernoulli.Sample(0.5);//In particular, when a==b==0, Beta distribution degradates to Bernoulli distribution.
+            }
+            else
+                do
+                {
+                    x = Gamma.SampleUnchecked(rnd, a, 1.0);
+                    y = Gamma.SampleUnchecked(rnd, b, 1.0);
+                } while (x == 0 && y == 0);//When a!=b, return value is not equally possible to be 0 or 1. Regenerate.
+            return x / (x + y);
         }
 
         internal static void SamplesUnchecked(System.Random rnd, double[] values, double a, double b)