From 1d3aa6734d0dd6ffff00c72c6bf1983c343f8e1d Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Sat, 4 Jan 2014 21:28:04 +0100
Subject: [PATCH] Distributions: Beta.InvCDF

---
 src/Numerics/Distributions/Beta.cs            | 33 ++++++++++++++++++-
 .../DistributionTests/Continuous/BetaTests.cs | 27 ++++++++-------
 2 files changed, 48 insertions(+), 12 deletions(-)

diff --git a/src/Numerics/Distributions/Beta.cs b/src/Numerics/Distributions/Beta.cs
index 38a4cd50..d7fa9013 100644
--- a/src/Numerics/Distributions/Beta.cs
+++ b/src/Numerics/Distributions/Beta.cs
@@ -32,6 +32,7 @@ using System;
 using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.RootFinding;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -326,6 +327,19 @@ namespace MathNet.Numerics.Distributions
             return CDF(_shapeA, _shapeB, x);
         }
 
+        /// <summary>
+        /// Computes the inverse of the cumulative distribution function (InvCDF) for the distribution
+        /// at the given probability. This is also known as the quantile or percent point function.
+        /// </summary>
+        /// <param name="p">The location at which to compute the inverse cumulative density.</param>
+        /// <returns>the inverse cumulative density at <paramref name="p"/>.</returns>
+        /// <seealso cref="InvCDF"/>
+        /// <remarks>WARNING: currently not an explicit implementation, hence slow and unreliable.</remarks>
+        public double InverseCumulativeDistribution(double p)
+        {
+            return InvCDF(_shapeA, _shapeB, p);
+        }
+
         /// <summary>
         /// Generates a sample from the Beta distribution.
         /// </summary>
@@ -458,7 +472,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         /// <param name="x">The location at which to compute the cumulative distribution function.</param>
         /// <param name="a">The α shape parameter of the Beta distribution. Range: α ≥ 0.</param>
-        /// <param name="b">The β shape parameter of the Beta distribution. Range: β ≥ 0.</param>>
+        /// <param name="b">The β shape parameter of the Beta distribution. Range: β ≥ 0.</param>
         /// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
         /// <seealso cref="CumulativeDistribution"/>
         public static double CDF(double a, double b, double x)
@@ -500,6 +514,23 @@ namespace MathNet.Numerics.Distributions
             return SpecialFunctions.BetaRegularized(a, b, x);
         }
 
+        /// <summary>
+        /// Computes the inverse of the cumulative distribution function (InvCDF) for the distribution
+        /// at the given probability. This is also known as the quantile or percent point function.
+        /// </summary>
+        /// <param name="p">The location at which to compute the inverse cumulative density.</param>
+        /// <param name="a">The α shape parameter of the Beta distribution. Range: α ≥ 0.</param>
+        /// <param name="b">The β shape parameter of the Beta distribution. Range: β ≥ 0.</param>
+        /// <returns>the inverse cumulative density at <paramref name="p"/>.</returns>
+        /// <seealso cref="InverseCumulativeDistribution"/>
+        /// <remarks>WARNING: currently not an explicit implementation, hence slow and unreliable.</remarks>
+        public static double InvCDF(double a, double b, double p)
+        {
+            if (a < 0.0 || b < 0.0 || p < 0.0 || p > 1.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
+
+            return Brent.FindRoot(x => SpecialFunctions.BetaRegularized(a, b, x) - p, 0.0, 1.0, accuracy: 1e-8);
+        }
+
         /// <summary>
         /// Generates a sample from the distribution.
         /// </summary>
diff --git a/src/UnitTests/DistributionTests/Continuous/BetaTests.cs b/src/UnitTests/DistributionTests/Continuous/BetaTests.cs
index 3e129758..e914b7a8 100644
--- a/src/UnitTests/DistributionTests/Continuous/BetaTests.cs
+++ b/src/UnitTests/DistributionTests/Continuous/BetaTests.cs
@@ -418,13 +418,6 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
             AssertHelpers.AlmostEqualRelative(pdfln, Beta.PDFLn(a, b, x), 13);
         }
 
-        /// <summary>
-        /// Validate cumulative distribution.
-        /// </summary>
-        /// <param name="a">Parameter A.</param>
-        /// <param name="b">Parameter B.</param>
-        /// <param name="x">Input value X.</param>
-        /// <param name="cdf">Cumulative distribution value.</param>
         [TestCase(0.0, 0.0, 0.0, 0.5)]
         [TestCase(0.0, 0.0, 0.5, 0.5)]
         [TestCase(0.0, 0.0, 1.0, 1.0)]
@@ -455,11 +448,23 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
         [TestCase(Double.PositiveInfinity, 0.0, 0.0, 0.0)]
         [TestCase(Double.PositiveInfinity, 0.0, 0.5, 0.0)]
         [TestCase(Double.PositiveInfinity, 0.0, 1.0, 1.0)]
-        public void ValidateCumulativeDistribution(double a, double b, double x, double cdf)
+        public void ValidateCumulativeDistribution(double a, double b, double x, double p)
         {
-            var n = new Beta(a, b);
-            AssertHelpers.AlmostEqualRelative(cdf, n.CumulativeDistribution(x), 13);
-            AssertHelpers.AlmostEqualRelative(cdf, Beta.CDF(a, b, x), 13);
+            var dist = new Beta(a, b);
+            Assert.That(dist.CumulativeDistribution(x), Is.EqualTo(p).Within(1e-13));
+            Assert.That(Beta.CDF(a, b, x), Is.EqualTo(p).Within(1e-13));
+        }
+
+        [TestCase(1.0, 1.0, 1.0, 1.0)]
+        [TestCase(9.0, 1.0, 0.0, 0.0)]
+        [TestCase(9.0, 1.0, 0.5, 0.001953125)]
+        [TestCase(9.0, 1.0, 1.0, 1.0)]
+        [TestCase(5.0, 100, 0.0, 0.0)]
+        public void ValidateInverseCumulativeDistribution(double a, double b, double x, double p)
+        {
+            var dist = new Beta(a, b);
+            Assert.That(dist.InverseCumulativeDistribution(p), Is.EqualTo(x).Within(1e-6));
+            Assert.That(Beta.InvCDF(a, b, p), Is.EqualTo(x).Within(1e-6));
         }
     }
 }