Statistics: Median direct implementation instead of R8-compatible 0.5-quantile

src/Numerics/Statistics/ArrayStatistics.cs

@@ -279,16 +279,17 @@ namespace MathNet.Numerics.Statistics
 
         /// <summary>
         /// Estimates the median value from the unsorted data array.
-        /// Approximately median-unbiased regardless of the sample distribution (R8).
         /// WARNING: Works inplace and can thus causes the data array to be reordered.
         /// </summary>
         /// <param name="data">Sample array, no sorting is assumed. Will be reordered.</param>
         public static double MedianInplace(double[] data)
         {
-            return QuantileInplace(data, 0.5d);
+            var k = data.Length/2;
+            return data.Length.IsOdd()
+                ? SelectInplace(data, k)
+                : (SelectInplace(data, k - 1) + SelectInplace(data, k))/2.0;
         }
 
-
         /// <summary>
         /// Estimates the p-Percentile value from the unsorted data array.
         /// If a non-integer Percentile is needed, use Quantile instead.

src/Numerics/Statistics/SortedArrayStatistics.cs

@@ -81,7 +81,12 @@ namespace MathNet.Numerics.Statistics
         /// <param name="data">Sample array, must be sorted ascendingly.</param>
         public static double Median(double[] data)
         {
-            return Quantile(data, 0.5d);
+            if (data.Length == 0) return double.NaN;
+
+            var k = data.Length/2;
+            return data.Length.IsOdd()
+                ? data[k]
+                : (data[k - 1] + data[k])/2.0;
         }
 
         /// <summary>

src/UnitTests/StatisticsTests/StatisticsTests.cs

@@ -912,6 +912,20 @@ namespace MathNet.Numerics.UnitTests.StatisticsTests
             Assert.AreEqual(21.578697, a.Variance(), 1e-5);
             Assert.AreEqual(21.578231, a.PopulationVariance(), 1e-5);
         }
+
+        [Test]
+        public void MedianIsRobustOnCloseInfinities()
+        {
+            Assert.That(Statistics.Median(new[] { 2.0, double.NegativeInfinity, double.PositiveInfinity }), Is.EqualTo(2.0));
+            Assert.That(Statistics.Median(new[] { 2.0, double.NegativeInfinity, 3.0, double.PositiveInfinity }), Is.EqualTo(2.5));
+            Assert.That(ArrayStatistics.MedianInplace(new[] { 2.0, double.NegativeInfinity, double.PositiveInfinity}), Is.EqualTo(2.0));
+            Assert.That(ArrayStatistics.MedianInplace(new[] { double.NegativeInfinity, 2.0, double.PositiveInfinity }), Is.EqualTo(2.0));
+            Assert.That(ArrayStatistics.MedianInplace(new[] { double.NegativeInfinity, double.PositiveInfinity, 2.0 }), Is.EqualTo(2.0));
+            Assert.That(ArrayStatistics.MedianInplace(new[] { double.NegativeInfinity, 2.0, 3.0, double.PositiveInfinity }), Is.EqualTo(2.5));
+            Assert.That(ArrayStatistics.MedianInplace(new[] { double.NegativeInfinity, 2.0, double.PositiveInfinity, 3.0,  }), Is.EqualTo(2.5));
+            Assert.That(SortedArrayStatistics.Median(new[] { double.NegativeInfinity, 2.0, double.PositiveInfinity }), Is.EqualTo(2.0));
+            Assert.That(SortedArrayStatistics.Median(new[] { double.NegativeInfinity, 2.0, 3.0, double.PositiveInfinity }), Is.EqualTo(2.5));
+        }
     }
 }