This is an automated email from the ASF dual-hosted git repository.
erans 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 5b79ca6 STATISTICS-25: Use normal distribution for large number of
degrees of freedom.
5b79ca6 is described below
commit 5b79ca6b24f05eb4b45a48531129abba4d0b683a
Author: Gilles Sadowski <[email protected]>
AuthorDate: Mon Oct 26 01:11:50 2020 +0100
STATISTICS-25: Use normal distribution for large number of degrees of
freedom.
Current workaround is only accurate to about 1e-9.
---
.../statistics/distribution/TDistribution.java | 22 ++++++++++-----
.../statistics/distribution/TDistributionTest.java | 33 ++++++++++++++++++++++
2 files changed, 48 insertions(+), 7 deletions(-)
diff --git
a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/TDistribution.java
b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/TDistribution.java
index 5f74ff5..6e69b33 100644
---
a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/TDistribution.java
+++
b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/TDistribution.java
@@ -17,6 +17,7 @@
package org.apache.commons.statistics.distribution;
import org.apache.commons.numbers.gamma.RegularizedBeta;
+import org.apache.commons.numbers.gamma.Erf;
import org.apache.commons.numbers.gamma.LogGamma;
/**
@@ -25,6 +26,10 @@ import org.apache.commons.numbers.gamma.LogGamma;
public class TDistribution extends AbstractContinuousDistribution {
/** 2. */
private static final double TWO = 2;
+ /** 1 / sqrt(2). */
+ private static final double ONE_OVER_SQRT_TWO = 1 / Math.sqrt(2);
+ /** Number of degrees of freedom above which to use the normal
distribution. */
+ private static final double DOF_THRESHOLD_NORMAL = 2.99e6;
/** The degrees of freedom. */
private final double degreesOfFreedom;
@@ -94,13 +99,16 @@ public class TDistribution extends
AbstractContinuousDistribution {
if (x == 0) {
return 0.5;
} else {
- final double t =
- RegularizedBeta.value(degreesOfFreedom / (degreesOfFreedom +
(x * x)),
- dofOver2,
- 0.5);
- return x < 0 ?
- 0.5 * t :
- 1 - 0.5 * t;
+ if (degreesOfFreedom > DOF_THRESHOLD_NORMAL) {
+ return 0.5 * (1 + Erf.value(x * ONE_OVER_SQRT_TWO));
+ } else {
+ final double a = 1 / (1 + x * x / degreesOfFreedom);
+ final double t = RegularizedBeta.value(a, dofOver2, 0.5);
+
+ return x < 0 ?
+ 0.5 * t :
+ 1 - 0.5 * t;
+ }
}
}
diff --git
a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/TDistributionTest.java
b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/TDistributionTest.java
index 36ec1a8..a48ff7d 100644
---
a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/TDistributionTest.java
+++
b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/TDistributionTest.java
@@ -164,6 +164,39 @@ class TDistributionTest extends
ContinuousDistributionAbstractTest {
return;
}
+ // See https://issues.apache.org/jira/browse/STATISTICS-25
+ @Test
+ public void testStatistics25() {
+ final double[] df = {1, 10, 1e2, 1e3, 1e4, 1e5,
+ 1e6, 2e6, 2.98e6, 2.99e6, 3e6, 4e6,
+ 1e7, 1e8, 1e9, 1e10};
+ // Generated from Python.
+ final double[] expected = {0.507956089912,
+ 0.509726595102,
+ 0.509947608093,
+ 0.509970024339,
+ 0.509972268782,
+ 0.509972493254,
+ 0.509972515701,
+ 0.509972516948,
+ 0.509972517358,
+ 0.509972517361,
+ 0.509972517364,
+ 0.509972517572,
+ 0.509972517946,
+ 0.50997251817,
+ 0.509972518193,
+ 0.509972518195};
+
+ final double c = 0.025;
+ final double tol = 1e-9;
+ for (int i = 0; i < df.length; i++) {
+ final TDistribution dist = new TDistribution(df[i]);
+ final double x = dist.cumulativeProbability(c);
+ Assertions.assertEquals(expected[i], x, tol);
+ }
+ }
+
private double[] makeNistResults(double[] args, int df) {
final TDistribution td = new TDistribution(df);
final double[] res = new double[args.length];
