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 5ad1a09f Compute survival function using summation
5ad1a09f is described below
commit 5ad1a09f7caf84a79cb000d3205629d551ec2309
Author: Alex Herbert <[email protected]>
AuthorDate: Fri Aug 28 17:37:43 2026 +0100
Compute survival function using summation
This change sums the terms excluded from the CDF so the CDF + SF raw
sums are the N-th harmonic number.
Requires updates to the high precision SF test. The expected results
have been independently computed in high precision. The previous values
from scipy.stats.zipfian had differences in the second significant digit
to these updated results.
Increases precision when x approaches n. The previous method computed
the n-th harmonic number minus the (x+1)-th harmonic number. These may
be very close and there is a loss of precision in the difference
compared to summing the upper terms to create the difference.
---
.../statistics/distribution/ZipfDistribution.java | 34 +++++++++-------------
.../distribution/ZipfDistributionTest.java | 28 +++++++++++++-----
.../statistics/distribution/test.zipf.3.properties | 2 ++
3 files changed, 37 insertions(+), 27 deletions(-)
diff --git
a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/ZipfDistribution.java
b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/ZipfDistribution.java
index 2f1d12fd..49d0ae0f 100644
---
a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/ZipfDistribution.java
+++
b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/ZipfDistribution.java
@@ -72,7 +72,7 @@ public final class ZipfDistribution extends
AbstractDiscreteDistribution {
double exponent) {
this.numberOfElements = numberOfElements;
this.exponent = exponent;
- this.nthHarmonic = generalizedHarmonic(numberOfElements, exponent);
+ this.nthHarmonic = generalizedHarmonic(1, numberOfElements, exponent);
logNthHarmonic = Math.log(nthHarmonic);
}
@@ -151,7 +151,7 @@ public final class ZipfDistribution extends
AbstractDiscreteDistribution {
return 1;
}
- return generalizedHarmonic(x, exponent) / nthHarmonic;
+ return generalizedHarmonic(1, x, exponent) / nthHarmonic;
}
/** {@inheritDoc} */
@@ -163,18 +163,9 @@ public final class ZipfDistribution extends
AbstractDiscreteDistribution {
return 0;
}
- // See http://www.math.wm.edu/~leemis/chart/UDR/PDFs/Zipf.pdf
- // S(x) = P(X > x) = ((x+1)^a Hn,a - (x+1)^a Hx+1,a + 1) / ((x+1)^a
Hn,a)
- // where a = exponent and Hx,a is the generalized harmonic for x with
exponent a.
- final double z = Math.pow(x + 1.0, exponent);
- // Compute generalizedHarmonic(x, exponent) and
generalizedHarmonic(x+1, exponent)
- final double hx = generalizedHarmonic(x, exponent);
- final double hx1 = hx + Math.pow(x + 1.0, -exponent);
- // Compute the survival function
- final double p = (z * (nthHarmonic - hx1) + 1) / (z * nthHarmonic);
- // May overflow for large exponent so validate the probability.
- // If this occurs revert to 1 - CDF(x), reusing the generalized
harmonic for x
- return p <= 1.0 ? p : 1.0 - hx / nthHarmonic;
+ // Compute summation of terms omitted in the CDF.
+ // The raw sums in CDF(x) + SF(x) = N-th harmonic
+ return generalizedHarmonic(x + 1, numberOfElements, exponent) /
nthHarmonic;
}
/**
@@ -249,26 +240,29 @@ public final class ZipfDistribution extends
AbstractDiscreteDistribution {
}
/**
- * Calculates the Nth generalized harmonic number. See
+ * Calculates the sum of terms of the
* <a href="https://mathworld.wolfram.com/HarmonicSeries.html";>Harmonic
* Series</a>.
*
* <pre>
* 1
- * sum ----- for k in [1, n]
+ * sum ----- for k in [from, to]
* k^m
* </pre>
*
+ * <p>When {@code from = 1} the result is the N-th harmonic number where
{@code N = to}.
+ *
* <p>Assumes {@code exponent > 0} to arrange the terms to sum from small
to large.
*
- * @param n Term in the series to calculate (must be larger than 0)
+ * @param from First term in the series to calculate.
+ * @param to Last term in the series to calculate.
* @param m Exponent (special case {@code m = 1} is the harmonic series).
- * @return the n<sup>th</sup> generalized harmonic number.
+ * @return the sum
*/
- private static double generalizedHarmonic(final int n, final double m) {
+ private static double generalizedHarmonic(final int from, final int to,
final double m) {
double value = 0;
// Sum small to large
- for (int k = n; k >= 1; k--) {
+ for (int k = to; k >= from; k--) {
value += Math.pow(k, -m);
}
return value;
diff --git
a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/ZipfDistributionTest.java
b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/ZipfDistributionTest.java
index a7f0c723..fcb349d5 100644
---
a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/ZipfDistributionTest.java
+++
b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/ZipfDistributionTest.java
@@ -87,20 +87,34 @@ class ZipfDistributionTest extends
BaseDiscreteDistributionTest {
}
static Stream<Arguments> testAdditionlSurvivalProbabilityHighPrecision() {
- // computed using scipy.stats (1.7.1) zipfian
+ // Computed from the generalized harmonic number and the upper
+ // series of terms using Matlab R2023a VPA, e.g.:
+ // vpa(symsum(1/k^10, k, 1, 60), 30)
+ // vpa(symsum(1/k^10, k, 58, 60), 30)
+ // vpa(symsum(1/k^10, k, 60, 60), 30)
+ // Generalized harmonic numbers, the upper summations are inlined below
+ final double k60e10 = 1.00099457512781807511565108861;
+ final double k60e505 = 1.00000000000000062803698427773;
+ final double k60e1005 = 1.0;
return Stream.of(
Arguments.of(60, 10,
new int[] {57, 59},
- new double[] {2.3189337454689757e-18, 1.6521739576668957e-18},
- DoubleTolerances.absolute(1e-25)),
+ new double[] {
+ 0.00000000000000000593155740928262062795573723513 / k60e10,
+ 0.0000000000000000016538171687920201866246676489 / k60e10},
+ DoubleTolerances.relative(1e-14)),
Arguments.of(60, 50.5,
new int[] {57, 59},
- new double[] {8.8488396450491320e-90, 1.5972093932264611e-90},
- DoubleTolerances.absolute(1e-95)),
+ new double[] {
+ 1.41783221158702775324465028711e-89 / k60e505,
+ 1.59720939322646230873883366414e-90 / k60e505},
+ DoubleTolerances.relative(1e-14)),
Arguments.of(60, 100.5,
new int[] {57, 59},
- new double[] {5.9632998443758656e-178,
1.9760564023408183e-179},
- DoubleTolerances.absolute(1e-185))
+ new double[] {
+ 7.23087851232732627244617172568e-178 / k60e1005,
+ 1.9760564023408181841715991846e-179 / k60e1005},
+ DoubleTolerances.relative(1e-14))
);
}
diff --git
a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zipf.3.properties
b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zipf.3.properties
index 265da1ef..47679093 100644
---
a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zipf.3.properties
+++
b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zipf.3.properties
@@ -45,3 +45,5 @@ sf.values = \
0.01185995655971061 , 0.008563148175057236 ,
\
0.0058478019355387234, 0.003576338667001895 ,
\
0.0016508814447150884, 0.
+
+sf.relative = 5e-14
