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..8eec4c14 100644
--- a/src/Numerics/Distributions/Beta.cs
+++ b/src/Numerics/Distributions/Beta.cs
@@ -394,7 +394,7 @@ namespace MathNet.Numerics.Distributions
         {
             var x = Gamma.SampleUnchecked(rnd, a, 1.0);
             var y = Gamma.SampleUnchecked(rnd, b, 1.0);
-            return x/(x + y);
+            return x == 0 ? 0 : x / (x + y);
         }
 
         internal static void SamplesUnchecked(System.Random rnd, double[] values, double a, double b)