You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2018/01/21 14:05:46 UTC
[11/16] commons-statistics git commit: STATISTICS-2: Migrate
"o.a.c.math4.distribution" from Commons Math.
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/NormalDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/NormalDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/NormalDistribution.java
new file mode 100644
index 0000000..632657f
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/NormalDistribution.java
@@ -0,0 +1,216 @@
+/*
+ * 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.Erfc;
+import org.apache.commons.numbers.gamma.InverseErf;
+import org.apache.commons.numbers.gamma.ErfDifference;
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
+import org.apache.commons.rng.sampling.distribution.GaussianSampler;
+import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler;
+
+/**
+ * Implementation of the <a href="http://en.wikipedia.org/wiki/Normal_distribution">normal (Gaussian) distribution</a>.
+ */
+public class NormalDistribution extends AbstractContinuousDistribution {
+ /** √(2) */
+ private static final double SQRT2 = Math.sqrt(2.0);
+ /** Mean of this distribution. */
+ private final double mean;
+ /** Standard deviation of this distribution. */
+ private final double standardDeviation;
+ /** The value of {@code log(sd) + 0.5*log(2*pi)} stored for faster computation. */
+ private final double logStandardDeviationPlusHalfLog2Pi;
+
+ /**
+ * Create a normal distribution with mean equal to zero and standard
+ * deviation equal to one.
+ */
+ public NormalDistribution() {
+ this(0, 1);
+ }
+
+ /**
+ * Creates a distribution.
+ *
+ * @param mean Mean for this distribution.
+ * @param sd Standard deviation for this distribution.
+ * @throws IllegalArgumentException if {@code sd <= 0}.
+ */
+ public NormalDistribution(double mean,
+ double sd) {
+ if (sd <= 0) {
+ throw new DistributionException(DistributionException.NEGATIVE, sd);
+ }
+
+ this.mean = mean;
+ standardDeviation = sd;
+ logStandardDeviationPlusHalfLog2Pi = Math.log(sd) + 0.5 * Math.log(2 * Math.PI);
+ }
+
+ /**
+ * Access the mean.
+ *
+ * @return the mean for this distribution.
+ */
+ public double getMean() {
+ return mean;
+ }
+
+ /**
+ * Access the standard deviation.
+ *
+ * @return the standard deviation for this distribution.
+ */
+ public double getStandardDeviation() {
+ return standardDeviation;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double density(double x) {
+ return Math.exp(logDensity(x));
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double logDensity(double x) {
+ final double x0 = x - mean;
+ final double x1 = x0 / standardDeviation;
+ return -0.5 * x1 * x1 - logStandardDeviationPlusHalfLog2Pi;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * If {@code x} is more than 40 standard deviations from the mean, 0 or 1
+ * is returned, as in these cases the actual value is within
+ * {@code Double.MIN_VALUE} of 0 or 1.
+ */
+ @Override
+ public double cumulativeProbability(double x) {
+ final double dev = x - mean;
+ if (Math.abs(dev) > 40 * standardDeviation) {
+ return dev < 0 ? 0.0d : 1.0d;
+ }
+ return 0.5 * Erfc.value(-dev / (standardDeviation * SQRT2));
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double inverseCumulativeProbability(final double p) {
+ if (p < 0 ||
+ p > 1) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
+ }
+ return mean + standardDeviation * SQRT2 * InverseErf.value(2 * p - 1);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double probability(double x0,
+ double x1) {
+ if (x0 > x1) {
+ throw new DistributionException(DistributionException.TOO_LARGE,
+ x0, x1);
+ }
+ final double denom = standardDeviation * SQRT2;
+ final double v0 = (x0 - mean) / denom;
+ final double v1 = (x1 - mean) / denom;
+ return 0.5 * ErfDifference.value(v0, v1);
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * For mean parameter {@code mu}, the mean is {@code mu}.
+ */
+ @Override
+ public double getNumericalMean() {
+ return getMean();
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * For standard deviation parameter {@code s}, the variance is {@code s^2}.
+ */
+ @Override
+ public double getNumericalVariance() {
+ final double s = getStandardDeviation();
+ return s * s;
+ }
+
+ /**
+ * {@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;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
+ return new ContinuousDistribution.Sampler() {
+ /** Gaussian distribution sampler. */
+ private final ContinuousSampler sampler =
+ new GaussianSampler(new ZigguratNormalizedGaussianSampler(rng),
+ mean, standardDeviation);
+
+ /** {@inheritDoc} */
+ @Override
+ public double sample() {
+ return sampler.sample();
+ }
+ };
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ParetoDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ParetoDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ParetoDistribution.java
new file mode 100644
index 0000000..2bbd42d
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ParetoDistribution.java
@@ -0,0 +1,225 @@
+/*
+ * 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.InverseTransformParetoSampler;
+
+/**
+ * Implementation of the <a href="http://en.wikipedia.org/wiki/Pareto_distribution">Pareto distribution</a>.
+ *
+ * <p>
+ * <strong>Parameters:</strong>
+ * The probability distribution function of {@code X} is given by (for {@code x >= k}):
+ * <pre>
+ * α * k^α / x^(α + 1)
+ * </pre>
+ * <ul>
+ * <li>{@code k} is the <em>scale</em> parameter: this is the minimum possible value of {@code X},</li>
+ * <li>{@code α} is the <em>shape</em> parameter: this is the Pareto index</li>
+ * </ul>
+ */
+public class ParetoDistribution extends AbstractContinuousDistribution {
+ /** The scale parameter of this distribution. */
+ private final double scale;
+ /** The shape parameter of this distribution. */
+ private final double shape;
+
+ /**
+ * Creates a Pareto distribution with a scale of {@code 1} and a shape of {@code 1}.
+ */
+ public ParetoDistribution() {
+ this(1, 1);
+ }
+
+ /**
+ * Creates a Pareto distribution.
+ *
+ * @param scale Scale parameter of this distribution.
+ * @param shape Shape parameter of this distribution.
+ * @throws IllegalArgumentException if {@code scale <= 0} or {@code shape <= 0}.
+ */
+ public ParetoDistribution(double scale,
+ double shape) {
+ if (scale <= 0) {
+ throw new DistributionException(DistributionException.NEGATIVE, scale);
+ }
+
+ if (shape <= 0) {
+ throw new DistributionException(DistributionException.NEGATIVE, shape);
+ }
+
+ this.scale = scale;
+ this.shape = shape;
+ }
+
+ /**
+ * Returns the scale parameter of this distribution.
+ *
+ * @return the scale parameter
+ */
+ public double getScale() {
+ return scale;
+ }
+
+ /**
+ * Returns the shape parameter of this distribution.
+ *
+ * @return the shape parameter
+ */
+ public double getShape() {
+ return shape;
+ }
+
+ /**
+ * {@inheritDoc}
+ * <p>
+ * For scale {@code k}, and shape {@code α} of this distribution, the PDF
+ * is given by
+ * <ul>
+ * <li>{@code 0} if {@code x < k},</li>
+ * <li>{@code α * k^α / x^(α + 1)} otherwise.</li>
+ * </ul>
+ */
+ @Override
+ public double density(double x) {
+ if (x < scale) {
+ return 0;
+ }
+ return Math.pow(scale, shape) / Math.pow(x, shape + 1) * shape;
+ }
+
+ /** {@inheritDoc}
+ *
+ * See documentation of {@link #density(double)} for computation details.
+ */
+ @Override
+ public double logDensity(double x) {
+ if (x < scale) {
+ return Double.NEGATIVE_INFINITY;
+ }
+ return Math.log(scale) * shape - Math.log(x) * (shape + 1) + Math.log(shape);
+ }
+
+ /**
+ * {@inheritDoc}
+ * <p>
+ * For scale {@code k}, and shape {@code α} of this distribution, the CDF is given by
+ * <ul>
+ * <li>{@code 0} if {@code x < k},</li>
+ * <li>{@code 1 - (k / x)^α} otherwise.</li>
+ * </ul>
+ */
+ @Override
+ public double cumulativeProbability(double x) {
+ if (x <= scale) {
+ return 0;
+ }
+ return 1 - Math.pow(scale / x, shape);
+ }
+
+ /**
+ * {@inheritDoc}
+ * <p>
+ * For scale {@code k} and shape {@code α}, the mean is given by
+ * <ul>
+ * <li>{@code ∞} if {@code α <= 1},</li>
+ * <li>{@code α * k / (α - 1)} otherwise.</li>
+ * </ul>
+ */
+ @Override
+ public double getNumericalMean() {
+ if (shape <= 1) {
+ return Double.POSITIVE_INFINITY;
+ }
+ return shape * scale / (shape - 1);
+ }
+
+ /**
+ * {@inheritDoc}
+ * <p>
+ * For scale {@code k} and shape {@code α}, the variance is given by
+ * <ul>
+ * <li>{@code ∞} if {@code 1 < α <= 2},</li>
+ * <li>{@code k^2 * α / ((α - 1)^2 * (α - 2))} otherwise.</li>
+ * </ul>
+ */
+ @Override
+ public double getNumericalVariance() {
+ if (shape <= 2) {
+ return Double.POSITIVE_INFINITY;
+ }
+ double s = shape - 1;
+ return scale * scale * shape / (s * s) / (shape - 2);
+ }
+
+ /**
+ * {@inheritDoc}
+ * <p>
+ * The lower bound of the support is equal to the scale parameter {@code k}.
+ *
+ * @return lower bound of the support
+ */
+ @Override
+ public double getSupportLowerBound() {
+ return scale;
+ }
+
+ /**
+ * {@inheritDoc}
+ * <p>
+ * 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}
+ * <p>
+ * 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() {
+ /**
+ * Pareto distribution sampler.
+ */
+ private final ContinuousSampler sampler =
+ new InverseTransformParetoSampler(rng, scale, shape);
+
+ /**{@inheritDoc} */
+ @Override
+ public double sample() {
+ return sampler.sample();
+ }
+ };
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/PascalDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/PascalDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/PascalDistribution.java
new file mode 100644
index 0000000..8bdf6b6
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/PascalDistribution.java
@@ -0,0 +1,211 @@
+/*
+ * 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.combinatorics.BinomialCoefficientDouble;
+import org.apache.commons.numbers.combinatorics.LogBinomialCoefficient;
+import org.apache.commons.numbers.gamma.RegularizedBeta;
+
+/**
+ * Implementation of the <a href="http://en.wikipedia.org/wiki/Negative_binomial_distribution">Pascal distribution.</a>
+ *
+ * The Pascal distribution is a special case of the Negative Binomial distribution
+ * where the number of successes parameter is an integer.
+ *
+ * There are various ways to express the probability mass and distribution
+ * functions for the Pascal distribution. The present implementation represents
+ * the distribution of the number of failures before {@code r} successes occur.
+ * This is the convention adopted in e.g.
+ * <a href="http://mathworld.wolfram.com/NegativeBinomialDistribution.html">MathWorld</a>,
+ * but <em>not</em> in
+ * <a href="http://en.wikipedia.org/wiki/Negative_binomial_distribution">Wikipedia</a>.
+ *
+ * For a random variable {@code X} whose values are distributed according to this
+ * distribution, the probability mass function is given by<br>
+ * {@code P(X = k) = C(k + r - 1, r - 1) * p^r * (1 - p)^k,}<br>
+ * where {@code r} is the number of successes, {@code p} is the probability of
+ * success, and {@code X} is the total number of failures. {@code C(n, k)} is
+ * the binomial coefficient ({@code n} choose {@code k}). The mean and variance
+ * of {@code X} are<br>
+ * {@code E(X) = (1 - p) * r / p, var(X) = (1 - p) * r / p^2.}<br>
+ * Finally, the cumulative distribution function is given by<br>
+ * {@code P(X <= k) = I(p, r, k + 1)},
+ * where I is the regularized incomplete Beta function.
+ */
+public class PascalDistribution extends AbstractDiscreteDistribution {
+ /** The number of successes. */
+ private final int numberOfSuccesses;
+ /** The probability of success. */
+ private final double probabilityOfSuccess;
+ /** The value of {@code log(p)}, where {@code p} is the probability of success,
+ * stored for faster computation. */
+ private final double logProbabilityOfSuccess;
+ /** The value of {@code log(1-p)}, where {@code p} is the probability of success,
+ * stored for faster computation. */
+ private final double log1mProbabilityOfSuccess;
+
+ /**
+ * Create a Pascal distribution with the given number of successes and
+ * probability of success.
+ *
+ * @param r Number of successes.
+ * @param p Probability of success.
+ * @throws IllegalArgumentException if {@code r <= 0} or {@code p < 0}
+ * or {@code p > 1}.
+ */
+ public PascalDistribution(int r,
+ double p) {
+ if (r <= 0) {
+ throw new DistributionException(DistributionException.NEGATIVE,
+ r);
+ }
+ if (p < 0 ||
+ p > 1) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
+ }
+
+ numberOfSuccesses = r;
+ probabilityOfSuccess = p;
+ logProbabilityOfSuccess = Math.log(p);
+ log1mProbabilityOfSuccess = Math.log1p(-p);
+ }
+
+ /**
+ * Access the number of successes for this distribution.
+ *
+ * @return the number of successes.
+ */
+ public int getNumberOfSuccesses() {
+ return numberOfSuccesses;
+ }
+
+ /**
+ * Access the probability of success for this distribution.
+ *
+ * @return the probability of success.
+ */
+ public double getProbabilityOfSuccess() {
+ return probabilityOfSuccess;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double probability(int x) {
+ double ret;
+ if (x < 0) {
+ ret = 0.0;
+ } else {
+ ret = BinomialCoefficientDouble.value(x +
+ numberOfSuccesses - 1, numberOfSuccesses - 1) *
+ Math.pow(probabilityOfSuccess, numberOfSuccesses) *
+ Math.pow(1.0 - probabilityOfSuccess, x);
+ }
+ return ret;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double logProbability(int x) {
+ double ret;
+ if (x < 0) {
+ ret = Double.NEGATIVE_INFINITY;
+ } else {
+ ret = LogBinomialCoefficient.value(x +
+ numberOfSuccesses - 1, numberOfSuccesses - 1) +
+ logProbabilityOfSuccess * numberOfSuccesses +
+ log1mProbabilityOfSuccess * x;
+ }
+ return ret;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double cumulativeProbability(int x) {
+ double ret;
+ if (x < 0) {
+ ret = 0.0;
+ } else {
+ ret = RegularizedBeta.value(probabilityOfSuccess,
+ numberOfSuccesses, x + 1.0);
+ }
+ return ret;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * For number of successes {@code r} and probability of success {@code p},
+ * the mean is {@code r * (1 - p) / p}.
+ */
+ @Override
+ public double getNumericalMean() {
+ final double p = getProbabilityOfSuccess();
+ final double r = getNumberOfSuccesses();
+ return (r * (1 - p)) / p;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * For number of successes {@code r} and probability of success {@code p},
+ * the variance is {@code r * (1 - p) / p^2}.
+ */
+ @Override
+ public double getNumericalVariance() {
+ final double p = getProbabilityOfSuccess();
+ final double r = getNumberOfSuccesses();
+ return r * (1 - p) / (p * p);
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * The lower bound of the support is always 0 no matter the parameters.
+ *
+ * @return lower bound of the support (always 0)
+ */
+ @Override
+ public int getSupportLowerBound() {
+ return 0;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * The upper bound of the support is always positive infinity no matter the
+ * parameters. Positive infinity is symbolized by {@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;
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/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
new file mode 100644
index 0000000..225b8f1
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/PoissonDistribution.java
@@ -0,0 +1,238 @@
+/*
+ * 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 π). */
+ 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/9c794a15/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
new file mode 100644
index 0000000..7bb847a
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/SaddlePointExpansion.java
@@ -0,0 +1,191 @@
+/*
+ * 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 π */
+ private static final double TWO_PI = 2 * Math.PI;
+ /** 1/2 * log(2 π). */
+ 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/9c794a15/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
new file mode 100644
index 0000000..ef34c1f
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TDistribution.java
@@ -0,0 +1,180 @@
+/*
+ * 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'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/9c794a15/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
new file mode 100644
index 0000000..6a94b62
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TriangularDistribution.java
@@ -0,0 +1,222 @@
+/*
+ * 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/9c794a15/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
new file mode 100644
index 0000000..2bb9a0b
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformContinuousDistribution.java
@@ -0,0 +1,168 @@
+/*
+ * 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/9c794a15/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
new file mode 100644
index 0000000..41df2bd
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformDiscreteDistribution.java
@@ -0,0 +1,159 @@
+/*
+ * 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/9c794a15/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
new file mode 100644
index 0000000..0e396d8
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/WeibullDistribution.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.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/9c794a15/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
new file mode 100644
index 0000000..f0d54a1
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ZipfDistribution.java
@@ -0,0 +1,236 @@
+/*
+ * 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/9c794a15/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
new file mode 100644
index 0000000..98315c6
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/package-info.java
@@ -0,0 +1,20 @@
+/*
+ * 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/9c794a15/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/AbstractContinuousDistributionTest.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/AbstractContinuousDistributionTest.java b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/AbstractContinuousDistributionTest.java
new file mode 100644
index 0000000..3e30763
--- /dev/null
+++ b/commons-statistics-distribution/src/test/java/org/apache/commons/statistics/distribution/AbstractContinuousDistributionTest.java
@@ -0,0 +1,209 @@
+/*
+ * 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.math3.analysis.UnivariateFunction;
+import org.apache.commons.math3.analysis.integration.RombergIntegrator;
+import org.apache.commons.math3.analysis.integration.UnivariateIntegrator;
+import org.junit.Assert;
+import org.junit.Test;
+
+/** Various tests related to MATH-699. */
+public class AbstractContinuousDistributionTest {
+
+ @Test
+ public void testContinuous() {
+ final double x0 = 0.0;
+ final double x1 = 1.0;
+ final double x2 = 2.0;
+ final double x3 = 3.0;
+ final double p12 = 0.5;
+ final AbstractContinuousDistribution distribution;
+ distribution = new AbstractContinuousDistribution() {
+ private static final long serialVersionUID = 1L;
+
+ @Override
+ public double cumulativeProbability(final double x) {
+ if (x < x0 ||
+ x > x3) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, x, x0, x3);
+ }
+ if (x <= x1) {
+ return p12 * (x - x0) / (x1 - x0);
+ } else if (x <= x2) {
+ return p12;
+ } else if (x <= x3) {
+ return p12 + (1.0 - p12) * (x - x2) / (x3 - x2);
+ }
+ return 0.0;
+ }
+
+ @Override
+ public double density(final double x) {
+ if (x < x0 ||
+ x > x3) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, x, x0, x3);
+ }
+ if (x <= x1) {
+ return p12 / (x1 - x0);
+ } else if (x <= x2) {
+ return 0.0;
+ } else if (x <= x3) {
+ return (1.0 - p12) / (x3 - x2);
+ }
+ return 0.0;
+ }
+
+ @Override
+ public double getNumericalMean() {
+ return ((x0 + x1) * p12 + (x2 + x3) * (1.0 - p12)) / 2.0;
+ }
+
+ @Override
+ public double getNumericalVariance() {
+ final double meanX = getNumericalMean();
+ final double meanX2;
+ meanX2 = ((x0 * x0 + x0 * x1 + x1 * x1) * p12 + (x2 * x2 + x2
+ * x3 + x3 * x3)
+ * (1.0 - p12)) / 3.0;
+ return meanX2 - meanX * meanX;
+ }
+
+ @Override
+ public double getSupportLowerBound() {
+ return x0;
+ }
+
+ @Override
+ public double getSupportUpperBound() {
+ return x3;
+ }
+
+ @Override
+ public boolean isSupportConnected() {
+ return false;
+ }
+
+ @Override
+ public double probability(final double x) {
+ throw new UnsupportedOperationException();
+ }
+ };
+ final double expected = x1;
+ final double actual = distribution.inverseCumulativeProbability(p12);
+ Assert.assertEquals("", expected, actual, 1e-8);
+ }
+
+ @Test
+ public void testDiscontinuous() {
+ final double x0 = 0.0;
+ final double x1 = 0.25;
+ final double x2 = 0.5;
+ final double x3 = 0.75;
+ final double x4 = 1.0;
+ final double p12 = 1.0 / 3.0;
+ final double p23 = 2.0 / 3.0;
+ final AbstractContinuousDistribution distribution;
+ distribution = new AbstractContinuousDistribution() {
+ @Override
+ public double cumulativeProbability(final double x) {
+ if (x < x0 ||
+ x > x4) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, x, x0, x4);
+ }
+ if (x <= x1) {
+ return p12 * (x - x0) / (x1 - x0);
+ } else if (x <= x2) {
+ return p12;
+ } else if (x <= x3) {
+ return p23;
+ } else {
+ return (1.0 - p23) * (x - x3) / (x4 - x3) + p23;
+ }
+ }
+
+ @Override
+ public double density(final double x) {
+ if (x < x0 ||
+ x > x4) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, x, x0, x4);
+ }
+ if (x <= x1) {
+ return p12 / (x1 - x0);
+ } else if (x <= x2) {
+ return 0.0;
+ } else if (x <= x3) {
+ return 0.0;
+ } else {
+ return (1.0 - p23) / (x4 - x3);
+ }
+ }
+
+ @Override
+ public double getNumericalMean() {
+ final UnivariateFunction f = new UnivariateFunction() {
+
+ @Override
+ public double value(final double x) {
+ return x * density(x);
+ }
+ };
+ final UnivariateIntegrator integrator = new RombergIntegrator();
+ return integrator.integrate(Integer.MAX_VALUE, f, x0, x4);
+ }
+
+ @Override
+ public double getNumericalVariance() {
+ final double meanX = getNumericalMean();
+ final UnivariateFunction f = new UnivariateFunction() {
+
+ @Override
+ public double value(final double x) {
+ return x * x * density(x);
+ }
+ };
+ final UnivariateIntegrator integrator = new RombergIntegrator();
+ final double meanX2 = integrator.integrate(Integer.MAX_VALUE,
+ f, x0, x4);
+ return meanX2 - meanX * meanX;
+ }
+
+ @Override
+ public double getSupportLowerBound() {
+ return x0;
+ }
+
+ @Override
+ public double getSupportUpperBound() {
+ return x4;
+ }
+
+ @Override
+ public boolean isSupportConnected() {
+ return false;
+ }
+
+ @Override
+ public double probability(final double x) {
+ throw new UnsupportedOperationException();
+ }
+ };
+ final double expected = x2;
+ final double actual = distribution.inverseCumulativeProbability(p23);
+ Assert.assertEquals("", expected, actual, 1e-8);
+ }
+}