diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/distribution/FDistribution.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/distribution/FDistribution.java | 341 |
1 files changed, 341 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/distribution/FDistribution.java b/src/main/java/org/apache/commons/math3/distribution/FDistribution.java new file mode 100644 index 0000000..3269f8d --- /dev/null +++ b/src/main/java/org/apache/commons/math3/distribution/FDistribution.java @@ -0,0 +1,341 @@ +/* + * 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.math3.distribution; + +import org.apache.commons.math3.exception.NotStrictlyPositiveException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.random.RandomGenerator; +import org.apache.commons.math3.random.Well19937c; +import org.apache.commons.math3.special.Beta; +import org.apache.commons.math3.util.FastMath; + +/** + * 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 AbstractRealDistribution { + /** + * 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; + + /** Cached numerical variance */ + private double numericalVariance = Double.NaN; + + /** Whether or not the numerical variance has been calculated */ + private boolean numericalVarianceIsCalculated = false; + + /** + * Creates an F distribution using the given degrees of freedom. + * + * <p><b>Note:</b> this constructor will implicitly create an instance of {@link Well19937c} as + * random generator to be used for sampling only (see {@link #sample()} and {@link + * #sample(int)}). In case no sampling is needed for the created distribution, it is advised to + * pass {@code null} as random generator via the appropriate constructors to avoid the + * additional initialisation overhead. + * + * @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); + } + + /** + * Creates an F distribution using the given degrees of freedom and inverse cumulative + * probability accuracy. + * + * <p><b>Note:</b> this constructor will implicitly create an instance of {@link Well19937c} as + * random generator to be used for sampling only (see {@link #sample()} and {@link + * #sample(int)}). In case no sampling is needed for the created distribution, it is advised to + * pass {@code null} as random generator via the appropriate constructors to avoid the + * additional initialisation overhead. + * + * @param numeratorDegreesOfFreedom Numerator degrees of freedom. + * @param denominatorDegreesOfFreedom Denominator degrees of freedom. + * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability + * estimates. + * @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0} or {@code + * denominatorDegreesOfFreedom <= 0}. + * @since 2.1 + */ + public FDistribution( + double numeratorDegreesOfFreedom, + double denominatorDegreesOfFreedom, + double inverseCumAccuracy) + throws NotStrictlyPositiveException { + this( + new Well19937c(), + numeratorDegreesOfFreedom, + denominatorDegreesOfFreedom, + inverseCumAccuracy); + } + + /** + * Creates an F distribution. + * + * @param rng Random number generator. + * @param numeratorDegreesOfFreedom Numerator degrees of freedom. + * @param denominatorDegreesOfFreedom Denominator degrees of freedom. + * @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0} or {@code + * denominatorDegreesOfFreedom <= 0}. + * @since 3.3 + */ + public FDistribution( + RandomGenerator rng, + double numeratorDegreesOfFreedom, + double denominatorDegreesOfFreedom) + throws NotStrictlyPositiveException { + this( + rng, + numeratorDegreesOfFreedom, + denominatorDegreesOfFreedom, + DEFAULT_INVERSE_ABSOLUTE_ACCURACY); + } + + /** + * Creates an F distribution. + * + * @param rng Random number generator. + * @param numeratorDegreesOfFreedom Numerator degrees of freedom. + * @param denominatorDegreesOfFreedom Denominator degrees of freedom. + * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability + * estimates. + * @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0} or {@code + * denominatorDegreesOfFreedom <= 0}. + * @since 3.1 + */ + public FDistribution( + RandomGenerator rng, + double numeratorDegreesOfFreedom, + double denominatorDegreesOfFreedom, + double inverseCumAccuracy) + throws NotStrictlyPositiveException { + super(rng); + + 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) { + return FastMath.exp(logDensity(x)); + } + + /** {@inheritDoc} * */ + @Override + public double logDensity(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 nhalf * logn + + nhalf * logx + - logx + + mhalf * logm + - nhalf * lognxm + - mhalf * lognxm + - Beta.logBeta(nhalf, mhalf); + } + + /** + * {@inheritDoc} + * + * <p>The implementation of this method is based on + * + * <ul> + * <li><a href="http://mathworld.wolfram.com/F-Distribution.html">F-Distribution</a>, equation + * (4). + * </ul> + */ + 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; + } + + /** + * 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} */ + @Override + protected double getSolverAbsoluteAccuracy() { + return solverAbsoluteAccuracy; + } + + /** + * {@inheritDoc} + * + * <p>For denominator degrees of freedom parameter {@code b}, the mean is + * + * <ul> + * <li>if {@code b > 2} then {@code b / (b - 2)}, + * <li>else undefined ({@code Double.NaN}). + * </ul> + */ + public double getNumericalMean() { + final double denominatorDF = getDenominatorDegreesOfFreedom(); + + if (denominatorDF > 2) { + return denominatorDF / (denominatorDF - 2); + } + + return Double.NaN; + } + + /** + * {@inheritDoc} + * + * <p>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>else undefined ({@code Double.NaN}). + * </ul> + */ + public double getNumericalVariance() { + if (!numericalVarianceIsCalculated) { + numericalVariance = calculateNumericalVariance(); + numericalVarianceIsCalculated = true; + } + return numericalVariance; + } + + /** + * used by {@link #getNumericalVariance()} + * + * @return the variance of this distribution + */ + 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} + * + * <p>The lower bound of the support is always 0 no matter the parameters. + * + * @return lower bound of the support (always 0) + */ + public double getSupportLowerBound() { + return 0; + } + + /** + * {@inheritDoc} + * + * <p>The upper bound of the support is always positive infinity no matter the parameters. + * + * @return upper bound of the support (always Double.POSITIVE_INFINITY) + */ + public double getSupportUpperBound() { + return Double.POSITIVE_INFINITY; + } + + /** {@inheritDoc} */ + public boolean isSupportLowerBoundInclusive() { + return false; + } + + /** {@inheritDoc} */ + public boolean isSupportUpperBoundInclusive() { + return false; + } + + /** + * {@inheritDoc} + * + * <p>The support of this distribution is connected. + * + * @return {@code true} + */ + public boolean isSupportConnected() { + return true; + } +} |