    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017
    018    package org.apache.commons.math3.distribution;
    019
    020    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    021    import org.apache.commons.math3.exception.util.LocalizedFormats;
    022    import org.apache.commons.math3.special.Beta;
    023    import org.apache.commons.math3.util.FastMath;
    024    import org.apache.commons.math3.random.RandomGenerator;
    025    import org.apache.commons.math3.random.Well19937c;
    026
    027    /**
    028     * Implementation of the F-distribution.
    029     *
    030     * @see <a href="http://en.wikipedia.org/wiki/F-distribution">F-distribution (Wikipedia)</a>
    031     * @see <a href="http://mathworld.wolfram.com/F-Distribution.html">F-distribution (MathWorld)</a>
    032     * @version $Id: FDistribution.java 1416643 2012-12-03 19:37:14Z tn $
    033     */
    034    public class FDistribution extends AbstractRealDistribution {
    035        /**
    036         * Default inverse cumulative probability accuracy.
    037         * @since 2.1
    038         */
    039        public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
    040        /** Serializable version identifier. */
    041        private static final long serialVersionUID = -8516354193418641566L;
    042        /** The numerator degrees of freedom. */
    043        private final double numeratorDegreesOfFreedom;
    044        /** The numerator degrees of freedom. */
    045        private final double denominatorDegreesOfFreedom;
    046        /** Inverse cumulative probability accuracy. */
    047        private final double solverAbsoluteAccuracy;
    048        /** Cached numerical variance */
    049        private double numericalVariance = Double.NaN;
    050        /** Whether or not the numerical variance has been calculated */
    051        private boolean numericalVarianceIsCalculated = false;
    052
    053        /**
    054         * Creates an F distribution using the given degrees of freedom.
    055         *
    056         * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
    057         * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
    058         * @throws NotStrictlyPositiveException if
    059         * {@code numeratorDegreesOfFreedom <= 0} or
    060         * {@code denominatorDegreesOfFreedom <= 0}.
    061         */
    062        public FDistribution(double numeratorDegreesOfFreedom,
    063                             double denominatorDegreesOfFreedom)
    064            throws NotStrictlyPositiveException {
    065            this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom,
    066                 DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    067        }
    068
    069        /**
    070         * Creates an F distribution using the given degrees of freedom
    071         * and inverse cumulative probability accuracy.
    072         *
    073         * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
    074         * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
    075         * @param inverseCumAccuracy the maximum absolute error in inverse
    076         * cumulative probability estimates.
    077         * @throws NotStrictlyPositiveException if
    078         * {@code numeratorDegreesOfFreedom <= 0} or
    079         * {@code denominatorDegreesOfFreedom <= 0}.
    080         * @since 2.1
    081         */
    082        public FDistribution(double numeratorDegreesOfFreedom,
    083                             double denominatorDegreesOfFreedom,
    084                             double inverseCumAccuracy)
    085            throws NotStrictlyPositiveException {
    086            this(new Well19937c(), numeratorDegreesOfFreedom,
    087                 denominatorDegreesOfFreedom, inverseCumAccuracy);
    088        }
    089
    090        /**
    091         * Creates an F distribution.
    092         *
    093         * @param rng Random number generator.
    094         * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
    095         * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
    096         * @param inverseCumAccuracy the maximum absolute error in inverse
    097         * cumulative probability estimates.
    098         * @throws NotStrictlyPositiveException if
    099         * {@code numeratorDegreesOfFreedom <= 0} or
    100         * {@code denominatorDegreesOfFreedom <= 0}.
    101         * @since 3.1
    102         */
    103        public FDistribution(RandomGenerator rng,
    104                             double numeratorDegreesOfFreedom,
    105                             double denominatorDegreesOfFreedom,
    106                             double inverseCumAccuracy)
    107            throws NotStrictlyPositiveException {
    108            super(rng);
    109
    110            if (numeratorDegreesOfFreedom <= 0) {
    111                throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
    112                                                       numeratorDegreesOfFreedom);
    113            }
    114            if (denominatorDegreesOfFreedom <= 0) {
    115                throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
    116                                                       denominatorDegreesOfFreedom);
    117            }
    118            this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom;
    119            this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom;
    120            solverAbsoluteAccuracy = inverseCumAccuracy;
    121        }
    122
    123        /**
    124         * {@inheritDoc}
    125         *
    126         * @since 2.1
    127         */
    128        public double density(double x) {
    129            final double nhalf = numeratorDegreesOfFreedom / 2;
    130            final double mhalf = denominatorDegreesOfFreedom / 2;
    131            final double logx = FastMath.log(x);
    132            final double logn = FastMath.log(numeratorDegreesOfFreedom);
    133            final double logm = FastMath.log(denominatorDegreesOfFreedom);
    134            final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x +
    135                                               denominatorDegreesOfFreedom);
    136            return FastMath.exp(nhalf * logn + nhalf * logx - logx +
    137                                mhalf * logm - nhalf * lognxm - mhalf * lognxm -
    138                                Beta.logBeta(nhalf, mhalf));
    139        }
    140
    141        /**
    142         * {@inheritDoc}
    143         *
    144         * The implementation of this method is based on
    145         * <ul>
    146         *  <li>
    147         *   <a href="http://mathworld.wolfram.com/F-Distribution.html">
    148         *   F-Distribution</a>, equation (4).
    149         *  </li>
    150         * </ul>
    151         */
    152        public double cumulativeProbability(double x)  {
    153            double ret;
    154            if (x <= 0) {
    155                ret = 0;
    156            } else {
    157                double n = numeratorDegreesOfFreedom;
    158                double m = denominatorDegreesOfFreedom;
    159
    160                ret = Beta.regularizedBeta((n * x) / (m + n * x),
    161                    0.5 * n,
    162                    0.5 * m);
    163            }
    164            return ret;
    165        }
    166
    167        /**
    168         * Access the numerator degrees of freedom.
    169         *
    170         * @return the numerator degrees of freedom.
    171         */
    172        public double getNumeratorDegreesOfFreedom() {
    173            return numeratorDegreesOfFreedom;
    174        }
    175
    176        /**
    177         * Access the denominator degrees of freedom.
    178         *
    179         * @return the denominator degrees of freedom.
    180         */
    181        public double getDenominatorDegreesOfFreedom() {
    182            return denominatorDegreesOfFreedom;
    183        }
    184
    185        /** {@inheritDoc} */
    186        @Override
    187        protected double getSolverAbsoluteAccuracy() {
    188            return solverAbsoluteAccuracy;
    189        }
    190
    191        /**
    192         * {@inheritDoc}
    193         *
    194         * For denominator degrees of freedom parameter {@code b}, the mean is
    195         * <ul>
    196         *  <li>if {@code b > 2} then {@code b / (b - 2)},</li>
    197         *  <li>else undefined ({@code Double.NaN}).
    198         * </ul>
    199         */
    200        public double getNumericalMean() {
    201            final double denominatorDF = getDenominatorDegreesOfFreedom();
    202
    203            if (denominatorDF > 2) {
    204                return denominatorDF / (denominatorDF - 2);
    205            }
    206
    207            return Double.NaN;
    208        }
    209
    210        /**
    211         * {@inheritDoc}
    212         *
    213         * For numerator degrees of freedom parameter {@code a} and denominator
    214         * degrees of freedom parameter {@code b}, the variance is
    215         * <ul>
    216         *  <li>
    217         *    if {@code b > 4} then
    218         *    {@code [2 * b^2 * (a + b - 2)] / [a * (b - 2)^2 * (b - 4)]},
    219         *  </li>
    220         *  <li>else undefined ({@code Double.NaN}).
    221         * </ul>
    222         */
    223        public double getNumericalVariance() {
    224            if (!numericalVarianceIsCalculated) {
    225                numericalVariance = calculateNumericalVariance();
    226                numericalVarianceIsCalculated = true;
    227            }
    228            return numericalVariance;
    229        }
    230
    231        /**
    232         * used by {@link #getNumericalVariance()}
    233         *
    234         * @return the variance of this distribution
    235         */
    236        protected double calculateNumericalVariance() {
    237            final double denominatorDF = getDenominatorDegreesOfFreedom();
    238
    239            if (denominatorDF > 4) {
    240                final double numeratorDF = getNumeratorDegreesOfFreedom();
    241                final double denomDFMinusTwo = denominatorDF - 2;
    242
    243                return ( 2 * (denominatorDF * denominatorDF) * (numeratorDF + denominatorDF - 2) ) /
    244                       ( (numeratorDF * (denomDFMinusTwo * denomDFMinusTwo) * (denominatorDF - 4)) );
    245            }
    246
    247            return Double.NaN;
    248        }
    249
    250        /**
    251         * {@inheritDoc}
    252         *
    253         * The lower bound of the support is always 0 no matter the parameters.
    254         *
    255         * @return lower bound of the support (always 0)
    256         */
    257        public double getSupportLowerBound() {
    258            return 0;
    259        }
    260
    261        /**
    262         * {@inheritDoc}
    263         *
    264         * The upper bound of the support is always positive infinity
    265         * no matter the parameters.
    266         *
    267         * @return upper bound of the support (always Double.POSITIVE_INFINITY)
    268         */
    269        public double getSupportUpperBound() {
    270            return Double.POSITIVE_INFINITY;
    271        }
    272
    273        /** {@inheritDoc} */
    274        public boolean isSupportLowerBoundInclusive() {
    275            return false;
    276        }
    277
    278        /** {@inheritDoc} */
    279        public boolean isSupportUpperBoundInclusive() {
    280            return false;
    281        }
    282
    283        /**
    284         * {@inheritDoc}
    285         *
    286         * The support of this distribution is connected.
    287         *
    288         * @return {@code true}
    289         */
    290        public boolean isSupportConnected() {
    291            return true;
    292        }
    293    }