Merge 2fab0f2 into 85b31f7

commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/BetaDistribution.java

@@ -40,6 +40,13 @@ public class BetaDistribution extends AbstractContinuousDistribution {
      */
     public BetaDistribution(double alpha,
                             double beta) {
+        if (alpha <= 0) {
+            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, alpha);
+        }
+        if (beta <= 0) {
+            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, beta);
+        }
+
         this.alpha = alpha;
         this.beta = beta;
         z = LogGamma.value(alpha) + LogGamma.value(beta) - LogGamma.value(alpha + beta);
@@ -63,28 +70,48 @@ public double getBeta() {
         return beta;
     }
 
-    /** {@inheritDoc} */
+    /** {@inheritDoc}
+     *
+     * <p>The density is not defined when {@code x = 0, alpha < 1}, or {@code x = 1, beta < 1}.
+     * In this case the limit of infinity is returned.
+     */
     @Override
     public double density(double x) {
         return Math.exp(logDensity(x));
     }
 
-    /** {@inheritDoc} **/
+    /** {@inheritDoc}
+     *
+     * <p>The density is not defined when {@code x = 0, alpha < 1}, or {@code x = 1, beta < 1}.
+     * In this case the limit of infinity is returned.
+     */
     @Override
     public double logDensity(double x) {
         if (x < 0 ||
             x > 1) {
             return Double.NEGATIVE_INFINITY;
         } else if (x == 0) {
             if (alpha < 1) {
-                throw new DistributionException(DistributionException.TOO_SMALL,
-                                                alpha, 1.0);
+                // Distribution is not valid when x=0, alpha<1
+                // due to a divide by zero error.
+                // Do not raise an exception and return the limit.
+                return Double.POSITIVE_INFINITY;
+            }
+            // Special case of cancellation: x^(a-1) (1-x)^(b-1) / B(a, b)
+            if (alpha == 1) {
+                return -z;
             }
             return Double.NEGATIVE_INFINITY;
         } else if (x == 1) {
             if (beta < 1) {
-                throw new DistributionException(DistributionException.TOO_SMALL,
-                                                beta, 1.0);
+                // Distribution is not valid when x=1, beta<1
+                // due to a divide by zero error.
+                // Do not raise an exception and return the limit.
+                return Double.POSITIVE_INFINITY;
+            }
+            // Special case of cancellation: x^(a-1) (1-x)^(b-1) / B(a, b)
+            if (beta == 1) {
+                return -z;
             }
             return Double.NEGATIVE_INFINITY;
         } else {

commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/GumbelDistribution.java

@@ -27,10 +27,10 @@ public class GumbelDistribution extends AbstractContinuousDistribution {
     /** &pi;<sup>2</sup>/6. */
     private static final double PI_SQUARED_OVER_SIX = Math.PI * Math.PI / 6;
     /**
-     * <a href="http://mathworld.wolfram.com/Euler-MascheroniConstantApproximations.html">
+     * <a href="https://en.wikipedia.org/wiki/Euler%27s_constant">
      * Approximation of Euler's constant</a>.
      */
-    private static final double EULER = Math.PI / (2 * Math.E);
+    private static final double EULER = 0.57721566490153286060;
     /** Location parameter. */
     private final double mu;
     /** Scale parameter. */

commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/LevyDistribution.java

@@ -39,9 +39,15 @@ public class LevyDistribution extends AbstractContinuousDistribution {
      *
      * @param mu Location parameter.
      * @param c Scale parameter.
+     * @throws IllegalArgumentException if {@code c <= 0}.
      */
     public LevyDistribution(final double mu,
                             final double c) {
+        if (c <= 0) {
+            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE,
+                                            c);
+        }
+
         this.mu = mu;
         this.c = c;
         this.halfC = 0.5 * c;

commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/LogisticDistribution.java

@@ -30,8 +30,8 @@ public class LogisticDistribution extends AbstractContinuousDistribution {
     private final double mu;
     /** Scale parameter. */
     private final double scale;
-    /** Inverse of "scale". */
-    private final double oneOverScale;
+    /** Logarithm of "scale". */
+    private final double logScale;
 
     /**
      * Creates a distribution.
@@ -49,7 +49,7 @@ public LogisticDistribution(double mu,
 
         this.mu = mu;
         this.scale = scale;
-        this.oneOverScale = 1 / scale;
+        this.logScale = Math.log(scale);
     }
 
     /**
@@ -78,9 +78,9 @@ public double density(double x) {
             return 0;
         }
 
-        final double z = oneOverScale * (x - mu);
+        final double z = (x - mu) / scale;
         final double v = Math.exp(-z);
-        return oneOverScale * v / ((1 + v) * (1 + v));
+        return v / ((1 + v) * (1 + v)) / scale;
     }
 
     /** {@inheritDoc} */
@@ -91,22 +91,22 @@ public double logDensity(double x) {
             return Double.NEGATIVE_INFINITY;
         }
 
-        final double z = oneOverScale * (x - mu);
+        final double z = (x - mu) / scale;
         final double v = Math.exp(-z);
-        return -Math.log(scale) - z - 2 * Math.log(1 + v);
+        return -z - 2 * Math.log1p(v) - logScale;
     }
 
     /** {@inheritDoc} */
     @Override
     public double cumulativeProbability(double x) {
-        final double z = oneOverScale * (x - mu);
+        final double z = (x - mu) / scale;
         return 1 / (1 + Math.exp(-z));
     }
 
     /** {@inheritDoc} */
     @Override
     public double survivalProbability(double x) {
-        final double z = oneOverScale * (x - mu);
+        final double z = (x - mu) / scale;
         return 1 / (1 + Math.exp(z));
     }
 

commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/ZipfDistribution.java

@@ -18,18 +18,19 @@
 package org.apache.commons.statistics.distribution;
 
 import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.sampling.distribution.DiscreteUniformSampler;
 import org.apache.commons.rng.sampling.distribution.RejectionInversionZipfSampler;
 
 /**
  * Implementation of the <a href="https://en.wikipedia.org/wiki/Zipf's_law">Zipf distribution</a>.
  * <p>
  * <strong>Parameters:</strong>
  * For a random variable {@code X} whose values are distributed according to this
- * distribution, the probability mass function is given by
+ * distribution, the probability mass function is given by:
  * <pre>
  *   P(X = k) = H(N,s) * 1 / k^s    for {@code k = 1,2,...,N}.
  * </pre>
- * {@code H(N,s)} is the normalizing constant
+ * <p>{@code H(N,s)} is the normalizing constant
  * which corresponds to the generalized harmonic number of order N of s.
  * <ul>
  * <li>{@code N} is the number of elements</li>
@@ -43,6 +44,8 @@ public class ZipfDistribution extends AbstractDiscreteDistribution {
     private final double exponent;
     /** Cached value of the nth generalized harmonic. */
     private final double nthHarmonic;
+    /** Cached value of the log of the nth generalized harmonic. */
+    private final double logNthHarmonic;
 
     /**
      * Creates a distribution.
@@ -58,14 +61,15 @@ public ZipfDistribution(int numberOfElements,
             throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE,
                                             numberOfElements);
         }
-        if (exponent <= 0) {
-            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE,
+        if (exponent < 0) {
+            throw new DistributionException(DistributionException.NEGATIVE,
                                             exponent);
         }
 
         this.numberOfElements = numberOfElements;
         this.exponent = exponent;
         this.nthHarmonic = generalizedHarmonic(numberOfElements, exponent);
+        logNthHarmonic = Math.log(nthHarmonic);
     }
 
     /**
@@ -103,7 +107,7 @@ public double logProbability(int x) {
             return Double.NEGATIVE_INFINITY;
         }
 
-        return -Math.log(x) * exponent - Math.log(nthHarmonic);
+        return -Math.log(x) * exponent - logNthHarmonic;
     }
 
     /** {@inheritDoc} */
@@ -118,6 +122,29 @@ public double cumulativeProbability(final int x) {
         return generalizedHarmonic(x, exponent) / nthHarmonic;
     }
 
+    /** {@inheritDoc} */
+    @Override
+    public double survivalProbability(int x) {
+        if (x <= 0) {
+            return 1;
+        } else if (x >= numberOfElements) {
+            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, exponent);
+        // Compute generalizedHarmonic(x, exponent) and generalizedHarmonic(x+1, exponent)
+        final double hx = generalizedHarmonic(x, exponent);
+        final double hx1 = hx + Math.pow(x + 1, -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;
+    }
+
     /**
      * {@inheritDoc}
      *
@@ -133,7 +160,7 @@ public double getMean() {
         final int N = getNumberOfElements();
         final double s = getExponent();
 
-        final double Hs1 = generalizedHarmonic(N, s - 1);
+        final double Hs1 = generalizedHarmonicAscendingSum(N, s - 1);
 
         return Hs1 / nthHarmonic;
     }
@@ -154,8 +181,8 @@ public double getVariance() {
         final int N = getNumberOfElements();
         final double s = getExponent();
 
-        final double Hs2 = generalizedHarmonic(N, s - 2);
-        final double Hs1 = generalizedHarmonic(N, s - 1);
+        final double Hs2 = generalizedHarmonicAscendingSum(N, s - 2);
+        final double Hs1 = generalizedHarmonicAscendingSum(N, s - 1);
         final double Hs = nthHarmonic;
 
         return (Hs2 / Hs) - ((Hs1 * Hs1) / (Hs * Hs));
@@ -166,14 +193,42 @@ public double getVariance() {
      * <a href="http://mathworld.wolfram.com/HarmonicSeries.html">Harmonic
      * Series</a>.
      *
+     * <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 1)
      * @param m Exponent (special case {@code m = 1} is the harmonic series).
      * @return the n<sup>th</sup> generalized harmonic number.
      */
     private static double generalizedHarmonic(final int n, final double m) {
         double value = 0;
-        for (int k = n; k > 0; --k) {
-            value += 1 / Math.pow(k, m);
+        // Sum small to large
+        for (int k = n; k >= 1; k--) {
+            value += Math.pow(k, -m);
+        }
+        return value;
+    }
+
+    /**
+     * Calculates the Nth generalized harmonic number.
+     *
+     * <p>Checks the value of the {@code exponent} to arrange the terms to sum from from small to large.
+     *
+     * @param n Term in the series to calculate (must be larger than 1)
+     * @param m Exponent (special case {@code m = 1} is the harmonic series).
+     * @return the n<sup>th</sup> generalized harmonic number.
+     */
+    private static double generalizedHarmonicAscendingSum(final int n, final double m) {
+        double value = 0;
+        // Sum small to large
+        // If m < 0 then sum ascending, otherwise descending
+        if (m < 0) {
+            for (int k = 1; k <= n; k++) {
+                value += Math.pow(k, -m);
+            }
+        } else {
+            for (int k = n; k >= 1; k--) {
+                value += Math.pow(k, -m);
+            }
         }
         return value;
     }
@@ -217,7 +272,11 @@ public boolean isSupportConnected() {
     /** {@inheritDoc} */
     @Override
     public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
-        // Zipf distribution sampler.
-        return new RejectionInversionZipfSampler(rng, numberOfElements, exponent)::sample;
+        if (exponent > 0) {
+            // Zipf distribution sampler.
+            return RejectionInversionZipfSampler.of(rng, numberOfElements, exponent)::sample;
+        }
+        // Special case when exponent is zero then the sample is uniform
+        return DiscreteUniformSampler.of(rng, 1, numberOfElements)::sample;
     }
 }