From b7ab86af1b54569afe907c21a5206e799926c857 Mon Sep 17 00:00:00 2001
From: Tom Minka <8955276+tminka@users.noreply.github.com>
Date: Fri, 7 Jun 2019 20:36:09 +0100
Subject: [PATCH] InnerQuantiles and OuterQuantiles throw on infinite
quantiles. InnerQuantiles.GetLogProb and GetQuantile are non-decreasing.
---
src/Runtime/Core/Maths/DenseVector.cs | 3 +-
src/Runtime/Distributions/InnerQuantiles.cs | 15 +++--
src/Runtime/Distributions/OuterQuantiles.cs | 22 ++++---
.../Distributions/QuantileEstimator.cs | 4 +-
test/TestApp/Program.cs | 2 -
test/Tests/QuantileTests.cs | 64 +++++++++++++++++--
6 files changed, 88 insertions(+), 22 deletions(-)
diff --git a/src/Runtime/Core/Maths/DenseVector.cs b/src/Runtime/Core/Maths/DenseVector.cs
index 7dc62bf5..eb40b7a9 100644
--- a/src/Runtime/Core/Maths/DenseVector.cs
+++ b/src/Runtime/Core/Maths/DenseVector.cs
@@ -2585,9 +2585,8 @@ namespace Microsoft.ML.Probabilistic.Math
// Therefore scale <= double.MaxValue/Y.Count/abs(Y[j]*US[j,k]).
double absY = System.Math.Abs(Y[j]);
double absU = System.Math.Abs(US[j, k]);
- bool Ygt1 = (absY > 1);
double thisUpperBound;
- if (Ygt1)
+ if (absY > 1)
{
thisUpperBound = limit / absY / absU;
}
diff --git a/src/Runtime/Distributions/InnerQuantiles.cs b/src/Runtime/Distributions/InnerQuantiles.cs
index ab9d688b..7fe7b2eb 100644
--- a/src/Runtime/Distributions/InnerQuantiles.cs
+++ b/src/Runtime/Distributions/InnerQuantiles.cs
@@ -48,6 +48,8 @@ namespace Microsoft.ML.Probabilistic.Distributions
{
this.quantiles[i] = canGetQuantile.GetQuantile((i + 1.0) / (quantileCount + 1.0));
}
+ OuterQuantiles.AssertFinite(quantiles, nameof(canGetQuantile));
+ OuterQuantiles.AssertNondecreasing(quantiles, nameof(canGetQuantile));
lowerGaussian = GetLowerGaussian(quantiles);
upperGaussian = GetUpperGaussian(quantiles);
}
@@ -72,16 +74,17 @@ namespace Microsoft.ML.Probabilistic.Distributions
return quantiles;
}
+ ///
public double GetProbLessThan(double x)
{
+ int n = quantiles.Length;
if (x < quantiles[0])
{
- return lowerGaussian.GetProbLessThan(x);
+ return Math.Min(lowerGaussian.GetProbLessThan(x), 1.0/(n+1));
}
- int n = quantiles.Length;
if (x > quantiles[n - 1])
{
- return upperGaussian.GetProbLessThan(x);
+ return Math.Max(upperGaussian.GetProbLessThan(x), n/(n+1.0));
}
return GetProbLessThan(x, quantiles);
}
@@ -156,7 +159,7 @@ namespace Microsoft.ML.Probabilistic.Distributions
double mean, stddev;
GetGaussianFromQuantiles(quantiles[0], p0, quantiles[i], p1, out mean, out stddev);
Gaussian result = Gaussian.FromMeanAndVariance(mean, stddev * stddev);
- if (!result.IsProper() || double.IsNaN(result.GetMean()) || double.IsInfinity(result.GetMean()))
+ if (result.IsPointMass || !result.IsProper() || double.IsNaN(result.GetMean()) || double.IsInfinity(result.GetMean()))
{
return Gaussian.PointMass(quantiles[0]);
}
@@ -184,6 +187,10 @@ namespace Microsoft.ML.Probabilistic.Distributions
double p1 = (i + 1.0) / (n + 1);
double mean, stddev;
GetGaussianFromQuantiles(quantiles[n - 1], p0, quantiles[i], p1, out mean, out stddev);
+ if (double.IsNaN(mean) || double.IsInfinity(mean))
+ {
+ return Gaussian.PointMass(quantiles[n - 1]);
+ }
Gaussian result = Gaussian.FromMeanAndVariance(mean, stddev * stddev);
if (!result.IsProper() || double.IsNaN(result.GetMean()) || double.IsInfinity(result.GetMean()))
{
diff --git a/src/Runtime/Distributions/OuterQuantiles.cs b/src/Runtime/Distributions/OuterQuantiles.cs
index 1cb47824..d1083a9b 100644
--- a/src/Runtime/Distributions/OuterQuantiles.cs
+++ b/src/Runtime/Distributions/OuterQuantiles.cs
@@ -32,8 +32,8 @@ namespace Microsoft.ML.Probabilistic.Distributions
{
for (int i = 0; i < array.Length; i++)
{
- if (double.IsInfinity(array[i])) throw new ArgumentOutOfRangeException(paramName, $"{paramName}[{i}] {array[i]}");
- if (double.IsNaN(array[i])) throw new ArgumentOutOfRangeException(paramName, $"{paramName}[{i}] {array[i]}");
+ if (double.IsInfinity(array[i])) throw new ArgumentOutOfRangeException(paramName, array[i], $"Array element is infinite: {paramName}[{i}]={array[i]}");
+ if (double.IsNaN(array[i])) throw new ArgumentOutOfRangeException(paramName, array[i], $"Array element is NaN: {paramName}[{i}]={array[i]}");
}
}
@@ -41,7 +41,7 @@ namespace Microsoft.ML.Probabilistic.Distributions
{
for (int i = 1; i < array.Length; i++)
{
- if (array[i] < array[i - 1]) throw new ArgumentException($"Array is not non-decreasing: {paramName}[{i}] {array[i]} < {paramName}[{i - 1}] {array[i - 1]}", paramName);
+ if (array[i] < array[i - 1]) throw new ArgumentException($"Decreasing array elements: {paramName}[{i}] {array[i]} < {paramName}[{i - 1}] {array[i - 1]}", paramName);
}
}
@@ -124,19 +124,23 @@ namespace Microsoft.ML.Probabilistic.Distributions
///
/// Returns the largest quantile such that ((quantile - lowerItem)/(upperItem - lowerItem) + lowerIndex)/(n-1) <= probability.
///
- ///
+ /// A number between 0 and 1.
///
- ///
- ///
- ///
+ /// Must be finite.
+ /// Must be finite and at least lowerItem.
+ /// Must be greater than 1
///
internal static double GetQuantile(double probability, double lowerIndex, double lowerItem, double upperItem, long n)
{
- if (n == 1) throw new ArgumentOutOfRangeException("n == 1");
+ if(probability < 0) throw new ArgumentOutOfRangeException(nameof(probability), probability, $"{nameof(probability)} < 0");
+ if (probability > 1) throw new ArgumentOutOfRangeException(nameof(probability), probability, $"{nameof(probability)} > 1");
+ if (n <= 1) throw new ArgumentOutOfRangeException(nameof(n), n, "n <= 1");
double pos = MMath.LargestDoubleProduct(n - 1, probability);
double frac = MMath.LargestDoubleSum(-lowerIndex, pos);
+ if (upperItem < lowerItem) throw new ArgumentOutOfRangeException(nameof(upperItem), upperItem, $"{nameof(upperItem)} ({upperItem}) < {nameof(lowerItem)} ({lowerItem})");
+ if (upperItem == lowerItem) return lowerItem;
+ // The above check ensures diff > 0
double diff = upperItem - lowerItem;
- if (diff == 0) return lowerItem;
double offset = MMath.LargestDoubleProduct(diff, frac);
return MMath.LargestDoubleSum(lowerItem, offset);
}
diff --git a/src/Runtime/Distributions/QuantileEstimator.cs b/src/Runtime/Distributions/QuantileEstimator.cs
index e54a3322..b54759f2 100644
--- a/src/Runtime/Distributions/QuantileEstimator.cs
+++ b/src/Runtime/Distributions/QuantileEstimator.cs
@@ -137,10 +137,12 @@ namespace Microsoft.ML.Probabilistic.Distributions
if (probability > 1.0) throw new ArgumentOutOfRangeException(nameof(probability), "probability > 1.0");
// compute the min and max of the retained items
double lowerBound = double.PositiveInfinity, upperBound = double.NegativeInfinity;
+ bool countGreaterThanZero = false;
for (int bufferIndex = 0; bufferIndex < buffers.Length; bufferIndex++)
{
double[] buffer = buffers[bufferIndex];
int count = countInBuffer[bufferIndex];
+ if (count > 0) countGreaterThanZero = true;
for (int i = 0; i < count; i++)
{
double item = buffer[i];
@@ -150,7 +152,7 @@ namespace Microsoft.ML.Probabilistic.Distributions
}
if (probability == 1.0) return MMath.NextDouble(upperBound);
if (probability == 0.0) return lowerBound;
- if (double.IsPositiveInfinity(lowerBound)) throw new Exception("QuantileEstimator has no data");
+ if (!countGreaterThanZero) throw new Exception("QuantileEstimator has no data");
if (lowerBound == upperBound) return upperBound;
// use bisection
while (true)
diff --git a/test/TestApp/Program.cs b/test/TestApp/Program.cs
index 059e9f77..c545039c 100644
--- a/test/TestApp/Program.cs
+++ b/test/TestApp/Program.cs
@@ -104,8 +104,6 @@ namespace TestApp
#endif
watch.Stop();
Console.WriteLine("elapsed time = {0}ms", watch.ElapsedMilliseconds);
-
- Console.ReadKey();
}
}
}
diff --git a/test/Tests/QuantileTests.cs b/test/Tests/QuantileTests.cs
index 212e2f9d..5e300d48 100644
--- a/test/Tests/QuantileTests.cs
+++ b/test/Tests/QuantileTests.cs
@@ -36,6 +36,30 @@ namespace Microsoft.ML.Probabilistic.Tests
});
}
+ [Fact]
+ public void InnerQuantiles_InfinityTest()
+ {
+ Assert.Throws(() =>
+ {
+ var inner = new InnerQuantiles(new double[] { double.PositiveInfinity });
+ });
+ Assert.Throws(() =>
+ {
+ var est = new QuantileEstimator(0.1);
+ est.Add(double.PositiveInfinity);
+ //est.Add(double.NegativeInfinity);
+ var inner = new InnerQuantiles(10, est);
+ });
+ }
+
+ [Fact]
+ public void QuantileEstimator_InfinityTest()
+ {
+ var est = new QuantileEstimator(0.1);
+ est.Add(double.PositiveInfinity);
+ Assert.Equal(double.PositiveInfinity, est.GetQuantile(0.1));
+ }
+
[Fact]
public void QuantileEstimator_SinglePointIsMedian()
{
@@ -63,6 +87,38 @@ namespace Microsoft.ML.Probabilistic.Tests
inner.GetQuantile((double)quantiles.Length / (quantiles.Length + 1));
}
+ [Fact]
+ public void InnerQuantiles_GetProbLessThan_IsIncreasing()
+ {
+ var quantiles = new double[] { -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.166666666666667, -7.0149625935162021, -6.741895261845384, -6.4688279301745659, -6.1957605985037265, -5.9226932668329084, -5.64962593516209, -5.3765586034912722, -5.1034912718204541, -4.8304239401496138, -4.5573566084787958, -4.2842892768079777, -4.01122194513716, -3.738154613466341, -3.4650872817955016, -3.1920199501246835, -2.9189526184538654, -2.6458852867830256, -2.3728179551122075, -2.0997506234413894, -1.8266832917705711, -1.553615960099753, -1.2805486284289136, -1.0074812967580955, -0.7344139650872773, -0.46134663341645915, -0.18827930174564103, 0.084788029925200167, 0.35785536159601822, 0.63092269326683625, 0.903990024937654, 1.1770573566084719, 1.4501246882793115, 1.7231920199501294, 1.996259351620969, 2.2693266832917867, 2.5423940149626052, 2.8154613466334251, 3.0885286783042427, 3.3615960099750803, 3.6346633416458989, 3.9077306733167165, 4.1807980049875351, 4.4538653366583523, 4.7269326683291935 };
+ var innerQuantiles = new InnerQuantiles(quantiles);
+
+ const double left = -7.1666666666666679;
+ const double right = -7.166666666666667;
+ double probLeft = innerQuantiles.GetProbLessThan(left);
+ double probRight = innerQuantiles.GetProbLessThan(right);
+
+ Assert.True(left < right, "It is given that left < right");
+ Assert.True(probLeft <= probRight, "CDF must be non-decreasing");
+ }
+
+ [Fact]
+ public void InnerQuantiles_GetQuantile_IsIncreasing()
+ {
+ var quantiles = new double[] { 0, 4.94065645841247E-324, 4.94065645841247E-324, 4.0000000000000009 };
+ var innerQuantiles = new InnerQuantiles(quantiles);
+
+ int n = 100;
+ double previousQuantile = double.MinValue;
+ for (int i = 0; i < n; i++)
+ {
+ double probability = (i + 1.0) / (n + 1);
+ double quantile = innerQuantiles.GetQuantile(probability);
+ Assert.True(quantile >= previousQuantile);
+ previousQuantile = quantile;
+ }
+ }
+
///
/// Test that QuantileEstimator can handle billions of items.
///
@@ -209,7 +265,7 @@ namespace Microsoft.ML.Probabilistic.Tests
Assert.Equal(outer.GetQuantile(0.3), middle);
Assert.Equal(outer.GetQuantile(0.5), middle);
Assert.Equal(outer.GetQuantile(0.7), middle);
- CheckGetQuantile(inner, inner, (int)Math.Ceiling(100.0/8), (int)Math.Floor(100.0*7/8));
+ CheckGetQuantile(inner, inner, (int)Math.Ceiling(100.0 / 8), (int)Math.Floor(100.0 * 7 / 8));
var est = new QuantileEstimator(0.01);
est.AddRange(x);
Assert.Equal(est.GetQuantile(0.3), middle);
@@ -239,7 +295,7 @@ namespace Microsoft.ML.Probabilistic.Tests
Assert.Equal(next, outer.GetQuantile(1.0));
CheckGetQuantile(outer, outer);
var inner = new InnerQuantiles(5, outer);
- CheckGetQuantile(inner, inner, (int)Math.Ceiling(100.0/6), (int)Math.Floor(100.0*5/6));
+ CheckGetQuantile(inner, inner, (int)Math.Ceiling(100.0 / 6), (int)Math.Floor(100.0 * 5 / 6));
var est = new QuantileEstimator(0.01);
est.Add(first, 2);
est.Add(second, 2);
@@ -276,7 +332,7 @@ namespace Microsoft.ML.Probabilistic.Tests
Assert.Equal(data[2], outer.GetQuantile(0.3));
CheckGetQuantile(outer, outer);
var inner = new InnerQuantiles(7, outer);
- CheckGetQuantile(inner, inner, (int)Math.Ceiling(100.0/8), (int)Math.Floor(100.0*7/8));
+ CheckGetQuantile(inner, inner, (int)Math.Ceiling(100.0 / 8), (int)Math.Floor(100.0 * 7 / 8));
}
internal void CheckGetQuantile(CanGetQuantile canGetQuantile, CanGetProbLessThan canGetProbLessThan, int minPercent = 0, int maxPercent = 100)
@@ -323,7 +379,7 @@ namespace Microsoft.ML.Probabilistic.Tests
[Fact]
public void QuantileEstimatorTest()
{
- foreach(double maximumError in new[] { 0.05, 0.01, 0.005, 0.001 })
+ foreach (double maximumError in new[] { 0.05, 0.01, 0.005, 0.001 })
{
int n = (int)(2.0 / maximumError);
QuantileEstimatorTester(maximumError, n);