commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject [4/6] commons-statistics git commit: Fixed file paths.
Date Mon, 22 Jan 2018 01:15:33 GMT
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/PoissonDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/PoissonDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/PoissonDistribution.java
deleted file mode 100644
index 225b8f1..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/PoissonDistribution.java
+++ /dev/null
@@ -1,238 +0,0 @@
-/*
- * 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.
- */
-package org.apache.commons.statistics.distribution;
-
-import org.apache.commons.numbers.gamma.RegularizedGamma;
-import org.apache.commons.rng.UniformRandomProvider;
-import org.apache.commons.rng.sampling.distribution.DiscreteSampler;
-import org.apache.commons.rng.sampling.distribution.PoissonSampler;
-
-/**
- * Implementation of the <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution</a>.
- */
-public class PoissonDistribution extends AbstractDiscreteDistribution {
-    /** ln(2 &pi;). */
-    private static final double LOG_TWO_PI = Math.log(2 * Math.PI);
-    /** Default maximum number of iterations. */
-    private static final int DEFAULT_MAX_ITERATIONS = 10000000;
-    /** Default convergence criterion. */
-    private static final double DEFAULT_EPSILON = 1e-12;
-    /** Distribution used to compute normal approximation. */
-    private final NormalDistribution normal;
-    /** Mean of the distribution. */
-    private final double mean;
-    /** Maximum number of iterations for cumulative probability. */
-    private final int maxIterations;
-    /** Convergence criterion for cumulative probability. */
-    private final double epsilon;
-
-    /**
-     * Creates a new Poisson distribution with specified mean.
-     *
-     * @param p the Poisson mean
-     * @throws IllegalArgumentException if {@code p <= 0}.
-     */
-    public PoissonDistribution(double p) {
-        this(p, DEFAULT_EPSILON, DEFAULT_MAX_ITERATIONS);
-    }
-
-    /**
-     * Creates a new Poisson distribution with specified mean, convergence
-     * criterion and maximum number of iterations.
-     *
-     * @param p Poisson mean.
-     * @param epsilon Convergence criterion for cumulative probabilities.
-     * @param maxIterations the maximum number of iterations for cumulative
-     * probabilities.
-     * @throws IllegalArgumentException if {@code p <= 0}.
-     */
-    public PoissonDistribution(double p,
-                               double epsilon,
-                               int maxIterations) {
-        if (p <= 0) {
-            throw new DistributionException(DistributionException.NEGATIVE, p);
-        }
-        mean = p;
-        this.epsilon = epsilon;
-        this.maxIterations = maxIterations;
-
-        normal = new NormalDistribution(p, Math.sqrt(p));
-    }
-
-    /**
-     * Creates a new Poisson distribution with the specified mean and
-     * convergence criterion.
-     *
-     * @param p Poisson mean.
-     * @param epsilon Convergence criterion for cumulative probabilities.
-     * @throws IllegalArgumentException if {@code p <= 0}.
-     */
-    public PoissonDistribution(double p,
-                               double epsilon) {
-        this(p, epsilon, DEFAULT_MAX_ITERATIONS);
-    }
-
-    /**
-     * Creates a new Poisson distribution with the specified mean and maximum
-     * number of iterations.
-     *
-     * @param p Poisson mean.
-     * @param maxIterations Maximum number of iterations for cumulative
-     * probabilities.
-     */
-    public PoissonDistribution(double p,
-                               int maxIterations) {
-        this(p, DEFAULT_EPSILON, maxIterations);
-    }
-
-    /**
-     * Get the mean for the distribution.
-     *
-     * @return the mean for the distribution.
-     */
-    public double getMean() {
-        return mean;
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double probability(int x) {
-        final double logProbability = logProbability(x);
-        return logProbability == Double.NEGATIVE_INFINITY ? 0 : Math.exp(logProbability);
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double logProbability(int x) {
-        double ret;
-        if (x < 0 || x == Integer.MAX_VALUE) {
-            ret = Double.NEGATIVE_INFINITY;
-        } else if (x == 0) {
-            ret = -mean;
-        } else {
-            ret = -SaddlePointExpansion.getStirlingError(x) -
-                  SaddlePointExpansion.getDeviancePart(x, mean) -
-                  0.5 * LOG_TWO_PI - 0.5 * Math.log(x);
-        }
-        return ret;
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double cumulativeProbability(int x) {
-        if (x < 0) {
-            return 0;
-        }
-        if (x == Integer.MAX_VALUE) {
-            return 1;
-        }
-        return RegularizedGamma.Q.value((double) x + 1, mean, epsilon,
-                                        maxIterations);
-    }
-
-    /**
-     * Calculates the Poisson distribution function using a normal
-     * approximation. The {@code N(mean, sqrt(mean))} distribution is used
-     * to approximate the Poisson distribution. The computation uses
-     * "half-correction" (evaluating the normal distribution function at
-     * {@code x + 0.5}).
-     *
-     * @param x Upper bound, inclusive.
-     * @return the distribution function value calculated using a normal
-     * approximation.
-     */
-    public double normalApproximateProbability(int x)  {
-        // Calculate the probability using half-correction.
-        return normal.cumulativeProbability(x + 0.5);
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * For mean parameter {@code p}, the mean is {@code p}.
-     */
-    @Override
-    public double getNumericalMean() {
-        return getMean();
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * For mean parameter {@code p}, the variance is {@code p}.
-     */
-    @Override
-    public double getNumericalVariance() {
-        return getMean();
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The lower bound of the support is always 0 no matter the mean parameter.
-     *
-     * @return lower bound of the support (always 0)
-     */
-    @Override
-    public int getSupportLowerBound() {
-        return 0;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The upper bound of the support is positive infinity,
-     * regardless of the parameter values. There is no integer infinity,
-     * so this method returns {@code Integer.MAX_VALUE}.
-     *
-     * @return upper bound of the support (always {@code Integer.MAX_VALUE} for
-     * positive infinity)
-     */
-    @Override
-    public int getSupportUpperBound() {
-        return Integer.MAX_VALUE;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The support of this distribution is connected.
-     *
-     * @return {@code true}
-     */
-    @Override
-    public boolean isSupportConnected() {
-        return true;
-    }
-
-    /**{@inheritDoc} */
-    @Override
-    public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
-        return new DiscreteDistribution.Sampler() {
-            /**
-             * Poisson distribution sampler.
-             */
-            private final DiscreteSampler sampler = new PoissonSampler(rng, mean);
-
-            /**{@inheritDoc} */
-            @Override
-            public int sample() {
-                return sampler.sample();
-            }
-        };
-    }
-}

http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/SaddlePointExpansion.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/SaddlePointExpansion.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/SaddlePointExpansion.java
deleted file mode 100644
index 7bb847a..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/SaddlePointExpansion.java
+++ /dev/null
@@ -1,191 +0,0 @@
-/*
- * 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.
- */
-package org.apache.commons.statistics.distribution;
-
-import org.apache.commons.numbers.gamma.LogGamma;
-
-/**
- * Utility class used by various distributions to accurately compute their
- * respective probability mass functions. The implementation for this class is
- * based on the Catherine Loader's
- * <a href="http://www.herine.net/stat/software/dbinom.html">dbinom</a> routines.
- *
- * This class is not intended to be called directly.
- *
- * @since 1.0
- */
-final class SaddlePointExpansion {
-    /** 2 &pi; */
-    private static final double TWO_PI = 2 * Math.PI;
-    /** 1/2 * log(2 &pi;). */
-    private static final double HALF_LOG_TWO_PI = 0.5 * Math.log(TWO_PI);
-
-    /** exact Stirling expansion error for certain values. */
-    private static final double[] EXACT_STIRLING_ERRORS = { 0.0, /* 0.0 */
-    0.1534264097200273452913848, /* 0.5 */
-    0.0810614667953272582196702, /* 1.0 */
-    0.0548141210519176538961390, /* 1.5 */
-    0.0413406959554092940938221, /* 2.0 */
-    0.03316287351993628748511048, /* 2.5 */
-    0.02767792568499833914878929, /* 3.0 */
-    0.02374616365629749597132920, /* 3.5 */
-    0.02079067210376509311152277, /* 4.0 */
-    0.01848845053267318523077934, /* 4.5 */
-    0.01664469118982119216319487, /* 5.0 */
-    0.01513497322191737887351255, /* 5.5 */
-    0.01387612882307074799874573, /* 6.0 */
-    0.01281046524292022692424986, /* 6.5 */
-    0.01189670994589177009505572, /* 7.0 */
-    0.01110455975820691732662991, /* 7.5 */
-    0.010411265261972096497478567, /* 8.0 */
-    0.009799416126158803298389475, /* 8.5 */
-    0.009255462182712732917728637, /* 9.0 */
-    0.008768700134139385462952823, /* 9.5 */
-    0.008330563433362871256469318, /* 10.0 */
-    0.007934114564314020547248100, /* 10.5 */
-    0.007573675487951840794972024, /* 11.0 */
-    0.007244554301320383179543912, /* 11.5 */
-    0.006942840107209529865664152, /* 12.0 */
-    0.006665247032707682442354394, /* 12.5 */
-    0.006408994188004207068439631, /* 13.0 */
-    0.006171712263039457647532867, /* 13.5 */
-    0.005951370112758847735624416, /* 14.0 */
-    0.005746216513010115682023589, /* 14.5 */
-    0.005554733551962801371038690 /* 15.0 */
-    };
-
-    /**
-     * Forbid construction.
-     */
-    private SaddlePointExpansion() {}
-
-    /**
-     * Compute the error of Stirling's series at the given value.
-     * <p>
-     * References:
-     * <ol>
-     * <li>Eric W. Weisstein. "Stirling's Series." From MathWorld--A Wolfram Web
-     * Resource. <a target="_blank"
-     * href="http://mathworld.wolfram.com/StirlingsSeries.html">
-     * http://mathworld.wolfram.com/StirlingsSeries.html</a></li>
-     * </ol>
-     * </p>
-     *
-     * @param z the value.
-     * @return the Striling's series error.
-     */
-    static double getStirlingError(double z) {
-        double ret;
-        if (z < 15.0) {
-            double z2 = 2.0 * z;
-            if (Math.floor(z2) == z2) {
-                ret = EXACT_STIRLING_ERRORS[(int) z2];
-            } else {
-                ret = LogGamma.value(z + 1.0) - (z + 0.5) * Math.log(z) +
-                      z - HALF_LOG_TWO_PI;
-            }
-        } else {
-            double z2 = z * z;
-            ret = (0.083333333333333333333 -
-                    (0.00277777777777777777778 -
-                            (0.00079365079365079365079365 -
-                                    (0.000595238095238095238095238 -
-                                            0.0008417508417508417508417508 /
-                                            z2) / z2) / z2) / z2) / z;
-        }
-        return ret;
-    }
-
-    /**
-     * A part of the deviance portion of the saddle point approximation.
-     * <p>
-     * References:
-     * <ol>
-     * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial
-     * Probabilities.". <a target="_blank"
-     * href="http://www.herine.net/stat/papers/dbinom.pdf">
-     * http://www.herine.net/stat/papers/dbinom.pdf</a></li>
-     * </ol>
-     * </p>
-     *
-     * @param x the x value.
-     * @param mu the average.
-     * @return a part of the deviance.
-     */
-    static double getDeviancePart(double x, double mu) {
-        double ret;
-        if (Math.abs(x - mu) < 0.1 * (x + mu)) {
-            double d = x - mu;
-            double v = d / (x + mu);
-            double s1 = v * d;
-            double s = Double.NaN;
-            double ej = 2.0 * x * v;
-            v *= v;
-            int j = 1;
-            while (s1 != s) {
-                s = s1;
-                ej *= v;
-                s1 = s + ej / ((j * 2) + 1);
-                ++j;
-            }
-            ret = s1;
-        } else {
-            if (x == 0) {
-                return mu;
-            }
-            ret = x * Math.log(x / mu) + mu - x;
-        }
-        return ret;
-    }
-
-    /**
-     * Compute the logarithm of the PMF for a binomial distribution
-     * using the saddle point expansion.
-     *
-     * @param x the value at which the probability is evaluated.
-     * @param n the number of trials.
-     * @param p the probability of success.
-     * @param q the probability of failure (1 - p).
-     * @return log(p(x)).
-     */
-    static double logBinomialProbability(int x, int n, double p, double q) {
-        double ret;
-        if (x == 0) {
-            if (p < 0.1) {
-                ret = -getDeviancePart(n, n * q) - n * p;
-            } else {
-                if (n == 0) {
-                    return 0;
-                }
-                ret = n * Math.log(q);
-            }
-        } else if (x == n) {
-            if (q < 0.1) {
-                ret = -getDeviancePart(n, n * p) - n * q;
-            } else {
-                ret = n * Math.log(p);
-            }
-        } else {
-            ret = getStirlingError(n) - getStirlingError(x) -
-                  getStirlingError(n - x) - getDeviancePart(x, n * p) -
-                  getDeviancePart(n - x, n * q);
-            final double f = (TWO_PI * x * (n - x)) / n;
-            ret = -0.5 * Math.log(f) + ret;
-        }
-        return ret;
-    }
-}

http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TDistribution.java
deleted file mode 100644
index ef34c1f..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TDistribution.java
+++ /dev/null
@@ -1,180 +0,0 @@
-/*
- * 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.
- */
-package org.apache.commons.statistics.distribution;
-
-import org.apache.commons.numbers.gamma.RegularizedBeta;
-import org.apache.commons.numbers.gamma.LogGamma;
-
-/**
- * Implementation of <a href='http://en.wikipedia.org/wiki/Student&apos;s_t-distribution'>Student's t-distribution</a>.
- */
-public class TDistribution extends AbstractContinuousDistribution {
-    /** The degrees of freedom. */
-    private final double degreesOfFreedom;
-    /** degreesOfFreedom / 2 */
-    private final double dofOver2;
-    /** Cached value. */
-    private final double factor;
-
-    /**
-     * Creates a distribution.
-     *
-     * @param degreesOfFreedom Degrees of freedom.
-     * @throws IllegalArgumentException if {@code degreesOfFreedom <= 0}
-     */
-    public TDistribution(double degreesOfFreedom) {
-        if (degreesOfFreedom <= 0) {
-            throw new DistributionException(DistributionException.NEGATIVE,
-                                            degreesOfFreedom);
-        }
-        this.degreesOfFreedom = degreesOfFreedom;
-
-        dofOver2 = 0.5 * degreesOfFreedom;
-        factor = LogGamma.value(dofOver2 + 0.5) -
-                 0.5 * (Math.log(Math.PI) + Math.log(degreesOfFreedom)) -
-                 LogGamma.value(dofOver2);
-    }
-
-    /**
-     * Access the degrees of freedom.
-     *
-     * @return the degrees of freedom.
-     */
-    public double getDegreesOfFreedom() {
-        return degreesOfFreedom;
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double density(double x) {
-        return Math.exp(logDensity(x));
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double logDensity(double x) {
-        final double nPlus1Over2 = dofOver2 + 0.5;
-        return factor - nPlus1Over2 * Math.log(1 + x * x / degreesOfFreedom);
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double cumulativeProbability(double x) {
-        double ret;
-        if (x == 0) {
-            ret = 0.5;
-        } else {
-            final double t =
-                RegularizedBeta.value(degreesOfFreedom / (degreesOfFreedom + (x * x)),
-                                      dofOver2,
-                                      0.5);
-            if (x < 0) {
-                ret = 0.5 * t;
-            } else {
-                ret = 1 - 0.5 * t;
-            }
-        }
-
-        return ret;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * For degrees of freedom parameter {@code df}, the mean is
-     * <ul>
-     *  <li>if {@code df > 1} then {@code 0},</li>
-     * <li>else undefined ({@code Double.NaN}).</li>
-     * </ul>
-     */
-    @Override
-    public double getNumericalMean() {
-        final double df = getDegreesOfFreedom();
-
-        if (df > 1) {
-            return 0;
-        }
-
-        return Double.NaN;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * For degrees of freedom parameter {@code df}, the variance is
-     * <ul>
-     *  <li>if {@code df > 2} then {@code df / (df - 2)},</li>
-     *  <li>if {@code 1 < df <= 2} then positive infinity
-     *  ({@code Double.POSITIVE_INFINITY}),</li>
-     *  <li>else undefined ({@code Double.NaN}).</li>
-     * </ul>
-     */
-    @Override
-    public double getNumericalVariance() {
-        final double df = getDegreesOfFreedom();
-
-        if (df > 2) {
-            return df / (df - 2);
-        }
-
-        if (df > 1 && df <= 2) {
-            return Double.POSITIVE_INFINITY;
-        }
-
-        return Double.NaN;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The lower bound of the support is always negative infinity no matter the
-     * parameters.
-     *
-     * @return lower bound of the support (always
-     * {@code Double.NEGATIVE_INFINITY})
-     */
-    @Override
-    public double getSupportLowerBound() {
-        return Double.NEGATIVE_INFINITY;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The upper bound of the support is always positive infinity no matter the
-     * parameters.
-     *
-     * @return upper bound of the support (always
-     * {@code Double.POSITIVE_INFINITY})
-     */
-    @Override
-    public double getSupportUpperBound() {
-        return Double.POSITIVE_INFINITY;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The support of this distribution is connected.
-     *
-     * @return {@code true}
-     */
-    @Override
-    public boolean isSupportConnected() {
-        return true;
-    }
-}

http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TriangularDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TriangularDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TriangularDistribution.java
deleted file mode 100644
index 6a94b62..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TriangularDistribution.java
+++ /dev/null
@@ -1,222 +0,0 @@
-/*
- * 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.
- */
-
-package org.apache.commons.statistics.distribution;
-
-/**
- * Implementation of the triangular real distribution.
- *
- * @see <a href="http://en.wikipedia.org/wiki/Triangular_distribution">
- * Triangular distribution (Wikipedia)</a>
- *
- * @since 3.0
- */
-public class TriangularDistribution extends AbstractContinuousDistribution {
-    /** Serializable version identifier. */
-    private static final long serialVersionUID = 20160311L;
-    /** Lower limit of this distribution (inclusive). */
-    private final double a;
-    /** Upper limit of this distribution (inclusive). */
-    private final double b;
-    /** Mode of this distribution. */
-    private final double c;
-
-    /**
-     * Creates a distribution.
-     *
-     * @param a Lower limit of this distribution (inclusive).
-     * @param b Upper limit of this distribution (inclusive).
-     * @param c Mode of this distribution.
-     * @throws IllegalArgumentException if {@code a >= b}, if {@code c > b}
-     * or if {@code c < a}.
-     */
-    public TriangularDistribution(double a,
-                                  double c,
-                                  double b) {
-        if (a >= b) {
-            throw new DistributionException(DistributionException.TOO_LARGE,
-                                            a, b);
-        }
-        if (c < a) {
-            throw new DistributionException(DistributionException.TOO_SMALL,
-                                            c, a);
-        }
-        if (c > b) {
-            throw new DistributionException(DistributionException.TOO_LARGE,
-                                            c, b);
-        }
-
-        this.a = a;
-        this.c = c;
-        this.b = b;
-    }
-
-    /**
-     * Gets the mode.
-     *
-     * @return the mode of the distribution.
-     */
-    public double getMode() {
-        return c;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * For lower limit {@code a}, upper limit {@code b} and mode {@code c}, the
-     * PDF is given by
-     * <ul>
-     * <li>{@code 2 * (x - a) / [(b - a) * (c - a)]} if {@code a <= x < c},</li>
-     * <li>{@code 2 / (b - a)} if {@code x = c},</li>
-     * <li>{@code 2 * (b - x) / [(b - a) * (b - c)]} if {@code c < x <= b},</li>
-     * <li>{@code 0} otherwise.
-     * </ul>
-     */
-    @Override
-    public double density(double x) {
-        if (x < a) {
-            return 0;
-        }
-        if (a <= x && x < c) {
-            double divident = 2 * (x - a);
-            double divisor = (b - a) * (c - a);
-            return divident / divisor;
-        }
-        if (x == c) {
-            return 2 / (b - a);
-        }
-        if (c < x && x <= b) {
-            double divident = 2 * (b - x);
-            double divisor = (b - a) * (b - c);
-            return divident / divisor;
-        }
-        return 0;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * For lower limit {@code a}, upper limit {@code b} and mode {@code c}, the
-     * CDF is given by
-     * <ul>
-     * <li>{@code 0} if {@code x < a},</li>
-     * <li>{@code (x - a)^2 / [(b - a) * (c - a)]} if {@code a <= x < c},</li>
-     * <li>{@code (c - a) / (b - a)} if {@code x = c},</li>
-     * <li>{@code 1 - (b - x)^2 / [(b - a) * (b - c)]} if {@code c < x <= b},</li>
-     * <li>{@code 1} if {@code x > b}.</li>
-     * </ul>
-     */
-    @Override
-    public double cumulativeProbability(double x)  {
-        if (x < a) {
-            return 0;
-        }
-        if (a <= x && x < c) {
-            double divident = (x - a) * (x - a);
-            double divisor = (b - a) * (c - a);
-            return divident / divisor;
-        }
-        if (x == c) {
-            return (c - a) / (b - a);
-        }
-        if (c < x && x <= b) {
-            double divident = (b - x) * (b - x);
-            double divisor = (b - a) * (b - c);
-            return 1 - (divident / divisor);
-        }
-        return 1;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * For lower limit {@code a}, upper limit {@code b}, and mode {@code c},
-     * the mean is {@code (a + b + c) / 3}.
-     */
-    @Override
-    public double getNumericalMean() {
-        return (a + b + c) / 3;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * For lower limit {@code a}, upper limit {@code b}, and mode {@code c},
-     * the variance is {@code (a^2 + b^2 + c^2 - a * b - a * c - b * c) / 18}.
-     */
-    @Override
-    public double getNumericalVariance() {
-        return (a * a + b * b + c * c - a * b - a * c - b * c) / 18;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The lower bound of the support is equal to the lower limit parameter
-     * {@code a} of the distribution.
-     *
-     * @return lower bound of the support
-     */
-    @Override
-    public double getSupportLowerBound() {
-        return a;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The upper bound of the support is equal to the upper limit parameter
-     * {@code b} of the distribution.
-     *
-     * @return upper bound of the support
-     */
-    @Override
-    public double getSupportUpperBound() {
-        return b;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The support of this distribution is connected.
-     *
-     * @return {@code true}
-     */
-    @Override
-    public boolean isSupportConnected() {
-        return true;
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double inverseCumulativeProbability(double p) {
-        if (p < 0 ||
-            p > 1) {
-            throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
-        }
-        if (p == 0) {
-            return a;
-        }
-        if (p == 1) {
-            return b;
-        }
-        if (p < (c - a) / (b - a)) {
-            return a + Math.sqrt(p * (b - a) * (c - a));
-        }
-        return b - Math.sqrt((1 - p) * (b - a) * (b - c));
-    }
-}

http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformContinuousDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformContinuousDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformContinuousDistribution.java
deleted file mode 100644
index 2bb9a0b..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformContinuousDistribution.java
+++ /dev/null
@@ -1,168 +0,0 @@
-/*
- * 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.
- */
-
-package org.apache.commons.statistics.distribution;
-
-import org.apache.commons.rng.UniformRandomProvider;
-import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
-import org.apache.commons.rng.sampling.distribution.ContinuousUniformSampler;
-
-/**
- * Implementation of the <a href="http://en.wikipedia.org/wiki/Uniform_distribution_(continuous)">uniform distribution</a>.
- */
-public class UniformContinuousDistribution extends AbstractContinuousDistribution {
-    private final double lower;
-    /** Upper bound of this distribution (exclusive). */
-    private final double upper;
-    
-    /**
-     * Create a standard uniform real distribution with lower bound (inclusive)
-     * equal to zero and upper bound (exclusive) equal to one.
-     */
-    public UniformContinuousDistribution() {
-        this(0, 1);
-    }
-
-    /**
-     * Creates a uniform distribution.
-     *
-     * @param lower Lower bound of this distribution (inclusive).
-     * @param upper Upper bound of this distribution (exclusive).
-     * @throws IllegalArgumentException if {@code lower >= upper}.
-     */
-    public UniformContinuousDistribution(double lower,
-                                         double upper) {
-        if (lower >= upper) {
-            throw new DistributionException(DistributionException.TOO_LARGE,
-                                            lower, upper);
-        }
-
-        this.lower = lower;
-        this.upper = upper;
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double density(double x) {
-        if (x < lower ||
-            x > upper) {
-            return 0;
-        }
-        return 1 / (upper - lower);
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double cumulativeProbability(double x)  {
-        if (x <= lower) {
-            return 0;
-        }
-        if (x >= upper) {
-            return 1;
-        }
-        return (x - lower) / (upper - lower);
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double inverseCumulativeProbability(final double p) {
-        if (p < 0 ||
-            p > 1) {
-            throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
-        }
-        return p * (upper - lower) + lower;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * For lower bound {@code lower} and upper bound {@code upper}, the mean is
-     * {@code 0.5 * (lower + upper)}.
-     */
-    @Override
-    public double getNumericalMean() {
-        return 0.5 * (lower + upper);
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * For lower bound {@code lower} and upper bound {@code upper}, the
-     * variance is {@code (upper - lower)^2 / 12}.
-     */
-    @Override
-    public double getNumericalVariance() {
-        double ul = upper - lower;
-        return ul * ul / 12;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The lower bound of the support is equal to the lower bound parameter
-     * of the distribution.
-     *
-     * @return lower bound of the support
-     */
-    @Override
-    public double getSupportLowerBound() {
-        return lower;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The upper bound of the support is equal to the upper bound parameter
-     * of the distribution.
-     *
-     * @return upper bound of the support
-     */
-    @Override
-    public double getSupportUpperBound() {
-        return upper;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The support of this distribution is connected.
-     *
-     * @return {@code true}
-     */
-    @Override
-    public boolean isSupportConnected() {
-        return true;
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
-        return new ContinuousDistribution.Sampler() {
-            /**
-             * Uniform distribution sampler.
-             */
-            private final ContinuousSampler sampler =
-                new ContinuousUniformSampler(rng, lower, upper);
-
-            /**{@inheritDoc} */
-            @Override
-            public double sample() {
-                return sampler.sample();
-            }
-        };
-    }
-}

http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformDiscreteDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformDiscreteDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformDiscreteDistribution.java
deleted file mode 100644
index 41df2bd..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformDiscreteDistribution.java
+++ /dev/null
@@ -1,159 +0,0 @@
-/*
- * 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.
- */
-
-package org.apache.commons.statistics.distribution;
-
-import org.apache.commons.rng.UniformRandomProvider;
-import org.apache.commons.rng.sampling.distribution.DiscreteSampler;
-import org.apache.commons.rng.sampling.distribution.DiscreteUniformSampler;
-
-/**
- * Implementation of the <a href="http://en.wikipedia.org/wiki/Uniform_distribution_(discrete)">
- * uniform integer distribution</a>.
- */
-public class UniformDiscreteDistribution extends AbstractDiscreteDistribution {
-    /** 1 / 12 **/
-    private static final double ONE_TWELFTH = 1 / 12d;
-    /** Lower bound (inclusive) of this distribution. */
-    private final int lower;
-    /** Upper bound (inclusive) of this distribution. */
-    private final int upper;
-    /** "upper" + "lower" (to avoid overflow). */
-    private final double upperPlusLower;
-    /** "upper" - "lower" (to avoid overflow). */
-    private final double upperMinusLower;
-
-    /**
-     * Creates a new uniform integer distribution using the given lower and
-     * upper bounds (both inclusive).
-     *
-     * @param lower Lower bound (inclusive) of this distribution.
-     * @param upper Upper bound (inclusive) of this distribution.
-     * @throws IllegalArgumentException if {@code lower > upper}.
-     */
-    public UniformDiscreteDistribution(int lower,
-                                       int upper) {
-        if (lower > upper) {
-            throw new DistributionException(DistributionException.TOO_LARGE,
-                                            lower, upper);
-        }
-        this.lower = lower;
-        this.upper = upper;
-        upperPlusLower = (double) upper + (double) lower;
-        upperMinusLower = (double) upper - (double) lower;
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double probability(int x) {
-        if (x < lower || x > upper) {
-            return 0;
-        }
-        return 1 / (upperMinusLower + 1);
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double cumulativeProbability(int x) {
-        if (x < lower) {
-            return 0;
-        }
-        if (x > upper) {
-            return 1;
-        }
-        return (x - lower + 1) / (upperMinusLower + 1);
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * For lower bound {@code lower} and upper bound {@code upper}, the mean is
-     * {@code 0.5 * (lower + upper)}.
-     */
-    @Override
-    public double getNumericalMean() {
-        return 0.5 * upperPlusLower;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * For lower bound {@code lower} and upper bound {@code upper}, and
-     * {@code n = upper - lower + 1}, the variance is {@code (n^2 - 1) / 12}.
-     */
-    @Override
-    public double getNumericalVariance() {
-        double n = upperMinusLower + 1;
-        return ONE_TWELFTH * (n * n - 1);
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The lower bound of the support is equal to the lower bound parameter
-     * of the distribution.
-     *
-     * @return lower bound of the support
-     */
-    @Override
-    public int getSupportLowerBound() {
-        return lower;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The upper bound of the support is equal to the upper bound parameter
-     * of the distribution.
-     *
-     * @return upper bound of the support
-     */
-    @Override
-    public int getSupportUpperBound() {
-        return upper;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The support of this distribution is connected.
-     *
-     * @return {@code true}
-     */
-    @Override
-    public boolean isSupportConnected() {
-        return true;
-    }
-
-    /**{@inheritDoc} */
-    @Override
-    public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
-        return new DiscreteDistribution.Sampler() {
-            /**
-             * Discrete uniform distribution sampler.
-             */
-            private final DiscreteSampler sampler =
-                new DiscreteUniformSampler(rng, lower, upper);
-
-            /**{@inheritDoc} */
-            @Override
-            public int sample() {
-                return sampler.sample();
-            }
-        };
-    }
-}

http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/WeibullDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/WeibullDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/WeibullDistribution.java
deleted file mode 100644
index 0e396d8..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/WeibullDistribution.java
+++ /dev/null
@@ -1,220 +0,0 @@
-/*
- * 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.
- */
-
-package org.apache.commons.statistics.distribution;
-
-import org.apache.commons.numbers.gamma.LogGamma;
-
-/**
- * Implementation of the Weibull distribution. This implementation uses the
- * two parameter form of the distribution defined by
- * <a href="http://mathworld.wolfram.com/WeibullDistribution.html">
- * Weibull Distribution</a>, equations (1) and (2).
- *
- * @see <a href="http://en.wikipedia.org/wiki/Weibull_distribution">Weibull distribution (Wikipedia)</a>
- * @see <a href="http://mathworld.wolfram.com/WeibullDistribution.html">Weibull distribution (MathWorld)</a>
- *
- * @since 1.1
- */
-public class WeibullDistribution extends AbstractContinuousDistribution {
-    /** The shape parameter. */
-    private final double shape;
-    /** The scale parameter. */
-    private final double scale;
-
-    /**
-     * Creates a distribution.
-     *
-     * @param alpha Shape parameter.
-     * @param beta Scale parameter.
-     * @throws IllegalArgumentException if {@code alpha <= 0} or {@code beta <= 0}.
-     */
-    public WeibullDistribution(double alpha,
-                               double beta) {
-        if (alpha <= 0) {
-            throw new DistributionException(DistributionException.NEGATIVE,
-                                            alpha);
-        }
-        if (beta <= 0) {
-            throw new DistributionException(DistributionException.NEGATIVE,
-                                            beta);
-        }
-        scale = beta;
-        shape = alpha;
-    }
-
-    /**
-     * Access the shape parameter, {@code alpha}.
-     *
-     * @return the shape parameter, {@code alpha}.
-     */
-    public double getShape() {
-        return shape;
-    }
-
-    /**
-     * Access the scale parameter, {@code beta}.
-     *
-     * @return the scale parameter, {@code beta}.
-     */
-    public double getScale() {
-        return scale;
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double density(double x) {
-        if (x < 0) {
-            return 0;
-        }
-
-        final double xscale = x / scale;
-        final double xscalepow = Math.pow(xscale, shape - 1);
-
-        /*
-         * Math.pow(x / scale, shape) =
-         * Math.pow(xscale, shape) =
-         * Math.pow(xscale, shape - 1) * xscale
-         */
-        final double xscalepowshape = xscalepow * xscale;
-
-        return (shape / scale) * xscalepow * Math.exp(-xscalepowshape);
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double logDensity(double x) {
-        if (x < 0) {
-            return Double.NEGATIVE_INFINITY;
-        }
-
-        final double xscale = x / scale;
-        final double logxscalepow = Math.log(xscale) * (shape - 1);
-
-        /*
-         * Math.pow(x / scale, shape) =
-         * Math.pow(xscale, shape) =
-         * Math.pow(xscale, shape - 1) * xscale
-         */
-        final double xscalepowshape = Math.exp(logxscalepow) * xscale;
-
-        return Math.log(shape / scale) + logxscalepow - xscalepowshape;
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double cumulativeProbability(double x) {
-        double ret;
-        if (x <= 0.0) {
-            ret = 0.0;
-        } else {
-            ret = 1.0 - Math.exp(-Math.pow(x / scale, shape));
-        }
-        return ret;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * Returns {@code 0} when {@code p == 0} and
-     * {@code Double.POSITIVE_INFINITY} when {@code p == 1}.
-     */
-    @Override
-    public double inverseCumulativeProbability(double p) {
-        double ret;
-        if (p < 0 ||
-            p > 1) {
-            throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
-        } else if (p == 0) {
-            ret = 0.0;
-        } else  if (p == 1) {
-            ret = Double.POSITIVE_INFINITY;
-        } else {
-            ret = scale * Math.pow(-Math.log1p(-p), 1.0 / shape);
-        }
-        return ret;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The mean is {@code scale * Gamma(1 + (1 / shape))}, where {@code Gamma()}
-     * is the Gamma-function.
-     */
-    @Override
-    public double getNumericalMean() {
-        final double sh = getShape();
-        final double sc = getScale();
-
-        return sc * Math.exp(LogGamma.value(1 + (1 / sh)));
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The variance is {@code scale^2 * Gamma(1 + (2 / shape)) - mean^2}
-     * where {@code Gamma()} is the Gamma-function.
-     */
-    @Override
-    public double getNumericalVariance() {
-        final double sh = getShape();
-        final double sc = getScale();
-        final double mn = getNumericalMean();
-
-        return (sc * sc) * Math.exp(LogGamma.value(1 + (2 / sh))) -
-               (mn * mn);
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The lower bound of the support is always 0 no matter the parameters.
-     *
-     * @return lower bound of the support (always 0)
-     */
-    @Override
-    public double getSupportLowerBound() {
-        return 0;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The upper bound of the support is always positive infinity
-     * no matter the parameters.
-     *
-     * @return upper bound of the support (always
-     * {@code Double.POSITIVE_INFINITY})
-     */
-    @Override
-    public double getSupportUpperBound() {
-        return Double.POSITIVE_INFINITY;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The support of this distribution is connected.
-     *
-     * @return {@code true}
-     */
-    @Override
-    public boolean isSupportConnected() {
-        return true;
-    }
-}
-

http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ZipfDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ZipfDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ZipfDistribution.java
deleted file mode 100644
index f0d54a1..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ZipfDistribution.java
+++ /dev/null
@@ -1,236 +0,0 @@
-/*
- * 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.
- */
-
-package org.apache.commons.statistics.distribution;
-
-import org.apache.commons.rng.UniformRandomProvider;
-import org.apache.commons.rng.sampling.distribution.DiscreteSampler;
-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
- * <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
- * which corresponds to the generalized harmonic number of order N of s.
- * <ul>
- * <li>{@code N} is the number of elements</li>
- * <li>{@code s} is the exponent</li>
- * </ul>
- */
-public class ZipfDistribution extends AbstractDiscreteDistribution {
-    /** Number of elements. */
-    private final int numberOfElements;
-    /** Exponent parameter of the distribution. */
-    private final double exponent;
-    /** Cached values of the nth generalized harmonic. */
-    private final double nthHarmonic;
-
-    /**
-     * Creates a distribution.
-     *
-     * @param numberOfElements Number of elements.
-     * @param exponent Exponent.
-     * @exception IllegalArgumentException if {@code numberOfElements <= 0}
-     * or {@code exponent <= 0}.
-     */
-    public ZipfDistribution(int numberOfElements,
-                            double exponent) {
-        if (numberOfElements <= 0) {
-            throw new DistributionException(DistributionException.NEGATIVE,
-                                            numberOfElements);
-        }
-        if (exponent <= 0) {
-            throw new DistributionException(DistributionException.NEGATIVE,
-                                            exponent);
-        }
-
-        this.numberOfElements = numberOfElements;
-        this.exponent = exponent;
-        this.nthHarmonic = generalizedHarmonic(numberOfElements, exponent);
-    }
-
-    /**
-     * Get the number of elements (e.g. corpus size) for the distribution.
-     *
-     * @return the number of elements
-     */
-    public int getNumberOfElements() {
-        return numberOfElements;
-    }
-
-    /**
-     * Get the exponent characterizing the distribution.
-     *
-     * @return the exponent
-     */
-    public double getExponent() {
-        return exponent;
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double probability(final int x) {
-        if (x <= 0 || x > numberOfElements) {
-            return 0;
-        }
-
-        return (1 / Math.pow(x, exponent)) / nthHarmonic;
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double logProbability(int x) {
-        if (x <= 0 || x > numberOfElements) {
-            return Double.NEGATIVE_INFINITY;
-        }
-
-        return -Math.log(x) * exponent - Math.log(nthHarmonic);
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public double cumulativeProbability(final int x) {
-        if (x <= 0) {
-            return 0;
-        } else if (x >= numberOfElements) {
-            return 1;
-        }
-
-        return generalizedHarmonic(x, exponent) / nthHarmonic;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * For number of elements {@code N} and exponent {@code s}, the mean is
-     * {@code Hs1 / Hs}, where
-     * <ul>
-     *  <li>{@code Hs1 = generalizedHarmonic(N, s - 1)},</li>
-     *  <li>{@code Hs = generalizedHarmonic(N, s)}.</li>
-     * </ul>
-     */
-    @Override
-    public double getNumericalMean() {
-        final int N = getNumberOfElements();
-        final double s = getExponent();
-
-        final double Hs1 = generalizedHarmonic(N, s - 1);
-        final double Hs = nthHarmonic;
-
-        return Hs1 / Hs;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * For number of elements {@code N} and exponent {@code s}, the mean is
-     * {@code (Hs2 / Hs) - (Hs1^2 / Hs^2)}, where
-     * <ul>
-     *  <li>{@code Hs2 = generalizedHarmonic(N, s - 2)},</li>
-     *  <li>{@code Hs1 = generalizedHarmonic(N, s - 1)},</li>
-     *  <li>{@code Hs = generalizedHarmonic(N, s)}.</li>
-     * </ul>
-     */
-    @Override
-    public double getNumericalVariance() {
-        final int N = getNumberOfElements();
-        final double s = getExponent();
-
-        final double Hs2 = generalizedHarmonic(N, s - 2);
-        final double Hs1 = generalizedHarmonic(N, s - 1);
-        final double Hs = nthHarmonic;
-
-        return (Hs2 / Hs) - ((Hs1 * Hs1) / (Hs * Hs));
-    }
-
-    /**
-     * Calculates the Nth generalized harmonic number. See
-     * <a href="http://mathworld.wolfram.com/HarmonicSeries.html">Harmonic
-     * Series</a>.
-     *
-     * @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 double generalizedHarmonic(final int n, final double m) {
-        double value = 0;
-        for (int k = n; k > 0; --k) {
-            value += 1 / Math.pow(k, m);
-        }
-        return value;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The lower bound of the support is always 1 no matter the parameters.
-     *
-     * @return lower bound of the support (always 1)
-     */
-    @Override
-    public int getSupportLowerBound() {
-        return 1;
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The upper bound of the support is the number of elements.
-     *
-     * @return upper bound of the support
-     */
-    @Override
-    public int getSupportUpperBound() {
-        return getNumberOfElements();
-    }
-
-    /**
-     * {@inheritDoc}
-     *
-     * The support of this distribution is connected.
-     *
-     * @return {@code true}
-     */
-    @Override
-    public boolean isSupportConnected() {
-        return true;
-    }
-
-    /**{@inheritDoc} */
-    @Override
-    public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
-        return new DiscreteDistribution.Sampler() {
-            /**
-             * Zipf distribution sampler.
-             */
-            private final DiscreteSampler sampler =
-                new RejectionInversionZipfSampler(rng, numberOfElements, exponent);
-
-            /**{@inheritDoc} */
-            @Override
-            public int sample() {
-                return sampler.sample();
-            }
-        };
-    }
-}

http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/package-info.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/package-info.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/package-info.java
deleted file mode 100644
index 98315c6..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/package-info.java
+++ /dev/null
@@ -1,20 +0,0 @@
-/*
- * 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.
- */
-/**
- * Implementations of common discrete and continuous distributions.
- */
-package org.apache.commons.statistics.distribution;

http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractContinuousDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractContinuousDistribution.java b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractContinuousDistribution.java
new file mode 100644
index 0000000..1d6b254
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractContinuousDistribution.java
@@ -0,0 +1,453 @@
+/*
+ * 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.
+ */
+package org.apache.commons.statistics.distribution;
+
+import java.util.function.DoubleUnaryOperator;
+import org.apache.commons.numbers.core.Precision;
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.sampling.distribution.InverseTransformContinuousSampler;
+import org.apache.commons.rng.sampling.distribution.ContinuousInverseCumulativeProbabilityFunction;
+import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
+
+/**
+ * Base class for probability distributions on the reals.
+ * Default implementations are provided for some of the methods
+ * that do not vary from distribution to distribution.
+ *
+ * This base class provides a default factory method for creating
+ * a {@link ContinuousDistribution.Sampler sampler instance} that uses the
+ * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">
+ * inversion method</a> for generating random samples that follow the
+ * distribution.
+ */
+public abstract class AbstractContinuousDistribution
+    implements ContinuousDistribution {
+    /**
+     * For a random variable {@code X} whose values are distributed according
+     * to this distribution, this method returns {@code P(x0 < X <= x1)}.
+     *
+     * @param x0 Lower bound (excluded).
+     * @param x1 Upper bound (included).
+     * @return the probability that a random variable with this distribution
+     * takes a value between {@code x0} and {@code x1}, excluding the lower
+     * and including the upper endpoint.
+     * @throws IllegalArgumentException if {@code x0 > x1}.
+     *
+     * The default implementation uses the identity
+     * {@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}
+     */
+    @Override
+    public double probability(double x0,
+                              double x1) {
+        if (x0 > x1) {
+            throw new DistributionException(DistributionException.TOO_LARGE, x0, x1);
+        }
+        return cumulativeProbability(x1) - cumulativeProbability(x0);
+    }
+
+    /**
+     * {@inheritDoc}
+     *
+     * The default implementation returns
+     * <ul>
+     * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
+     * <li>{@link #getSupportUpperBound()} for {@code p = 1}.</li>
+     * </ul>
+     */
+    @Override
+    public double inverseCumulativeProbability(final double p) {
+        /*
+         * IMPLEMENTATION NOTES
+         * --------------------
+         * Where applicable, use is made of the one-sided Chebyshev inequality
+         * to bracket the root. This inequality states that
+         * P(X - mu >= k * sig) <= 1 / (1 + k^2),
+         * mu: mean, sig: standard deviation. Equivalently
+         * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
+         * F(mu + k * sig) >= k^2 / (1 + k^2).
+         *
+         * For k = sqrt(p / (1 - p)), we find
+         * F(mu + k * sig) >= p,
+         * and (mu + k * sig) is an upper-bound for the root.
+         *
+         * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
+         * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
+         * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
+         * P(X <= mu - k * sig) <= 1 / (1 + k^2),
+         * F(mu - k * sig) <= 1 / (1 + k^2).
+         *
+         * For k = sqrt((1 - p) / p), we find
+         * F(mu - k * sig) <= p,
+         * and (mu - k * sig) is a lower-bound for the root.
+         *
+         * In cases where the Chebyshev inequality does not apply, geometric
+         * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
+         * the root.
+         */
+        if (p < 0 ||
+            p > 1) {
+            throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
+        }
+
+        double lowerBound = getSupportLowerBound();
+        if (p == 0) {
+            return lowerBound;
+        }
+
+        double upperBound = getSupportUpperBound();
+        if (p == 1) {
+            return upperBound;
+        }
+
+        final double mu = getNumericalMean();
+        final double sig = Math.sqrt(getNumericalVariance());
+        final boolean chebyshevApplies;
+        chebyshevApplies = !(Double.isInfinite(mu) ||
+                             Double.isNaN(mu) ||
+                             Double.isInfinite(sig) ||
+                             Double.isNaN(sig));
+
+        if (lowerBound == Double.NEGATIVE_INFINITY) {
+            if (chebyshevApplies) {
+                lowerBound = mu - sig * Math.sqrt((1 - p) / p);
+            } else {
+                lowerBound = -1;
+                while (cumulativeProbability(lowerBound) >= p) {
+                    lowerBound *= 2;
+                }
+            }
+        }
+
+        if (upperBound == Double.POSITIVE_INFINITY) {
+            if (chebyshevApplies) {
+                upperBound = mu + sig * Math.sqrt(p / (1 - p));
+            } else {
+                upperBound = 1;
+                while (cumulativeProbability(upperBound) < p) {
+                    upperBound *= 2;
+                }
+            }
+        }
+
+        // XXX Values copied from defaults in class
+        // "o.a.c.math4.analysis.solvers.BaseAbstractUnivariateSolver"
+        final double solverRelativeAccuracy = 1e-14;
+        final double solverAbsoluteAccuracy = 1e-9;
+        final double solverFunctionValueAccuracy = 1e-15;
+
+        double x = new BrentSolver(solverRelativeAccuracy,
+                                   solverAbsoluteAccuracy,
+                                   solverFunctionValueAccuracy)
+            .solve((arg) -> cumulativeProbability(arg) - p,
+                   lowerBound,
+                   0.5 * (lowerBound + upperBound),
+                   upperBound);
+
+        if (!isSupportConnected()) {
+            /* Test for plateau. */
+            final double dx = solverAbsoluteAccuracy;
+            if (x - dx >= getSupportLowerBound()) {
+                double px = cumulativeProbability(x);
+                if (cumulativeProbability(x - dx) == px) {
+                    upperBound = x;
+                    while (upperBound - lowerBound > dx) {
+                        final double midPoint = 0.5 * (lowerBound + upperBound);
+                        if (cumulativeProbability(midPoint) < px) {
+                            lowerBound = midPoint;
+                        } else {
+                            upperBound = midPoint;
+                        }
+                    }
+                    return upperBound;
+                }
+            }
+        }
+        return x;
+    }
+
+    /**
+     * {@inheritDoc}
+     *
+     * @return zero.
+     */
+    @Override
+    public double probability(double x) {
+        return 0;
+    }
+
+    /**
+     * {@inheritDoc}
+     *
+     * The default implementation computes the logarithm of {@code density(x)}.
+     */
+    @Override
+    public double logDensity(double x) {
+        return Math.log(density(x));
+    }
+
+    /**
+     * Utility function for allocating an array and filling it with {@code n}
+     * samples generated by the given {@code sampler}.
+     *
+     * @param n Number of samples.
+     * @param sampler Sampler.
+     * @return an array of size {@code n}.
+     */
+    public static double[] sample(int n,
+                                  ContinuousDistribution.Sampler sampler) {
+        final double[] samples = new double[n];
+        for (int i = 0; i < n; i++) {
+            samples[i] = sampler.sample();
+        }
+        return samples;
+    }
+
+    /**{@inheritDoc} */
+    @Override
+    public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
+        return new ContinuousDistribution.Sampler() {
+            /**
+             * Inversion method distribution sampler.
+             */
+            private final ContinuousSampler sampler =
+                new InverseTransformContinuousSampler(rng, createICPF());
+
+            /** {@inheritDoc} */
+            @Override
+            public double sample() {
+                return sampler.sample();
+            }
+        };
+    }
+
+    /**
+     * @return an instance for use by {@link #createSampler(UniformRandomProvider)}
+     */
+    private ContinuousInverseCumulativeProbabilityFunction createICPF() {
+        return new ContinuousInverseCumulativeProbabilityFunction() {
+            /** {@inheritDoc} */
+            @Override
+            public double inverseCumulativeProbability(double p) {
+                return AbstractContinuousDistribution.this.inverseCumulativeProbability(p);
+            }
+        };
+    }
+
+    /**
+     * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
+     * Brent algorithm</a> for finding zeros of real univariate functions.
+     * The function should be continuous but not necessarily smooth.
+     * The {@code solve} method returns a zero {@code x} of the function {@code f}
+     * in the given interval {@code [a, b]} to within a tolerance
+     * {@code 2 eps abs(x) + t} where {@code eps} is the relative accuracy and
+     * {@code t} is the absolute accuracy.
+     * <p>The given interval must bracket the root.</p>
+     * <p>
+     *  The reference implementation is given in chapter 4 of
+     *  <blockquote>
+     *   <b>Algorithms for Minimization Without Derivatives</b>,
+     *   <em>Richard P. Brent</em>,
+     *   Dover, 2002
+     *  </blockquote>
+     *
+     * Used by {@link #inverseCumulativeProbability(double)}.
+     */
+    private static class BrentSolver {
+        /** Relative accuracy. */
+        private final double relativeAccuracy;
+        /** Absolutee accuracy. */
+        private final double absoluteAccuracy;
+        /** Function accuracy. */
+        private final double functionValueAccuracy;
+
+        /**
+         * Construct a solver.
+         *
+         * @param relativeAccuracy Relative accuracy.
+         * @param absoluteAccuracy Absolute accuracy.
+         * @param functionValueAccuracy Function value accuracy.
+         */
+        BrentSolver(double relativeAccuracy,
+                    double absoluteAccuracy,
+                    double functionValueAccuracy) {
+            this.relativeAccuracy = relativeAccuracy;
+            this.absoluteAccuracy = absoluteAccuracy;
+            this.functionValueAccuracy = functionValueAccuracy;
+        }
+
+        /**
+         * @param func Function to solve.
+         * @param min Lower bound.
+         * @param init Initial guess.
+         * @param max Upper bound.
+         * @return the root.
+         */
+        double solve(DoubleUnaryOperator func,
+                     double min,
+                     double initial,
+                     double max) {
+            if (min > max) {
+                throw new DistributionException(DistributionException.TOO_LARGE, min, max);
+            }
+            if (initial < min ||
+                initial > max) {
+                throw new DistributionException(DistributionException.OUT_OF_RANGE, initial, min, max);
+            }
+
+            // Return the initial guess if it is good enough.
+            double yInitial = func.applyAsDouble(initial);
+            if (Math.abs(yInitial) <= functionValueAccuracy) {
+                return initial;
+            }
+
+            // Return the first endpoint if it is good enough.
+            double yMin = func.applyAsDouble(min);
+            if (Math.abs(yMin) <= functionValueAccuracy) {
+                return min;
+            }
+
+            // Reduce interval if min and initial bracket the root.
+            if (yInitial * yMin < 0) {
+                return brent(func, min, initial, yMin, yInitial);
+            }
+
+            // Return the second endpoint if it is good enough.
+            double yMax = func.applyAsDouble(max);
+            if (Math.abs(yMax) <= functionValueAccuracy) {
+                return max;
+            }
+
+            // Reduce interval if initial and max bracket the root.
+            if (yInitial * yMax < 0) {
+                return brent(func, initial, max, yInitial, yMax);
+            }
+
+            throw new DistributionException(DistributionException.BRACKETING, min, yMin, max, yMax);
+        }
+
+        /**
+         * Search for a zero inside the provided interval.
+         * This implementation is based on the algorithm described at page 58 of
+         * the book
+         * <blockquote>
+         *  <b>Algorithms for Minimization Without Derivatives</b>,
+         *  <it>Richard P. Brent</it>,
+         *  Dover 0-486-41998-3
+         * </blockquote>
+         *
+         * @param func Function to solve.
+         * @param lo Lower bound of the search interval.
+         * @param hi Higher bound of the search interval.
+         * @param fLo Function value at the lower bound of the search interval.
+         * @param fHi Function value at the higher bound of the search interval.
+         * @return the value where the function is zero.
+         */
+        private double brent(DoubleUnaryOperator func,
+                             double lo, double hi,
+                             double fLo, double fHi) {
+            double a = lo;
+            double fa = fLo;
+            double b = hi;
+            double fb = fHi;
+            double c = a;
+            double fc = fa;
+            double d = b - a;
+            double e = d;
+
+            final double t = absoluteAccuracy;
+            final double eps = relativeAccuracy;
+
+            while (true) {
+                if (Math.abs(fc) < Math.abs(fb)) {
+                    a = b;
+                    b = c;
+                    c = a;
+                    fa = fb;
+                    fb = fc;
+                    fc = fa;
+                }
+
+                final double tol = 2 * eps * Math.abs(b) + t;
+                final double m = 0.5 * (c - b);
+
+                if (Math.abs(m) <= tol ||
+                    Precision.equals(fb, 0))  {
+                    return b;
+                }
+                if (Math.abs(e) < tol ||
+                    Math.abs(fa) <= Math.abs(fb)) {
+                    // Force bisection.
+                    d = m;
+                    e = d;
+                } else {
+                    double s = fb / fa;
+                    double p;
+                    double q;
+                    // The equality test (a == c) is intentional,
+                    // it is part of the original Brent's method and
+                    // it should NOT be replaced by proximity test.
+                    if (a == c) {
+                        // Linear interpolation.
+                        p = 2 * m * s;
+                        q = 1 - s;
+                    } else {
+                        // Inverse quadratic interpolation.
+                        q = fa / fc;
+                        final double r = fb / fc;
+                        p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
+                        q = (q - 1) * (r - 1) * (s - 1);
+                    }
+                    if (p > 0) {
+                        q = -q;
+                    } else {
+                        p = -p;
+                    }
+                    s = e;
+                    e = d;
+                    if (p >= 1.5 * m * q - Math.abs(tol * q) ||
+                        p >= Math.abs(0.5 * s * q)) {
+                        // Inverse quadratic interpolation gives a value
+                        // in the wrong direction, or progress is slow.
+                        // Fall back to bisection.
+                        d = m;
+                        e = d;
+                    } else {
+                        d = p / q;
+                    }
+                }
+                a = b;
+                fa = fb;
+
+                if (Math.abs(d) > tol) {
+                    b += d;
+                } else if (m > 0) {
+                    b += tol;
+                } else {
+                    b -= tol;
+                }
+                fb = func.applyAsDouble(b);
+                if ((fb > 0 && fc > 0) ||
+                    (fb <= 0 && fc <= 0)) {
+                    c = a;
+                    fc = fa;
+                    d = b - a;
+                    e = d;
+                }
+            }
+        }
+    }
+}

http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractDiscreteDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractDiscreteDistribution.java b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractDiscreteDistribution.java
new file mode 100644
index 0000000..faef96c
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractDiscreteDistribution.java
@@ -0,0 +1,220 @@
+/*
+ * 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.
+ */
+package org.apache.commons.statistics.distribution;
+
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.sampling.distribution.InverseTransformDiscreteSampler;
+import org.apache.commons.rng.sampling.distribution.DiscreteInverseCumulativeProbabilityFunction;
+import org.apache.commons.rng.sampling.distribution.DiscreteSampler;
+
+/**
+ * Base class for integer-valued discrete distributions.  Default
+ * implementations are provided for some of the methods that do not vary
+ * from distribution to distribution.
+ */
+public abstract class AbstractDiscreteDistribution
+    implements DiscreteDistribution {
+    /**
+     * {@inheritDoc}
+     *
+     * The default implementation uses the identity
+     * {@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}
+     */
+    @Override
+    public double probability(int x0,
+                              int x1) {
+        if (x1 < x0) {
+            throw new DistributionException(DistributionException.TOO_SMALL,
+                                            x1, x0);
+        }
+        return cumulativeProbability(x1) - cumulativeProbability(x0);
+    }
+
+    /**
+     * {@inheritDoc}
+     *
+     * The default implementation returns
+     * <ul>
+     * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
+     * <li>{@link #getSupportUpperBound()} for {@code p = 1}, and</li>
+     * <li>{@link #solveInverseCumulativeProbability(double, int, int)} for
+     *     {@code 0 < p < 1}.</li>
+     * </ul>
+     */
+    @Override
+    public int inverseCumulativeProbability(final double p) {
+        if (p < 0 ||
+            p > 1) {
+            throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
+        }
+
+        int lower = getSupportLowerBound();
+        if (p == 0.0) {
+            return lower;
+        }
+        if (lower == Integer.MIN_VALUE) {
+            if (checkedCumulativeProbability(lower) >= p) {
+                return lower;
+            }
+        } else {
+            lower -= 1; // this ensures cumulativeProbability(lower) < p, which
+                        // is important for the solving step
+        }
+
+        int upper = getSupportUpperBound();
+        if (p == 1.0) {
+            return upper;
+        }
+
+        // use the one-sided Chebyshev inequality to narrow the bracket
+        // cf. AbstractRealDistribution.inverseCumulativeProbability(double)
+        final double mu = getNumericalMean();
+        final double sigma = Math.sqrt(getNumericalVariance());
+        final boolean chebyshevApplies = !(Double.isInfinite(mu) ||
+                                           Double.isNaN(mu) ||
+                                           Double.isInfinite(sigma) ||
+                                           Double.isNaN(sigma) ||
+                                           sigma == 0.0);
+        if (chebyshevApplies) {
+            double k = Math.sqrt((1.0 - p) / p);
+            double tmp = mu - k * sigma;
+            if (tmp > lower) {
+                lower = ((int) Math.ceil(tmp)) - 1;
+            }
+            k = 1.0 / k;
+            tmp = mu + k * sigma;
+            if (tmp < upper) {
+                upper = ((int) Math.ceil(tmp)) - 1;
+            }
+        }
+
+        return solveInverseCumulativeProbability(p, lower, upper);
+    }
+
+    /**
+     * This is a utility function used by {@link
+     * #inverseCumulativeProbability(double)}. It assumes {@code 0 < p < 1} and
+     * that the inverse cumulative probability lies in the bracket {@code
+     * (lower, upper]}. The implementation does simple bisection to find the
+     * smallest {@code p}-quantile {@code inf{x in Z | P(X <= x) >= p}}.
+     *
+     * @param p Cumulative probability.
+     * @param lower Value satisfying {@code cumulativeProbability(lower) < p}.
+     * @param upper Value satisfying {@code p <= cumulativeProbability(upper)}.
+     * @return the smallest {@code p}-quantile of this distribution.
+     */
+    private int solveInverseCumulativeProbability(final double p,
+                                                  int lower,
+                                                  int upper) {
+        while (lower + 1 < upper) {
+            int xm = (lower + upper) / 2;
+            if (xm < lower || xm > upper) {
+                /*
+                 * Overflow.
+                 * There will never be an overflow in both calculation methods
+                 * for xm at the same time
+                 */
+                xm = lower + (upper - lower) / 2;
+            }
+
+            double pm = checkedCumulativeProbability(xm);
+            if (pm >= p) {
+                upper = xm;
+            } else {
+                lower = xm;
+            }
+        }
+        return upper;
+    }
+
+    /**
+     * Computes the cumulative probability function and checks for {@code NaN}
+     * values returned. Throws {@code MathInternalError} if the value is
+     * {@code NaN}. Rethrows any exception encountered evaluating the cumulative
+     * probability function. Throws {@code MathInternalError} if the cumulative
+     * probability function returns {@code NaN}.
+     *
+     * @param argument Input value.
+     * @return the cumulative probability.
+     * @throws IllegalStateException if the cumulative probability is {@code NaN}.
+     */
+    private double checkedCumulativeProbability(int argument) {
+        final double result = cumulativeProbability(argument);
+        if (Double.isNaN(result)) {
+            throw new IllegalStateException("Internal error");
+        }
+        return result;
+    }
+
+    /**
+     * {@inheritDoc}
+     *
+     * The default implementation simply computes the logarithm of {@code probability(x)}.
+     */
+    @Override
+    public double logProbability(int x) {
+        return Math.log(probability(x));
+    }
+
+    /**
+     * Utility function for allocating an array and filling it with {@code n}
+     * samples generated by the given {@code sampler}.
+     *
+     * @param n Number of samples.
+     * @param sampler Sampler.
+     * @return an array of size {@code n}.
+     */
+    public static int[] sample(int n,
+                               DiscreteDistribution.Sampler sampler) {
+        final int[] samples = new int[n];
+        for (int i = 0; i < n; i++) {
+            samples[i] = sampler.sample();
+        }
+        return samples;
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
+        return new DiscreteDistribution.Sampler() {
+            /**
+             * Inversion method distribution sampler.
+             */
+            private final DiscreteSampler sampler =
+                new InverseTransformDiscreteSampler(rng, createICPF());
+
+            /** {@inheritDoc} */
+            @Override
+            public int sample() {
+                return sampler.sample();
+            }
+        };
+    }
+
+    /**
+     * @return an instance for use by {@link #createSampler(UniformRandomProvider)}.
+     */
+    private DiscreteInverseCumulativeProbabilityFunction createICPF() {
+        return new DiscreteInverseCumulativeProbabilityFunction() {
+            /** {@inheritDoc} */
+            @Override
+            public int inverseCumulativeProbability(double p) {
+                return AbstractDiscreteDistribution.this.inverseCumulativeProbability(p);
+            }
+        };
+    }
+}


Mime
View raw message