http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DiscreteDistribution.java ---------------------------------------------------------------------- diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DiscreteDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DiscreteDistribution.java new file mode 100644 index 0000000..1cd4a35 --- /dev/null +++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DiscreteDistribution.java @@ -0,0 +1,163 @@ +/* + * 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; + +/** + * Interface for distributions on the integers. + */ +public interface DiscreteDistribution { + + /** + * For a random variable {@code X} whose values are distributed according to + * this distribution, this method returns {@code log(P(X = x))}, where + * {@code log} is the natural logarithm. In other words, this method + * represents the logarithm of the probability mass function (PMF) for the + * distribution. Note that due to the floating point precision and + * under/overflow issues, this method will for some distributions be more + * precise and faster than computing the logarithm of + * {@link #probability(int)}. + * + * @param x the point at which the PMF is evaluated + * @return the logarithm of the value of the probability mass function at {@code x} + */ + double logProbability(int x); + + /** + * For a random variable {@code X} whose values are distributed according + * to this distribution, this method returns {@code P(X = x)}. In other + * words, this method represents the probability mass function (PMF) + * for the distribution. + * + * @param x the point at which the PMF is evaluated + * @return the value of the probability mass function at {@code x} + */ + double probability(int x); + + /** + * For a random variable {@code X} whose values are distributed according + * to this distribution, this method returns {@code P(x0 < X <= x1)}. + * + * @param x0 the exclusive lower bound + * @param x1 the inclusive upper bound + * @return the probability that a random variable with this distribution + * will take a value between {@code x0} and {@code x1}, + * excluding the lower and including the upper endpoint + * @throws IllegalArgumentException if {@code x0 > x1} + */ + double probability(int x0, int x1); + + /** + * For a random variable {@code X} whose values are distributed according + * to this distribution, this method returns {@code P(X <= x)}. In other + * words, this method represents the (cumulative) distribution function + * (CDF) for this distribution. + * + * @param x the point at which the CDF is evaluated + * @return the probability that a random variable with this + * distribution takes a value less than or equal to {@code x} + */ + double cumulativeProbability(int x); + + /** + * Computes the quantile function of this distribution. + * For a random variable {@code X} distributed according to this distribution, + * the returned value is + * <ul> + * <li>{@code inf{x in Z | P(X<=x) >= p}} for {@code 0 < p <= 1},</li> + * <li>{@code inf{x in Z | P(X<=x) > 0}} for {@code p = 0}.</li> + * </ul> + * If the result exceeds the range of the data type {@code int}, + * then {@code Integer.MIN_VALUE} or {@code Integer.MAX_VALUE} is returned. + * + * @param p the cumulative probability + * @return the smallest {@code p}-quantile of this distribution + * (largest 0-quantile for {@code p = 0}) + * @throws IllegalArgumentException if {@code p < 0} or {@code p > 1} + */ + int inverseCumulativeProbability(double p); + + /** + * Use this method to get the numerical value of the mean of this + * distribution. + * + * @return the mean or {@code Double.NaN} if it is not defined + */ + double getNumericalMean(); + + /** + * Use this method to get the numerical value of the variance of this + * distribution. + * + * @return the variance (possibly {@code Double.POSITIVE_INFINITY} or + * {@code Double.NaN} if it is not defined) + */ + double getNumericalVariance(); + + /** + * Access the lower bound of the support. This method must return the same + * value as {@code inverseCumulativeProbability(0)}. In other words, this + * method must return + * <p>{@code inf {x in Z | P(X <= x) > 0}}.</p> + * + * @return lower bound of the support ({@code Integer.MIN_VALUE} + * for negative infinity) + */ + int getSupportLowerBound(); + + /** + * Access the upper bound of the support. This method must return the same + * value as {@code inverseCumulativeProbability(1)}. In other words, this + * method must return + * <p>{@code inf {x in R | P(X <= x) = 1}}.</p> + * + * @return upper bound of the support ({@code Integer.MAX_VALUE} + * for positive infinity) + */ + int getSupportUpperBound(); + + /** + * Use this method to get information about whether the support is + * connected, i.e. whether all integers between the lower and upper bound of + * the support are included in the support. + * + * @return whether the support is connected or not + */ + boolean isSupportConnected(); + + /** + * Creates a sampler. + * + * @param rng Generator of uniformly distributed numbers. + * @return a sampler that produces random numbers according this + * distribution. + */ + Sampler createSampler(UniformRandomProvider rng); + + /** + * Sampling functionality. + */ + interface Sampler { + /** + * Generates a random value sampled from this distribution. + * + * @return a random value. + */ + int sample(); + } +}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DistributionException.java ---------------------------------------------------------------------- diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DistributionException.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DistributionException.java new file mode 100644 index 0000000..31d7adc --- /dev/null +++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DistributionException.java @@ -0,0 +1,61 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.statistics.distribution; + +import java.text.MessageFormat; + +/** + * Package private exception class with constants for frequently used messages. + */ +class DistributionException extends IllegalArgumentException { + /** Error message for "too large" condition. */ + static final String TOO_LARGE = "{0} > {1}"; + /** Error message for "too small" condition. */ + static final String TOO_SMALL = "{0} < {1}"; + /** Error message for "out of range" condition. */ + static final String OUT_OF_RANGE = "Number {0} is out of range [{1}, {2}]"; + /** Error message for "out of range" condition. */ + static final String NEGATIVE = "Number {0} is negative"; + /** Error message for "mismatch" condition. */ + static final String MISMATCH = "Expected {1} but was {0}"; + /** Error message for "failed bracketing" condition. */ + static final String BRACKETING = "No bracketing: f({0})={1}, f({2})={3}"; + + /** Serializable version identifier. */ + private static final long serialVersionUID = 20180119L; + + /** Arguments for formatting the message. */ + private Object[] formatArguments; + + /** + * Create an exception where the message is constructed by applying + * the {@code format()} method from {@code java.text.MessageFormat}. + * + * @param message the exception message with replaceable parameters + * @param formatArguments the arguments for formatting the message + */ + DistributionException(String message, Object... formatArguments) { + super(message); + this.formatArguments = formatArguments; + } + + /** {@inheritDoc} */ + @Override + public String getMessage() { + return MessageFormat.format(super.getMessage(), formatArguments); + } +} http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ExponentialDistribution.java ---------------------------------------------------------------------- diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ExponentialDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ExponentialDistribution.java new file mode 100644 index 0000000..091c6a4 --- /dev/null +++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ExponentialDistribution.java @@ -0,0 +1,197 @@ +/* + * 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.AhrensDieterExponentialSampler; + +/** + * Implementation of the <a href="http://en.wikipedia.org/wiki/Exponential_distribution";>exponential distribution</a>. + */ +public class ExponentialDistribution extends AbstractContinuousDistribution { + /** The mean of this distribution. */ + private final double mean; + /** The logarithm of the mean, stored to reduce computing time. */ + private final double logMean; + + /** + * Creates a distribution. + * + * @param mean Mean of this distribution. + * @throws IllegalArgumentException if {@code mean <= 0}. + */ + public ExponentialDistribution(double mean) { + if (mean <= 0) { + throw new DistributionException(DistributionException.NEGATIVE, mean); + } + this.mean = mean; + logMean = Math.log(mean); + } + + /** + * Access the mean. + * + * @return the mean. + */ + public double getMean() { + return mean; + } + + /** {@inheritDoc} */ + @Override + public double density(double x) { + final double logDensity = logDensity(x); + return logDensity == Double.NEGATIVE_INFINITY ? 0 : Math.exp(logDensity); + } + + /** {@inheritDoc} **/ + @Override + public double logDensity(double x) { + if (x < 0) { + return Double.NEGATIVE_INFINITY; + } + return -x / mean - logMean; + } + + /** + * {@inheritDoc} + * + * The implementation of this method is based on: + * <ul> + * <li> + * <a href="http://mathworld.wolfram.com/ExponentialDistribution.html";> + * Exponential Distribution</a>, equation (1).</li> + * </ul> + */ + @Override + public double cumulativeProbability(double x) { + double ret; + if (x <= 0) { + ret = 0; + } else { + ret = 1 - Math.exp(-x / mean); + } + 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 == 1) { + ret = Double.POSITIVE_INFINITY; + } else { + ret = -mean * Math.log(1 - p); + } + + return ret; + } + + /** + * {@inheritDoc} + * + * For mean parameter {@code k}, the mean is {@code k}. + */ + @Override + public double getNumericalMean() { + return getMean(); + } + + /** + * {@inheritDoc} + * + * For mean parameter {@code k}, the variance is {@code k^2}. + */ + @Override + public double getNumericalVariance() { + final double m = getMean(); + return m * m; + } + + /** + * {@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 double getSupportLowerBound() { + return 0; + } + + /** + * {@inheritDoc} + * + * The upper bound of the support is always positive infinity + * no matter the mean parameter. + * + * @return upper bound of the support (always 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} + * + * <p>Sampling algorithm uses the + * <a href="http://www.jesus.ox.ac.uk/~clifford/a5/chap1/node5.html";> + * inversion method</a> to generate exponentially distributed + * random values from uniform deviates. + * </p> + */ + @Override + public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) { + return new ContinuousDistribution.Sampler() { + /** + * Exponential distribution sampler. + */ + private final ContinuousSampler sampler = + new AhrensDieterExponentialSampler(rng, mean); + + /**{@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/FDistribution.java ---------------------------------------------------------------------- diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/FDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/FDistribution.java new file mode 100644 index 0000000..a8b9890 --- /dev/null +++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/FDistribution.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.numbers.gamma.LogBeta; +import org.apache.commons.numbers.gamma.RegularizedBeta; + +/** + * Implementation of the F-distribution. + * + * @see <a href="http://en.wikipedia.org/wiki/F-distribution";>F-distribution (Wikipedia)</a> + * @see <a href="http://mathworld.wolfram.com/F-Distribution.html";>F-distribution (MathWorld)</a> + */ +public class FDistribution extends AbstractContinuousDistribution { + /** The numerator degrees of freedom. */ + private final double numeratorDegreesOfFreedom; + /** The numerator degrees of freedom. */ + private final double denominatorDegreesOfFreedom; + + /** + * Creates a distribution. + * + * @param numeratorDegreesOfFreedom Numerator degrees of freedom. + * @param denominatorDegreesOfFreedom Denominator degrees of freedom. + * @throws IllegalArgumentException if {@code numeratorDegreesOfFreedom <= 0} or + * {@code denominatorDegreesOfFreedom <= 0}. + */ + public FDistribution(double numeratorDegreesOfFreedom, + double denominatorDegreesOfFreedom) { + if (numeratorDegreesOfFreedom <= 0) { + throw new DistributionException(DistributionException.NEGATIVE, + numeratorDegreesOfFreedom); + } + if (denominatorDegreesOfFreedom <= 0) { + throw new DistributionException(DistributionException.NEGATIVE, + denominatorDegreesOfFreedom); + } + this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom; + this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom; + } + + /** + * {@inheritDoc} + */ + @Override + public double density(double x) { + return Math.exp(logDensity(x)); + } + + /** {@inheritDoc} **/ + @Override + public double logDensity(double x) { + final double nhalf = numeratorDegreesOfFreedom / 2; + final double mhalf = denominatorDegreesOfFreedom / 2; + final double logx = Math.log(x); + final double logn = Math.log(numeratorDegreesOfFreedom); + final double logm = Math.log(denominatorDegreesOfFreedom); + final double lognxm = Math.log(numeratorDegreesOfFreedom * x + + denominatorDegreesOfFreedom); + return nhalf * logn + nhalf * logx - logx + + mhalf * logm - nhalf * lognxm - mhalf * lognxm - + LogBeta.value(nhalf, mhalf); + } + + /** + * {@inheritDoc} + * + * The implementation of this method is based on + * <ul> + * <li> + * <a href="http://mathworld.wolfram.com/F-Distribution.html";> + * F-Distribution</a>, equation (4). + * </li> + * </ul> + */ + @Override + public double cumulativeProbability(double x) { + double ret; + if (x <= 0) { + ret = 0; + } else { + double n = numeratorDegreesOfFreedom; + double m = denominatorDegreesOfFreedom; + + ret = RegularizedBeta.value((n * x) / (m + n * x), + 0.5 * n, + 0.5 * m); + } + return ret; + } + + /** + * Access the numerator degrees of freedom. + * + * @return the numerator degrees of freedom. + */ + public double getNumeratorDegreesOfFreedom() { + return numeratorDegreesOfFreedom; + } + + /** + * Access the denominator degrees of freedom. + * + * @return the denominator degrees of freedom. + */ + public double getDenominatorDegreesOfFreedom() { + return denominatorDegreesOfFreedom; + } + + /** + * {@inheritDoc} + * + * For denominator degrees of freedom parameter {@code b}, the mean is + * <ul> + * <li>if {@code b > 2} then {@code b / (b - 2)},</li> + * <li>else undefined ({@code Double.NaN}). + * </ul> + */ + @Override + public double getNumericalMean() { + final double denominatorDF = getDenominatorDegreesOfFreedom(); + + if (denominatorDF > 2) { + return denominatorDF / (denominatorDF - 2); + } + + return Double.NaN; + } + + /** + * {@inheritDoc} + * + * For numerator degrees of freedom parameter {@code a} and denominator + * degrees of freedom parameter {@code b}, the variance is + * <ul> + * <li> + * if {@code b > 4} then + * {@code [2 * b^2 * (a + b - 2)] / [a * (b - 2)^2 * (b - 4)]}, + * </li> + * <li>else undefined ({@code Double.NaN}). + * </ul> + */ + @Override + public double getNumericalVariance() { + final double denominatorDF = getDenominatorDegreesOfFreedom(); + + if (denominatorDF > 4) { + final double numeratorDF = getNumeratorDegreesOfFreedom(); + final double denomDFMinusTwo = denominatorDF - 2; + + return (2 * (denominatorDF * denominatorDF) * (numeratorDF + denominatorDF - 2)) / + ((numeratorDF * (denomDFMinusTwo * denomDFMinusTwo) * (denominatorDF - 4))); + } + + return Double.NaN; + } + + /** + * {@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 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/GammaDistribution.java ---------------------------------------------------------------------- diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GammaDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GammaDistribution.java new file mode 100644 index 0000000..c13db5b --- /dev/null +++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GammaDistribution.java @@ -0,0 +1,354 @@ +/* + * 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.LanczosApproximation; +import org.apache.commons.numbers.gamma.RegularizedGamma; +import org.apache.commons.rng.UniformRandomProvider; +import org.apache.commons.rng.sampling.distribution.ContinuousSampler; +import org.apache.commons.rng.sampling.distribution.AhrensDieterMarsagliaTsangGammaSampler; + +/** + * Implementation of the <a href="http://en.wikipedia.org/wiki/Gamma_distribution";>Gamma distribution</a>. + */ +public class GammaDistribution extends AbstractContinuousDistribution { + /** Lanczos constant. */ + private static final double LANCZOS_G = LanczosApproximation.g(); + /** The shape parameter. */ + private final double shape; + /** The scale parameter. */ + private final double scale; + /** + * The constant value of {@code shape + g + 0.5}, where {@code g} is the + * Lanczos constant {@link LanczosApproximation#g()}. + */ + private final double shiftedShape; + /** + * The constant value of + * {@code shape / scale * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape)}, + * where {@code L(shape)} is the Lanczos approximation returned by + * {@link LanczosApproximation#value(double)}. This prefactor is used in + * {@link #density(double)}, when no overflow occurs with the natural + * calculation. + */ + private final double densityPrefactor1; + /** + * The constant value of + * {@code log(shape / scale * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape))}, + * where {@code L(shape)} is the Lanczos approximation returned by + * {@link LanczosApproximation#value(double)}. This prefactor is used in + * {@link #logDensity(double)}, when no overflow occurs with the natural + * calculation. + */ + private final double logDensityPrefactor1; + /** + * The constant value of + * {@code shape * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape)}, + * where {@code L(shape)} is the Lanczos approximation returned by + * {@link LanczosApproximation#value(double)}. This prefactor is used in + * {@link #density(double)}, when overflow occurs with the natural + * calculation. + */ + private final double densityPrefactor2; + /** + * The constant value of + * {@code log(shape * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape))}, + * where {@code L(shape)} is the Lanczos approximation returned by + * {@link LanczosApproximation#value(double)}. This prefactor is used in + * {@link #logDensity(double)}, when overflow occurs with the natural + * calculation. + */ + private final double logDensityPrefactor2; + /** + * Lower bound on {@code y = x / scale} for the selection of the computation + * method in {@link #density(double)}. For {@code y <= minY}, the natural + * calculation overflows. + */ + private final double minY; + /** + * Upper bound on {@code log(y)} ({@code y = x / scale}) for the selection + * of the computation method in {@link #density(double)}. For + * {@code log(y) >= maxLogY}, the natural calculation overflows. + */ + private final double maxLogY; + + /** + * Creates a distribution. + * + * @param shape the shape parameter + * @param scale the scale parameter + * @throws IllegalArgumentException if {@code shape <= 0} or + * {@code scale <= 0}. + */ + public GammaDistribution(double shape, + double scale) { + if (shape <= 0) { + throw new DistributionException(DistributionException.NEGATIVE, shape); + } + if (scale <= 0) { + throw new DistributionException(DistributionException.NEGATIVE, scale); + } + + this.shape = shape; + this.scale = scale; + this.shiftedShape = shape + LANCZOS_G + 0.5; + final double aux = Math.E / (2.0 * Math.PI * shiftedShape); + this.densityPrefactor2 = shape * Math.sqrt(aux) / LanczosApproximation.value(shape); + this.logDensityPrefactor2 = Math.log(shape) + 0.5 * Math.log(aux) - + Math.log(LanczosApproximation.value(shape)); + this.densityPrefactor1 = this.densityPrefactor2 / scale * + Math.pow(shiftedShape, -shape) * // XXX FastMath vs Math + Math.exp(shape + LANCZOS_G); + this.logDensityPrefactor1 = this.logDensityPrefactor2 - Math.log(scale) - + Math.log(shiftedShape) * shape + + shape + LANCZOS_G; + this.minY = shape + LANCZOS_G - Math.log(Double.MAX_VALUE); + this.maxLogY = Math.log(Double.MAX_VALUE) / (shape - 1.0); + } + + /** + * Returns the shape parameter of {@code this} distribution. + * + * @return the shape parameter + */ + public double getShape() { + return shape; + } + + /** + * Returns the scale parameter of {@code this} distribution. + * + * @return the scale parameter + */ + public double getScale() { + return scale; + } + + /** {@inheritDoc} */ + @Override + public double density(double x) { + /* The present method must return the value of + * + * 1 x a - x + * ---------- (-) exp(---) + * x Gamma(a) b b + * + * where a is the shape parameter, and b the scale parameter. + * Substituting the Lanczos approximation of Gamma(a) leads to the + * following expression of the density + * + * a e 1 y a + * - sqrt(------------------) ---- (-----------) exp(a - y + g), + * x 2 pi (a + g + 0.5) L(a) a + g + 0.5 + * + * where y = x / b. The above formula is the "natural" computation, which + * is implemented when no overflow is likely to occur. If overflow occurs + * with the natural computation, the following identity is used. It is + * based on the BOOST library + * http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_gamma/igamma.html + * Formula (15) needs adaptations, which are detailed below. + * + * y a + * (-----------) exp(a - y + g) + * a + g + 0.5 + * y - a - g - 0.5 y (g + 0.5) + * = exp(a log1pm(---------------) - ----------- + g), + * a + g + 0.5 a + g + 0.5 + * + * where log1pm(z) = log(1 + z) - z. Therefore, the value to be + * returned is + * + * a e 1 + * - sqrt(------------------) ---- + * x 2 pi (a + g + 0.5) L(a) + * y - a - g - 0.5 y (g + 0.5) + * * exp(a log1pm(---------------) - ----------- + g). + * a + g + 0.5 a + g + 0.5 + */ + if (x < 0) { + return 0; + } + final double y = x / scale; + if ((y <= minY) || (Math.log(y) >= maxLogY)) { + /* + * Overflow. + */ + final double aux1 = (y - shiftedShape) / shiftedShape; + final double aux2 = shape * (Math.log1p(aux1) - aux1); // XXX FastMath vs Math + final double aux3 = -y * (LANCZOS_G + 0.5) / shiftedShape + LANCZOS_G + aux2; + return densityPrefactor2 / x * Math.exp(aux3); + } + /* + * Natural calculation. + */ + return densityPrefactor1 * Math.exp(-y) * Math.pow(y, shape - 1); + } + + /** {@inheritDoc} **/ + @Override + public double logDensity(double x) { + /* + * see the comment in {@link #density(double)} for computation details + */ + if (x < 0) { + return Double.NEGATIVE_INFINITY; + } + final double y = x / scale; + if ((y <= minY) || (Math.log(y) >= maxLogY)) { + /* + * Overflow. + */ + final double aux1 = (y - shiftedShape) / shiftedShape; + final double aux2 = shape * (Math.log1p(aux1) - aux1); + final double aux3 = -y * (LANCZOS_G + 0.5) / shiftedShape + LANCZOS_G + aux2; + return logDensityPrefactor2 - Math.log(x) + aux3; + } + /* + * Natural calculation. + */ + return logDensityPrefactor1 - y + Math.log(y) * (shape - 1); + } + + /** + * {@inheritDoc} + * + * The implementation of this method is based on: + * <ul> + * <li> + * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html";> + * Chi-Squared Distribution</a>, equation (9). + * </li> + * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>. + * Belmont, CA: Duxbury Press. + * </li> + * </ul> + */ + @Override + public double cumulativeProbability(double x) { + double ret; + + if (x <= 0) { + ret = 0; + } else { + ret = RegularizedGamma.P.value(shape, x / scale); + } + + return ret; + } + + /** + * {@inheritDoc} + * + * For shape parameter {@code alpha} and scale parameter {@code beta}, the + * mean is {@code alpha * beta}. + */ + @Override + public double getNumericalMean() { + return shape * scale; + } + + /** + * {@inheritDoc} + * + * For shape parameter {@code alpha} and scale parameter {@code beta}, the + * variance is {@code alpha * beta^2}. + * + * @return {@inheritDoc} + */ + @Override + public double getNumericalVariance() { + return shape * scale * scale; + } + + /** + * {@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 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} + * + * <p> + * Sampling algorithms: + * <ul> + * <li> + * For {@code 0 < shape < 1}: + * <blockquote> + * Ahrens, J. H. and Dieter, U., + * <i>Computer methods for sampling from gamma, beta, Poisson and binomial distributions,</i> + * Computing, 12, 223-246, 1974. + * </blockquote> + * </li> + * <li> + * For {@code shape >= 1}: + * <blockquote> + * Marsaglia and Tsang, <i>A Simple Method for Generating + * Gamma Variables.</i> ACM Transactions on Mathematical Software, + * Volume 26 Issue 3, September, 2000. + * </blockquote> + * </li> + * </ul> + */ + @Override + public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) { + return new ContinuousDistribution.Sampler() { + /** + * Gamma distribution sampler. + */ + private final ContinuousSampler sampler = + new AhrensDieterMarsagliaTsangGammaSampler(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/GeometricDistribution.java ---------------------------------------------------------------------- diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GeometricDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GeometricDistribution.java new file mode 100644 index 0000000..fba2580 --- /dev/null +++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GeometricDistribution.java @@ -0,0 +1,160 @@ +/* + * 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 <a href="http://en.wikipedia.org/wiki/Geometric_distribution";>geometric distribution</a>. + */ +public class GeometricDistribution extends AbstractDiscreteDistribution { + /** The probability of success. */ + private final double probabilityOfSuccess; + /** {@code log(p)} where p is the probability of success. */ + private final double logProbabilityOfSuccess; + /** {@code log(1 - p)} where p is the probability of success. */ + private final double log1mProbabilityOfSuccess; + + /** + * Creates a geometric distribution. + * + * @param p Probability of success. + * @throws IllegalArgumentException if {@code p <= 0} or {@code p > 1}. + */ + public GeometricDistribution(double p) { + if (p <= 0 || p > 1) { + throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1); + } + + probabilityOfSuccess = p; + logProbabilityOfSuccess = Math.log(p); + log1mProbabilityOfSuccess = Math.log1p(-p); + } + + /** + * 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) { + if (x < 0) { + return 0.0; + } else { + return Math.exp(log1mProbabilityOfSuccess * x) * probabilityOfSuccess; + } + } + + /** {@inheritDoc} */ + @Override + public double logProbability(int x) { + if (x < 0) { + return Double.NEGATIVE_INFINITY; + } else { + return x * log1mProbabilityOfSuccess + logProbabilityOfSuccess; + } + } + + /** {@inheritDoc} */ + @Override + public double cumulativeProbability(int x) { + if (x < 0) { + return 0.0; + } else { + return -Math.expm1(log1mProbabilityOfSuccess * (x + 1)); + } + } + + /** + * {@inheritDoc} + * + * For probability parameter {@code p}, the mean is {@code (1 - p) / p}. + */ + @Override + public double getNumericalMean() { + return (1 - probabilityOfSuccess) / probabilityOfSuccess; + } + + /** + * {@inheritDoc} + * + * For probability parameter {@code p}, the variance is + * {@code (1 - p) / (p * p)}. + */ + @Override + public double getNumericalVariance() { + return (1 - probabilityOfSuccess) / (probabilityOfSuccess * probabilityOfSuccess); + } + + /** + * {@inheritDoc} + * + * The lower bound of the support is always 0. + * + * @return lower bound of the support (always 0) + */ + @Override + public int getSupportLowerBound() { + return 0; + } + + /** + * {@inheritDoc} + * + * The upper bound of the support is infinite (which we approximate as + * {@code Integer.MAX_VALUE}). + * + * @return upper bound of the support (always Integer.MAX_VALUE) + */ + @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 int inverseCumulativeProbability(double p) { + if (p < 0 || + p > 1) { + throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1); + } + if (p == 1) { + return Integer.MAX_VALUE; + } + if (p == 0) { + return 0; + } + return Math.max(0, (int) Math.ceil(Math.log1p(-p)/log1mProbabilityOfSuccess-1)); + } +} http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GumbelDistribution.java ---------------------------------------------------------------------- diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GumbelDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GumbelDistribution.java new file mode 100644 index 0000000..8898b5e --- /dev/null +++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GumbelDistribution.java @@ -0,0 +1,128 @@ +/* + * 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; + +/** + * This class implements the <a href="http://en.wikipedia.org/wiki/Gumbel_distribution";>Gumbel distribution</a>. + */ +public class GumbelDistribution extends AbstractContinuousDistribution { + /** π<sup>2</sup>/6. */ + private static final double PI_SQUARED_OVER_SIX = Math.PI * Math.PI / 6; + /** + * <a href="http://mathworld.wolfram.com/Euler-MascheroniConstantApproximations.html";> + * Approximation of Euler's constant</a>. + */ + private static final double EULER = Math.PI / (2 * Math.E); + /** Location parameter. */ + private final double mu; + /** Scale parameter. */ + private final double beta; + + /** + * Creates a distribution. + * + * @param mu location parameter + * @param beta scale parameter (must be positive) + * @throws IllegalArgumenException if {@code beta <= 0} + */ + public GumbelDistribution(double mu, + double beta) { + if (beta <= 0) { + throw new DistributionException(DistributionException.NEGATIVE, beta); + } + + this.beta = beta; + this.mu = mu; + } + + /** + * Gets the location parameter. + * + * @return the location parameter. + */ + public double getLocation() { + return mu; + } + + /** + * Gets the scale parameter. + * + * @return the scale parameter. + */ + public double getScale() { + return beta; + } + + /** {@inheritDoc} */ + @Override + public double density(double x) { + final double z = (x - mu) / beta; + final double t = Math.exp(-z); + return Math.exp(-z - t) / beta; + } + + /** {@inheritDoc} */ + @Override + public double cumulativeProbability(double x) { + final double z = (x - mu) / beta; + return Math.exp(-Math.exp(-z)); + } + + /** {@inheritDoc} */ + @Override + public double inverseCumulativeProbability(double p) { + if (p < 0 || p > 1) { + throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1); + } else if (p == 0) { + return Double.NEGATIVE_INFINITY; + } else if (p == 1) { + return Double.POSITIVE_INFINITY; + } + return mu - Math.log(-Math.log(p)) * beta; + } + + /** {@inheritDoc} */ + @Override + public double getNumericalMean() { + return mu + EULER * beta; + } + + /** {@inheritDoc} */ + @Override + public double getNumericalVariance() { + return PI_SQUARED_OVER_SIX * beta * beta; + } + + /** {@inheritDoc} */ + @Override + public double getSupportLowerBound() { + return Double.NEGATIVE_INFINITY; + } + + /** {@inheritDoc} */ + @Override + public double getSupportUpperBound() { + return Double.POSITIVE_INFINITY; + } + + /** {@inheritDoc} */ + @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/HypergeometricDistribution.java ---------------------------------------------------------------------- diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/HypergeometricDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/HypergeometricDistribution.java new file mode 100644 index 0000000..732a253 --- /dev/null +++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/HypergeometricDistribution.java @@ -0,0 +1,293 @@ +/* + * 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 <a href="http://en.wikipedia.org/wiki/Hypergeometric_distribution";>hypergeometric distribution</a>. + */ +public class HypergeometricDistribution extends AbstractDiscreteDistribution { + /** The number of successes in the population. */ + private final int numberOfSuccesses; + /** The population size. */ + private final int populationSize; + /** The sample size. */ + private final int sampleSize; + + /** + * Creates a new hypergeometric distribution. + * + * @param populationSize Population size. + * @param numberOfSuccesses Number of successes in the population. + * @param sampleSize Sample size. + * @throws IllegalArgumentException if {@code numberOfSuccesses < 0}, or + * {@code populationSize <= 0} or {@code numberOfSuccesses > populationSize}, + * or {@code sampleSize > populationSize}. + */ + public HypergeometricDistribution(int populationSize, + int numberOfSuccesses, + int sampleSize) { + if (populationSize <= 0) { + throw new DistributionException(DistributionException.NEGATIVE, + populationSize); + } + if (numberOfSuccesses < 0) { + throw new DistributionException(DistributionException.NEGATIVE, + numberOfSuccesses); + } + if (sampleSize < 0) { + throw new DistributionException(DistributionException.NEGATIVE, + sampleSize); + } + + if (numberOfSuccesses > populationSize) { + throw new DistributionException(DistributionException.TOO_LARGE, + numberOfSuccesses, populationSize); + } + if (sampleSize > populationSize) { + throw new DistributionException(DistributionException.TOO_LARGE, + sampleSize, populationSize); + } + + this.numberOfSuccesses = numberOfSuccesses; + this.populationSize = populationSize; + this.sampleSize = sampleSize; + } + + /** {@inheritDoc} */ + @Override + public double cumulativeProbability(int x) { + double ret; + + int[] domain = getDomain(populationSize, numberOfSuccesses, sampleSize); + if (x < domain[0]) { + ret = 0.0; + } else if (x >= domain[1]) { + ret = 1.0; + } else { + ret = innerCumulativeProbability(domain[0], x, 1); + } + + return ret; + } + + /** + * Return the domain for the given hypergeometric distribution parameters. + * + * @param n Population size. + * @param m Number of successes in the population. + * @param k Sample size. + * @return a two element array containing the lower and upper bounds of the + * hypergeometric distribution. + */ + private int[] getDomain(int n, int m, int k) { + return new int[] { getLowerDomain(n, m, k), getUpperDomain(m, k) }; + } + + /** + * Return the lowest domain value for the given hypergeometric distribution + * parameters. + * + * @param n Population size. + * @param m Number of successes in the population. + * @param k Sample size. + * @return the lowest domain value of the hypergeometric distribution. + */ + private int getLowerDomain(int n, int m, int k) { + return Math.max(0, m - (n - k)); + } + + /** + * Access the number of successes. + * + * @return the number of successes. + */ + public int getNumberOfSuccesses() { + return numberOfSuccesses; + } + + /** + * Access the population size. + * + * @return the population size. + */ + public int getPopulationSize() { + return populationSize; + } + + /** + * Access the sample size. + * + * @return the sample size. + */ + public int getSampleSize() { + return sampleSize; + } + + /** + * Return the highest domain value for the given hypergeometric distribution + * parameters. + * + * @param m Number of successes in the population. + * @param k Sample size. + * @return the highest domain value of the hypergeometric distribution. + */ + private int getUpperDomain(int m, int k) { + return Math.min(k, m); + } + + /** {@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; + + int[] domain = getDomain(populationSize, numberOfSuccesses, sampleSize); + if (x < domain[0] || x > domain[1]) { + ret = Double.NEGATIVE_INFINITY; + } else { + double p = (double) sampleSize / (double) populationSize; + double q = (double) (populationSize - sampleSize) / (double) populationSize; + double p1 = SaddlePointExpansion.logBinomialProbability(x, + numberOfSuccesses, p, q); + double p2 = + SaddlePointExpansion.logBinomialProbability(sampleSize - x, + populationSize - numberOfSuccesses, p, q); + double p3 = + SaddlePointExpansion.logBinomialProbability(sampleSize, populationSize, p, q); + ret = p1 + p2 - p3; + } + + return ret; + } + + /** + * For this distribution, {@code X}, this method returns {@code P(X >= x)}. + * + * @param x Value at which the CDF is evaluated. + * @return the upper tail CDF for this distribution. + * @since 1.1 + */ + public double upperCumulativeProbability(int x) { + double ret; + + final int[] domain = getDomain(populationSize, numberOfSuccesses, sampleSize); + if (x <= domain[0]) { + ret = 1.0; + } else if (x > domain[1]) { + ret = 0.0; + } else { + ret = innerCumulativeProbability(domain[1], x, -1); + } + + return ret; + } + + /** + * For this distribution, {@code X}, this method returns + * {@code P(x0 <= X <= x1)}. + * This probability is computed by summing the point probabilities for the + * values {@code x0, x0 + 1, x0 + 2, ..., x1}, in the order directed by + * {@code dx}. + * + * @param x0 Inclusive lower bound. + * @param x1 Inclusive upper bound. + * @param dx Direction of summation (1 indicates summing from x0 to x1, and + * 0 indicates summing from x1 to x0). + * @return {@code P(x0 <= X <= x1)}. + */ + private double innerCumulativeProbability(int x0, int x1, int dx) { + double ret = probability(x0); + while (x0 != x1) { + x0 += dx; + ret += probability(x0); + } + return ret; + } + + /** + * {@inheritDoc} + * + * For population size {@code N}, number of successes {@code m}, and sample + * size {@code n}, the mean is {@code n * m / N}. + */ + @Override + public double getNumericalMean() { + return getSampleSize() * (getNumberOfSuccesses() / (double) getPopulationSize()); + } + + /** + * {@inheritDoc} + * + * For population size {@code N}, number of successes {@code m}, and sample + * size {@code n}, the variance is + * {@code (n * m * (N - n) * (N - m)) / (N^2 * (N - 1))}. + */ + @Override + public double getNumericalVariance() { + final double N = getPopulationSize(); + final double m = getNumberOfSuccesses(); + final double n = getSampleSize(); + return (n * m * (N - n) * (N - m)) / (N * N * (N - 1)); + } + + /** + * {@inheritDoc} + * + * For population size {@code N}, number of successes {@code m}, and sample + * size {@code n}, the lower bound of the support is + * {@code max(0, n + m - N)}. + * + * @return lower bound of the support + */ + @Override + public int getSupportLowerBound() { + return Math.max(0, + getSampleSize() + getNumberOfSuccesses() - getPopulationSize()); + } + + /** + * {@inheritDoc} + * + * For number of successes {@code m} and sample size {@code n}, the upper + * bound of the support is {@code min(m, n)}. + * + * @return upper bound of the support + */ + @Override + public int getSupportUpperBound() { + return Math.min(getNumberOfSuccesses(), getSampleSize()); + } + + /** + * {@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/LaplaceDistribution.java ---------------------------------------------------------------------- diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LaplaceDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LaplaceDistribution.java new file mode 100644 index 0000000..0d1a8bf --- /dev/null +++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LaplaceDistribution.java @@ -0,0 +1,132 @@ +/* + * 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; + +/** + * This class implements the Laplace distribution. + * + * @see <a href="http://en.wikipedia.org/wiki/Laplace_distribution";>Laplace distribution (Wikipedia)</a> + * + * @since 3.4 + */ +public class LaplaceDistribution extends AbstractContinuousDistribution { + + /** Serializable version identifier. */ + private static final long serialVersionUID = 20160311L; + + /** The location parameter. */ + private final double mu; + /** The scale parameter. */ + private final double beta; + + /** + * Creates a distribution. + * + * @param mu location parameter + * @param beta scale parameter (must be positive) + * @throws IllegalArgumentException if {@code beta <= 0} + */ + public LaplaceDistribution(double mu, + double beta) { + if (beta <= 0.0) { + throw new DistributionException(DistributionException.NEGATIVE, beta); + } + + this.mu = mu; + this.beta = beta; + } + + /** + * Access the location parameter, {@code mu}. + * + * @return the location parameter. + */ + public double getLocation() { + return mu; + } + + /** + * Access the scale parameter, {@code beta}. + * + * @return the scale parameter. + */ + public double getScale() { + return beta; + } + + /** {@inheritDoc} */ + @Override + public double density(double x) { + return Math.exp(-Math.abs(x - mu) / beta) / (2.0 * beta); + } + + /** {@inheritDoc} */ + @Override + public double cumulativeProbability(double x) { + if (x <= mu) { + return Math.exp((x - mu) / beta) / 2.0; + } else { + return 1.0 - Math.exp((mu - x) / beta) / 2.0; + } + } + + /** {@inheritDoc} */ + @Override + public double inverseCumulativeProbability(double p) { + if (p < 0 || + p > 1) { + throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1); + } else if (p == 0) { + return Double.NEGATIVE_INFINITY; + } else if (p == 1) { + return Double.POSITIVE_INFINITY; + } + double x = (p > 0.5) ? -Math.log(2.0 - 2.0 * p) : Math.log(2.0 * p); + return mu + beta * x; + } + + /** {@inheritDoc} */ + @Override + public double getNumericalMean() { + return mu; + } + + /** {@inheritDoc} */ + @Override + public double getNumericalVariance() { + return 2.0 * beta * beta; + } + + /** {@inheritDoc} */ + @Override + public double getSupportLowerBound() { + return Double.NEGATIVE_INFINITY; + } + + /** {@inheritDoc} */ + @Override + public double getSupportUpperBound() { + return Double.POSITIVE_INFINITY; + } + + /** {@inheritDoc} */ + @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/LevyDistribution.java ---------------------------------------------------------------------- diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LevyDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LevyDistribution.java new file mode 100644 index 0000000..d16da8d --- /dev/null +++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LevyDistribution.java @@ -0,0 +1,161 @@ +/* + * 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.InverseErfc; + +/** + * This class implements the <a href="http://en.wikipedia.org/wiki/L%C3%A9vy_distribution";> + * Lévy distribution</a>. + */ +public class LevyDistribution extends AbstractContinuousDistribution { + /** Location parameter. */ + private final double mu; + /** Scale parameter. */ + private final double c; + /** Half of c (for calculations). */ + private final double halfC; + + /** + * Creates a distribution. + * + * @param mu location + * @param c scale parameter + */ + public LevyDistribution(final double mu, + final double c) { + this.mu = mu; + this.c = c; + this.halfC = 0.5 * c; + } + + /** {@inheritDoc} + * <p> + * From Wikipedia: The probability density function of the Lévy distribution + * over the domain is + * </p> + * <div style="white-space: pre"><code> + * f(x; μ, c) = √(c / 2π) * e<sup>-c / 2 (x - μ)</sup> / (x - μ)<sup>3/2</sup> + * </code></div> + * <p> + * For this distribution, {@code X}, this method returns {@code P(X < x)}. + * If {@code x} is less than location parameter μ, {@code Double.NaN} is + * returned, as in these cases the distribution is not defined. + * </p> + */ + @Override + public double density(final double x) { + if (x < mu) { + return Double.NaN; + } + + final double delta = x - mu; + final double f = halfC / delta; + return Math.sqrt(f / Math.PI) * Math.exp(-f) /delta; + } + + /** {@inheritDoc} + * + * See documentation of {@link #density(double)} for computation details. + */ + @Override + public double logDensity(double x) { + if (x < mu) { + return Double.NaN; + } + + final double delta = x - mu; + final double f = halfC / delta; + return 0.5 * Math.log(f / Math.PI) - f - Math.log(delta); + } + + /** {@inheritDoc} + * <p> + * From Wikipedia: the cumulative distribution function is + * </p> + * <pre> + * f(x; u, c) = erfc (√ (c / 2 (x - u ))) + * </pre> + */ + @Override + public double cumulativeProbability(final double x) { + if (x < mu) { + return Double.NaN; + } + return Erfc.value(Math.sqrt(halfC / (x - mu))); + } + + /** {@inheritDoc} */ + @Override + public double inverseCumulativeProbability(final double p) { + if (p < 0 || + p > 1) { + throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1); + } + final double t = InverseErfc.value(p); + return mu + halfC / (t * t); + } + + /** + * Gets the scale parameter of the distribution. + * + * @return scale parameter of the distribution + */ + public double getScale() { + return c; + } + + /** + * Gets the location parameter of the distribution. + * + * @return location parameter of the distribution + */ + public double getLocation() { + return mu; + } + + /** {@inheritDoc} */ + @Override + public double getNumericalMean() { + return Double.POSITIVE_INFINITY; + } + + /** {@inheritDoc} */ + @Override + public double getNumericalVariance() { + return Double.POSITIVE_INFINITY; + } + + /** {@inheritDoc} */ + @Override + public double getSupportLowerBound() { + return mu; + } + + /** {@inheritDoc} */ + @Override + public double getSupportUpperBound() { + return Double.POSITIVE_INFINITY; + } + + /** {@inheritDoc} */ + @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/LogNormalDistribution.java ---------------------------------------------------------------------- diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LogNormalDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LogNormalDistribution.java new file mode 100644 index 0000000..25bdd33 --- /dev/null +++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LogNormalDistribution.java @@ -0,0 +1,266 @@ +/* + * 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.Erf; +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.LogNormalSampler; +import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler; + +/** + * Implementation of the <a href="http://en.wikipedia.org/wiki/Log-normal_distribution";>log-normal distribution</a>. + * + * <p> + * <strong>Parameters:</strong> + * {@code X} is log-normally distributed if its natural logarithm {@code log(X)} + * is normally distributed. The probability distribution function of {@code X} + * is given by (for {@code x > 0}) + * </p> + * <p> + * {@code exp(-0.5 * ((ln(x) - m) / s)^2) / (s * sqrt(2 * pi) * x)} + * </p> + * <ul> + * <li>{@code m} is the <em>scale</em> parameter: this is the mean of the + * normally distributed natural logarithm of this distribution,</li> + * <li>{@code s} is the <em>shape</em> parameter: this is the standard + * deviation of the normally distributed natural logarithm of this + * distribution. + * </ul> + */ +public class LogNormalDistribution extends AbstractContinuousDistribution { + /** √(2 π) */ + private static final double SQRT2PI = Math.sqrt(2 * Math.PI); + /** √(2) */ + private static final double SQRT2 = Math.sqrt(2); + /** The scale parameter of this distribution. */ + private final double scale; + /** The shape parameter of this distribution. */ + private final double shape; + /** The value of {@code log(shape) + 0.5 * log(2*PI)} stored for faster computation. */ + private final double logShapePlusHalfLog2Pi; + + /** + * Creates a log-normal distribution, where the mean and standard deviation + * of the {@link NormalDistribution normally distributed} natural + * logarithm of the log-normal distribution are equal to zero and one + * respectively. In other words, the scale of the returned distribution is + * {@code 0}, while its shape is {@code 1}. + */ + public LogNormalDistribution() { + this(0, 1); + } + + /** + * Creates a log-normal distribution. + * + * @param scale Scale parameter of this distribution. + * @param shape Shape parameter of this distribution. + * @throws IllegalArgumentException if {@code shape <= 0}. + */ + public LogNormalDistribution(double scale, + double shape) { + if (shape <= 0) { + throw new DistributionException(DistributionException.NEGATIVE, shape); + } + + this.scale = scale; + this.shape = shape; + this.logShapePlusHalfLog2Pi = Math.log(shape) + 0.5 * Math.log(2 * Math.PI); + } + + /** + * 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} + * + * For scale {@code m}, and shape {@code s} of this distribution, the PDF + * is given by + * <ul> + * <li>{@code 0} if {@code x <= 0},</li> + * <li>{@code exp(-0.5 * ((ln(x) - m) / s)^2) / (s * sqrt(2 * pi) * x)} + * otherwise.</li> + * </ul> + */ + @Override + public double density(double x) { + if (x <= 0) { + return 0; + } + final double x0 = Math.log(x) - scale; + final double x1 = x0 / shape; + return Math.exp(-0.5 * x1 * x1) / (shape * SQRT2PI * x); + } + + /** {@inheritDoc} + * + * See documentation of {@link #density(double)} for computation details. + */ + @Override + public double logDensity(double x) { + if (x <= 0) { + return Double.NEGATIVE_INFINITY; + } + final double logX = Math.log(x); + final double x0 = logX - scale; + final double x1 = x0 / shape; + return -0.5 * x1 * x1 - (logShapePlusHalfLog2Pi + logX); + } + + /** + * {@inheritDoc} + * + * For scale {@code m}, and shape {@code s} of this distribution, the CDF + * is given by + * <ul> + * <li>{@code 0} if {@code x <= 0},</li> + * <li>{@code 0} if {@code ln(x) - m < 0} and {@code m - ln(x) > 40 * s}, as + * in these cases the actual value is within {@code Double.MIN_VALUE} of 0, + * <li>{@code 1} if {@code ln(x) - m >= 0} and {@code ln(x) - m > 40 * s}, + * as in these cases the actual value is within {@code Double.MIN_VALUE} of + * 1,</li> + * <li>{@code 0.5 + 0.5 * erf((ln(x) - m) / (s * sqrt(2))} otherwise.</li> + * </ul> + */ + @Override + public double cumulativeProbability(double x) { + if (x <= 0) { + return 0; + } + final double dev = Math.log(x) - scale; + if (Math.abs(dev) > 40 * shape) { + return dev < 0 ? 0.0d : 1.0d; + } + return 0.5 + 0.5 * Erf.value(dev / (shape * SQRT2)); + } + + /** {@inheritDoc} */ + @Override + public double probability(double x0, + double x1) { + if (x0 > x1) { + throw new DistributionException(DistributionException.TOO_LARGE, + x0, x1); + } + if (x0 <= 0 || x1 <= 0) { + return super.probability(x0, x1); + } + final double denom = shape * SQRT2; + final double v0 = (Math.log(x0) - scale) / denom; + final double v1 = (Math.log(x1) - scale) / denom; + return 0.5 * ErfDifference.value(v0, v1); + } + + /** + * {@inheritDoc} + * + * For scale {@code m} and shape {@code s}, the mean is + * {@code exp(m + s^2 / 2)}. + */ + @Override + public double getNumericalMean() { + double s = shape; + return Math.exp(scale + (s * s / 2)); + } + + /** + * {@inheritDoc} + * + * For scale {@code m} and shape {@code s}, the variance is + * {@code (exp(s^2) - 1) * exp(2 * m + s^2)}. + */ + @Override + public double getNumericalVariance() { + final double s = shape; + final double ss = s * s; + return (Math.expm1(ss)) * Math.exp(2 * scale + ss); + } + + /** + * {@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; + } + + /** {@inheritDoc} */ + @Override + public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) { + return new ContinuousDistribution.Sampler() { + /** + * Log normal distribution sampler. + */ + private final ContinuousSampler sampler = + new LogNormalSampler(new ZigguratNormalizedGaussianSampler(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/LogisticDistribution.java ---------------------------------------------------------------------- diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LogisticDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LogisticDistribution.java new file mode 100644 index 0000000..28a6657 --- /dev/null +++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LogisticDistribution.java @@ -0,0 +1,128 @@ +/* + * 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 <a href="http://en.wikipedia.org/wiki/Logistic_distribution";>Logistic distribution</a>. + */ +public class LogisticDistribution extends AbstractContinuousDistribution { + /** π<sup>2</sup>/3. */ + private static final double PI_SQUARED_OVER_THREE = Math.PI * Math.PI / 3; + /** Location parameter. */ + private final double mu; + /** Scale parameter. */ + private final double scale; + /** Inverse of "scale". */ + private final double oneOverScale; + + /** + * Creates a distribution. + * + * @param mu Location parameter. + * @param scale Scale parameter (must be positive). + * @throws IllegalArgumentException if {@code scale <= 0}. + */ + public LogisticDistribution(double mu, + double scale) { + if (scale <= 0) { + throw new DistributionException(DistributionException.NEGATIVE, + scale); + } + + this.mu = mu; + this.scale = scale; + this.oneOverScale = 1 / scale; + } + + /** + * Gets the location parameter. + * + * @return the location parameter. + */ + public double getLocation() { + return mu; + } + + /** + * Gets the scale parameter. + * + * @return the scale parameter. + */ + public double getScale() { + return scale; + } + + /** {@inheritDoc} */ + @Override + public double density(double x) { + final double z = oneOverScale * (x - mu); + final double v = Math.exp(-z); + return oneOverScale * v / ((1 + v) * (1 + v)); + } + + /** {@inheritDoc} */ + @Override + public double cumulativeProbability(double x) { + final double z = oneOverScale * (x - mu); + return 1 / (1 + Math.exp(-z)); + } + + /** {@inheritDoc} */ + @Override + public double inverseCumulativeProbability(double p) { + if (p < 0 || + p > 1) { + throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1); + } else if (p == 0) { + return 0; + } else if (p == 1) { + return Double.POSITIVE_INFINITY; + } else { + return scale * Math.log(p / (1 - p)) + mu; + } + } + + /** {@inheritDoc} */ + @Override + public double getNumericalMean() { + return mu; + } + + /** {@inheritDoc} */ + @Override + public double getNumericalVariance() { + return oneOverScale * oneOverScale * PI_SQUARED_OVER_THREE; + } + + /** {@inheritDoc} */ + @Override + public double getSupportLowerBound() { + return Double.NEGATIVE_INFINITY; + } + + /** {@inheritDoc} */ + @Override + public double getSupportUpperBound() { + return Double.POSITIVE_INFINITY; + } + + /** {@inheritDoc} */ + @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/NakagamiDistribution.java ---------------------------------------------------------------------- diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/NakagamiDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/NakagamiDistribution.java new file mode 100644 index 0000000..9bf7d2f --- /dev/null +++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/NakagamiDistribution.java @@ -0,0 +1,117 @@ +/* + * 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.Gamma; +import org.apache.commons.numbers.gamma.RegularizedGamma; + +/** + * This class implements the <a href="http://en.wikipedia.org/wiki/Nakagami_distribution";>Nakagami distribution</a>. + */ +public class NakagamiDistribution extends AbstractContinuousDistribution { + /** The shape parameter. */ + private final double mu; + /** The scale parameter. */ + private final double omega; + + /** + * Creates a distribution. + * + * @param mu shape parameter + * @param omega scale parameter (must be positive) + * @throws IllegalArgumentException if {@code mu < 0.5} or if + * {@code omega <= 0}. + */ + public NakagamiDistribution(double mu, + double omega) { + if (mu < 0.5) { + throw new DistributionException(DistributionException.TOO_SMALL, mu, 0.5); + } + if (omega <= 0) { + throw new DistributionException(DistributionException.NEGATIVE, omega); + } + + this.mu = mu; + this.omega = omega; + } + + /** + * Access the shape parameter, {@code mu}. + * + * @return the shape parameter. + */ + public double getShape() { + return mu; + } + + /** + * Access the scale parameter, {@code omega}. + * + * @return the scale parameter. + */ + public double getScale() { + return omega; + } + + /** {@inheritDoc} */ + @Override + public double density(double x) { + if (x <= 0) { + return 0.0; + } + return 2.0 * Math.pow(mu, mu) / (Gamma.value(mu) * Math.pow(omega, mu)) * + Math.pow(x, 2 * mu - 1) * Math.exp(-mu * x * x / omega); + } + + /** {@inheritDoc} */ + @Override + public double cumulativeProbability(double x) { + return RegularizedGamma.P.value(mu, mu * x * x / omega); + } + + /** {@inheritDoc} */ + @Override + public double getNumericalMean() { + return Gamma.value(mu + 0.5) / Gamma.value(mu) * Math.sqrt(omega / mu); + } + + /** {@inheritDoc} */ + @Override + public double getNumericalVariance() { + double v = Gamma.value(mu + 0.5) / Gamma.value(mu); + return omega * (1 - 1 / mu * v * v); + } + + /** {@inheritDoc} */ + @Override + public double getSupportLowerBound() { + return 0; + } + + /** {@inheritDoc} */ + @Override + public double getSupportUpperBound() { + return Double.POSITIVE_INFINITY; + } + + /** {@inheritDoc} */ + @Override + public boolean isSupportConnected() { + return true; + } + +}
