This is an automated email from the ASF dual-hosted git repository.
asf-gitbox-commits 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 d2e3525f STATISTICS-101: Add a zeta distribution
d2e3525f is described below
commit d2e3525f90ac880971b75d5d11b538bcca108c4b
Author: Alex Herbert <[email protected]>
AuthorDate: Mon Sep 14 19:51:44 2026 +0100
STATISTICS-101: Add a zeta distribution
---
.../statistics/distribution/ZetaDistribution.java | 445 +++++++++++++++++++++
.../distribution/ZetaDistributionTest.java | 263 ++++++++++++
.../statistics/distribution/test.zeta.1.properties | 48 +++
.../statistics/distribution/test.zeta.2.properties | 92 +++++
.../statistics/distribution/test.zeta.3.properties | 40 ++
.../statistics/distribution/test.zeta.4.properties | 48 +++
src/changes/changes.xml | 4 +
7 files changed, 940 insertions(+)
diff --git
a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/ZetaDistribution.java
b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/ZetaDistribution.java
new file mode 100644
index 00000000..337748fe
--- /dev/null
+++
b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/ZetaDistribution.java
@@ -0,0 +1,445 @@
+/*
+ * 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
+ *
+ * https://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.
+ */
+
+package org.apache.commons.statistics.distribution;
+
+import java.util.function.ToDoubleFunction;
+import org.apache.commons.rng.UniformRandomProvider;
+
+/**
+ * Implementation of the zeta distribution.
+ *
+ * <p>The probability mass function of \( X \) is:
+ *
+ * <p>\[ f(k; s) = \frac{k^-s}{\zeta(s)} \]
+ *
+ * <p>for \( s \gt 1 \) the exponent characterizing the distribution,
+ * \( k \in \{1, 2, \dots, \infty \} \), and
+ * \( \zeta(s) \) is the Riemann zeta function.
+ *
+ * <p>\[ \sum_{k=1}^\infty \frac{1}{k^s} \]
+ *
+ * <p><strong>Implementation note</strong>
+ *
+ * <p>The zeta distribution is the large \( N \) limit of the {@link
ZipfDistribution}.
+ * Note that the zeta distribution has an upper limit of positive infinity.
This
+ * implementation is clipped to {@link Integer#MAX_VALUE} to implement the
+ * {@link DiscreteDistribution} interface. As \( s \to 1 \) the survival
probability
+ * function ({@code sf}) at {@code x} = 2<sup>31</sup> - 1 will become
significantly larger
+ * than 0 as the distribution is truncated:
+ *
+ * <table border=1>
+ * <caption><b>Survival function values</b></caption>
+ * <tr><th>s</th><th>x</th><th>sf(x; s)</th>
+ * <tr><td>1.20889...</td> <td>2147483647</td> <td>0.01</td>
+ * <tr><td>1.10440...</td> <td>2147483647</td> <td>0.1</td>
+ * <tr><td>1.03141...</td> <td>2147483647</td> <td>0.5</td>
+ * <tr><td>1.00477...</td> <td>2147483647</td> <td>0.9</td>
+ * </table>
+ *
+ * <p>The sum of the power series can be computed using the <a
+ * href="https://en.wikipedia.org/wiki/Hurwitz_zeta_function";>Hurwitz zeta
function</a>:
+ *
+ * <p>\[ \zeta(s, a) = \sum_{n=0}^\infty \frac{1}{(n+a)^s} \]
+ *
+ * <p>The survival probability function is implemented using \( \zeta(s, a) \)
with efficient
+ * evaluation time and no loss of precision. The sum of the power series in
cumulative probability
+ * functions uses a direct sum of the probability mass function when the
number of terms is
+ * practical. In all other cases the implementation uses a difference of zeta
functions:
+ *
+ * <p>\[ \begin{aligned}
+ * \sum_{k=a}^b \frac{1}{k^s} &= \sum_{n=0}^\infty \frac{1}{(n+a)^s}
- \sum_{n=0}^\infty \frac{1}{(n+b+1)^s} \\
+ * &= \zeta(s, a) - \zeta(s, b+1)
\end{aligned} \]
+ *
+ * <p>This may incur significant cancellation in the two zeta terms. The
+ * {@link #cumulativeProbability(int) cdf(x)} method should be considered
imprecise as the
+ * value approaches 0. Precision loss is limited to {@code b}-bits for a value
\( > 2^{-b+1} \).
+ * The {@link #probability(int, int) probability(a, b)} method should be
considered potentially
+ * imprecise for any range {@code b - a >= 10}.
+ *
+ * <p>The following table shows the value of {@code x} for different {@code s}
where the survival
+ * function is 0.875. Computing {@code cdf(x) = 1 - sf(x)} is expected to lose
2-bits of precision.
+ * A sum of the power series requires an increasingly large run-time cost for
a similar precision
+ * result.
+ *
+ * <table border=1>
+ * <caption><b>Survival function values</b></caption>
+ * <tr><th>s</th><th>x</th><th>sf(x; s)</th>
+ * <tr><td>1.04565...</td> <td>10</td> <td>0.875</td>
+ * <tr><td>1.02575...</td> <td>100</td> <td>0.875</td>
+ * <tr><td>1.01784...</td> <td>1000</td> <td>0.875</td>
+ * <tr><td>1.01364...</td> <td>10000</td> <td>0.875</td>
+ * <tr><td>1.01104...</td> <td>100000</td> <td>0.875</td>
+ * <tr><td>1.00927...</td> <td>1000000</td> <td>0.875</td>
+ * <tr><td>1.00799...</td> <td>10000000</td> <td>0.875</td>
+ * </table>
+ *
+ * @see <a href="https://en.wikipedia.org/wiki/Zeta_distribution";>Zeta
distribution (Wikipedia)</a>
+ * @see <a href="https://en.wikipedia.org/wiki/Riemann_zeta_function";>Riemann
zeta function (Wikipedia)</a>
+ * @see <a href="https://en.wikipedia.org/wiki/Hurwitz_zeta_function";>Hurwitz
zeta function (Wikipedia)</a>
+ */
+public final class ZetaDistribution extends AbstractDiscreteDistribution {
+ /** Minimum number of terms required to use the Hurwitz zeta function for
cumulative
+ * probability functions. Below this level a regular sum of the terms is
used.
+ * Note the evaluation of the Hurwitz zeta function requires multiple
calls to
+ * {@link Math#pow(double, double)}. If the number of probability terms is
low it is
+ * more efficient to sum the terms directly. */
+ private static final int MIN_TERMS = 10;
+ /** 2 as a {@code double}. */
+ private static final double TWO = 2;
+ /** 3 as a {@code double}. */
+ private static final double THREE = 3;
+
+ /** Exponent parameter of the distribution. */
+ private final double exponent;
+ /** zeta(s, 1) where s is the exponent of the distribution. */
+ private final double zeta1;
+ /** Cached value of the log(zeta(s, 1)). */
+ private final double logZeta1;
+ /** Cached value of the mean. */
+ private double mean = Double.NaN;
+ /** Cached value of the variance. */
+ private double variance = Double.NaN;
+
+ /** Create an instance.
+ * @param exponent Exponent (s).
+ */
+ private ZetaDistribution(double exponent) {
+ this.exponent = exponent;
+ this.zeta1 = HurwitzZeta.value(exponent, 1);
+ logZeta1 = Math.log(zeta1);
+ }
+
+ /**
+ * Creates a zeta distribution.
+ *
+ * @param exponent Exponent.
+ * @return the distribution
+ * @exception IllegalArgumentException if {@code exponent <= 1} or is
{@code NaN}.
+ */
+ public static ZetaDistribution of(double exponent) {
+ if (!(exponent > 1)) {
+ // <= 1 or nan
+ throw new DistributionException(DistributionException.NEGATIVE,
+ exponent);
+ }
+ return new ZetaDistribution(exponent);
+ }
+
+ /**
+ * Gets the exponent parameter of this distribution.
+ *
+ * @return the exponent.
+ */
+ public double getExponent() {
+ return exponent;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double probability(final int x) {
+ if (x <= 0) {
+ return 0;
+ }
+ return Math.pow(x, -exponent) / zeta1;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * <p><strong>Implementation note</strong>
+ *
+ * <p>This is implemented using a direct sum of the probability mass
function when the
+ * size {@code x1 - x0} is practical. In all other cases the
implementation uses the
+ * survival probability function using the identity
+ * {@code P(x0 < X <= x1) = P(X > x0) - P(X > x1)}. This will suffer loss
of precision
+ * when the two survival functions are of a similar magnitude.
+ *
+ * @see #survivalProbability(int)
+ */
+ @Override
+ public double probability(int x0, int x1) {
+ if (x0 > x1) {
+ throw new
DistributionException(DistributionException.INVALID_RANGE_LOW_GT_HIGH, x0, x1);
+ }
+ if (x0 == x1 || x1 < 1) {
+ // (x0, x1] does not overlap [1, infinity]
+ return 0;
+ }
+ // If the lower range is outside the bounds use the cumulative
probability.
+ // Note: This cannot compare x1 to the unlimited upper bound to use the
+ // survival probability.
+ if (x0 < 1) {
+ return cumulativeProbability(x1);
+ }
+ return genHarmonic(x0 + 1, x1) / zeta1;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double logProbability(int x) {
+ if (x <= 0) {
+ return Double.NEGATIVE_INFINITY;
+ }
+ return -Math.log(x) * exponent - logZeta1;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * <p><strong>Implementation note</strong>
+ *
+ * <p>This is implemented using a direct sum of the probability mass
function when the
+ * size {@code x} is practical. In all other cases the implementation uses
+ * {@code 1 - survivalProbability(x)}. This will suffer loss of precision
when the
+ * survival function is above 0.5.
+ *
+ * <p>When the exponent {@code s} of the distribution approaches 1 a
significant part
+ * of the probability mass is truncated at 2<sup>31</sup> - 1, the CDF is
small for
+ * any {@code x}, and this function loses many bits of precision.
+ *
+ * <p>A high-precision method to compute the CDF is to sum the power terms
directly. This
+ * can be done from small to large terms; the normalizing constant
+ * \( 1 / \zeta(s) \) is provided using the probability mass function at
{@code x=1}.
+ *
+ * <pre>{@code
+ * double s = ...;
+ * ZetaDistribution dist = ZetaDistribution.of(s);
+ * double cdfX = IntStream.range(0, x)
+ * .mapToDouble(y -> Math.pow(x - y, -s))
+ * .sum() * dist.probability(1);
+ * }</pre>
+ *
+ * <p>Direct summation may require a high run-time cost and this should
only be used
+ * when the argument {@code x} has been bounded to a reasonable range.
+ *
+ * @see #survivalProbability(int)
+ * @see #probability(int)
+ */
+ @Override
+ public double cumulativeProbability(int x) {
+ if (x <= 0) {
+ return 0;
+ }
+ return genHarmonic(1, x) / zeta1;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double survivalProbability(int x) {
+ if (x <= 0) {
+ return 1;
+ }
+ // Add 1.0 to support integer max value without overflow
+ return HurwitzZeta.value(exponent, x + 1.0) / zeta1;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * <p>For \( s > 2 \) the mean is:
+ *
+ * <p>\[ \frac{\zeta(s - 1)}{\zeta(s)} \]
+ *
+ * <p>Otherwise the mean is infinity.
+ */
+ @Override
+ public double getMean() {
+ double m = mean;
+ if (Double.isNaN(m)) {
+ if (exponent <= TWO) {
+ m = Double.POSITIVE_INFINITY;
+ } else {
+ m = HurwitzZeta.value(exponent - 1, 1) / zeta1;
+ }
+ mean = m;
+ }
+ return m;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * <p>For \( s > 3 \) the variance is:
+ *
+ * <p>\[ \frac{\zeta(s)\zeta(s-2) - \zeta(s-1)^2}{\zeta(s)^2} \]
+ *
+ * <p>Otherwise the variance is infinity.
+ */
+ @Override
+ public double getVariance() {
+ double v = variance;
+ if (Double.isNaN(v)) {
+ if (exponent <= THREE) {
+ v = Double.POSITIVE_INFINITY;
+ } else {
+ final double s2 = HurwitzZeta.value(exponent - 2, 1);
+ final double s1 = HurwitzZeta.value(exponent - 1, 1);
+ final double s = zeta1;
+ v = (s * s2 - s1 * s1) / (s * s);
+ }
+ variance = v;
+ }
+ return v;
+ }
+
+ /**
+ * Calculates the sum of terms of the generalized
+ * <a href="https://mathworld.wolfram.com/HarmonicSeries.html";>Harmonic
+ * Series</a>.
+ *
+ * <pre>
+ * 1
+ * sum ----- for k in [from, to]
+ * k^m
+ * </pre>
+ *
+ * @param a First term in the series to calculate.
+ * @param b Last term in the series to calculate.
+ * @return the sum
+ */
+ private double genHarmonic(int a, int b) {
+ final double s = exponent;
+ // Entirely define the sum of the series using the zeta function
+ // unless the number of terms is small.
+ // Number of terms is b - a + 1 so use >= not >.
+ if (b - a >= MIN_TERMS) {
+ // a == 1 when used in the CDF
+ final double z1 = a == 1 ? zeta1 : HurwitzZeta.value(s, a);
+ // Add 1.0 to support integer max value without overflow
+ final double z2 = HurwitzZeta.value(s, b + 1.0);
+ return z1 - z2;
+ }
+ return ZipfDistribution.generalizedHarmonic(a, b, s);
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * <p>The lower bound of the support is always 1.
+ *
+ * @return 1.
+ */
+ @Override
+ public int getSupportLowerBound() {
+ return 1;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * <p>The upper bound of the support is always positive infinity.
+ *
+ * @return {@link Integer#MAX_VALUE}
+ */
+ @Override
+ public int getSupportUpperBound() {
+ return Integer.MAX_VALUE;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public DiscreteDistribution.Sampler createSampler(final
UniformRandomProvider rng) {
+ return new ZetaSampler(rng, exponent);
+ }
+
+ /**
+ * Sample from a zeta distribution.
+ * <ol>
+ * <li>Devroye, L (2015)
+ * Non-uniform random variate generation.
+ * Springer New York, NY. pp 550-552.</li>
+ * </ol>
+ */
+ private static final class ZetaSampler implements
DiscreteDistribution.Sampler {
+ /**
+ * The threshold to bias the extreme sample to 1 or infinity. Change
the
+ * extreme sample of the zeta distribution using the midpoint of the
support
+ * domain, i.e. x = 2^31 / 2; cdf(x; a) = sf(x; a) ~ 0.5.
+ */
+ private static final double THRESHOLD = 1.0324376395045163;
+
+ /** Source of randomness. */
+ private final UniformRandomProvider rng;
+ /** a - 1. */
+ private final double am1;
+ /** Reciprocal of (a - 1) = 1 / (a - 1). */
+ private final double ram1;
+ /** b = 2^(a-1). This constants is {@code (b-1) / b}. */
+ private final double bm1Db;
+ /** Function to compute u in [0, 1]. */
+ private final ToDoubleFunction<UniformRandomProvider> nextU;
+
+ /**
+ * Create an instance.
+ *
+ * @param rng Source of randomness.
+ * @param a Exponent of the zeta distribution ({@code a > 1}).
+ */
+ ZetaSampler(UniformRandomProvider rng, double a) {
+ this.rng = rng;
+ am1 = a - 1;
+ ram1 = 1 / am1;
+ final double b = Math.pow(2, am1);
+ bm1Db = b == Double.POSITIVE_INFINITY ? 1 : (b - 1) / b;
+ // Note:
+ // u in [0, 1]
+ // u == 0 : x == inf
+ // u == 1 : x == 1
+ // When a -> 1 then bias to infinity; otherwise bias to 1.
+ nextU = a <= THRESHOLD ?
+ // u in [0, 1)
+ UniformRandomProvider::nextDouble :
+ // u in (0, 1]
+ g -> 1.0 - g.nextDouble();
+ }
+
+ @Override
+ public int sample() {
+ double u;
+ double v;
+ double x;
+ double t;
+ for (;;) {
+ // Generate iid uniform [0, 1] random variate U, V.
+ u = nextU.applyAsDouble(rng);
+ v = rng.nextDouble();
+ // X = floor ( U^{-1/(a-1)} ) , X in [1, inf]
+ x = Math.floor(Math.pow(u, -ram1));
+ t = Math.pow(1 + 1 / x, am1);
+
+ // Until:
+ // T-1 T
+ // VX --- <= -
+ // b-1 b
+
+ // If (a-1) -> inf then t & b -> inf; b >= t
+ // Avoid inf / inf = NaN and accept.
+ // Large a will mostly sample X=1.
+
+ // v * x * (t - 1) / (b - 1) <= t / b
+ // Rearrange terms to ratios of similar magnitude and guard
infinity:
+ // v * x <= (t / (t - 1)) * ((b - 1) / b)
+ final double tDtm1 = t == Double.POSITIVE_INFINITY ? 1 : t /
(t - 1);
+ if (v * x <= tDtm1 * bm1Db) {
+ // Truncates x >= 2^31 to integer max
+ return (int) x;
+ }
+ }
+ }
+ }
+}
diff --git
a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/ZetaDistributionTest.java
b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/ZetaDistributionTest.java
new file mode 100644
index 00000000..00ca8d9d
--- /dev/null
+++
b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/ZetaDistributionTest.java
@@ -0,0 +1,263 @@
+/*
+ * 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
+ *
+ * https://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.
+ */
+
+package org.apache.commons.statistics.distribution;
+
+import java.util.Arrays;
+import java.util.function.DoubleUnaryOperator;
+import java.util.stream.IntStream;
+import java.util.stream.Stream;
+import org.apache.commons.numbers.core.DD;
+import org.apache.commons.numbers.rootfinder.BrentSolver;
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.simple.RandomSource;
+import org.apache.commons.statistics.distribution.DiscreteDistribution.Sampler;
+import org.junit.jupiter.api.Assertions;
+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;
+import org.junit.jupiter.params.provider.ValueSource;
+
+/**
+ * Test cases for {@link ZetaDistribution}.
+ * Extends {@link BaseDiscreteDistributionTest}. See javadoc of that class for
details.
+ */
+class ZetaDistributionTest extends BaseDiscreteDistributionTest {
+ @Override
+ DiscreteDistribution makeDistribution(Object... parameters) {
+ final double e = (Double) parameters[0];
+ return ZetaDistribution.of(e);
+ }
+
+ @Override
+ Object[][] makeInvalidParameters() {
+ return new Object[][] {
+ {-0.1},
+ {1.0},
+ {Double.NaN},
+ };
+ }
+
+ @Override
+ String[] getParameterNames() {
+ return new String[] {"Exponent"};
+ }
+
+ @Override
+ protected double getRelativeTolerance() {
+ return 1e-14;
+ }
+
+ //-------------------- Additional test cases
-------------------------------
+
+ @ParameterizedTest
+ @CsvSource({
+ // Generated using scipy 1.17.1 using scipy.stats.zipf.stats(s)
+ "1.512, Infinity, Infinity, 1e-15",
+ "2.512, 1.919627298627266, Infinity, 1e-15",
+ "3.512, 1.1879319717598735, 0.8692042723730669, 1e-15",
+ })
+ void testAdditionalMoments(double s, double mean, double variance, double
eps) {
+ final DoubleTolerance tolerance = createRelTolerance(eps);
+ final ZetaDistribution dist = ZetaDistribution.of(s);
+ testMoments(dist, mean, variance, tolerance);
+ // Run twice to check the cached values
+ testMoments(dist, mean, variance, tolerance);
+ }
+
+ /**
+ * Test the suggested code in the javadoc for computing the cumulative
probability in
+ * high precision.
+ */
+ @ParameterizedTest
+ @MethodSource
+ void testJavadocCumulativeProbabilityHighPrecision(double s, int x, double
expected) {
+ final DoubleTolerance tolerance = DoubleTolerances.ulps(1);
+ final ZetaDistribution dist = ZetaDistribution.of(s);
+ // Method suggested in the javadoc
+ double cdfX = IntStream.range(0, x)
+ .mapToDouble(y -> Math.pow(x - y, -s))
+ .sum() * dist.probability(1);
+ TestUtils.assertEquals(expected, cdfX, tolerance);
+ }
+
+ static Stream<Arguments> testJavadocCumulativeProbabilityHighPrecision() {
+ return Stream.of(
+ // Computed from scipy.stats (1.17.1) zipf(s)
+ // Note: large x will have a long runtime so are avoided
+ Arguments.of(2.1, 1, 0.6409366767538137),
+ Arguments.of(2.1, 10, 0.9561735629717311),
+ Arguments.of(2.1, 100, 0.9963437521110778),
+ Arguments.of(1.1, 10, 0.2532163286767671),
+ Arguments.of(1.1, 100, 0.40418015518943295),
+ Arguments.of(1.01, 100, 0.05054241012268188),
+ Arguments.of(1.00001, 100, 5.187242037932912e-05),
+ Arguments.of(1.00000001, 100, 5.1873773506313715e-08),
+ Arguments.of(1.0000000000000002, 100, 1.1518291915012786e-15),
+ Arguments.of(1.0000000000000002, 1000, 1.6621084199087341e-15),
+ Arguments.of(1.0000000000000002, 10000, 2.1732851154353238e-15),
+ Arguments.of(1.0000000000000002, 10001, 2.1733073176755921e-15)
+ );
+ }
+
+ /**
+ * Suggested method to compute the CDF in [1, x].
+ */
+ @ParameterizedTest
+ @ValueSource(doubles = {1.3, 2.4})
+ void testCumulativeProbabilityHighPrecisionSinglePass(double s) {
+ final int x = 5;
+ final double zeta1 = HurwitzZeta.value(s, 1);
+ final double[] expected = IntStream.rangeClosed(1, x)
+ .mapToDouble(i -> ZipfDistribution.generalizedHarmonic(1, i, s) /
zeta1)
+ .toArray();
+ final double[] cdf = new double[x];
+ final double p1 = ZetaDistribution.of(s).probability(1);
+ cdf[0] = p1;
+ DD sum = DD.ONE;
+ for (int k = 2; k <= x; k++) {
+ sum = sum.add(Math.pow(k, -s));
+ cdf[k - 1] = sum.doubleValue() * p1;
+ }
+ final DoubleTolerance tolerance = DoubleTolerances.ulps(1);
+ for (int i = 0; i < x; i++) {
+ TestUtils.assertEquals(expected[i], cdf[i], tolerance);
+ }
+ }
+
+ /**
+ * Test inversion correctly returns the lower or upper bound with high
precision p.
+ */
+ @ParameterizedTest
+ @CsvSource({
+ // sf(1, 2, 3) = 8.673617380119933e-19, 2.3589825628243265e-29,
7.523175374682315e-37
+ "60, 1e-16, 1",
+ "60, 1e-20, 2",
+ "60, 1e-28, 2",
+ "60, 1e-29, 3",
+ "60, 1e-38, 4",
+ // sf(1, 2, 3) = 6.223015277861142e-61, 3.764861949599026e-96,
3.8725919148493183e-121
+ "200, 1e-60, 1",
+ "200, 1e-95, 2",
+ "200, 1e-96, 3",
+ "200, 1e-122, 4",
+ // sf(2^31-3, 2^31-2, 2^31-1] = 2.830881171681671e-10,
2.830881170363439e-10, 2.8308811690452073e-10
+ "2, 2.830881171e-10, 2147483646",
+ "2, 2.830881170e-10, 2147483647",
+ "2, 1e-10, 2147483647",
+ "2, 1e-100, 2147483647",
+ // sf(2^31-3, 2^31-2, 2^31-1] = 9.019557827587404e-20,
9.019557819187286e-20, 9.019557810787168e-20
+ "3, 9.01955782e-20, 2147483646",
+ "3, 9.01955781e-20, 2147483647",
+ "3, 1e-20, 2147483647",
+ "3, 1e-200, 2147483647",
+ })
+ void testAdditionalInverseSurvivalFunction(double s, double p, int x) {
+ final ZetaDistribution dist = ZetaDistribution.of(s);
+ Assertions.assertEquals(x, dist.inverseSurvivalProbability(p));
+ }
+
+ @ParameterizedTest
+ @CsvSource({
+ // Tiny s. cdf(2^31 - 1) == 4.8993649719501934e-15
+ "1.0000000000000002, 2147483647",
+ // Large s. sf(1) == 8.673617380119933e-19
+ "60, 1",
+ // Very large s. 2^(s-1) == Infinity
+ "1025, 1",
+ })
+ void testSamplingExtremeS(double s, int x) {
+ final UniformRandomProvider rng =
RandomSource.XO_SHI_RO_256_PP.create(123456789L);
+ final Sampler sampler = ZetaDistribution.of(s).createSampler(rng);
+ final int n = 10;
+ final int[] expected = new int[n];
+ Arrays.fill(expected, x);
+ final int[] sample = TestUtils.sample(n, sampler);
+ Assertions.assertArrayEquals(expected, sample);
+ }
+
+ /**
+ * This is added as the sampler can only be tested when the quantiles of
the distribution
+ * are spread within the range [1, 2^31). The test resources do not test
the sampler
+ * with many values of s.
+ *
+ * <p>The sensitivity of the quantiles to s can be observed using the pmf
at x=1:
+ * <pre>
+ * s pmf(1)
+ * 1.0625 0.06030727407685079
+ * 1.125 0.11646539687154238
+ * 1.25 0.21762256021191748
+ * 1.5 0.3827933839994266
+ * 2.25 0.6848321282518275
+ * </pre>
+ */
+ @ParameterizedTest
+ @ValueSource(doubles = {1.11, 1.12, 1.13, 1.14, 1.15, 1.2, 1.25, 1.3,
1.35})
+ void testAdditionalSampling(double s) {
+ testSampling(ZetaDistribution.of(s));
+ }
+
+ /**
+ * Test the value of the exponent s for critical points in the
distribution.
+ * These points are used in the main ZetaDistribution either in the code or
+ * the documentation. Some values are recorded in the STATISTICS-101 issue.
+ *
+ * @param p the desired survival function probability
+ * @param x the value at which to evaluate the survival function
+ * @param exponent the expected exponent
+ */
+ @ParameterizedTest
+ @CsvSource({
+ // Threshold to switch the sampler extreme value bias from x=1 to x=inf
+ "0.5, 1073741824, 1.0324376395045163",
+ // Used in the class javadoc to describe truncation of the distribution
+ "0.01, 2147483647, 1.2088900037546617",
+ "0.10, 2147483647, 1.1044011399570164",
+ "0.50, 2147483647, 1.0314183630697709",
+ "0.90, 2147483647, 1.0047751511822856",
+ // Test threshold for significant cancellation in the cumulative
probability.
+ // Each p value is 1 - 2^-b where (b-1) is the number of bits of
precision
+ // lost by cancellation in the CDF.
+ "0.5, 10, 1.238388603270506", // 0-bits
+ "0.75, 10, 1.0985236087748393", // 1-bit
+ "0.875, 10, 1.0456558205227693", // 2-bits
+ "0.9375, 10, 1.0220499988943923", // 3-bits
+ "0.96875, 10, 1.0108432918659283", // 4-bits
+ // Number of terms to sum to avoid a 2-bit loss of precision for
various s
+ "0.875, 100, 1.0257535317691207",
+ "0.875, 1000, 1.0178427184412626",
+ "0.875, 10000, 1.0136446856290058",
+ "0.875, 100000, 1.0110455910807687",
+ "0.875, 1000000, 1.0092782593613394",
+ "0.875, 10000000, 1.0079984961118713",
+ })
+ void testExponent(double p, int x, double exponent) {
+ // Search survival function so that sf(x; s) ~ p
+ final DoubleUnaryOperator fun =
+ s -> HurwitzZeta.value(s, x + 1.0) / HurwitzZeta.value(s, 1) - p;
+ final BrentSolver solver = new BrentSolver(0x1.0p-53,
Double.MIN_VALUE, Double.MIN_VALUE);
+ // Bracket using the expected value
+ final double lower = exponent * 0.5;
+ final double upper = exponent * 2;
+ final double s = solver.findRoot(fun, lower, upper);
+ final ZetaDistribution dist = ZetaDistribution.of(s);
+ final double q = dist.survivalProbability(x);
+ TestUtils.assertEquals(p, q, DoubleTolerances.relative(5e-15));
+ TestUtils.assertEquals(exponent, s, DoubleTolerances.relative(1e-15));
+ }
+}
diff --git
a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zeta.1.properties
b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zeta.1.properties
new file mode 100644
index 00000000..3ad90a2d
--- /dev/null
+++
b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zeta.1.properties
@@ -0,0 +1,48 @@
+# 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 = 6.6
+mean = 1.0120796439951505
+variance = 0.014616754077960703
+lower = 1
+upper = 2147483647
+cdf.points = -1, 0, 1, 2, 3, 4, 5, 10, 100, 1000
+# Reference values are from scipy.stats (1.17.1) zipf(6.6)
+cdf.values = \
+ 0.,\
+ 0. , 0.9889621791207941, 0.9991570495802089,\
+ 0.9998587947015799, 0.9999638901095526, 0.9999879879239209,\
+ 0.9999996670903828, 0.9999999999989163, 0.9999999999999997
+pmf.values = \
+ 0,\
+ 0.0000000000000000e+00, 9.8896217912079409e-01,\
+ 1.0194870459414751e-02, 7.0174512137096204e-04,\
+ 1.0509540797267701e-04, 2.4097814368311538e-05,\
+ 2.4841606790096716e-07, 6.2399295032940255e-14,\
+ 1.5673994252900586e-20
+# The scipy zipf sf is entirely defined by the cdf
+# Remove points that are inaccurate:
+# x=100 : 1.0836886943366153e-12
+# x=1000 : 3.3306690738754696e-16
+sf.points = -1, 0, 1, 2, 3, 4, 5, 10
+sf.values = \
+ 1,\
+ 1.0000000000000000e+00, 1.1037820879205906e-02,\
+ 8.4295041979109975e-04, 1.4120529842009955e-04,\
+ 3.6109890447444393e-05, 1.2012076079148670e-05,\
+ 3.3290961720666701e-07
+
+moments.relative = 5e-14
+sf.relative = 5e-10
diff --git
a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zeta.2.properties
b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zeta.2.properties
new file mode 100644
index 00000000..c44305ac
--- /dev/null
+++
b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zeta.2.properties
@@ -0,0 +1,92 @@
+# 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 = 6.6
+mean = 1.0120796439951505
+variance = 0.014616754077960871
+lower = 1
+upper = 2147483647
+cdf.points = 0, 1, 2, 3, 4, 5, 10, 100, 1000,\
+ 10000, 100000, 1000000, 10000000, 100000000, 1000000000, 2147483647
+
+# Reference values are from a custom Python implementation using mpmath based
on an extended
+# precision zeta function:
+#from mpmath import zeta, mp
+#import math
+## At east twice the 17 digits of a double
+#mp.dps = 50
+#class zetad:
+# def __init__(self, s):
+# self._s = s
+# self._z1 = zeta(s)
+# def pmf(self, x):
+# try:
+# i = iter(x)
+# return [self._pmf(_) for _ in i]
+# except TypeError:
+# return self._pmf(x)
+# def _pmf(self, x):
+# if x <= 0:
+# return 0
+# return float(mp.mpf(x)**-self._s / self._z1)
+# def cdf(self, x):
+# try:
+# i = iter(x)
+# return [self._cdf(_) for _ in i]
+# except TypeError:
+# return self._cdf(x)
+# def _cdf(self, x):
+# if x <= 0:
+# return 0
+# return float((self._z1 - zeta(self._s, x+1)) / self._z1)
+# def sf(self, x):
+# try:
+# i = iter(x)
+# return [self._sf(_) for _ in i]
+# except TypeError:
+# return self._sf(x)
+# def _sf(self, x):
+# if x <= 0:
+# return 1
+# return float(zeta(self._s, x+1) / self._z1)
+# def moments(self):
+# s = self._s
+# m, v = math.inf, math.inf
+# if s > 2:
+# s1 = zeta(s - 1)
+# s0 = self._z1
+# m = s1 / s0
+# if v > 3:
+# s2 = zeta(s - 2)
+# v = (s0 * s2 - s1 * s1) / (s0 * s0)
+ return float(m), float(v)
+
+cdf.values = \
+ 0, 0.9889621791207942, 0.9991570495802089, 0.9998587947015799,
0.9999638901095526,\
+ 0.9999879879239209, 0.9999996670903828, 0.9999999999989165, 1.0, 1.0, 1.0,
1.0, 1.0,\
+ 1.0, 1.0, 1.0
+pmf.values = \
+ 0, 0.9889621791207942, 0.010194870459414753, 0.0007017451213709621,
0.00010509540797267703,\
+ 2.4097814368311538e-05, 2.4841606790096716e-07, 6.239929503294026e-14,
1.567399425290059e-20,\
+ 3.9371293491420155e-27, 9.889621791207982e-34, 2.484160679009682e-40,
6.23992950329405e-47,\
+ 1.5673994252900653e-53, 3.937129349142032e-60, 2.5377133927894665e-62
+sf.values = \
+ 1, 0.011037820879205844, 0.0008429504197910912, 0.0001412052984201291,
3.6109890447452063e-05,\
+ 1.2012076079140522e-05, 3.329096171741983e-07, 1.0834166368179869e-12,
2.7910991687218175e-18,\
+ 7.028619775335426e-24, 1.7659544437221128e-29, 4.435988791727557e-35,
1.1142728135917826e-40,\
+ 2.7989274667908607e-46, 7.030588103782268e-52, 9.731603579790481e-54
+
+moments.relative = 5e-14
+logpmf.relative = 2e-14
diff --git
a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zeta.3.properties
b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zeta.3.properties
new file mode 100644
index 00000000..9c755f73
--- /dev/null
+++
b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zeta.3.properties
@@ -0,0 +1,40 @@
+# 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.
+
+# Note: This parameterization allows reasonable quartiles to test the sampler
+parameters = 1.1
+mean = Infinity
+variance = Infinity
+lower = 1
+upper = 2147483647
+cdf.points = 0, 1, 2, 3, 4, 5, 10, 100, 1000, 10000
+# Reference values are from scipy.stats (1.17.1) zipf(1.1)
+# Note: The cdf is defined by sum(pdf) so x should be limited in size.
+cdf.values = \
+ 0. , 0.09447823411029742, 0.13855388881382025,\
+ 0.16677007985446593, 0.18733209983545154, 0.20341871833692377,\
+ 0.2532163286767671 , 0.40418015518943295, 0.5265108233845649,\
+ 0.6238772559832313
+pmf.values = \
+ 0.0000000000000000e+00, 9.4478234110297418e-02,\
+ 4.4075654703522822e-02, 2.8216191040645679e-02,\
+ 2.0562019980985623e-02, 1.6086618501472229e-02,\
+ 7.5046728920699938e-03, 5.9611735705411247e-04,\
+ 4.7351284791729730e-05, 3.7612462460541389e-06
+sf.values = \
+ 1. , 0.9055217658897026, 0.8614461111861798,\
+ 0.8332299201455341, 0.8126679001645485, 0.7965812816630762,\
+ 0.746783671323233 , 0.595819844810567 , 0.4734891766154351,\
+ 0.3761227440167687
diff --git
a/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zeta.4.properties
b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zeta.4.properties
new file mode 100644
index 00000000..d762e128
--- /dev/null
+++
b/commons-statistics-distribution/src/test/resources/org/apache/commons/statistics/distribution/test.zeta.4.properties
@@ -0,0 +1,48 @@
+# 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.
+
+# s -> 1 leads to most of the probability mass above 2^31
+parameters = 1.000001
+mean = Infinity
+variance = Infinity
+lower = 1
+upper = 2147483647
+cdf.points = 0, 1, 2, 3, 4, 5, 10, 100, 1000,\
+ 10000, 100000, 1000000, 10000000, 100000000, 1000000000, 2147483647
+
+# Reference values are from a custom Python implementation using mpmath based
on an extended
+# precision zeta function (see test.zeta2.properties).
+# Summation of the pdf to convergence for the cdf is not practical.
+
+cdf.values = \
+ 0, 9.99999422702329e-07, 1.4999987874802234e-06, 1.833331562177316e-06,
2.083331071279748e-06,\
+ 2.2833306339330765e-06, 2.928963870908051e-06, 5.187363969024896e-06,
7.485442750095092e-06,\
+ 9.787558043007328e-06, 1.2090072949574676e-05, 1.43926230531156e-05,
1.669517190477966e-05,\
+ 1.8997715859628422e-05, 2.1300254553179146e-05, 2.2064534741011326e-05
+pmf.values = \
+ 0, 9.99999422702329e-07, 4.999993647778944e-07, 3.333327746970927e-07,
2.499995091024323e-07,\
+ 1.9999956265332818e-07, 9.999971201212163e-08, 9.999948175454056e-09,
9.999925149748968e-10,\
+ 9.999902124096898e-11, 9.999879098497847e-12, 9.999856072951814e-13,
9.9998330474588e-14,\
+ 9.999810022018803e-15, 9.999786996631824e-16, 4.656510128866161e-16
+sf.values = \
+ 1, 0.9999990000005773, 0.9999985000012125, 0.9999981666684378,
0.9999979166689287,\
+ 0.9999977166693661, 0.9999970710361291, 0.999994812636031,
0.9999925145572499,\
+ 0.999990212441957, 0.9999879099270504, 0.9999856073769469,
0.9999833048280952,\
+ 0.9999810022841403, 0.9999786997454468, 0.999977935465259
+
+# cdf is defined using 1 - sf. This loses precision when sf -> 1.
+cdf.relative = 1e-10
+# cdf(x1) - cdf(x0) != probability(x0, x1)
+consistency.disable = true
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index d74f096f..3ef56471 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -53,6 +53,10 @@ If the output is not quite correct, check for invisible
trailing spaces!
</properties>
<body>
<release version="1.4" date="TBD" description="Adds new features (requires
Java 8).">
+ <action dev="aherbert" type="add" due-to="Alex Herbert"
issue="STATISTICS-101">
+ "ZetaDistribution": Add a zeta distribution. This is the large N limit
+ of the ZipfDistribution.
+ </action>
<action dev="aherbert" type="update" due-to="Alex Herbert"
issue="STATISTICS-100">
"ZipfDistribution": Use the Hurwitz zeta function when parameters are
suitable
for a significant performance improvment when the number of elements
is large.
