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 481e922 STATISTICS-25: Specialise t-distribution for infinite degrees
freedom
481e922 is described below
commit 481e922f44e72c4b5c4d4d991326077ad87a51d7
Author: Alex Herbert <[email protected]>
AuthorDate: Fri Oct 22 20:02:53 2021 +0100
STATISTICS-25: Specialise t-distribution for infinite degrees freedom
If the variance of the distribution matches the standard normal
distribution then delegate to a standard normal distribution.
---
.../statistics/distribution/TDistribution.java | 327 +++++++++++++--------
.../statistics/distribution/test.t.3.properties | 35 +++
.../statistics/distribution/test.t.4.properties | 36 +++
3 files changed, 279 insertions(+), 119 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 bda9b40..11dce96 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,51 +17,221 @@
package org.apache.commons.statistics.distribution;
import org.apache.commons.numbers.gamma.RegularizedBeta;
-import org.apache.commons.numbers.gamma.Erf;
+import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.numbers.gamma.LogGamma;
/**
* Implementation of <a
href='http://en.wikipedia.org/wiki/Student's_t-distribution'>Student's
t-distribution</a>.
*/
-public final 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;
+public abstract class TDistribution extends AbstractContinuousDistribution {
+ /** A standard normal distribution used for calculations.
+ * This is immutable and thread-safe and can be used across instances. */
+ static final NormalDistribution STANDARD_NORMAL = NormalDistribution.of(0,
1);
/** The degrees of freedom. */
private final double degreesOfFreedom;
- /** degreesOfFreedom / 2. */
- private final double dofOver2;
- /** Cached value. */
- private final double factor;
- /** Cached value. */
- private final double mean;
- /** Cached value. */
- private final double variance;
+
+ /**
+ * Specialisation of the T-distribution used when there are infinite
degrees of freedom.
+ * In this case the distribution matches a normal distribution. This is
used when the
+ * variance is not different from 1.0.
+ *
+ * <p>This delegates all methods to the standard normal distribution.
Instances are
+ * allowed to provide access to the degrees of freedom used during
construction.
+ */
+ private static class NormalTDistribution extends TDistribution {
+ /**
+ * @param degreesOfFreedom Degrees of freedom.
+ */
+ NormalTDistribution(double degreesOfFreedom) {
+ super(degreesOfFreedom);
+ }
+
+ @Override
+ public double density(double x) {
+ return STANDARD_NORMAL.density(x);
+ }
+
+ @Override
+ public double probability(double x0, double x1) {
+ return STANDARD_NORMAL.probability(x0, x1);
+ }
+
+ @Override
+ public double logDensity(double x) {
+ return STANDARD_NORMAL.logDensity(x);
+ }
+
+ @Override
+ public double cumulativeProbability(double x) {
+ return STANDARD_NORMAL.cumulativeProbability(x);
+ }
+
+ @Override
+ public double inverseCumulativeProbability(double p) {
+ return STANDARD_NORMAL.inverseCumulativeProbability(p);
+ }
+
+ @Override
+ public double getMean() {
+ return 0;
+ }
+
+ @Override
+ public double getVariance() {
+ return 1.0;
+ }
+
+ @Override
+ public Sampler createSampler(UniformRandomProvider rng) {
+ return STANDARD_NORMAL.createSampler(rng);
+ }
+ }
+
+ /**
+ * Implementation of Student's T-distribution.
+ */
+ private static class StudentsTDistribution extends TDistribution {
+ /** 2. */
+ private static final double TWO = 2;
+ /** Number of degrees of freedom above which to use the normal
distribution.
+ * This is used to check the CDF when the degrees of freedom is large.
*/
+ private static final double DOF_THRESHOLD_NORMAL = 2.99e6;
+
+ /** degreesOfFreedom / 2. */
+ private final double dofOver2;
+ /** Cached value. */
+ private final double factor;
+ /** Cached value. */
+ private final double mean;
+ /** Cached value. */
+ private final double variance;
+
+ /**
+ * @param degreesOfFreedom Degrees of freedom.
+ * @param variance Precomputed variance
+ */
+ StudentsTDistribution(double degreesOfFreedom, double variance) {
+ super(degreesOfFreedom);
+
+ dofOver2 = 0.5 * degreesOfFreedom;
+ factor = LogGamma.value(dofOver2 + 0.5) -
+ 0.5 * (Math.log(Math.PI) + Math.log(degreesOfFreedom)) -
+ LogGamma.value(dofOver2);
+ this.variance = variance;
+ mean = degreesOfFreedom > 1 ? 0 : Double.NaN;
+ }
+
+ /**
+ * @param degreesOfFreedom Degrees of freedom.
+ * @return the variance
+ */
+ static double computeVariance(double degreesOfFreedom) {
+ if (degreesOfFreedom == Double.POSITIVE_INFINITY) {
+ return 1;
+ } else if (degreesOfFreedom > TWO) {
+ return degreesOfFreedom / (degreesOfFreedom - 2);
+ } else if (degreesOfFreedom > 1) {
+ return Double.POSITIVE_INFINITY;
+ } else {
+ return Double.NaN;
+ }
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double density(double x) {
+ return Math.exp(logDensity(x));
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double logDensity(double x) {
+ final double nPlus1Over2 = dofOver2 + 0.5;
+ return factor - nPlus1Over2 * Math.log1p(x * x /
getDegreesOfFreedom());
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double cumulativeProbability(double x) {
+ if (x == 0) {
+ return 0.5;
+ }
+ final double df = getDegreesOfFreedom();
+ if (df > DOF_THRESHOLD_NORMAL) {
+ return STANDARD_NORMAL.cumulativeProbability(x);
+ }
+ final double x2 = x * x;
+ // z = 1 / (1 + x^2/df)
+ // Simplify by multiplication by df
+ final double z = df / (df + x2);
+
+ // The RegularizedBeta has the complement:
+ // I(z, a, b) = 1 - I(1 - z, a, b)
+ // This is used when z > (a + 1) / (2 + b + a).
+ // Detect this condition and directly use the complement.
+ if (z > (dofOver2 + 1) / (2.5 + dofOver2)) {
+ // zc = 1 - z; pc = 1 - p
+ final double zc = x2 / (df + x2);
+ final double pc = RegularizedBeta.value(zc, 0.5, dofOver2);
+
+ return x < 0 ?
+ // 0.5 * p == 0.5 * (1 - pc) = 0.5 - 0.5 * pc
+ 0.5 - 0.5 * pc :
+ // 1 - 0.5 * p == 1 - 0.5 * (1 - pc) = 0.5 + 0.5 * pc
+ 0.5 + 0.5 * pc;
+ }
+
+ final double p = RegularizedBeta.value(z, dofOver2, 0.5);
+
+ return x < 0 ?
+ 0.5 * p :
+ 1 - 0.5 * p;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * <p>For degrees of freedom parameter {@code df}, the mean is
+ * <ul>
+ * <li>zero if {@code df > 1}, and</li>
+ * <li>undefined ({@code Double.NaN}) otherwise.</li>
+ * </ul>
+ */
+ @Override
+ public double getMean() {
+ return mean;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * <p>For degrees of freedom parameter {@code df}, the variance is
+ * <ul>
+ * <li>{@code df / (df - 2)} if {@code df > 2},</li>
+ * <li>infinite ({@code Double.POSITIVE_INFINITY}) if {@code 1 < df
<= 2}, and</li>
+ * <li>undefined ({@code Double.NaN}) otherwise.</li>
+ * </ul>
+ */
+ @Override
+ public double getVariance() {
+ return variance;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ protected double getMedian() {
+ // Overridden for the probability(double, double) method.
+ // This is intentionally not a public method.
+ return 0;
+ }
+ }
/**
* @param degreesOfFreedom Degrees of freedom.
*/
private TDistribution(double degreesOfFreedom) {
this.degreesOfFreedom = degreesOfFreedom;
-
- dofOver2 = 0.5 * degreesOfFreedom;
- factor = LogGamma.value(dofOver2 + 0.5) -
- 0.5 * (Math.log(Math.PI) + Math.log(degreesOfFreedom)) -
- LogGamma.value(dofOver2);
- if (degreesOfFreedom > TWO) {
- mean = 0;
- variance = degreesOfFreedom / (degreesOfFreedom - 2);
- } else if (degreesOfFreedom > 1) {
- mean = 0;
- variance = Double.POSITIVE_INFINITY;
- } else {
- mean = Double.NaN;
- variance = Double.NaN;
- }
}
/**
@@ -76,7 +246,13 @@ public final class TDistribution extends
AbstractContinuousDistribution {
throw new
DistributionException(DistributionException.NOT_STRICTLY_POSITIVE,
degreesOfFreedom);
}
- return new TDistribution(degreesOfFreedom);
+ // If the variance converges to 1 use a NormalDistribution.
+ // Occurs at 2^55 = 3.60e16
+ final double var =
StudentsTDistribution.computeVariance(degreesOfFreedom);
+ if (var == 1) {
+ return new NormalTDistribution(degreesOfFreedom);
+ }
+ return new StudentsTDistribution(degreesOfFreedom, var);
}
/**
@@ -90,56 +266,6 @@ public final class TDistribution extends
AbstractContinuousDistribution {
/** {@inheritDoc} */
@Override
- public double density(double x) {
- return Math.exp(logDensity(x));
- }
-
- /** {@inheritDoc} */
- @Override
- public double logDensity(double x) {
- final double nPlus1Over2 = dofOver2 + 0.5;
- return factor - nPlus1Over2 * Math.log1p(x * x / degreesOfFreedom);
- }
-
- /** {@inheritDoc} */
- @Override
- public double cumulativeProbability(double x) {
- if (x == 0) {
- return 0.5;
- }
- if (degreesOfFreedom > DOF_THRESHOLD_NORMAL) {
- return 0.5 * (1 + Erf.value(x * ONE_OVER_SQRT_TWO));
- }
- final double x2 = x * x;
- // z = 1 / (1 + x^2/df)
- // Simplify by multiplication by df
- final double z = degreesOfFreedom / (degreesOfFreedom + x2);
-
- // The RegularizedBeta has the complement:
- // I(z, a, b) = 1 - I(1 - z, a, b)
- // This is used when z > (a + 1) / (2 + b + a).
- // Detect this condition and directly use the complement.
- if (z > (dofOver2 + 1) / (2.5 + dofOver2)) {
- // zc = 1 - z; pc = 1 - p
- final double zc = x2 / (degreesOfFreedom + x2);
- final double pc = RegularizedBeta.value(zc, 0.5, dofOver2);
-
- return x < 0 ?
- // 0.5 * p == 0.5 * (1 - pc) = 0.5 - 0.5 * pc
- 0.5 - 0.5 * pc :
- // 1 - 0.5 * p == 1 - 0.5 * (1 - pc) = 0.5 + 0.5 * pc
- 0.5 + 0.5 * pc;
- }
-
- final double p = RegularizedBeta.value(z, dofOver2, 0.5);
-
- return x < 0 ?
- 0.5 * p :
- 1 - 0.5 * p;
- }
-
- /** {@inheritDoc} */
- @Override
public double survivalProbability(double x) {
// Exploit symmetry
return cumulativeProbability(-x);
@@ -148,35 +274,6 @@ public final class TDistribution extends
AbstractContinuousDistribution {
/**
* {@inheritDoc}
*
- * <p>For degrees of freedom parameter {@code df}, the mean is
- * <ul>
- * <li>zero if {@code df > 1}, and</li>
- * <li>undefined ({@code Double.NaN}) otherwise.</li>
- * </ul>
- */
- @Override
- public double getMean() {
- return mean;
- }
-
- /**
- * {@inheritDoc}
- *
- * <p>For degrees of freedom parameter {@code df}, the variance is
- * <ul>
- * <li>{@code df / (df - 2)} if {@code df > 2},</li>
- * <li>infinite ({@code Double.POSITIVE_INFINITY}) if {@code 1 < df <=
2}, and</li>
- * <li>undefined ({@code Double.NaN}) otherwise.</li>
- * </ul>
- */
- @Override
- public double getVariance() {
- return variance;
- }
-
- /**
- * {@inheritDoc}
- *
* <p>The lower bound of the support is always negative infinity..
*
* @return lower bound of the support (always
@@ -200,14 +297,6 @@ public final class TDistribution extends
AbstractContinuousDistribution {
return Double.POSITIVE_INFINITY;
}
- /** {@inheritDoc} */
- @Override
- protected double getMedian() {
- // Overridden for the probability(double, double) method.
- // This is intentionally not a public method.
- return 0;
- }
-
/**
* {@inheritDoc}
*
diff --git
a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.t.3.properties
b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.t.3.properties
new file mode 100644
index 0000000..e40171e
--- /dev/null
+++
b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.t.3.properties
@@ -0,0 +1,35 @@
+# Licensed to the Apache Software Foundation (ASF) under one or more
+# contributor license agreements. See the NOTICE file distributed with
+# this work for additional information regarding copyright ownership.
+# The ASF licenses this file to You under the Apache License, Version 2.0
+# (the "License"); you may not use this file except in compliance with
+# the License. You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+parameters = 2.0
+# Limited by cdf inverse mapping
+tolerance.relative = 1e-11
+# Computed using scipy stats
+mean = 0
+variance = Infinity
+cdf.points = \
+ -2, -1, 0, 1, 2, 3, 4, 5
+cdf.values = \
+ 0.09175170953613696, 0.21132486540518713, 0.5 ,\
+ 0.78867513459481287, 0.90824829046386302, 0.95226701686664539,\
+ 0.97140452079103168, 0.98112522432468807
+pdf.values = \
+ 0.06804138174397717, 0.19245008972987523, 0.35355339059327379,\
+ 0.19245008972987523, 0.06804138174397717, 0.02741012223434215,\
+ 0.0130945700219731 , 0.00712778110110649
+sf.values = \
+0.90824829046386302, 0.78867513459481287, 0.5 ,\
+ 0.21132486540518713, 0.09175170953613696, 0.04773298313335456,\
+ 0.02859547920896831, 0.01887477567531186
diff --git
a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.t.4.properties
b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.t.4.properties
new file mode 100644
index 0000000..a44fda5
--- /dev/null
+++
b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.t.4.properties
@@ -0,0 +1,36 @@
+# Licensed to the Apache Software Foundation (ASF) under one or more
+# contributor license agreements. See the NOTICE file distributed with
+# this work for additional information regarding copyright ownership.
+# The ASF licenses this file to You under the Apache License, Version 2.0
+# (the "License"); you may not use this file except in compliance with
+# the License. You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# Special case for large degrees of freedom is the standard normal distribution
+parameters = Infinity
+mean = 0
+variance = 1
+# normal(0, 1) computed using R
+cdf.points = \
+ -2, -1, 0, 1, 2, 3, 4, 5
+cdf.values = \
+ 0.022750131948179212055 0.158655253931457046468 0.500000000000000000000 \
+ 0.841344746068542925777 0.977249868051820791415 0.998650101968369896532 \
+ 0.999968328758166880021 0.999999713348428076465
+pdf.values = \
+ 5.3990966513188062836e-02 2.4197072451914336533e-01 \
+ 3.9894228040143270286e-01 2.4197072451914336533e-01 \
+ 5.3990966513188062836e-02 4.4318484119380075273e-03 \
+ 1.3383022576488536764e-04 1.4867195147342976779e-06
+sf.values = \
+ 9.7724986805182079141e-01 8.4134474606854292578e-01 \
+ 5.0000000000000000000e-01 1.5865525393145704647e-01 \
+ 2.2750131948179212055e-02 1.3498980316300945772e-03 \
+ 3.1671241833119924327e-05 2.8665157187919391185e-07
