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); + } +}
