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 8f1a59a63fc919a7ca644f07c4d8175d91bba0c7 Author: Alex Herbert <[email protected]> AuthorDate: Fri Aug 28 16:01:43 2026 +0100 Cache the values used for the mean and variance. Avoids expensive repeat computation of Order(N) for the moments. This benefits the inverse cumulative probability functions which use the mean and variance. --- .../statistics/distribution/ZipfDistribution.java | 55 +++++++++++++++++----- .../distribution/ZipfDistributionTest.java | 21 ++++++++- 2 files changed, 63 insertions(+), 13 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 8d92e1aa..d31b9b72 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 @@ -35,7 +35,7 @@ import org.apache.commons.rng.sampling.distribution.RejectionInversionZipfSample * generalized harmonic number</a> of order N of s. * * <p><strong>Note:</strong> The generalized harmonic number \( H_{N,s} \) is computed - * by direct summation of \( N \) terms. Construction of the distribution, and each + * by direct summation of \( N \) terms. Construction of the distribution, and the first * call to {@link #getMean()} or {@link #getVariance()}, is \( O(N) \); each call to * {@link #cumulativeProbability(int) cumulativeProbability(x)} or * {@link #survivalProbability(int) survivalProbability(x)} is \( O(x) \) (the @@ -59,6 +59,10 @@ public final class ZipfDistribution extends AbstractDiscreteDistribution { private final double nthHarmonic; /** Cached value of the log of the nth generalized harmonic. */ private final double logNthHarmonic; + /** Cached value of the nth generalized harmonic using (exponent - 1). */ + private double nthHarmonicM1 = Double.NaN; + /** Cached value of the nth generalized harmonic using (exponent - 2). */ + private double nthHarmonicM2 = Double.NaN; /** Create an instance. * @param numberOfElements Number of elements. @@ -186,11 +190,7 @@ public final class ZipfDistribution extends AbstractDiscreteDistribution { */ @Override public double getMean() { - final int N = getNumberOfElements(); - final double s = getExponent(); - - final double Hs1 = generalizedHarmonicAscendingSum(N, s - 1); - + final double Hs1 = nthHarmonicExpMinus1(); return Hs1 / nthHarmonic; } @@ -207,14 +207,45 @@ public final class ZipfDistribution extends AbstractDiscreteDistribution { */ @Override public double getVariance() { - final int N = getNumberOfElements(); - final double s = getExponent(); - - final double Hs2 = generalizedHarmonicAscendingSum(N, s - 2); - final double Hs1 = generalizedHarmonicAscendingSum(N, s - 1); + final double Hs2 = nthHarmonicExpMinus2(); + final double Hs1 = nthHarmonicExpMinus1(); final double Hs = nthHarmonic; + // (Hs2 / Hs) - ((Hs1 * Hs1) / (Hs * Hs)) + // Values are ascending magnitude: Hs < Hs1 < Hs2. + // (small * large) - (mid * mid) with a common denominator: + return (Hs2 * Hs - Hs1 * Hs1) / (Hs * Hs); + } + + /** + * Compute the N-th harmonic number using the {@code exponent - 1}. + * This is cached to avoid repeat expensive computation across all N + * for the mean and variance. + * + * @return the number + */ + private double nthHarmonicExpMinus1() { + double h = nthHarmonicM1; + if (Double.isNaN(h)) { + h = generalizedHarmonicAscendingSum(getNumberOfElements(), getExponent() - 1); + nthHarmonicM1 = h; + } + return h; + } - return (Hs2 / Hs) - ((Hs1 * Hs1) / (Hs * Hs)); + /** + * Compute the N-th harmonic number using the {@code exponent - 2}. + * This is cached to avoid repeat expensive computation across all N + * for the variance. + * + * @return the number + */ + private double nthHarmonicExpMinus2() { + double h = nthHarmonicM2; + if (Double.isNaN(h)) { + h = generalizedHarmonicAscendingSum(getNumberOfElements(), getExponent() - 2); + nthHarmonicM2 = h; + } + return h; } /** 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 403f165e..a7f0c723 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 @@ -24,13 +24,14 @@ import org.junit.jupiter.api.Assertions; import org.junit.jupiter.api.Test; import org.junit.jupiter.params.ParameterizedTest; import org.junit.jupiter.params.provider.Arguments; +import org.junit.jupiter.params.provider.CsvSource; import org.junit.jupiter.params.provider.MethodSource; /** * Test cases for {@link ZipfDistribution}. * Extends {@link BaseDiscreteDistributionTest}. See javadoc of that class for details. */ -class ZipfDistributionTest extends BaseDiscreteDistributionTest { +class ZipfDistributionTest extends BaseDiscreteDistributionTest { @Override DiscreteDistribution makeDistribution(Object... parameters) { final int n = (Integer) parameters[0]; @@ -61,6 +62,24 @@ class ZipfDistributionTest extends BaseDiscreteDistributionTest { //-------------------- Additional test cases ------------------------------- + /** + * Test additional moments. + */ + @ParameterizedTest + @CsvSource({ + // Generated using scipy 1.16.3 using scipy.stats.zipfian.stats(exp, n) + "150, 0.512, 52.707637767916495, 1966.9356468021338", + "73, 1.67, 4.937625767687036, 87.76033876340095", + "999, 2.1, 3.5725516349635846, 343.7153292773371", + }) + void testAdditionalMoments(int n, double exp, double mean, double variance) { + final DoubleTolerance tolerance = createRelTolerance(1e-14); + final ZipfDistribution dist = ZipfDistribution.of(n, exp); + testMoments(dist, mean, variance, tolerance); + // Run twice to check the cached N-th harmonic numbers + testMoments(dist, mean, variance, tolerance); + } + @ParameterizedTest @MethodSource void testAdditionlSurvivalProbabilityHighPrecision(int n, double e, int[] x, double[] expected, DoubleTolerance tol) {
