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
The following commit(s) were added to refs/heads/master by this push:
new 2535c20a STATISTICS-102: Fix the Cauchy distribution quantile function
in the tails (#140)
2535c20a is described below
commit 2535c20ad85ce756c0d396151d83e7fc4e2c904c
Author: 陈志谦 <[email protected]>
AuthorDate: Fri Nov 6 17:36:48 2026 +0800
STATISTICS-102: Fix the Cauchy distribution quantile function in the tails
(#140)
CauchyDistribution.inverseCumulativeProbability(p) computes
location + scale * tan(pi * (p - 0.5)). For a probability close to 0 or 1
the
subtraction of one half cancels the significant digits of p, so the tangent
argument is rounded to +/- pi / 2 and the result saturates at about
+/- 1.633e16 * scale whatever the probability:
CauchyDistribution.of(0.5, 2).inverseCumulativeProbability(1e-20)
// -3.266247870639074e16, expected -6.366197723675814e19
CauchyDistribution.of(0.5, 2).inverseCumulativeProbability(1 - 1e-16)
// 3.957875932190438e15, expected 5.73416113922266e15
The loss of accuracy starts in the middle of the tail, not just at the
saturation: the value at p = 1e-5 is out by 5e-12 relative (the value at
p = 1e-10 by 5e-7). The function is not even monotonic (successive
probabilities below 1e-18 can return the same value) and the median is off
by
a rounding of pi (2.1e-16 relative for the scale 2, location 1.2 test).
Use the equivalent forms tan(pi (p - 1/2)) = -cot(pi p) and
cot(pi (1 - p)) for the lower and upper tail. The cotangent is computed with
1 / tan(x) or cos(x) / sin(x) to avoid the relative error of sin(x) for
small
x, and is exactly zero at the median.
Assisted-By: DeepSeek dsv41
Signed-off-by: simpleqt <[email protected]>
---
.../distribution/CauchyDistribution.java | 45 ++++++++++++-
.../distribution/CauchyDistributionTest.java | 74 ++++++++++++++++++++++
src/changes/changes.xml | 7 ++
3 files changed, 124 insertions(+), 2 deletions(-)
diff --git
a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/CauchyDistribution.java
b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/CauchyDistribution.java
index b412db41..30670ae1 100644
---
a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/CauchyDistribution.java
+++
b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/CauchyDistribution.java
@@ -34,6 +34,8 @@ import
org.apache.commons.rng.sampling.distribution.StableSampler;
* @see <a href="https://mathworld.wolfram.com/CauchyDistribution.html";>Cauchy
distribution (MathWorld)</a>
*/
public final class CauchyDistribution extends AbstractContinuousDistribution {
+ /** The value of one half. */
+ private static final double HALF = 0.5;
/** The location of this distribution. */
private final double location;
/** The scale of this distribution. */
@@ -136,7 +138,7 @@ public final class CauchyDistribution extends
AbstractContinuousDistribution {
} else if (p == 1) {
return Double.POSITIVE_INFINITY;
}
- return location + scale * Math.tan(Math.PI * (p - 0.5));
+ return location + scale * tailTangent(p);
}
/**
@@ -153,7 +155,46 @@ public final class CauchyDistribution extends
AbstractContinuousDistribution {
} else if (p == 0) {
return Double.POSITIVE_INFINITY;
}
- return location - scale * Math.tan(Math.PI * (p - 0.5));
+ return location - scale * tailTangent(p);
+ }
+
+ /**
+ * Compute the value of \( \tan(\pi (p - 1/2)) \) for a probability \( p
\) strictly within (0, 1).
+ *
+ * <p>The subtraction of one half from \( p \) loses the significant
digits of a probability
+ * close to 0 or 1, so the tails are computed using the equivalent forms
+ * \( -\cot(\pi p) \) and \( \cot(\pi (1 - p)) \) respectively.
+ *
+ * @param p Probability.
+ * @return tan(pi (p - 1/2))
+ */
+ private static double tailTangent(double p) {
+ // The tangent is exactly zero at the median; this also avoids the
rounding
+ // of pi (p - 1/2) to a non-zero value.
+ if (p == HALF) {
+ return 0;
+ }
+ // Avoid the cancellation of the significant digits of p in (p - 1/2).
+ // Assume p != 1/2; the two branches agree at the rounding of pi p.
+ return p < HALF ?
+ -cotPi(p) :
+ cotPi(1 - p);
+ }
+
+ /**
+ * Compute the cotangent of {@code pi x} for a value {@code x} in (0, 0.5].
+ *
+ * <p>For small \( x \), \( \sin(\pi x) \) has a relative error of
+ * \( (\pi x)^2 / 6 \) so the division by it is avoided.
+ *
+ * @param x Value.
+ * @return cot(pi x)
+ */
+ private static double cotPi(double x) {
+ final double px = Math.PI * x;
+ return px < HALF ?
+ 1 / Math.tan(px) :
+ Math.cos(px) / Math.sin(px);
}
/**
diff --git
a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/CauchyDistributionTest.java
b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/CauchyDistributionTest.java
index f2a894e6..833112cd 100644
---
a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/CauchyDistributionTest.java
+++
b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/CauchyDistributionTest.java
@@ -17,6 +17,9 @@
package org.apache.commons.statistics.distribution;
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.Test;
+
/**
* Test cases for {@link CauchyDistribution}.
* Extends {@link BaseContinuousDistributionTest}. See javadoc of that class
for details.
@@ -48,4 +51,75 @@ class CauchyDistributionTest extends
BaseContinuousDistributionTest {
protected double getRelativeTolerance() {
return 5e-15;
}
+
+ //-------------------- Additional test cases
-------------------------------
+
+ /**
+ * Test the quantile function in the tails of the distribution.
+ *
+ * <p>The standard computation \( x_0 + \gamma \tan(\pi (p - 1/2)) \)
loses every
+ * significant digit of a probability close to 0 or 1: the tangent
argument is
+ * rounded to \( \pm \pi / 2 \) and the result saturates at
+ * \( \pm 1.633 \times 10^{16} \gamma \). The expected values are computed
using
+ * {@code scipy.stats.cauchy(loc=x0, scale=gamma)}.
+ */
+ @Test
+ void testInverseProbabilityTails() {
+ // Position and scale with location 0.5 and scale 2
+ final CauchyDistribution dist = CauchyDistribution.of(0.5, 2);
+ final double[] p = {
+ 0.5, 1e-1, 1e-5, 1e-10, 1e-20, 1e-100, 1e-300,
+ 1 - 1e-7, 1 - 1e-12, 1 - 1e-16, 1 - 2 * 0x1.0p-53, 0x1.0p-53,
0x1.0p-60
+ };
+ final double[] icdf = {
+ 0.5, -5.655367074350508, -63661.47721581417, -6366197723.175813,
+ -6.366197723675814e19, -6.366197723675814e99,
-6.366197723675814e299,
+ 6366198.22702649, 636633855804.0593, 5734161139222660.0,
2867080569611330.0,
+ -5734161139222658.0, -7.339726258205004e17
+ };
+ for (int i = 0; i < p.length; i++) {
+ final double pi = p[i];
+ final double expected = icdf[i];
+ final double actual = dist.inverseCumulativeProbability(pi);
+ final double tolerance = Math.abs(expected) * 4e-16 + 1e-15;
+ Assertions.assertEquals(expected, actual, tolerance, () -> "icdf("
+ pi + ")");
+ // The survival function is the reflection of the cumulative value
in the location
+ Assertions.assertEquals(2 * dist.getMedian() - expected,
dist.inverseSurvivalProbability(pi),
+ tolerance, () -> "isf(" + pi + ")");
+ }
+ }
+
+ /**
+ * Test the quantile function in the tails is monotonic and antisymmetric.
+ *
+ * <p>The tail values can not be checked with a round trip through the
cumulative
+ * probability: \( 1/2 + \arctan(x) / \pi \) cancels by \( \pi / 2 \) so
the relative
+ * accuracy of {@code cdf(icdf(p))} is \( \mathrm{ulp}(\pi/2) / p \) and
is zero for
+ * the smallest probabilities.
+ */
+ @Test
+ void testInverseProbabilityTailsAreMonotonic() {
+ // The quantile function of the standard distribution is cot(pi p)
+ final CauchyDistribution dist = CauchyDistribution.of(0, 1);
+ double previous = dist.inverseCumulativeProbability(0.5);
+ Assertions.assertEquals(0, previous, "icdf(0.5)");
+ for (int i = -1; i >= -300; i--) {
+ final double p = Math.pow(10, i);
+ final double x = dist.inverseCumulativeProbability(p);
+ // Strictly decreasing towards the lower tail and finite
+ Assertions.assertTrue(x < previous, () -> "icdf(" + p + ")");
+ Assertions.assertTrue(x < 0, () -> "icdf(" + p + ")");
+ // The survival value is the reflection in the support centre
+ Assertions.assertEquals(-x, dist.inverseSurvivalProbability(p),
Math.abs(x) * 1e-15, () -> "isf(" + p + ")");
+ previous = x;
+ }
+ // And towards the upper tail
+ previous = dist.inverseCumulativeProbability(0.5);
+ for (int i = 1; i <= 16; i++) {
+ final double p = 1 - Math.pow(10, -i);
+ final double x = dist.inverseCumulativeProbability(p);
+ Assertions.assertTrue(x > 0 && x > previous, () -> "icdf(" + p +
")");
+ previous = x;
+ }
+ }
}
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 35284f03..24645a30 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -106,6 +106,13 @@ If the output is not quite correct, check for invisible
trailing spaces!
<action dev="aherbert" type="fix" due-to="Security scan, Alex Herbert"
issue="STATISTICS-95">
"NormalConfidenceInterval": Validate mean is finite and variance is
positive finite.
</action>
+ <action dev="aherbert" type="fix" due-to="simpleqt"
issue="STATISTICS-102">
+ "CauchyDistribution": Correct the inverse cumulative probability
function.
+ Subtraction of one half from a probability close to 0 or 1 discarded
the
+ significant digits, so the result saturated in the tails and was not
+ monotonic. The tails are evaluated using the equivalent cotangent forms
+ of the tangent.
+ </action>
</release>
<release version="1.3" date="2026-05-01" description="Adds new features
(requires Java 8).">
