MMath.Median fix. (#43)

Gamma supports GetQuantile when Shape=1.

src/Runtime/Core/Maths/SpecialFunctions.cs

@@ -3527,7 +3527,7 @@ else if (m < 20.0 - 60.0/11.0 * s) {
             if (length % 2 == 0)
             {
                 // average the two middle elements
-                return (a[middle - 1] + a[middle]) / 2;
+                return Average(a[middle - 1], a[middle]);
             }
             else
             {

src/Runtime/Distributions/Beta.cs

@@ -34,7 +34,7 @@ namespace Microsoft.ML.Probabilistic.Distributions
                          SettableToRatio<Beta>, SettableToPower<Beta>, SettableToWeightedSum<Beta>, CanGetLogAverageOf<Beta>, CanGetLogAverageOfPower<Beta>,
                          CanGetAverageLog<Beta>, CanGetLogNormalizer,
                          Sampleable<double>, CanGetMean<double>, CanGetVariance<double>, CanGetMeanAndVarianceOut<double, double>, CanSetMeanAndVariance<double, double>,
-                         CanGetMode<double>
+                         CanGetMode<double>, CanGetProbLessThan
     {
         /// <summary>
         /// True count

src/Runtime/Distributions/Gamma.cs

@@ -38,7 +38,8 @@ namespace Microsoft.ML.Probabilistic.Distributions
                           Sampleable<double>, CanGetMean<double>, CanGetVariance<double>,
                           CanGetMeanAndVarianceOut<double, double>, CanSetMeanAndVariance<double, double>,
                           CanGetLogAverageOf<Gamma>, CanGetLogAverageOfPower<Gamma>,
-                          CanGetAverageLog<Gamma>, CanGetLogNormalizer, CanGetMode<double>
+                          CanGetAverageLog<Gamma>, CanGetLogNormalizer, CanGetMode<double>,
+                          CanGetProbLessThan, CanGetQuantile
     {
         /// <summary>
         /// Rate parameter for the distribution
@@ -455,6 +456,34 @@ namespace Microsoft.ML.Probabilistic.Distributions
                 return MMath.GammaLower(Shape, x * Rate);
         }
 
+        /// <summary>
+        /// Returns the value x such that GetProbLessThan(x) == probability.
+        /// </summary>
+        /// <param name="probability">A real number in [0,1].</param>
+        /// <returns></returns>
+        public double GetQuantile(double probability)
+        {
+            if (probability < 0) throw new ArgumentOutOfRangeException("probability < 0");
+            if (probability > 1) throw new ArgumentOutOfRangeException("probability > 1");
+            if (this.IsPointMass)
+            {
+                return (probability == 1.0) ? MMath.NextDouble(this.Point) : this.Point;
+            }
+            else if (!IsProper())
+            {
+                throw new ImproperDistributionException(this);
+            }
+            else if (Shape == 1)
+            {
+                // cdf is 1 - exp(-x*rate)
+                return -Math.Log(1 - probability) / Rate;
+            }
+            else
+            {
+                throw new NotImplementedException("Shape != 1 is not implemented");
+            }
+        }
+
         /// <summary>
         /// Asks whether the instance is a point mass
         /// </summary>

src/Runtime/Distributions/Gaussian.cs

@@ -273,6 +273,8 @@ namespace Microsoft.ML.Probabilistic.Distributions
         /// <returns></returns>
         public double GetQuantile(double probability)
         {
+            if (probability < 0) throw new ArgumentOutOfRangeException("probability < 0");
+            if (probability > 1) throw new ArgumentOutOfRangeException("probability > 1");
             if (this.IsPointMass)
             {
                 return (probability == 1.0) ? MMath.NextDouble(this.Point) : this.Point;

test/Tests/Core/SpecialFunctionsTest.cs

@@ -23599,5 +23599,13 @@ exp(x*x/4)*pcfu(0.5+n,-x)
             };
             CheckFunctionValues("ulp", MMath.Ulp, Ulp_pairs);
         }
+
+        [Fact]
+        public void MedianTest()
+        {
+            Assert.Equal(double.MaxValue, MMath.Median(new[] { double.MaxValue }));
+            Assert.Equal(double.MaxValue, MMath.Median(new[] { double.MaxValue, double.MaxValue }));
+            Assert.Equal(3, MMath.Median(new[] { double.MaxValue, 3, double.MinValue }));
+        }
     }
 }

test/Tests/Distributions/DistributionTests.cs

@@ -1399,6 +1399,17 @@ namespace Microsoft.ML.Probabilistic.Tests
             Assert.True(MMath.AbsDiff((r * denominator).GetMean(), numerator.GetMean(), 1e-8) < 1e-10);
         }
 
+        [Fact]
+        public void GammaGetProbLessThanTest()
+        {
+            // exponential distribution with density exp(-x/m)/m and cdf 1-exp(-x/m)
+            double m = 2.3;
+            Gamma g = new Gamma(1.0, m);
+            double median = -m * System.Math.Log(0.5);
+            Assert.Equal(0.5, g.GetProbLessThan(median), 1e-4);
+            Assert.Equal(median, g.GetQuantile(0.5));
+        }
+
         [Fact]
         public void GammaModeTest()
         {