From 97f15820f9051396d7c8d1ee52d4d1dfde4687e5 Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Wed, 15 Jan 2014 19:03:29 +0100
Subject: [PATCH] Distributions: FisherSnedecor.InvCDF

---
 src/Numerics/Distributions/FisherSnedecor.cs  | 35 ++++++++++++++++++-
 .../Continuous/FisherSnedecorTests.cs         | 30 +++++++++++-----
 2 files changed, 56 insertions(+), 9 deletions(-)

diff --git a/src/Numerics/Distributions/FisherSnedecor.cs b/src/Numerics/Distributions/FisherSnedecor.cs
index 4fe5947e..e7cc83b1 100644
--- a/src/Numerics/Distributions/FisherSnedecor.cs
+++ b/src/Numerics/Distributions/FisherSnedecor.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
 {
@@ -259,6 +260,19 @@ namespace MathNet.Numerics.Distributions
         {
             return SpecialFunctions.BetaRegularized(_freedom1/2.0, _freedom2/2.0, _freedom1*x/((_freedom1*x) + _freedom2));
         }
+        
+        /// <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(_freedom1, _freedom2, p);
+        }
 
         /// <summary>
         /// Generates a sample from the <c>FisherSnedecor</c> distribution.
@@ -333,7 +347,26 @@ namespace MathNet.Numerics.Distributions
         {
             if (d1 <= 0.0 || d2 <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
 
-            return SpecialFunctions.BetaRegularized(d1/2.0, d2/2.0, d1*x/((d1*x) + d2));
+            return SpecialFunctions.BetaRegularized(d1/2.0, d2/2.0, d1*x/(d1*x + d2));
+        }
+
+        /// <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="d1">The first degree of freedom (d1) of the distribution. Range: d1 > 0.</param>
+        /// <param name="d2">The second degree of freedom (d2) of the distribution. Range: d2 > 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 d1, double d2, double p)
+        {
+            if (d1 <= 0.0 || d2 <= 0.0) throw new ArgumentOutOfRangeException(Resources.InvalidDistributionParameters);
+
+            return Brent.FindRoot(
+                x => SpecialFunctions.BetaRegularized(d1/2.0, d2/2.0, d1*x/(d1*x + d2)) - p,
+                0, 1000, accuracy: 1e-8);
         }
 
         /// <summary>
diff --git a/src/UnitTests/DistributionTests/Continuous/FisherSnedecorTests.cs b/src/UnitTests/DistributionTests/Continuous/FisherSnedecorTests.cs
index a6002799..5f39557d 100644
--- a/src/UnitTests/DistributionTests/Continuous/FisherSnedecorTests.cs
+++ b/src/UnitTests/DistributionTests/Continuous/FisherSnedecorTests.cs
@@ -425,12 +425,6 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
             ied.Take(5).ToArray();
         }
 
-        /// <summary>
-        /// Validate cumulative distribution.
-        /// </summary>
-        /// <param name="d1">Degrees of freedom 1</param>
-        /// <param name="d2">Degrees of freedom 2</param>
-        /// <param name="x">Input X value</param>
         [TestCase(0.1, 0.1, 1.0)]
         [TestCase(1.0, 0.1, 1.0)]
         [TestCase(10.0, 0.1, 1.0)]
@@ -447,8 +441,28 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous
         {
             var n = new FisherSnedecor(d1, d2);
             double expected = SpecialFunctions.BetaRegularized(d1/2.0, d2/2.0, d1*x/(d2 + (x*d1)));
-            Assert.AreEqual(expected, n.CumulativeDistribution(x));
-            Assert.AreEqual(expected, FisherSnedecor.CDF(d1, d2, x));
+            Assert.That(n.CumulativeDistribution(x), Is.EqualTo(expected));
+            Assert.That(FisherSnedecor.CDF(d1, d2, x), Is.EqualTo(expected));
+        }
+
+        [TestCase(0.1, 0.1, 1.0)]
+        [TestCase(1.0, 0.1, 1.0)]
+        [TestCase(10.0, 0.1, 1.0)]
+        [TestCase(0.1, 1.0, 1.0)]
+        [TestCase(1.0, 1.0, 1.0)]
+        [TestCase(10.0, 1.0, 1.0)]
+        [TestCase(0.1, 0.1, 10.0)]
+        [TestCase(1.0, 0.1, 10.0)]
+        [TestCase(10.0, 0.1, 10.0)]
+        [TestCase(0.1, 1.0, 10.0)]
+        [TestCase(1.0, 1.0, 10.0)]
+        [TestCase(10.0, 1.0, 10.0)]
+        public void ValidateInverseCumulativeDistribution(double d1, double d2, double x)
+        {
+            var n = new FisherSnedecor(d1, d2);
+            double p = SpecialFunctions.BetaRegularized(d1/2.0, d2/2.0, d1*x/(d2 + (x*d1)));
+            Assert.That(n.InverseCumulativeDistribution(p), Is.EqualTo(x).Within(1e-8));
+            Assert.That(FisherSnedecor.InvCDF(d1, d2, p), Is.EqualTo(x).Within(1e-8));
         }
     }
 }