You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2021/12/16 12:17:47 UTC
[commons-statistics] 01/05: Ensure default inverse probability is robust to distribution truncation
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-statistics.git
commit 5bc2f6c514976212713a790f041ab508800538db
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 16 10:53:22 2021 +0000
Ensure default inverse probability is robust to distribution truncation
This was discovered as the F-distribution with numerator DF=1 and
denominator DF=2 computed a survival probability:
sf(infinity) = 0.0
sf(max_value) = 5.56e-309
If the inverse SF was called with a p-value lower than 5.56e-309 (e.g.
Double.MIN_VALUE) the required value of x is outside the finite range of
a double.
---
.../AbstractContinuousDistribution.java | 32 +++++-
.../AbstractContinuousDistributionTest.java | 109 +++++++++++++++++++++
2 files changed, 138 insertions(+), 3 deletions(-)
diff --git a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractContinuousDistribution.java b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractContinuousDistribution.java
index 4d8b5f8..5e26812 100644
--- a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractContinuousDistribution.java
+++ b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractContinuousDistribution.java
@@ -269,15 +269,39 @@ abstract class AbstractContinuousDistribution
}
}
+ // Here the bracket [lower, upper] uses finite values. If the support
+ // is infinite the bracket can truncate the distribution and the target
+ // probability can be outside the range of [lower, upper].
+ if (upperBound == Double.MAX_VALUE) {
+ if (complement) {
+ if (survivalProbability(upperBound) > q) {
+ return Double.POSITIVE_INFINITY;
+ }
+ } else if (cumulativeProbability(upperBound) < p) {
+ return Double.POSITIVE_INFINITY;
+ }
+ }
+ if (lowerBound == -Double.MAX_VALUE) {
+ if (complement) {
+ if (survivalProbability(lowerBound) < q) {
+ return Double.NEGATIVE_INFINITY;
+ }
+ } else if (cumulativeProbability(lowerBound) > p) {
+ return Double.NEGATIVE_INFINITY;
+ }
+ }
+
final DoubleUnaryOperator fun = complement ?
arg -> survivalProbability(arg) - q :
arg -> cumulativeProbability(arg) - p;
+ // Note the initial value is robust to overflow.
+ // Do not use 0.5 * (lowerBound + upperBound).
final double x = new BrentSolver(SOLVER_RELATIVE_ACCURACY,
SOLVER_ABSOLUTE_ACCURACY,
SOLVER_FUNCTION_VALUE_ACCURACY)
.findRoot(fun,
lowerBound,
- 0.5 * (lowerBound + upperBound),
+ lowerBound + 0.5 * (upperBound - lowerBound),
upperBound);
if (!isSupportConnected()) {
@@ -331,8 +355,10 @@ abstract class AbstractContinuousDistribution
* @return the infimum y
*/
private double searchPlateau(boolean complement, double lower, final double x) {
- /* Test for plateau. Lower the value x if the probability is the same. */
- final double dx = SOLVER_ABSOLUTE_ACCURACY;
+ // Test for plateau. Lower the value x if the probability is the same.
+ // Ensure the step is robust to the solver accuracy being less
+ // than 1 ulp of x (e.g. dx=0 will infinite loop)
+ final double dx = Math.max(SOLVER_ABSOLUTE_ACCURACY, Math.ulp(x));
if (x - dx >= lower) {
final DoubleUnaryOperator fun = complement ?
this::survivalProbability :
diff --git a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/AbstractContinuousDistributionTest.java b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/AbstractContinuousDistributionTest.java
index 317f95b..e72698a 100644
--- a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/AbstractContinuousDistributionTest.java
+++ b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/AbstractContinuousDistributionTest.java
@@ -371,4 +371,113 @@ class AbstractContinuousDistributionTest {
Assertions.assertEquals(10 - x, distribution.inverseSurvivalProbability(0.25), tolerance, "Inverse SF");
Assertions.assertEquals(Double.POSITIVE_INFINITY, distribution.inverseSurvivalProbability(0), "Inverse CDF");
}
+
+ /**
+ * Create a distribution near positive infinity so that it is truncated by MAX_VALUE.
+ */
+ @Test
+ void testTruncatedDistributionAtPositiveInfinity() {
+ final double mean = Double.MAX_VALUE;
+ final double width = Double.MAX_VALUE / 2;
+ final CentredUniformDistribution dist = new CentredUniformDistribution(mean, width);
+ final double x = mean - width / 2;
+ Assertions.assertEquals(0, dist.cumulativeProbability(x));
+ Assertions.assertTrue(dist.cumulativeProbability(Math.nextUp(x)) > 0);
+ Assertions.assertEquals(0.25, dist.cumulativeProbability(mean - width / 4));
+ Assertions.assertEquals(0.5, dist.cumulativeProbability(mean));
+
+ // Truncated
+ Assertions.assertEquals(1.0, dist.cumulativeProbability(Math.nextUp(mean)));
+
+ // Inversion should be robust to return the upper infinite bound
+ Assertions.assertEquals(mean, dist.inverseCumulativeProbability(0.5));
+ Assertions.assertEquals(mean, dist.inverseSurvivalProbability(0.5));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, dist.inverseCumulativeProbability(0.75));
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, dist.inverseSurvivalProbability(0.25));
+ }
+
+ /**
+ * Create a distribution near negative infinity so that it is truncated by -MAX_VALUE.
+ */
+ @Test
+ void testTruncatedDistributionAtNegativeInfinity() {
+ final double mean = -Double.MAX_VALUE;
+ final double width = Double.MAX_VALUE / 2;
+ final CentredUniformDistribution dist = new CentredUniformDistribution(mean, width);
+ final double x = mean + width / 2;
+ Assertions.assertEquals(1, dist.cumulativeProbability(x));
+ Assertions.assertTrue(dist.cumulativeProbability(Math.nextDown(x)) < 1);
+ Assertions.assertEquals(0.75, dist.cumulativeProbability(mean + width / 4));
+ Assertions.assertEquals(0.5, dist.cumulativeProbability(mean));
+
+ // Truncated
+ Assertions.assertEquals(0.0, dist.cumulativeProbability(Math.nextDown(mean)));
+
+ // Inversion should be robust to return the lower infinite bound
+ Assertions.assertEquals(mean, dist.inverseCumulativeProbability(0.5));
+ Assertions.assertEquals(mean, dist.inverseSurvivalProbability(0.5));
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, dist.inverseCumulativeProbability(0.25));
+ Assertions.assertEquals(Double.NEGATIVE_INFINITY, dist.inverseSurvivalProbability(0.75));
+ }
+
+ /**
+ * Uniform distribution described with a centre and width.
+ * This can be use to place the centre of the distribution at the limit of a finite double
+ * value so that the distribution is truncated.
+ */
+ static final class CentredUniformDistribution extends AbstractContinuousDistribution {
+ private final double mean;
+ private final double width;
+ private final double density;
+ private final double lo;
+ private final double hi;
+
+ /**
+ * @param mean Mean
+ * @param width Width
+ */
+ CentredUniformDistribution(double mean, double width) {
+ this.mean = mean;
+ this.width = width;
+ density = 1 / width;
+ lo = mean - width / 2;
+ hi = mean + width / 2;
+ }
+
+ @Override
+ public double density(double x) {
+ return density;
+ }
+
+ @Override
+ public double cumulativeProbability(double x) {
+ if (x <= lo) {
+ return 0;
+ }
+ if (x >= hi) {
+ return 1;
+ }
+ return 0.5 + density * (x - mean);
+ }
+
+ @Override
+ public double getMean() {
+ return mean;
+ }
+
+ @Override
+ public double getVariance() {
+ return width * width / 12;
+ }
+
+ @Override
+ public double getSupportLowerBound() {
+ return lo;
+ }
+
+ @Override
+ public double getSupportUpperBound() {
+ return hi;
+ }
+ }
}