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 2022/01/26 18:10:53 UTC
[commons-statistics] 02/02: Update F distribution with the beta function
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 09cd1168cf0cdedaa4aa6932a59e088aff169fbd
Author: aherbert <ah...@apache.org>
AuthorDate: Wed Jan 26 18:05:11 2022 +0000
Update F distribution with the beta function
Added second PDF form to the class javadoc.
Use the beta function following the Boost documentation for the F
distribution.
The CDF and SF can use the regularized beta or its complement. This is
chosen depending on the x argument.
The PDF can use the beta derivative using the same logic for choosing
the value passed to the derivative function.
Updated the log PDF to reduce the number of terms.
Added SF values to test files.
Large degrees of freedom data (n=100, m=100) requires a more accurate
source implementation.
---
.../statistics/distribution/FDistribution.java | 91 ++++++++----
.../statistics/distribution/FDistributionTest.java | 2 +-
.../statistics/distribution/test.f.1.properties | 19 ++-
.../statistics/distribution/test.f.2.properties | 43 ++++--
.../statistics/distribution/test.f.3.properties | 43 ++++--
.../statistics/distribution/test.f.4.properties | 43 ++++--
.../statistics/distribution/test.f.5.properties | 15 ++
.../statistics/distribution/test.f.6.properties | 155 +++++++++++++++++----
8 files changed, 308 insertions(+), 103 deletions(-)
diff --git a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/FDistribution.java b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/FDistribution.java
index e6d3bc9..1a9ca94 100644
--- a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/FDistribution.java
+++ b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/FDistribution.java
@@ -25,9 +25,11 @@ import org.apache.commons.numbers.gamma.RegularizedBeta;
*
* <p>The probability density function of \( X \) is:
*
- * <p>\[ f(x; d_1, d_2) = \frac{1}{\operatorname{B}\left(\frac{d_1}{2},\frac{d_2}{2}\right)} \left(\frac{d_1}{d_2}\right)^{d_1/2} x^{d_1/2 - 1} \left(1+\frac{d_1}{d_2} \, x \right)^{-(d_1+d_2)/2} \]
+ * <p>\[ \begin{aligned}
+ * f(x; n, m) &= \frac{1}{\operatorname{B}\left(\frac{n}{2},\frac{m}{2}\right)} \left(\frac{n}{m}\right)^{n/2} x^{n/2 - 1} \left(1+\frac{n}{m} \, x \right)^{-(n+m)/2} \\
+ * &= \frac{n^{\frac n 2} m^{\frac m 2} x^{\frac{n}{2}-1} }{ (nx+m)^{\frac{(n+m)}{2}} \operatorname{B}\left(\frac{n}{2},\frac{m}{2}\right)} \end{aligned} \]
*
- * <p>for \( d_1, d_2 > 0 \) the degrees of freedom, \( \operatorname{B}(d_1 / 2, d_2 / 2) \) is the beta function,
+ * <p>for \( n, m > 0 \) the degrees of freedom, \( \operatorname{B}(a, b) \) is the beta function,
* and \( x \in [0, \infty) \).
*
* @see <a href="https://en.wikipedia.org/wiki/F-distribution">F-distribution (Wikipedia)</a>
@@ -47,12 +49,10 @@ public final class FDistribution extends AbstractContinuousDistribution {
private final double numeratorDegreesOfFreedom;
/** The denominator degrees of freedom. */
private final double denominatorDegreesOfFreedom;
- /** n/2 * log(n) with n = numerator DF. */
- private final double nhalfLogn;
- /** m/2 * log(m) with m = denominator DF. */
- private final double mhalfLogm;
+ /** n/2 * log(n) + m/2 * log(m) with n = numerator DF and m = denominator DF. */
+ private final double nHalfLogNmHalfLogM;
/** LogBeta(n/2, n/2) with n = numerator DF. */
- private final double logBetaNhalfNhalf;
+ private final double logBetaNhalfMhalf;
/** Cached value for inverse probability function. */
private final double mean;
/** Cached value for inverse probability function. */
@@ -68,9 +68,9 @@ public final class FDistribution extends AbstractContinuousDistribution {
this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom;
final double nhalf = numeratorDegreesOfFreedom / 2;
final double mhalf = denominatorDegreesOfFreedom / 2;
- nhalfLogn = nhalf * Math.log(numeratorDegreesOfFreedom);
- mhalfLogm = mhalf * Math.log(denominatorDegreesOfFreedom);
- logBetaNhalfNhalf = LogBeta.value(nhalf, mhalf);
+ nHalfLogNmHalfLogM = nhalf * Math.log(numeratorDegreesOfFreedom) +
+ mhalf * Math.log(denominatorDegreesOfFreedom);
+ logBetaNhalfMhalf = LogBeta.value(nhalf, mhalf);
if (denominatorDegreesOfFreedom > MIN_DENOMINATOR_DF_FOR_MEAN) {
mean = denominatorDegreesOfFreedom / (denominatorDegreesOfFreedom - 2);
@@ -141,7 +141,32 @@ public final class FDistribution extends AbstractContinuousDistribution {
*/
@Override
public double density(double x) {
- return Math.exp(logDensity(x));
+ if (x <= SUPPORT_LO ||
+ x >= SUPPORT_HI) {
+ // Special case x=0:
+ // PDF reduces to the term x^(df1 / 2 - 1) multiplied by a scaling factor
+ if (x == SUPPORT_LO && numeratorDegreesOfFreedom <= 2) {
+ return numeratorDegreesOfFreedom == 2 ?
+ 1 :
+ Double.POSITIVE_INFINITY;
+ }
+ return 0;
+ }
+
+ // Keep the z argument to the regularized beta well away from 1 to avoid rounding error.
+ // See: https://www.boost.org/doc/libs/1_78_0/libs/math/doc/html/math_toolkit/dist_ref/dists/f_dist.html
+
+ final double n = numeratorDegreesOfFreedom;
+ final double m = denominatorDegreesOfFreedom;
+ final double nx = n * x;
+ final double z = m + nx;
+ final double y = n * m / (z * z);
+ if (nx > m) {
+ return y * RegularizedBeta.derivative(m / z,
+ m / 2, n / 2);
+ }
+ return y * RegularizedBeta.derivative(nx / z,
+ n / 2, m / 2);
}
/** {@inheritDoc}
@@ -168,14 +193,12 @@ public final class FDistribution extends AbstractContinuousDistribution {
return Double.NEGATIVE_INFINITY;
}
- final double nhalf = numeratorDegreesOfFreedom / 2;
- final double mhalf = denominatorDegreesOfFreedom / 2;
- final double logx = Math.log(x);
- final double lognxm = Math.log(numeratorDegreesOfFreedom * x +
- denominatorDegreesOfFreedom);
- return nhalfLogn + nhalf * logx - logx +
- mhalfLogm - nhalf * lognxm - mhalf * lognxm -
- logBetaNhalfNhalf;
+ final double n = numeratorDegreesOfFreedom;
+ final double m = denominatorDegreesOfFreedom;
+
+ // This may suffer cancellation effects due to summation of opposing terms.
+ return nHalfLogNmHalfLogM + (n / 2 - 1) * Math.log(x) - logBetaNhalfMhalf -
+ ((n + m) / 2) * Math.log(n * x + m);
}
/**
@@ -200,9 +223,18 @@ public final class FDistribution extends AbstractContinuousDistribution {
final double n = numeratorDegreesOfFreedom;
final double m = denominatorDegreesOfFreedom;
- return RegularizedBeta.value((n * x) / (m + n * x),
- 0.5 * n,
- 0.5 * m);
+ // Keep the z argument to the regularized beta well away from 1 to avoid rounding error.
+ // See: https://www.boost.org/doc/libs/1_78_0/libs/math/doc/html/math_toolkit/dist_ref/dists/f_dist.html
+ // Note the logic in the Boost documentation for pdf and cdf is contradictory.
+ // This order will keep the argument far from 1.
+
+ final double nx = n * x;
+ if (nx > m) {
+ return RegularizedBeta.complement(m / (m + nx),
+ m / 2, n / 2);
+ }
+ return RegularizedBeta.value(nx / (m + nx),
+ n / 2, m / 2);
}
/** {@inheritDoc} */
@@ -217,12 +249,15 @@ public final class FDistribution extends AbstractContinuousDistribution {
final double n = numeratorDegreesOfFreedom;
final double m = denominatorDegreesOfFreedom;
- // Compute the complement of the regularized beta function
- // with direct computation of 1 - z:
- // 1 - I(z, a, b) = I(1 - z, b, a)
- return RegularizedBeta.value(m / (m + n * x),
- 0.5 * m,
- 0.5 * n);
+ // Do the opposite of the CDF
+
+ final double nx = n * x;
+ if (nx > m) {
+ return RegularizedBeta.value(m / (m + nx),
+ m / 2, n / 2);
+ }
+ return RegularizedBeta.complement(nx / (m + nx),
+ n / 2, m / 2);
}
/**
diff --git a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/FDistributionTest.java b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/FDistributionTest.java
index 44c57d8..f095114 100644
--- a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/FDistributionTest.java
+++ b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/FDistributionTest.java
@@ -48,7 +48,7 @@ class FDistributionTest extends BaseContinuousDistributionTest {
@Override
protected double getRelativeTolerance() {
- return 1e-14;
+ return 8e-15;
}
//-------------------- Additional test cases -------------------------------
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.1.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.1.properties
index 7a93d9c..40a69e2 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.1.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.1.properties
@@ -24,13 +24,20 @@ cdf.points = \
20.802663959456705101 8.745895256019915465 5.987565126046930253 \
4.387374187406129167 3.107511666638931302 Inf
cdf.values = \
- 0, 0.001, 0.01, 0.025, 0.05, 0.1, 0.999,\
- 0.990, 0.975, 0.950, 0.900, 1
+ 0.0000000000000000000000 0.0010000000000000024061 0.0100000000000000088818 \
+ 0.0249999999999999840405 0.0500000000000000860423 0.1000000000000000749401 \
+ 0.9989999999999999991118 0.9899999999999999911182 0.9749999999999999777955 \
+ 0.9499999999999999555911 0.9000000000000000222045 1.0000000000000000000000
pdf.values = \
- 0.00000000000000000000 0.06891565767071070048 0.23673565319275807761 \
- 0.36407413194146687196 0.48157078964913180297 0.59588047999337290239 \
- 0.00013344391565753327 0.00286681303402549050 0.00969192007503559549 \
- 0.02428838614718745875 0.06054913146582809741 0.00000000000000000000
+ 0.0000000000000000000000 0.0689156576707107004776 0.2367356531927580776120 \
+ 0.3640741319414668719645 0.4815707896491318029675 0.5958804799933729023920 \
+ 0.0001334439156575332650 0.0028668130340254904982 0.0096919200750355954943 \
+ 0.0242883861471874587523 0.0605491314658280974093 0.0000000000000000000000
+sf.values = \
+ 1.000000000000000000000 0.998999999999999999112 0.989999999999999991118 \
+ 0.975000000000000088818 0.949999999999999955591 0.899999999999999911182 \
+ 0.001000000000000001105 0.010000000000000010617 0.025000000000000001388 \
+ 0.050000000000000044409 0.100000000000000005551 0.000000000000000000000
# Computed using WolframAlpha
cdf.hp.points = 1e-7, 4e-8, 9e-8
cdf.hp.values = 1.578691625481747e-17, 1.597523916857153e-18, 1.2131195257872846e-17
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.2.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.2.properties
index 2aaaff9..d2ae231 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.2.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.2.properties
@@ -25,20 +25,20 @@ cdf.points = \
3.0, 3.25, 3.5, 3.75, \
4.0, 4.25, 4.5, 4.75, 5.0
cdf.values = \
- 0. , 0.07053456158585982849,\
- 0.0995037190209989153 , 0.12156613477096617215,\
- 0.14002800840280096861, 0.15617376188860607189,\
- 0.21821789023599236224, 0.33333333333333331483,\
- 0.4472135954999579277 , 0.52223296786709350048,\
- 0.57735026918962573106, 0.62017367294604230921,\
- 0.65465367070797708671, 0.68313005106397317601,\
- 0.70710678118654757274, 0.72760687510899890729,\
- 0.74535599249992989801, 0.76088591025268215162,\
- 0.77459666924148340428, 0.78679579246944320037,\
- 0.79772403521746559907, 0.80757285308724824358,\
- 0.81649658092772603446, 0.82462112512353213933,\
- 0.83205029433784372106, 0.83887049280786107897,\
- 0.84515425472851657407
+ 0. , 0.07053456158585984237,\
+ 0.09950371902099894306, 0.12156613477096619991,\
+ 0.14002800840280099637, 0.15617376188860609965,\
+ 0.21821789023599241775, 0.33333333333333337034,\
+ 0.44721359549995803873, 0.5222329678670936115 ,\
+ 0.57735026918962584208, 0.62017367294604242023,\
+ 0.65465367070797719773, 0.68313005106397328703,\
+ 0.70710678118654768376, 0.72760687510899901831,\
+ 0.74535599249993000903, 0.76088591025268237367,\
+ 0.77459666924148362632, 0.78679579246944342241,\
+ 0.79772403521746582111, 0.80757285308724846562,\
+ 0.8164965809277262565 , 0.82462112512353236138,\
+ 0.8320502943378439431 , 0.83887049280786130101,\
+ 0.84515425472851679611
pdf.values = \
inf, 3.50918216845073827059,\
2.46296334210393341735, 1.9961598484559310851 ,\
@@ -54,3 +54,18 @@ pdf.values = \
0.03402069087198857783, 0.03104456000465061899,\
0.0284461639089860982 , 0.02616360211486505952,\
0.02414726442081477353
+sf.values = \
+ 1. , 0.92946543841414108744,\
+ 0.90049628097900102919, 0.87843386522903432745,\
+ 0.8599719915971988371 , 0.84382623811139423342,\
+ 0.78178210976400741572, 0.66666666666666607455,\
+ 0.55278640450004201679, 0.47776703213290649952,\
+ 0.42264973081037426894, 0.37982632705395780182,\
+ 0.34534632929202285778, 0.31686994893602682399,\
+ 0.29289321881345248277, 0.27239312489100114822,\
+ 0.25464400750007010199, 0.23911408974731787613,\
+ 0.22540333075851665123, 0.2132042075305568829 ,\
+ 0.20227596478253440093, 0.19242714691275175642,\
+ 0.1835034190722739933 , 0.17537887487646791618,\
+ 0.16794970566215636221, 0.16112950719213892103,\
+ 0.15484574527148342593
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.3.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.3.properties
index 1f7fa43..26bebc8 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.3.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.3.properties
@@ -25,20 +25,20 @@ cdf.points = \
3.0, 3.25, 3.5, 3.75, \
4.0, 4.25, 4.5, 4.75, 5.0
cdf.values = \
- 0. , 0.0098524570233256923 ,\
- 0.01941932430907984328, 0.02871413764273582012,\
- 0.03774955135062372374, 0.04653741075440769037,\
- 0.08712907082472313991, 0.18350341907227396554,\
- 0.29289321881345242726, 0.36754446796632406214,\
- 0.42264973081037415792, 0.46547751617515120692,\
- 0.5 , 0.52859547920896821083,\
- 0.55278640450004190576, 0.57359856728877900434,\
- 0.59175170953613676073, 0.60776772972363179992,\
- 0.62203552699077246935, 0.63485162832988906167,\
- 0.64644660940672582505, 0.65700282971498213946,\
- 0.66666666666666596353, 0.6755571577384745785 ,\
- 0.68377223398316155922, 0.6913933000758156755 ,\
- 0.69848865542223581571
+ 0. , 0.00985245702332569404,\
+ 0.01941932430907984675, 0.02871413764273582359,\
+ 0.03774955135062373068, 0.04653741075440769731,\
+ 0.08712907082472316767, 0.1835034190722739933 ,\
+ 0.29289321881345248277, 0.36754446796632411765,\
+ 0.42264973081037426894, 0.46547751617515131795,\
+ 0.5 , 0.52859547920896832185,\
+ 0.55278640450004201679, 0.57359856728877911536,\
+ 0.59175170953613676073, 0.60776772972363191094,\
+ 0.62203552699077258037, 0.63485162832988906167,\
+ 0.64644660940672593608, 0.65700282971498225049,\
+ 0.66666666666666607455, 0.67555715773847468952,\
+ 0.68377223398316178127, 0.69139330007581589754,\
+ 0.69848865542223603775
# Scipy computes pdf(x=0) as nan. This has been changed to the limit of 1
# computed by R and matlab.
pdf.values = \
@@ -56,3 +56,18 @@ pdf.values = \
0.03703703703703703498, 0.03415187813279210033,\
0.03162277660168379828, 0.02939111427849369282,\
0.02741012223434214148
+sf.values = \
+ 1. , 0.99014754297667428862,\
+ 0.98058067569092011162, 0.97128586235726432907,\
+ 0.96225044864937647748, 0.95346258924559257331,\
+ 0.91287092917527712377, 0.8164965809277262565 ,\
+ 0.70710678118654768376, 0.63245553203367599338,\
+ 0.57735026918962584208, 0.5345224838248489041 ,\
+ 0.50000000000000011102, 0.47140452079103178917,\
+ 0.44721359549995803873, 0.42640143271122099566,\
+ 0.40824829046386312825, 0.39223227027636820008,\
+ 0.37796447300922730861, 0.36514837167011077179,\
+ 0.35355339059327384188, 0.34299717028501774951,\
+ 0.33333333333333337034, 0.32444284226152514394,\
+ 0.31622776601683799669, 0.30860669992418388041,\
+ 0.30151134457776368469
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.4.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.4.properties
index 9b66af5..92479b1 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.4.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.4.properties
@@ -25,20 +25,20 @@ cdf.points = \
3.0, 3.25, 3.5, 3.75, \
4.0, 4.25, 4.5, 4.75, 5.0
cdf.values = \
- 0.00000000000000000000e+00, 9.29052717957205116076e-05,\
- 4.94825147927420373913e-04, 1.28567736456418856006e-03,\
- 2.49182929403110511174e-03, 4.11522633744856002058e-03,\
- 1.78885438199983225205e-02, 9.17416675955684129962e-02,\
- 2.30048145833311790120e-01, 3.43485618042780493919e-01,\
- 4.31201150371692265573e-01, 4.99534136695643560255e-01,\
- 5.53788770758154358376e-01, 5.97721629744720872601e-01,\
- 6.33938145260609098308e-01, 6.64265367526923222741e-01,\
- 6.90009431395109507079e-01, 7.12123553165658562669e-01,\
- 7.31317294952380514417e-01, 7.48128571959779331557e-01,\
- 7.62971987806307816449e-01, 7.76171921004541243150e-01,\
- 7.87985610946770664853e-01, 7.98619560106298975732e-01,\
- 8.08241385334160011844e-01, 8.16988515143581617295e-01,\
- 8.24974664479918184945e-01
+ 0.00000000000000000000e+00, 9.29052717957204845026e-05,\
+ 4.94825147927420265492e-04, 1.28567736456418834322e-03,\
+ 2.49182929403110467806e-03, 4.11522633744855915322e-03,\
+ 1.78885438199983190510e-02, 9.17416675955683991184e-02,\
+ 2.30048145833311734609e-01, 3.43485618042780438408e-01,\
+ 4.31201150371692154550e-01, 4.99534136695643449233e-01,\
+ 5.53788770758154247353e-01, 5.97721629744720761579e-01,\
+ 6.33938145260608987286e-01, 6.64265367526923111718e-01,\
+ 6.90009431395109396057e-01, 7.12123553165658451647e-01,\
+ 7.31317294952380403394e-01, 7.48128571959779220535e-01,\
+ 7.62971987806307594404e-01, 7.76171921004541021105e-01,\
+ 7.87985610946770442808e-01, 7.98619560106298753688e-01,\
+ 8.08241385334159789799e-01, 8.16988515143581395250e-01,\
+ 8.24974664479917962900e-01
pdf.values = \
0. , 0.02265982238920011618,\
0.05890775570564530522, 0.09966491198172006127,\
@@ -54,3 +54,18 @@ pdf.values = \
0.04477190971288462495, 0.04041085693137497908,\
0.03665493811039276933, 0.03339758876417302236,\
0.03055461720295993625
+sf.values = \
+ 1. , 0.99990709472820427095,\
+ 0.9995051748520725754 , 0.99871432263543580277,\
+ 0.99750817070596886627, 0.99588477366255145906,\
+ 0.98211145618000172952, 0.90825833240443165639,\
+ 0.7699518541666883209 , 0.65651438195721956159,\
+ 0.56879884962830784545, 0.50046586330435649526,\
+ 0.44621122924184580816, 0.40227837025527923842,\
+ 0.36606185473939101271, 0.3357346324730769993 ,\
+ 0.30999056860489049292, 0.28787644683434149284,\
+ 0.26868270504761954109, 0.25187142804022072395,\
+ 0.23702801219369248886, 0.22382807899545911767,\
+ 0.21201438905322947392, 0.20138043989370113529,\
+ 0.19175861466584026571, 0.18301148485641868802,\
+ 0.17502533552008206486
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.5.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.5.properties
index 3c4125c..f627461 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.5.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.5.properties
@@ -54,3 +54,18 @@ pdf.values = \
0.04246189204400971362, 0.0390782759592947504 ,\
0.03612040151810168298, 0.03351702510403412744,\
0.03121149973172263753
+sf.values = \
+ 1. , 0.99999841044682435953,\
+ 0.99996588598985725049, 0.99982068296046555655,\
+ 0.99946266639724357095, 0.99880653316997969338,\
+ 0.98988044026456623037, 0.92661196522925959318,\
+ 0.8123301291303970384 , 0.72493972426305530732,\
+ 0.65910686769794013529, 0.60788777958062267803,\
+ 0.56676396727561950684, 0.53287485882577778629,\
+ 0.5043524956168800033 , 0.47993058333457078168,\
+ 0.45872048292518796675, 0.4400794424630952717 ,\
+ 0.42353027808445919966, 0.40871066227075381061,\
+ 0.39534005574662101079, 0.38319751238636967638,\
+ 0.37210639425702707328, 0.3619236021123256597 ,\
+ 0.35253183118944514396, 0.34383389995642488213,\
+ 0.33574852826886469881\
diff --git a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.6.properties b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.6.properties
index ed0ee31..693fa98 100644
--- a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.6.properties
+++ b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.f.6.properties
@@ -14,18 +14,16 @@
# limitations under the License.
parameters = 100.0 100.0
-# Limited by logpdf values
-tolerance.relative = 6e-13
+# Limited by pdf values
+tolerance.relative = 3e-13
+
# Computed using scipy.stats f
mean = 1.0204081632653061
variance = 0.04295085381091212
lower = 0
cdf.points = \
- 0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.1, 0.25, 0.5, 0.75, \
- 1.0, 1.25, 1.5, 1.75, \
- 2.0, 2.25, 2.5, 2.75, \
- 3.0, 3.25, 3.5, 3.75, \
- 4.0, 4.25, 4.5, 4.75, 5.0
+ 0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.1, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.25, 4.5, 4.75, 5.0
+
cdf.values = \
0.00000000000000000000e+00, 1.90195008441456752357e-72,\
8.15317293024812868281e-58, 1.99842701059155630730e-49,\
@@ -41,11 +39,116 @@ cdf.values = \
9.99999999986713850042e-01, 9.99999999997967625731e-01,\
9.99999999999670818873e-01, 9.99999999999943711693e-01,\
9.99999999999989785948e-01
+
+# TODO:
+# PDF values are inconsistent
+# Reference values from Mathematica, Matlab, Scipy only agree to approximately ~1e-13.
+# scipy is computed using exp(logpdf). The method of Mathematica and Matlab are unknown.
+# Note that the PDF approaches 1 for x=0.75 and x=1.25. This will be computed with
+# cancellation effects using logs (as the log pdf -> 0 with positive and negative terms in
+# a summation)
+
+# Wolfram Mathematica - pass at 8e-13 (fails on the logpdf, no log pdf is available)
+# These match closely with the current implementation.
+# Note these floating-point numbers are inexact and the full decimal representation
+# was obtained from e.g. new BigDecimal(0.1).toString().
+# 0.0 \
+# 9.3251653263638622991e-69 \
+# 1.9599569387331353308e-54 \
+# 3.1406021899708862397e-46 \
+# 1.5827035031843255814e-40 \
+# 3.4071733011342913426e-36 \
+# 1.8303131613029907954e-23 \
+# 1.621261398867281e-09 \
+# 0.0110204220067404 \
+# 0.946217456596435 \
+# 1.9897309346794690375, \
+# 0.8553281090792538 \
+# 0.172291854230817 \
+# 0.02384453848700535 \
+# 0.0027551055016848700217 \
+# 0.000295393685948072 \
+# 0.00003119248021908925 \
+# 3.357137609156691e-6 \
+# 3.7560923959053191116e-7 \
+# 4.417591985387405e-8 \
+# 5.493848943473314e-9 \
+# 7.243777042409837e-10 \
+# 1.0132883742919304427e-10 \
+# 1.502869083243862e-11 \
+# 2.360142057328525e-12 \
+# 3.91738558683421e-13 \
+# 6.8580254516416926559e-14
+
+# Matlab - pass at 7e-13 (fails on the logpdf, no log pdf is available from Matlab)
+# 0 \
+# 9.3251653263634157e-69 \
+# 1.959956938733149e-54 \
+# 3.1406021899708079e-46 \
+# 1.5827035031842939e-40 \
+# 3.4071733011341871e-36 \
+# 1.8303131613029334e-23 \
+# 1.6212613988670735e-09 \
+# 0.0110204220067393 \
+# 0.94621745659646195 \
+# 1.9897309346794265 \
+# 0.85532810907907153 \
+# 0.17229185423079005 \
+# 0.023844538487001965 \
+# 0.0027551055016848541 \
+# 0.00029539368594804681 \
+# 3.1192480219083929e-05 \
+# 3.357137609155976e-06 \
+# 3.7560923959052268e-07 \
+# 4.4175919853867784e-08 \
+# 5.4938489434729238e-09 \
+# 7.2437770424093207e-10 \
+# 1.0132883742919134e-10 \
+# 1.5028690832436914e-11 \
+# 2.3601420573281563e-12 \
+# 3.9173855868335976e-13 \
+# 6.8580254516415638e-14
+
+# scipy.stats - pass at 2e-12
+# 0.00000000000000000000e+00, 9.32516532636367921761e-69,\
+# 1.95995693873326028191e-54, 3.14060218997071882346e-46,\
+# 1.58270350318433878314e-40, 3.40717330113428396147e-36,\
+# 1.83031316130290750821e-23, 1.62126139886700438187e-09,\
+# 1.10204220067397704041e-02, 9.46217456596381234490e-01,\
+# 1.98973093467951134272e+00, 8.55328109079107945512e-01,\
+# 1.72291854230807173209e-01, 2.38445384870012849765e-02,\
+# 2.75510550168477563390e-03, 2.95393685948038410387e-04,\
+# 3.11924802190892479052e-05, 3.35713760915630972549e-06,\
+# 3.75609239590597380614e-07, 4.41759198538690344890e-08,\
+# 5.49384894347331426951e-09, 7.24377704240983564821e-10,\
+# 1.01328837429192775725e-10, 1.50286908324386235308e-11,\
+# 2.36014205732825684207e-12, 3.91738558683376469201e-13,\
+# 6.85802545164195382838e-14
+
+# R - Pass at 1e-12
+# 0.0000000000000000000e+00 9.3251653263639139462e-69 \
+# 1.9599569387331759554e-54 3.1406021899709451856e-46 \
+# 1.5827035031843542807e-40 3.4071733011343073482e-36 \
+# 1.8303131613029665768e-23 1.6212613988670885475e-09 \
+# 1.1020422006739544890e-02 9.4621745659648370808e-01 \
+# 1.9897309346794691542e+00 8.5532810907908718434e-01 \
+# 1.7229185423079382278e-01 2.3844538487002037846e-02 \
+# 2.7551055016848810184e-03 2.9539368594805125818e-04 \
+# 3.1192480219084748466e-05 3.3571376091560772149e-06 \
+# 3.7560923959053067677e-07 4.4175919853868194073e-08 \
+# 5.4938489434729833973e-09 7.2437770424094313637e-10 \
+# 1.0132883742919338319e-10 1.5028690832437327830e-11 \
+# 2.3601420573281809095e-12 3.9173855868336505912e-13 \
+# 6.8580254516417531422e-14
+
+# This is a hybrid set of data using Mathematica and the previous values
+# from a different version of scipy.stats. This data allows the lowest test
+# tolerance of 3e-13 to be used.
pdf.values = \
- 0.00000000000000000000e+00, 9.32516532636421044555e-69,\
- 1.95995693873337184787e-54, 3.14060218997089734622e-46,\
- 1.58270350318442891391e-40, 3.40717330113447773692e-36,\
- 1.83031316130301153946e-23, 1.62126139886709640572e-09,\
+ 0.0, 9.3251653263638622991e-69,\
+ 1.9599569387331353308e-54, 3.1406021899708862397e-46,\
+ 1.5827035031843255814e-40, 3.4071733011342913426e-36,\
+ 1.8303131613029907954e-23, 1.621261398867281e-09,\
1.10204220067403966393e-02, 9.46217456596435080307e-01,\
1.98973093467962458547e+00, 8.55328109079156573280e-01,\
1.72291854230816970928e-01, 2.38445384870026415303e-02,\
@@ -56,19 +159,19 @@ pdf.values = \
1.01328837429192775725e-10, 1.50286908324386235308e-11,\
2.36014205732825684207e-12, 3.91738558683376469201e-13,\
6.85802545164195382838e-14
-# 1-CDF is not accurate as the CDF -> 1.0 so include the explicit SF values
+
sf.values = \
- 1.0000000000000000e+00, 9.9999999999999989e-01,\
- 9.9999999999999989e-01, 9.9999999999999989e-01,\
- 9.9999999999999989e-01, 9.9999999999999989e-01,\
- 9.9999999999999989e-01, 9.9999999998671385e-01,\
- 9.9969086313739297e-01, 9.2397996563976537e-01,\
- 5.0000000000000000e-01, 1.3311560227797117e-01,\
- 2.1930442130085225e-02, 2.7750470811179150e-03,\
- 3.0913686260705529e-04, 3.2868343740149095e-05,\
- 3.4975350236030054e-06, 3.8297320241308051e-07,\
- 4.3848470459334639e-08, 5.2963327957952971e-09,\
- 6.7791444025696802e-10, 9.2110771275951864e-11,\
- 1.3286199631860416e-11, 2.0324013825740018e-12,\
- 3.2917389024730059e-13, 5.6335109943945702e-14,\
- 1.0165268915392895e-14
+ 1.00000000000000000000e+00, 9.99999999999999888978e-01,\
+ 9.99999999999999888978e-01, 9.99999999999999888978e-01,\
+ 9.99999999999999888978e-01, 9.99999999999999888978e-01,\
+ 9.99999999999999888978e-01, 9.99999999986713850042e-01,\
+ 9.99690863137392971005e-01, 9.23979965639765365992e-01,\
+ 5.00000000000000000000e-01, 1.33115602277971173395e-01,\
+ 2.19304421300852252219e-02, 2.77504708111791500644e-03,\
+ 3.09136862607055286790e-04, 3.28683437401490946344e-05,\
+ 3.49753502360300536433e-06, 3.82973202413080513972e-07,\
+ 4.38484704593346391169e-08, 5.29633279579529713456e-09,\
+ 6.77914440256968018369e-10, 9.21107712759518643957e-11,\
+ 1.32861996318604157343e-11, 2.03240138257400182105e-12,\
+ 3.29173890247300592167e-13, 5.63351099439457020525e-14,\
+ 1.01652689153928951395e-14