From 32354a1039ceccd4a6a71baea3329b27dff7a124 Mon Sep 17 00:00:00 2001 From: Sebastien Brisard Date: Sat, 26 Nov 2011 06:17:49 +0000 Subject: [PATCH] - Merged ExponentialDistribution and ExponentialDistributionImpl (MATH-711). - Merged FDistribution and FDistributionImpl (MATH-711). - Merged GammaDistribution and GammaDistributionImpl (MATH-711). git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1206399 13f79535-47bb-0310-9956-ffa450edef68 --- .../distribution/ChiSquaredDistribution.java | 428 +++++++++--------- .../distribution/ExponentialDistribution.java | 232 +++++++++- .../ExponentialDistributionImpl.java | 279 ------------ .../math/distribution/FDistribution.java | 268 ++++++++++- .../math/distribution/FDistributionImpl.java | 318 ------------- .../math/distribution/GammaDistribution.java | 259 ++++++++++- .../distribution/GammaDistributionImpl.java | 297 ------------ .../commons/math/random/RandomDataImpl.java | 8 +- .../math/stat/inference/OneWayAnovaImpl.java | 3 +- .../ExponentialDistributionTest.java | 16 +- .../math/distribution/FDistributionTest.java | 26 +- .../distribution/GammaDistributionTest.java | 22 +- .../commons/math/random/RandomDataTest.java | 16 +- 13 files changed, 977 insertions(+), 1195 deletions(-) delete mode 100644 src/main/java/org/apache/commons/math/distribution/ExponentialDistributionImpl.java delete mode 100644 src/main/java/org/apache/commons/math/distribution/FDistributionImpl.java delete mode 100644 src/main/java/org/apache/commons/math/distribution/GammaDistributionImpl.java diff --git a/src/main/java/org/apache/commons/math/distribution/ChiSquaredDistribution.java b/src/main/java/org/apache/commons/math/distribution/ChiSquaredDistribution.java index 15c3cc728..efcd955f5 100644 --- a/src/main/java/org/apache/commons/math/distribution/ChiSquaredDistribution.java +++ b/src/main/java/org/apache/commons/math/distribution/ChiSquaredDistribution.java @@ -1,214 +1,214 @@ -/* - * 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.math.distribution; - -import java.io.Serializable; - - -/** - * Implementation of the chi-squared distribution. - * - * @see Chi-squared distribution (Wikipedia) - * @see Chi-squared Distribution (MathWorld) - * @version $Id$ - */ -public class ChiSquaredDistribution - extends AbstractContinuousDistribution - implements Serializable { - /** - * Default inverse cumulative probability accuracy - * @since 2.1 - */ - public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; - /** Serializable version identifier */ - private static final long serialVersionUID = -8352658048349159782L; - /** Internal Gamma distribution. */ - private final GammaDistribution gamma; - /** Inverse cumulative probability accuracy */ - private final double solverAbsoluteAccuracy; - - /** - * Create a Chi-Squared distribution with the given degrees of freedom. - * - * @param degreesOfFreedom Degrees of freedom. - */ - public ChiSquaredDistribution(double degreesOfFreedom) { - this(degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); - } - - /** - * Create a Chi-Squared distribution with the given degrees of freedom and - * inverse cumulative probability accuracy. - * - * @param degreesOfFreedom Degrees of freedom. - * @param inverseCumAccuracy the maximum absolute error in inverse - * cumulative probability estimates (defaults to - * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). - * @since 2.1 - */ - public ChiSquaredDistribution(double degreesOfFreedom, - double inverseCumAccuracy) { - gamma = new GammaDistributionImpl(degreesOfFreedom / 2, 2); - solverAbsoluteAccuracy = inverseCumAccuracy; - } - - /** - * Access the number of degrees of freedom. - * - * @return the degrees of freedom. - */ - public double getDegreesOfFreedom() { - return gamma.getAlpha() * 2.0; - } - - /** {@inheritDoc} */ - public double density(double x) { - return gamma.density(x); - } - - /** {@inheritDoc} */ - public double cumulativeProbability(double x) { - return gamma.cumulativeProbability(x); - } - - /** - * {@inheritDoc} - * - * Returns {@code 0} when {@code p == 0} and - * {@code Double.POSITIVE_INFINITY} when {@code p == 1}. - */ - @Override - public double inverseCumulativeProbability(final double p) { - if (p == 0) { - return 0d; - } - if (p == 1) { - return Double.POSITIVE_INFINITY; - } - return super.inverseCumulativeProbability(p); - } - - /** {@inheritDoc} */ - @Override - protected double getDomainLowerBound(double p) { - return Double.MIN_VALUE * gamma.getBeta(); - } - - /** {@inheritDoc} */ - @Override - protected double getDomainUpperBound(double p) { - // NOTE: chi squared is skewed to the left - // NOTE: therefore, P(X < μ) > .5 - - double ret; - - if (p < .5) { - // use mean - ret = getDegreesOfFreedom(); - } else { - // use max - ret = Double.MAX_VALUE; - } - - return ret; - } - - /** {@inheritDoc} */ - @Override - protected double getInitialDomain(double p) { - // NOTE: chi squared is skewed to the left - // NOTE: therefore, P(X < μ) > 0.5 - - double ret; - - if (p < 0.5) { - // use 1/2 mean - ret = getDegreesOfFreedom() * 0.5; - } else { - // use mean - ret = getDegreesOfFreedom(); - } - - return ret; - } - - /** {@inheritDoc} */ - @Override - protected double getSolverAbsoluteAccuracy() { - return solverAbsoluteAccuracy; - } - - /** - * {@inheritDoc} - * - * The lower bound of the support is always 0 no matter the - * degrees of freedom. - * - * @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 - * degrees of freedom. - * - * @return upper bound of the support (always Double.POSITIVE_INFINITY) - */ - @Override - public double getSupportUpperBound() { - return Double.POSITIVE_INFINITY; - } - - /** - * {@inheritDoc} - * - * For {@code k} degrees of freedom, the mean is {@code k}. - */ - @Override - protected double calculateNumericalMean() { - return getDegreesOfFreedom(); - } - - /** - * {@inheritDoc} - * - * For {@code k} degrees of freedom, the variance is {@code 2 * k}. - * - * @return {@inheritDoc} - */ - @Override - protected double calculateNumericalVariance() { - return 2*getDegreesOfFreedom(); - } - - /** {@inheritDoc} */ - @Override - public boolean isSupportLowerBoundInclusive() { - return true; - } - - /** {@inheritDoc} */ - @Override - public boolean isSupportUpperBoundInclusive() { - return false; - } -} +/* + * 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.math.distribution; + +import java.io.Serializable; + + +/** + * Implementation of the chi-squared distribution. + * + * @see Chi-squared distribution (Wikipedia) + * @see Chi-squared Distribution (MathWorld) + * @version $Id: ChiSquaredDistribution.java 1206060 2011-11-25 05:16:56Z celestin $ + */ +public class ChiSquaredDistribution + extends AbstractContinuousDistribution + implements Serializable { + /** + * Default inverse cumulative probability accuracy + * @since 2.1 + */ + public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; + /** Serializable version identifier */ + private static final long serialVersionUID = -8352658048349159782L; + /** Internal Gamma distribution. */ + private final GammaDistribution gamma; + /** Inverse cumulative probability accuracy */ + private final double solverAbsoluteAccuracy; + + /** + * Create a Chi-Squared distribution with the given degrees of freedom. + * + * @param degreesOfFreedom Degrees of freedom. + */ + public ChiSquaredDistribution(double degreesOfFreedom) { + this(degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); + } + + /** + * Create a Chi-Squared distribution with the given degrees of freedom and + * inverse cumulative probability accuracy. + * + * @param degreesOfFreedom Degrees of freedom. + * @param inverseCumAccuracy the maximum absolute error in inverse + * cumulative probability estimates (defaults to + * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). + * @since 2.1 + */ + public ChiSquaredDistribution(double degreesOfFreedom, + double inverseCumAccuracy) { + gamma = new GammaDistribution(degreesOfFreedom / 2, 2); + solverAbsoluteAccuracy = inverseCumAccuracy; + } + + /** + * Access the number of degrees of freedom. + * + * @return the degrees of freedom. + */ + public double getDegreesOfFreedom() { + return gamma.getAlpha() * 2.0; + } + + /** {@inheritDoc} */ + public double density(double x) { + return gamma.density(x); + } + + /** {@inheritDoc} */ + public double cumulativeProbability(double x) { + return gamma.cumulativeProbability(x); + } + + /** + * {@inheritDoc} + * + * Returns {@code 0} when {@code p == 0} and + * {@code Double.POSITIVE_INFINITY} when {@code p == 1}. + */ + @Override + public double inverseCumulativeProbability(final double p) { + if (p == 0) { + return 0d; + } + if (p == 1) { + return Double.POSITIVE_INFINITY; + } + return super.inverseCumulativeProbability(p); + } + + /** {@inheritDoc} */ + @Override + protected double getDomainLowerBound(double p) { + return Double.MIN_VALUE * gamma.getBeta(); + } + + /** {@inheritDoc} */ + @Override + protected double getDomainUpperBound(double p) { + // NOTE: chi squared is skewed to the left + // NOTE: therefore, P(X < μ) > .5 + + double ret; + + if (p < .5) { + // use mean + ret = getDegreesOfFreedom(); + } else { + // use max + ret = Double.MAX_VALUE; + } + + return ret; + } + + /** {@inheritDoc} */ + @Override + protected double getInitialDomain(double p) { + // NOTE: chi squared is skewed to the left + // NOTE: therefore, P(X < μ) > 0.5 + + double ret; + + if (p < 0.5) { + // use 1/2 mean + ret = getDegreesOfFreedom() * 0.5; + } else { + // use mean + ret = getDegreesOfFreedom(); + } + + return ret; + } + + /** {@inheritDoc} */ + @Override + protected double getSolverAbsoluteAccuracy() { + return solverAbsoluteAccuracy; + } + + /** + * {@inheritDoc} + * + * The lower bound of the support is always 0 no matter the + * degrees of freedom. + * + * @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 + * degrees of freedom. + * + * @return upper bound of the support (always Double.POSITIVE_INFINITY) + */ + @Override + public double getSupportUpperBound() { + return Double.POSITIVE_INFINITY; + } + + /** + * {@inheritDoc} + * + * For {@code k} degrees of freedom, the mean is {@code k}. + */ + @Override + protected double calculateNumericalMean() { + return getDegreesOfFreedom(); + } + + /** + * {@inheritDoc} + * + * For {@code k} degrees of freedom, the variance is {@code 2 * k}. + * + * @return {@inheritDoc} + */ + @Override + protected double calculateNumericalVariance() { + return 2*getDegreesOfFreedom(); + } + + /** {@inheritDoc} */ + @Override + public boolean isSupportLowerBoundInclusive() { + return true; + } + + /** {@inheritDoc} */ + @Override + public boolean isSupportUpperBoundInclusive() { + return false; + } +} diff --git a/src/main/java/org/apache/commons/math/distribution/ExponentialDistribution.java b/src/main/java/org/apache/commons/math/distribution/ExponentialDistribution.java index 60da38277..934421935 100644 --- a/src/main/java/org/apache/commons/math/distribution/ExponentialDistribution.java +++ b/src/main/java/org/apache/commons/math/distribution/ExponentialDistribution.java @@ -16,24 +16,234 @@ */ package org.apache.commons.math.distribution; +import java.io.Serializable; + +import org.apache.commons.math.exception.NotStrictlyPositiveException; +import org.apache.commons.math.exception.OutOfRangeException; +import org.apache.commons.math.exception.util.LocalizedFormats; +import org.apache.commons.math.util.FastMath; + /** - * The Exponential Distribution. - * - *

- * References: - *

- *

+ * Implementation of the exponential distribution. * + * @see Exponential distribution (Wikipedia) + * @see Exponential distribution (MathWorld) * @version $Id$ */ -public interface ExponentialDistribution extends ContinuousDistribution { +public class ExponentialDistribution extends AbstractContinuousDistribution + implements Serializable { + /** + * Default inverse cumulative probability accuracy. + * @since 2.1 + */ + public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; + /** Serializable version identifier */ + private static final long serialVersionUID = 2401296428283614780L; + /** The mean of this distribution. */ + private final double mean; + /** Inverse cumulative probability accuracy. */ + private final double solverAbsoluteAccuracy; + + /** + * Create a exponential distribution with the given mean. + * @param mean mean of this distribution. + */ + public ExponentialDistribution(double mean) { + this(mean, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); + } + + /** + * Create a exponential distribution with the given mean. + * + * @param mean Mean of this distribution. + * @param inverseCumAccuracy Maximum absolute error in inverse + * cumulative probability estimates (defaults to + * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). + * @throws NotStrictlyPositiveException if {@code mean <= 0}. + * @since 2.1 + */ + public ExponentialDistribution(double mean, double inverseCumAccuracy) + throws NotStrictlyPositiveException{ + if (mean <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean); + } + this.mean = mean; + solverAbsoluteAccuracy = inverseCumAccuracy; + } + /** * Access the mean. * * @return the mean. */ - double getMean(); + public double getMean() { + return mean; + } + + /** {@inheritDoc} */ + public double density(double x) { + if (x < 0) { + return 0; + } + return FastMath.exp(-x / mean) / mean; + } + + /** + * {@inheritDoc} + * + * The implementation of this method is based on: + * + */ + public double cumulativeProbability(double x) { + double ret; + if (x <= 0.0) { + ret = 0.0; + } else { + ret = 1.0 - FastMath.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) throws OutOfRangeException { + double ret; + + if (p < 0.0 || p > 1.0) { + throw new OutOfRangeException(p, 0.0, 1.0); + } else if (p == 1.0) { + ret = Double.POSITIVE_INFINITY; + } else { + ret = -mean * FastMath.log(1.0 - p); + } + + return ret; + } + + /** + * {@inheritDoc} + * + *

Algorithm Description: this implementation uses the + * + * Inversion Method to generate exponentially distributed random values + * from uniform deviates.

+ * + * @return a random value. + * @since 2.2 + */ + @Override + public double sample() { + return randomData.nextExponential(mean); + } + + /** {@inheritDoc} */ + @Override + protected double getDomainLowerBound(double p) { + return 0; + } + + /** {@inheritDoc} */ + @Override + protected double getDomainUpperBound(double p) { + // NOTE: exponential is skewed to the left + // NOTE: therefore, P(X < μ) > .5 + + if (p < 0.5) { + // use mean + return mean; + } else { + // use max + return Double.MAX_VALUE; + } + } + + /** {@inheritDoc} */ + @Override + protected double getInitialDomain(double p) { + // TODO: try to improve on this estimate + // TODO: what should really happen here is not derive from + // AbstractContinuousDistribution + // TODO: because the inverse cumulative distribution is simple. + // Exponential is skewed to the left, therefore, P(X < μ) > .5 + if (p < 0.5) { + // use 1/2 mean + return mean * 0.5; + } else { + // use mean + return mean; + } + } + + /** {@inheritDoc} */ + @Override + protected double getSolverAbsoluteAccuracy() { + return solverAbsoluteAccuracy; + } + + /** + * {@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} + * + * For mean parameter {@code k}, the mean is {@code k}. + */ + @Override + protected double calculateNumericalMean() { + return getMean(); + } + + /** + * {@inheritDoc} + * + * For mean parameter {@code k}, the variance is {@code k^2}. + */ + @Override + protected double calculateNumericalVariance() { + final double m = getMean(); + return m * m; + } + + /** {@inheritDoc} */ + @Override + public boolean isSupportLowerBoundInclusive() { + return true; + } + + /** {@inheritDoc} */ + @Override + public boolean isSupportUpperBoundInclusive() { + return false; + } } diff --git a/src/main/java/org/apache/commons/math/distribution/ExponentialDistributionImpl.java b/src/main/java/org/apache/commons/math/distribution/ExponentialDistributionImpl.java deleted file mode 100644 index 7106d8971..000000000 --- a/src/main/java/org/apache/commons/math/distribution/ExponentialDistributionImpl.java +++ /dev/null @@ -1,279 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math.distribution; - -import java.io.Serializable; - -import org.apache.commons.math.exception.NotStrictlyPositiveException; -import org.apache.commons.math.exception.OutOfRangeException; -import org.apache.commons.math.exception.util.LocalizedFormats; -import org.apache.commons.math.util.FastMath; - -/** - * The default implementation of {@link ExponentialDistribution}. - * - * @version $Id$ - */ -public class ExponentialDistributionImpl extends AbstractContinuousDistribution - implements ExponentialDistribution, Serializable { - /** - * Default inverse cumulative probability accuracy. - * @since 2.1 - */ - public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; - /** Serializable version identifier */ - private static final long serialVersionUID = 2401296428283614780L; - /** The mean of this distribution. */ - private final double mean; - /** Inverse cumulative probability accuracy. */ - private final double solverAbsoluteAccuracy; - - /** - * Create a exponential distribution with the given mean. - * @param mean mean of this distribution. - */ - public ExponentialDistributionImpl(double mean) { - this(mean, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); - } - - /** - * Create a exponential distribution with the given mean. - * - * @param mean Mean of this distribution. - * @param inverseCumAccuracy Maximum absolute error in inverse - * cumulative probability estimates (defaults to - * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). - * @throws NotStrictlyPositiveException if {@code mean <= 0}. - * @since 2.1 - */ - public ExponentialDistributionImpl(double mean, double inverseCumAccuracy) { - if (mean <= 0) { - throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean); - } - this.mean = mean; - solverAbsoluteAccuracy = inverseCumAccuracy; - } - - /** - * {@inheritDoc} - */ - public double getMean() { - return mean; - } - - /** - * {@inheritDoc} - */ - public double density(double x) { - if (x < 0) { - return 0; - } - return FastMath.exp(-x / mean) / mean; - } - - /** - * {@inheritDoc} - * - * The implementation of this method is based on: - * - */ - public double cumulativeProbability(double x) { - double ret; - if (x <= 0.0) { - ret = 0.0; - } else { - ret = 1.0 - FastMath.exp(-x / mean); - } - return ret; - } - - /** - * {@inheritDoc} - * - * It will return {@code 0} when {@code p = 0} and - * {@code Double.POSITIVE_INFINITY} when {@code p = 1}. - */ - @Override - public double inverseCumulativeProbability(double p) throws OutOfRangeException { - double ret; - - if (p < 0.0 || p > 1.0) { - throw new OutOfRangeException(p, 0.0, 1.0); - } else if (p == 1.0) { - ret = Double.POSITIVE_INFINITY; - } else { - ret = -mean * FastMath.log(1.0 - p); - } - - return ret; - } - - /** - * Generates a random value sampled from this distribution. - * - *

Algorithm Description: Uses the Inversion - * Method to generate exponentially distributed random values from - * uniform deviates.

- * - * @return a random value. - * @since 2.2 - */ - @Override - public double sample() { - return randomData.nextExponential(mean); - } - - /** - * Access the domain value lower bound, based on {@code p}, used to - * bracket a CDF root. - * - * @param p Desired probability for the critical value. - * @return the domain value lower bound, i.e. {@code P(X < 'lower bound') < p}. - */ - @Override - protected double getDomainLowerBound(double p) { - return 0; - } - - /** - * Access the domain value upper bound, based on {@code p}, used to - * bracket a CDF root. - * - * @param p Desired probability for the critical value. - * @return the domain value upper bound, i.e. {@code P(X < 'upper bound') > p}. - */ - @Override - protected double getDomainUpperBound(double p) { - // NOTE: exponential is skewed to the left - // NOTE: therefore, P(X < μ) > .5 - - if (p < 0.5) { - // use mean - return mean; - } else { - // use max - return Double.MAX_VALUE; - } - } - - /** - * Access the initial domain value, based on {@code p}, used to - * bracket a CDF root. - * - * @param p Desired probability for the critical value. - * @return the initial domain value. - */ - @Override - protected double getInitialDomain(double p) { - // TODO: try to improve on this estimate - // TODO: what should really happen here is not derive from AbstractContinuousDistribution - // TODO: because the inverse cumulative distribution is simple. - // Exponential is skewed to the left, therefore, P(X < μ) > .5 - if (p < 0.5) { - // use 1/2 mean - return mean * 0.5; - } else { - // use mean - return mean; - } - } - - /** - * Return the absolute accuracy setting of the solver used to estimate - * inverse cumulative probabilities. - * - * @return the solver absolute accuracy. - * @since 2.1 - */ - @Override - protected double getSolverAbsoluteAccuracy() { - return solverAbsoluteAccuracy; - } - - /** - * {@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} - * - * For mean parameter k, the mean is - * k - * - * @return {@inheritDoc} - */ - @Override - protected double calculateNumericalMean() { - return getMean(); - } - - /** - * {@inheritDoc} - * - * For mean parameter k, the variance is - * k^2 - * - * @return {@inheritDoc} - */ - @Override - protected double calculateNumericalVariance() { - final double m = getMean(); - return m * m; - } - - /** - * {@inheritDoc} - */ - @Override - public boolean isSupportLowerBoundInclusive() { - return true; - } - - /** - * {@inheritDoc} - */ - @Override - public boolean isSupportUpperBoundInclusive() { - return false; - } -} diff --git a/src/main/java/org/apache/commons/math/distribution/FDistribution.java b/src/main/java/org/apache/commons/math/distribution/FDistribution.java index f9c897bd3..18985fe91 100644 --- a/src/main/java/org/apache/commons/math/distribution/FDistribution.java +++ b/src/main/java/org/apache/commons/math/distribution/FDistribution.java @@ -14,33 +14,277 @@ * See the License for the specific language governing permissions and * limitations under the License. */ + package org.apache.commons.math.distribution; +import java.io.Serializable; + +import org.apache.commons.math.exception.NotStrictlyPositiveException; +import org.apache.commons.math.exception.OutOfRangeException; +import org.apache.commons.math.exception.util.LocalizedFormats; +import org.apache.commons.math.special.Beta; +import org.apache.commons.math.util.FastMath; + /** - * F-Distribution. - * - *

- * References: - *

- *

+ * Implementation of the F-distribution. * + * @see F-distribution (Wikipedia) + * @see F-distribution (MathWorld) * @version $Id$ */ -public interface FDistribution extends ContinuousDistribution { +public class FDistribution + extends AbstractContinuousDistribution + implements Serializable { + /** + * Default inverse cumulative probability accuracy. + * @since 2.1 + */ + public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; + /** Serializable version identifier. */ + private static final long serialVersionUID = -8516354193418641566L; + /** The numerator degrees of freedom. */ + private final double numeratorDegreesOfFreedom; + /** The numerator degrees of freedom. */ + private final double denominatorDegreesOfFreedom; + /** Inverse cumulative probability accuracy. */ + private final double solverAbsoluteAccuracy; + + /** + * Create a F distribution using the given degrees of freedom. + * @param numeratorDegreesOfFreedom Numerator degrees of freedom. + * @param denominatorDegreesOfFreedom Denominator degrees of freedom. + * @throws NotStrictlyPositiveException if + * {@code numeratorDegreesOfFreedom <= 0} or + * {@code denominatorDegreesOfFreedom <= 0}. + */ + public FDistribution(double numeratorDegreesOfFreedom, + double denominatorDegreesOfFreedom) + throws NotStrictlyPositiveException { + this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom, + DEFAULT_INVERSE_ABSOLUTE_ACCURACY); + } + + /** + * Create an F distribution using the given degrees of freedom + * and inverse cumulative probability accuracy. + * @param numeratorDegreesOfFreedom Numerator degrees of freedom. + * @param denominatorDegreesOfFreedom Denominator degrees of freedom. + * @param inverseCumAccuracy the maximum absolute error in inverse + * cumulative probability estimates. + * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}) + * @throws NotStrictlyPositiveException if + * {@code numeratorDegreesOfFreedom <= 0} or + * {@code denominatorDegreesOfFreedom <= 0}. + * @since 2.1 + */ + public FDistribution(double numeratorDegreesOfFreedom, + double denominatorDegreesOfFreedom, + double inverseCumAccuracy) + throws NotStrictlyPositiveException { + if (numeratorDegreesOfFreedom <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM, + numeratorDegreesOfFreedom); + } + if (denominatorDegreesOfFreedom <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM, + denominatorDegreesOfFreedom); + } + this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom; + this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom; + solverAbsoluteAccuracy = inverseCumAccuracy; + } + + /** + * {@inheritDoc} + * + * @since 2.1 + */ + public double density(double x) { + final double nhalf = numeratorDegreesOfFreedom / 2; + final double mhalf = denominatorDegreesOfFreedom / 2; + final double logx = FastMath.log(x); + final double logn = FastMath.log(numeratorDegreesOfFreedom); + final double logm = FastMath.log(denominatorDegreesOfFreedom); + final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x + + denominatorDegreesOfFreedom); + return FastMath.exp(nhalf * logn + nhalf * logx - logx + + mhalf * logm - nhalf * lognxm - mhalf * lognxm - + Beta.logBeta(nhalf, mhalf)); + } + + /** + * {@inheritDoc} + * + * The implementation of this method is based on + * + */ + public double cumulativeProbability(double x) { + double ret; + if (x <= 0) { + ret = 0; + } else { + double n = numeratorDegreesOfFreedom; + double m = denominatorDegreesOfFreedom; + + ret = Beta.regularizedBeta((n * x) / (m + n * x), + 0.5 * n, + 0.5 * m); + } + return ret; + } + + /** + * {@inheritDoc} + * + * Returns {@code 0} when {@code p == 0} and + * {@code Double.POSITIVE_INFINITY} when {@code p == 1}. + */ + @Override + public double inverseCumulativeProbability(final double p) throws OutOfRangeException { + if (p == 0) { + return 0; + } + if (p == 1) { + return Double.POSITIVE_INFINITY; + } + return super.inverseCumulativeProbability(p); + } + + /** {@inheritDoc} */ + @Override + protected double getDomainLowerBound(double p) { + return 0; + } + + /** {@inheritDoc} */ + @Override + protected double getDomainUpperBound(double p) { + return Double.MAX_VALUE; + } + + /** {@inheritDoc} */ + @Override + protected double getInitialDomain(double p) { + double ret = 1; + double d = denominatorDegreesOfFreedom; + if (d > 2) { + // use mean + ret = d / (d - 2); + } + return ret; + } + /** * Access the numerator degrees of freedom. * * @return the numerator degrees of freedom. */ - double getNumeratorDegreesOfFreedom(); + public double getNumeratorDegreesOfFreedom() { + return numeratorDegreesOfFreedom; + } /** * Access the denominator degrees of freedom. * * @return the denominator degrees of freedom. */ - double getDenominatorDegreesOfFreedom(); + public double getDenominatorDegreesOfFreedom() { + return denominatorDegreesOfFreedom; + } + + /** {@inheritDoc} */ + @Override + protected double getSolverAbsoluteAccuracy() { + return solverAbsoluteAccuracy; + } + + /** + * {@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} + * + * For denominator degrees of freedom parameter {@code b}, the mean is + * + */ + @Override + protected double calculateNumericalMean() { + 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 + * + */ + @Override + protected double calculateNumericalVariance() { + 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} */ + @Override + public boolean isSupportLowerBoundInclusive() { + return true; + } + + /** {@inheritDoc} */ + @Override + public boolean isSupportUpperBoundInclusive() { + return false; + } } diff --git a/src/main/java/org/apache/commons/math/distribution/FDistributionImpl.java b/src/main/java/org/apache/commons/math/distribution/FDistributionImpl.java deleted file mode 100644 index cd4689265..000000000 --- a/src/main/java/org/apache/commons/math/distribution/FDistributionImpl.java +++ /dev/null @@ -1,318 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math.distribution; - -import java.io.Serializable; - -import org.apache.commons.math.exception.NotStrictlyPositiveException; -import org.apache.commons.math.exception.OutOfRangeException; -import org.apache.commons.math.exception.util.LocalizedFormats; -import org.apache.commons.math.special.Beta; -import org.apache.commons.math.util.FastMath; - -/** - * Default implementation of - * {@link org.apache.commons.math.distribution.FDistribution}. - * - * @version $Id$ - */ -public class FDistributionImpl - extends AbstractContinuousDistribution - implements FDistribution, Serializable { - /** - * Default inverse cumulative probability accuracy. - * @since 2.1 - */ - public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; - /** Serializable version identifier. */ - private static final long serialVersionUID = -8516354193418641566L; - /** The numerator degrees of freedom. */ - private final double numeratorDegreesOfFreedom; - /** The numerator degrees of freedom. */ - private final double denominatorDegreesOfFreedom; - /** Inverse cumulative probability accuracy. */ - private final double solverAbsoluteAccuracy; - - /** - * Create a F distribution using the given degrees of freedom. - * @param numeratorDegreesOfFreedom Numerator degrees of freedom. - * @param denominatorDegreesOfFreedom Denominator degrees of freedom. - * @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0} - * or {@code denominatorDegreesOfFreedom <= 0}. - */ - public FDistributionImpl(double numeratorDegreesOfFreedom, - double denominatorDegreesOfFreedom) { - this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom, - DEFAULT_INVERSE_ABSOLUTE_ACCURACY); - } - - /** - * Create an F distribution using the given degrees of freedom - * and inverse cumulative probability accuracy. - * @param numeratorDegreesOfFreedom Numerator degrees of freedom. - * @param denominatorDegreesOfFreedom Denominator degrees of freedom. - * @param inverseCumAccuracy the maximum absolute error in inverse - * cumulative probability estimates. - * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}) - * @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0} - * or {@code denominatorDegreesOfFreedom <= 0}. - * @since 2.1 - */ - public FDistributionImpl(double numeratorDegreesOfFreedom, - double denominatorDegreesOfFreedom, - double inverseCumAccuracy) { - if (numeratorDegreesOfFreedom <= 0) { - throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM, - numeratorDegreesOfFreedom); - } - if (denominatorDegreesOfFreedom <= 0) { - throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM, - denominatorDegreesOfFreedom); - } - this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom; - this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom; - solverAbsoluteAccuracy = inverseCumAccuracy; - } - - /** - * {@inheritDoc} - * - * @since 2.1 - */ - public double density(double x) { - final double nhalf = numeratorDegreesOfFreedom / 2; - final double mhalf = denominatorDegreesOfFreedom / 2; - final double logx = FastMath.log(x); - final double logn = FastMath.log(numeratorDegreesOfFreedom); - final double logm = FastMath.log(denominatorDegreesOfFreedom); - final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x + - denominatorDegreesOfFreedom); - return FastMath.exp(nhalf * logn + nhalf * logx - logx + - mhalf * logm - nhalf * lognxm - mhalf * lognxm - - Beta.logBeta(nhalf, mhalf)); - } - - /** - * {@inheritDoc} - * - * The implementation of this method is based on - * - */ - public double cumulativeProbability(double x) { - double ret; - if (x <= 0) { - ret = 0; - } else { - double n = numeratorDegreesOfFreedom; - double m = denominatorDegreesOfFreedom; - - ret = Beta.regularizedBeta((n * x) / (m + n * x), - 0.5 * n, - 0.5 * m); - } - return ret; - } - - /** - * {@inheritDoc} - * - * It will return {@code 0} when {@code p = 0} and - * {@code Double.POSITIVE_INFINITY} when {@code p = 1}. - */ - @Override - public double inverseCumulativeProbability(final double p) throws OutOfRangeException { - if (p == 0) { - return 0; - } - if (p == 1) { - return Double.POSITIVE_INFINITY; - } - return super.inverseCumulativeProbability(p); - } - - /** - * Access the domain value lower bound, based on {@code p}, used to - * bracket a CDF root. This method is used by - * {@link #inverseCumulativeProbability(double)} to find critical values. - * - * @param p Desired probability for the critical value. - * @return the domain value lower bound, i.e. {@code P(X < 'lower bound') < p}. - */ - @Override - protected double getDomainLowerBound(double p) { - return 0; - } - - /** - * Access the domain value upper bound, based on {@code p}, used to - * bracket a CDF root. This method is used by - * {@link #inverseCumulativeProbability(double)} to find critical values. - * - * @param p Desired probability for the critical value. - * @return the domain value upper bound, i.e. {@code P(X < 'upper bound') > p}. - */ - @Override - protected double getDomainUpperBound(double p) { - return Double.MAX_VALUE; - } - - /** - * Access the initial domain value, based on {@code p}, used to - * bracket a CDF root. This method is used by - * {@link #inverseCumulativeProbability(double)} to find critical values. - * - * @param p Desired probability for the critical value. - * @return the initial domain value. - */ - @Override - protected double getInitialDomain(double p) { - double ret = 1; - double d = denominatorDegreesOfFreedom; - if (d > 2) { - // use mean - ret = d / (d - 2); - } - return ret; - } - - /** - * {@inheritDoc} - */ - public double getNumeratorDegreesOfFreedom() { - return numeratorDegreesOfFreedom; - } - - /** - * {@inheritDoc} - */ - public double getDenominatorDegreesOfFreedom() { - return denominatorDegreesOfFreedom; - } - - /** - * Return the absolute accuracy setting of the solver used to estimate - * inverse cumulative probabilities. - * - * @return the solver absolute accuracy - * @since 2.1 - */ - @Override - protected double getSolverAbsoluteAccuracy() { - return solverAbsoluteAccuracy; - } - - /** - * {@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} - * - * For denominator degrees of freedom parameter b, - * the mean is - * - * - * @return {@inheritDoc} - */ - @Override - protected double calculateNumericalMean() { - final double denominatorDF = getDenominatorDegreesOfFreedom(); - - if (denominatorDF > 2) { - return denominatorDF / (denominatorDF - 2); - } - - return Double.NaN; - } - - /** - * {@inheritDoc} - * - * For numerator degrees of freedom parameter a - * and denominator degrees of freedom parameter b, - * the variance is - * - * - * @return {@inheritDoc} - */ - @Override - protected double calculateNumericalVariance() { - 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} - */ - @Override - public boolean isSupportLowerBoundInclusive() { - return true; - } - - /** - * {@inheritDoc} - */ - @Override - public boolean isSupportUpperBoundInclusive() { - return false; - } -} diff --git a/src/main/java/org/apache/commons/math/distribution/GammaDistribution.java b/src/main/java/org/apache/commons/math/distribution/GammaDistribution.java index ab2e7a206..e7bcf1932 100644 --- a/src/main/java/org/apache/commons/math/distribution/GammaDistribution.java +++ b/src/main/java/org/apache/commons/math/distribution/GammaDistribution.java @@ -16,31 +16,254 @@ */ package org.apache.commons.math.distribution; +import java.io.Serializable; + +import org.apache.commons.math.exception.NotStrictlyPositiveException; +import org.apache.commons.math.exception.util.LocalizedFormats; +import org.apache.commons.math.special.Gamma; +import org.apache.commons.math.util.FastMath; + /** - * The Gamma Distribution. - * - *

- * References: - *

- *

+ * Implementation of the Gamma distribution. * + * @see Gamma distribution (Wikipedia) + * @see Gamma distribution (MathWorld) * @version $Id$ */ -public interface GammaDistribution extends ContinuousDistribution { +public class GammaDistribution extends AbstractContinuousDistribution + implements Serializable { /** - * Access the alpha shape parameter. - * - * @return alpha. + * Default inverse cumulative probability accuracy. + * @since 2.1 */ - double getAlpha(); + public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; + /** Serializable version identifier. */ + private static final long serialVersionUID = -3239549463135430361L; + /** The shape parameter. */ + private final double alpha; + /** The scale parameter. */ + private final double beta; + /** Inverse cumulative probability accuracy. */ + private final double solverAbsoluteAccuracy; /** - * Access the beta scale parameter. - * - * @return beta. + * Create a new gamma distribution with the given {@code alpha} and + * {@code beta} values. + * @param alpha the shape parameter. + * @param beta the scale parameter. */ - double getBeta(); + public GammaDistribution(double alpha, double beta) { + this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); + } + + /** + * Create a new gamma distribution with the given {@code alpha} and + * {@code beta} values. + * + * @param alpha Shape parameter. + * @param beta Scale parameter. + * @param inverseCumAccuracy Maximum absolute error in inverse + * cumulative probability estimates (defaults to + * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). + * @throws NotStrictlyPositiveException if {@code alpha <= 0} or + * {@code beta <= 0}. + * @since 2.1 + */ + public GammaDistribution(double alpha, double beta, double inverseCumAccuracy) + throws NotStrictlyPositiveException { + if (alpha <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.ALPHA, alpha); + } + if (beta <= 0) { + throw new NotStrictlyPositiveException(LocalizedFormats.BETA, beta); + } + + this.alpha = alpha; + this.beta = beta; + solverAbsoluteAccuracy = inverseCumAccuracy; + } + + /** + * {@inheritDoc} + * + * The implementation of this method is based on: + * + */ + public double cumulativeProbability(double x) { + double ret; + + if (x <= 0) { + ret = 0; + } else { + ret = Gamma.regularizedGammaP(alpha, x / beta); + } + + return ret; + } + + /** + * {@inheritDoc} + * + * Returns {@code 0} when {@code p == 0} and + * {@code Double.POSITIVE_INFINITY} when {@code p == 1}. + */ + @Override + public double inverseCumulativeProbability(final double p) { + if (p == 0) { + return 0; + } + if (p == 1) { + return Double.POSITIVE_INFINITY; + } + return super.inverseCumulativeProbability(p); + } + + /** + * Access the {@code alpha} shape parameter. + * + * @return {@code alpha}. + */ + public double getAlpha() { + return alpha; + } + + /** + * Access the {@code beta} scale parameter. + * + * @return {@code beta}. + */ + public double getBeta() { + return beta; + } + + /** {@inheritDoc} */ + public double density(double x) { + if (x < 0) { + return 0; + } + return FastMath.pow(x / beta, alpha - 1) / beta * + FastMath.exp(-x / beta) / FastMath.exp(Gamma.logGamma(alpha)); + } + + /** {@inheritDoc} */ + @Override + protected double getDomainLowerBound(double p) { + // TODO: try to improve on this estimate + return Double.MIN_VALUE; + } + + /** {@inheritDoc} */ + @Override + protected double getDomainUpperBound(double p) { + // TODO: try to improve on this estimate + // NOTE: gamma is skewed to the left + // NOTE: therefore, P(X < μ) > .5 + + double ret; + + if (p < 0.5) { + // use mean + ret = alpha * beta; + } else { + // use max value + ret = Double.MAX_VALUE; + } + + return ret; + } + + /** {@inheritDoc} */ + @Override + protected double getInitialDomain(double p) { + // TODO: try to improve on this estimate + // Gamma is skewed to the left, therefore, P(X < μ) > .5 + + double ret; + + if (p < 0.5) { + // use 1/2 mean + ret = alpha * beta * 0.5; + } else { + // use mean + ret = alpha * beta; + } + + return ret; + } + + /** {@inheritDoc} */ + @Override + protected double getSolverAbsoluteAccuracy() { + return solverAbsoluteAccuracy; + } + + /** + * {@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} + * + * For shape parameter {@code alpha} and scale parameter {@code beta}, the + * mean is {@code alpha * beta}. + */ + @Override + protected double calculateNumericalMean() { + return getAlpha() * getBeta(); + } + + /** + * {@inheritDoc} + * + * For shape parameter {@code alpha} and scale parameter {@code beta}, the + * variance is {@code alpha * beta^2}. + * + * @return {@inheritDoc} + */ + @Override + protected double calculateNumericalVariance() { + final double b = getBeta(); + return getAlpha() * b * b; + } + + /** {@inheritDoc} */ + @Override + public boolean isSupportLowerBoundInclusive() { + return true; + } + + /** {@inheritDoc} */ + @Override + public boolean isSupportUpperBoundInclusive() { + return false; + } } diff --git a/src/main/java/org/apache/commons/math/distribution/GammaDistributionImpl.java b/src/main/java/org/apache/commons/math/distribution/GammaDistributionImpl.java deleted file mode 100644 index 013e5de05..000000000 --- a/src/main/java/org/apache/commons/math/distribution/GammaDistributionImpl.java +++ /dev/null @@ -1,297 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math.distribution; - -import java.io.Serializable; - -import org.apache.commons.math.exception.NotStrictlyPositiveException; -import org.apache.commons.math.exception.util.LocalizedFormats; -import org.apache.commons.math.special.Gamma; -import org.apache.commons.math.util.FastMath; - -/** - * The default implementation of {@link GammaDistribution}. - * - * @version $Id$ - */ -public class GammaDistributionImpl extends AbstractContinuousDistribution - implements GammaDistribution, Serializable { - /** - * Default inverse cumulative probability accuracy. - * @since 2.1 - */ - public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; - /** Serializable version identifier. */ - private static final long serialVersionUID = -3239549463135430361L; - /** The shape parameter. */ - private final double alpha; - /** The scale parameter. */ - private final double beta; - /** Inverse cumulative probability accuracy. */ - private final double solverAbsoluteAccuracy; - - /** - * Create a new gamma distribution with the given alpha and beta values. - * @param alpha the shape parameter. - * @param beta the scale parameter. - */ - public GammaDistributionImpl(double alpha, double beta) { - this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); - } - - /** - * Create a new gamma distribution with the given alpha and beta values. - * - * @param alpha Shape parameter. - * @param beta Scale parameter. - * @param inverseCumAccuracy Maximum absolute error in inverse - * cumulative probability estimates (defaults to - * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). - * @throws NotStrictlyPositiveException if {@code alpha <= 0} or - * {@code beta <= 0}. - * @since 2.1 - */ - public GammaDistributionImpl(double alpha, double beta, double inverseCumAccuracy) { - if (alpha <= 0) { - throw new NotStrictlyPositiveException(LocalizedFormats.ALPHA, alpha); - } - if (beta <= 0) { - throw new NotStrictlyPositiveException(LocalizedFormats.BETA, beta); - } - - this.alpha = alpha; - this.beta = beta; - solverAbsoluteAccuracy = inverseCumAccuracy; - } - - /** - * {@inheritDoc} - * - * The implementation of this method is based on: - * - */ - public double cumulativeProbability(double x) { - double ret; - - if (x <= 0) { - ret = 0; - } else { - ret = Gamma.regularizedGammaP(alpha, x / beta); - } - - return ret; - } - - /** - * {@inheritDoc} - * - * It will return {@code 0} when {@cod p = 0} and - * {@code Double.POSITIVE_INFINITY} when {@code p = 1}. - */ - @Override - public double inverseCumulativeProbability(final double p) { - if (p == 0) { - return 0; - } - if (p == 1) { - return Double.POSITIVE_INFINITY; - } - return super.inverseCumulativeProbability(p); - } - - /** - * {@inheritDoc} - */ - public double getAlpha() { - return alpha; - } - - /** - * {@inheritDoc} - */ - public double getBeta() { - return beta; - } - - /** - * {@inheritDoc} - */ - public double density(double x) { - if (x < 0) { - return 0; - } - return FastMath.pow(x / beta, alpha - 1) / beta * - FastMath.exp(-x / beta) / FastMath.exp(Gamma.logGamma(alpha)); - } - - /** - * Access the domain value lower bound, based on {@code p}, used to - * bracket a CDF root. This method is used by - * {@link #inverseCumulativeProbability(double)} to find critical values. - * - * @param p Desired probability for the critical value. - * @return the domain value lower bound, i.e. {@code P(X < 'lower bound') < p}. - */ - @Override - protected double getDomainLowerBound(double p) { - // TODO: try to improve on this estimate - return Double.MIN_VALUE; - } - - /** - * Access the domain value upper bound, based on {@code p}, used to - * bracket a CDF root. This method is used by - * {@link #inverseCumulativeProbability(double)} to find critical values. - * - * @param p Desired probability for the critical value. - * @return the domain value upper bound, i.e. {@code P(X < 'upper bound') > p}. - */ - @Override - protected double getDomainUpperBound(double p) { - // TODO: try to improve on this estimate - // NOTE: gamma is skewed to the left - // NOTE: therefore, P(X < μ) > .5 - - double ret; - - if (p < 0.5) { - // use mean - ret = alpha * beta; - } else { - // use max value - ret = Double.MAX_VALUE; - } - - return ret; - } - - /** - * Access the initial domain value, based on {@code p}, used to - * bracket a CDF root. This method is used by - * {@link #inverseCumulativeProbability(double)} to find critical values. - * - * @param p Desired probability for the critical value. - * @return the initial domain value. - */ - @Override - protected double getInitialDomain(double p) { - // TODO: try to improve on this estimate - // Gamma is skewed to the left, therefore, P(X < μ) > .5 - - double ret; - - if (p < 0.5) { - // use 1/2 mean - ret = alpha * beta * 0.5; - } else { - // use mean - ret = alpha * beta; - } - - return ret; - } - - /** - * Return the absolute accuracy setting of the solver used to estimate - * inverse cumulative probabilities. - * - * @return the solver absolute accuracy. - * @since 2.1 - */ - @Override - protected double getSolverAbsoluteAccuracy() { - return solverAbsoluteAccuracy; - } - - /** - * {@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} - * - * For shape parameter alpha and scale - * parameter beta, the mean is - * alpha * beta - * - * @return {@inheritDoc} - */ - @Override - protected double calculateNumericalMean() { - return getAlpha() * getBeta(); - } - - /** - * {@inheritDoc} - * - * For shape parameter alpha and scale - * parameter beta, the variance is - * alpha * beta^2 - * - * @return {@inheritDoc} - */ - @Override - protected double calculateNumericalVariance() { - final double b = getBeta(); - return getAlpha() * b * b; - } - - /** - * {@inheritDoc} - */ - @Override - public boolean isSupportLowerBoundInclusive() { - return true; - } - - /** - * {@inheritDoc} - */ - @Override - public boolean isSupportUpperBoundInclusive() { - return false; - } -} diff --git a/src/main/java/org/apache/commons/math/random/RandomDataImpl.java b/src/main/java/org/apache/commons/math/random/RandomDataImpl.java index 34f941a22..ba6c7db39 100644 --- a/src/main/java/org/apache/commons/math/random/RandomDataImpl.java +++ b/src/main/java/org/apache/commons/math/random/RandomDataImpl.java @@ -29,7 +29,7 @@ import org.apache.commons.math.distribution.BinomialDistribution; import org.apache.commons.math.distribution.CauchyDistribution; import org.apache.commons.math.distribution.ChiSquaredDistribution; import org.apache.commons.math.distribution.ContinuousDistribution; -import org.apache.commons.math.distribution.FDistributionImpl; +import org.apache.commons.math.distribution.FDistribution; import org.apache.commons.math.distribution.HypergeometricDistributionImpl; import org.apache.commons.math.distribution.IntegerDistribution; import org.apache.commons.math.distribution.PascalDistributionImpl; @@ -654,7 +654,7 @@ public class RandomDataImpl implements RandomData, Serializable { } /** - * Generates a random value from the {@link FDistributionImpl F Distribution}. + * Generates a random value from the {@link FDistribution F Distribution}. * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion} * to generate random values. * @@ -664,12 +664,12 @@ public class RandomDataImpl implements RandomData, Serializable { * @since 2.2 */ public double nextF(double numeratorDf, double denominatorDf) { - return nextInversionDeviate(new FDistributionImpl(numeratorDf, denominatorDf)); + return nextInversionDeviate(new FDistribution(numeratorDf, denominatorDf)); } /** *

Generates a random value from the - * {@link org.apache.commons.math.distribution.GammaDistributionImpl Gamma Distribution}.

+ * {@link org.apache.commons.math.distribution.GammaDistribution Gamma Distribution}.

* *

This implementation uses the following algorithms:

* diff --git a/src/main/java/org/apache/commons/math/stat/inference/OneWayAnovaImpl.java b/src/main/java/org/apache/commons/math/stat/inference/OneWayAnovaImpl.java index 2d686eab1..5dadc8046 100644 --- a/src/main/java/org/apache/commons/math/stat/inference/OneWayAnovaImpl.java +++ b/src/main/java/org/apache/commons/math/stat/inference/OneWayAnovaImpl.java @@ -21,7 +21,6 @@ import java.util.Collection; import org.apache.commons.math.MathException; import org.apache.commons.math.MathRuntimeException; import org.apache.commons.math.distribution.FDistribution; -import org.apache.commons.math.distribution.FDistributionImpl; import org.apache.commons.math.exception.util.LocalizedFormats; import org.apache.commons.math.stat.descriptive.summary.Sum; import org.apache.commons.math.stat.descriptive.summary.SumOfSquares; @@ -84,7 +83,7 @@ public class OneWayAnovaImpl implements OneWayAnova { public double anovaPValue(Collection categoryData) throws IllegalArgumentException, MathException { AnovaStats a = anovaStats(categoryData); - FDistribution fdist = new FDistributionImpl(a.dfbg, a.dfwg); + FDistribution fdist = new FDistribution(a.dfbg, a.dfwg); return 1.0 - fdist.cumulativeProbability(a.F); } diff --git a/src/test/java/org/apache/commons/math/distribution/ExponentialDistributionTest.java b/src/test/java/org/apache/commons/math/distribution/ExponentialDistributionTest.java index 3ccb1ebe8..7c80247fb 100644 --- a/src/test/java/org/apache/commons/math/distribution/ExponentialDistributionTest.java +++ b/src/test/java/org/apache/commons/math/distribution/ExponentialDistributionTest.java @@ -43,7 +43,7 @@ public class ExponentialDistributionTest extends ContinuousDistributionAbstractT /** Creates the default continuous distribution instance to use in tests. */ @Override public ExponentialDistribution makeDistribution() { - return new ExponentialDistributionImpl(5.0); + return new ExponentialDistribution(5.0); } /** Creates the default cumulative probability distribution test input values */ @@ -92,14 +92,14 @@ public class ExponentialDistributionTest extends ContinuousDistributionAbstractT @Test public void testDensity() { - ExponentialDistribution d1 = new ExponentialDistributionImpl(1); + ExponentialDistribution d1 = new ExponentialDistribution(1); Assert.assertTrue(Precision.equals(0.0, d1.density(-1e-9), 1)); Assert.assertTrue(Precision.equals(1.0, d1.density(0.0), 1)); Assert.assertTrue(Precision.equals(0.0, d1.density(1000.0), 1)); Assert.assertTrue(Precision.equals(FastMath.exp(-1), d1.density(1.0), 1)); Assert.assertTrue(Precision.equals(FastMath.exp(-2), d1.density(2.0), 1)); - ExponentialDistribution d2 = new ExponentialDistributionImpl(3); + ExponentialDistribution d2 = new ExponentialDistribution(3); Assert.assertTrue(Precision.equals(1/3.0, d2.density(0.0), 1)); // computed using print(dexp(1, rate=1/3), digits=10) in R 2.5 Assert.assertEquals(0.2388437702, d2.density(1.0), 1e-8); @@ -116,19 +116,19 @@ public class ExponentialDistributionTest extends ContinuousDistributionAbstractT @Test(expected=NotStrictlyPositiveException.class) public void testPreconditions() { - new ExponentialDistributionImpl(0); + new ExponentialDistribution(0); } @Test public void testMoments() { final double tol = 1e-9; ExponentialDistribution dist; - - dist = new ExponentialDistributionImpl(11d); + + dist = new ExponentialDistribution(11d); Assert.assertEquals(dist.getNumericalMean(), 11d, tol); Assert.assertEquals(dist.getNumericalVariance(), 11d * 11d, tol); - - dist = new ExponentialDistributionImpl(10.5d); + + dist = new ExponentialDistribution(10.5d); Assert.assertEquals(dist.getNumericalMean(), 10.5d, tol); Assert.assertEquals(dist.getNumericalVariance(), 10.5d * 10.5d, tol); } diff --git a/src/test/java/org/apache/commons/math/distribution/FDistributionTest.java b/src/test/java/org/apache/commons/math/distribution/FDistributionTest.java index 6f35ad832..d60a9bf04 100644 --- a/src/test/java/org/apache/commons/math/distribution/FDistributionTest.java +++ b/src/test/java/org/apache/commons/math/distribution/FDistributionTest.java @@ -34,7 +34,7 @@ public class FDistributionTest extends ContinuousDistributionAbstractTest { /** Creates the default continuous distribution instance to use in tests. */ @Override public FDistribution makeDistribution() { - return new FDistributionImpl(5.0, 6.0); + return new FDistribution(5.0, 6.0); } /** Creates the default cumulative probability distribution test input values */ @@ -91,13 +91,13 @@ public class FDistributionTest extends ContinuousDistributionAbstractTest { @Test public void testPreconditions() { try { - new FDistributionImpl(0, 1); + new FDistribution(0, 1); Assert.fail("Expecting NotStrictlyPositiveException for df = 0"); } catch (NotStrictlyPositiveException ex) { // Expected. } try { - new FDistributionImpl(1, 0); + new FDistribution(1, 0); Assert.fail("Expecting NotStrictlyPositiveException for df = 0"); } catch (NotStrictlyPositiveException ex) { // Expected. @@ -106,7 +106,7 @@ public class FDistributionTest extends ContinuousDistributionAbstractTest { @Test public void testLargeDegreesOfFreedom() throws Exception { - FDistributionImpl fd = new FDistributionImpl(100000, 100000); + FDistribution fd = new FDistribution(100000, 100000); double p = fd.cumulativeProbability(.999); double x = fd.inverseCumulativeProbability(p); Assert.assertEquals(.999, x, 1.0e-5); @@ -114,12 +114,12 @@ public class FDistributionTest extends ContinuousDistributionAbstractTest { @Test public void testSmallDegreesOfFreedom() throws Exception { - FDistributionImpl fd = new FDistributionImpl(1, 1); + FDistribution fd = new FDistribution(1, 1); double p = fd.cumulativeProbability(0.975); double x = fd.inverseCumulativeProbability(p); Assert.assertEquals(0.975, x, 1.0e-5); - fd = new FDistributionImpl(1, 2); + fd = new FDistribution(1, 2); p = fd.cumulativeProbability(0.975); x = fd.inverseCumulativeProbability(p); Assert.assertEquals(0.975, x, 1.0e-5); @@ -129,17 +129,17 @@ public class FDistributionTest extends ContinuousDistributionAbstractTest { public void testMoments() { final double tol = 1e-9; FDistribution dist; - - dist = new FDistributionImpl(1, 2); + + dist = new FDistribution(1, 2); Assert.assertTrue(Double.isNaN(dist.getNumericalMean())); Assert.assertTrue(Double.isNaN(dist.getNumericalVariance())); - - dist = new FDistributionImpl(1, 3); + + dist = new FDistribution(1, 3); Assert.assertEquals(dist.getNumericalMean(), 3d / (3d - 2d), tol); Assert.assertTrue(Double.isNaN(dist.getNumericalVariance())); - - dist = new FDistributionImpl(1, 5); + + dist = new FDistribution(1, 5); Assert.assertEquals(dist.getNumericalMean(), 5d / (5d - 2d), tol); - Assert.assertEquals(dist.getNumericalVariance(), (2d * 5d * 5d * 4d) / 9d, tol); + Assert.assertEquals(dist.getNumericalVariance(), (2d * 5d * 5d * 4d) / 9d, tol); } } diff --git a/src/test/java/org/apache/commons/math/distribution/GammaDistributionTest.java b/src/test/java/org/apache/commons/math/distribution/GammaDistributionTest.java index f2ba4c1b8..e38d1073f 100644 --- a/src/test/java/org/apache/commons/math/distribution/GammaDistributionTest.java +++ b/src/test/java/org/apache/commons/math/distribution/GammaDistributionTest.java @@ -35,7 +35,7 @@ public class GammaDistributionTest extends ContinuousDistributionAbstractTest { /** Creates the default continuous distribution instance to use in tests. */ @Override public GammaDistribution makeDistribution() { - return new GammaDistributionImpl(4d, 2d); + return new GammaDistribution(4d, 2d); } /** Creates the default cumulative probability distribution test input values */ @@ -77,13 +77,13 @@ public class GammaDistributionTest extends ContinuousDistributionAbstractTest { @Test public void testPreconditions() { try { - new GammaDistributionImpl(0, 1); + new GammaDistribution(0, 1); Assert.fail("Expecting NotStrictlyPositiveException for alpha = 0"); } catch (NotStrictlyPositiveException ex) { // Expected. } try { - new GammaDistributionImpl(1, 0); + new GammaDistribution(1, 0); Assert.fail("Expecting NotStrictlyPositiveException for alpha = 0"); } catch (NotStrictlyPositiveException ex) { // Expected. @@ -108,13 +108,13 @@ public class GammaDistributionTest extends ContinuousDistributionAbstractTest { } private void testProbability(double x, double a, double b, double expected) throws Exception { - GammaDistribution distribution = new GammaDistributionImpl( a, b ); + GammaDistribution distribution = new GammaDistribution( a, b ); double actual = distribution.cumulativeProbability(x); Assert.assertEquals("probability for " + x, expected, actual, 10e-4); } private void testValue(double expected, double a, double b, double p) throws Exception { - GammaDistribution distribution = new GammaDistributionImpl( a, b ); + GammaDistribution distribution = new GammaDistribution( a, b ); double actual = distribution.inverseCumulativeProbability(p); Assert.assertEquals("critical value for " + p, expected, actual, 10e-4); } @@ -141,7 +141,7 @@ public class GammaDistributionTest extends ContinuousDistributionAbstractTest { } private void checkDensity(double alpha, double rate, double[] x, double[] expected) { - GammaDistribution d = new GammaDistributionImpl(alpha, 1 / rate); + GammaDistribution d = new GammaDistribution(alpha, 1 / rate); for (int i = 0; i < x.length; i++) { Assert.assertEquals(expected[i], d.density(x[i]), 1e-5); } @@ -158,12 +158,12 @@ public class GammaDistributionTest extends ContinuousDistributionAbstractTest { public void testMoments() { final double tol = 1e-9; GammaDistribution dist; - - dist = new GammaDistributionImpl(1, 2); + + dist = new GammaDistribution(1, 2); Assert.assertEquals(dist.getNumericalMean(), 2, tol); - Assert.assertEquals(dist.getNumericalVariance(), 4, tol); - - dist = new GammaDistributionImpl(1.1, 4.2); + Assert.assertEquals(dist.getNumericalVariance(), 4, tol); + + dist = new GammaDistribution(1.1, 4.2); Assert.assertEquals(dist.getNumericalMean(), 1.1d * 4.2d, tol); Assert.assertEquals(dist.getNumericalVariance(), 1.1d * 4.2d * 4.2d, tol); } diff --git a/src/test/java/org/apache/commons/math/random/RandomDataTest.java b/src/test/java/org/apache/commons/math/random/RandomDataTest.java index 9859e1d36..67ee45c64 100644 --- a/src/test/java/org/apache/commons/math/random/RandomDataTest.java +++ b/src/test/java/org/apache/commons/math/random/RandomDataTest.java @@ -30,9 +30,9 @@ import org.apache.commons.math.distribution.BinomialDistribution; import org.apache.commons.math.distribution.BinomialDistributionTest; import org.apache.commons.math.distribution.CauchyDistribution; import org.apache.commons.math.distribution.ChiSquaredDistribution; -import org.apache.commons.math.distribution.ExponentialDistributionImpl; -import org.apache.commons.math.distribution.FDistributionImpl; -import org.apache.commons.math.distribution.GammaDistributionImpl; +import org.apache.commons.math.distribution.ExponentialDistribution; +import org.apache.commons.math.distribution.FDistribution; +import org.apache.commons.math.distribution.GammaDistribution; import org.apache.commons.math.distribution.HypergeometricDistributionImpl; import org.apache.commons.math.distribution.HypergeometricDistributionTest; import org.apache.commons.math.distribution.PascalDistributionImpl; @@ -616,7 +616,7 @@ public class RandomDataTest { long[] counts; // Mean 1 - quartiles = TestUtils.getDistributionQuartiles(new ExponentialDistributionImpl(1)); + quartiles = TestUtils.getDistributionQuartiles(new ExponentialDistribution(1)); counts = new long[4]; randomData.reSeed(1000); for (int i = 0; i < 1000; i++) { @@ -626,7 +626,7 @@ public class RandomDataTest { TestUtils.assertChiSquareAccept(expected, counts, 0.001); // Mean 5 - quartiles = TestUtils.getDistributionQuartiles(new ExponentialDistributionImpl(5)); + quartiles = TestUtils.getDistributionQuartiles(new ExponentialDistribution(5)); counts = new long[4]; randomData.reSeed(1000); for (int i = 0; i < 1000; i++) { @@ -896,7 +896,7 @@ public class RandomDataTest { @Test public void testNextF() throws Exception { - double[] quartiles = TestUtils.getDistributionQuartiles(new FDistributionImpl(12, 5)); + double[] quartiles = TestUtils.getDistributionQuartiles(new FDistribution(12, 5)); long[] counts = new long[4]; randomData.reSeed(1000); for (int i = 0; i < 1000; i++) { @@ -912,7 +912,7 @@ public class RandomDataTest { long[] counts; // Tests shape > 1, one case in the rejection sampling - quartiles = TestUtils.getDistributionQuartiles(new GammaDistributionImpl(4, 2)); + quartiles = TestUtils.getDistributionQuartiles(new GammaDistribution(4, 2)); counts = new long[4]; randomData.reSeed(1000); for (int i = 0; i < 1000; i++) { @@ -922,7 +922,7 @@ public class RandomDataTest { TestUtils.assertChiSquareAccept(expected, counts, 0.001); // Tests shape <= 1, another case in the rejection sampling - quartiles = TestUtils.getDistributionQuartiles(new GammaDistributionImpl(0.3, 3)); + quartiles = TestUtils.getDistributionQuartiles(new GammaDistribution(0.3, 3)); counts = new long[4]; randomData.reSeed(1000); for (int i = 0; i < 1000; i++) {