    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    package org.apache.commons.math.distribution;
    018
    019    import java.io.Serializable;
    020
    021    import org.apache.commons.math.MathException;
    022    import org.apache.commons.math.MathRuntimeException;
    023    import org.apache.commons.math.exception.util.LocalizedFormats;
    024    import org.apache.commons.math.special.Beta;
    025    import org.apache.commons.math.util.FastMath;
    026
    027    /**
    028     * Default implementation of
    029     * {@link org.apache.commons.math.distribution.FDistribution}.
    030     *
    031     * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $
    032     */
    033    public class FDistributionImpl
    034        extends AbstractContinuousDistribution
    035        implements FDistribution, Serializable  {
    036
    037        /**
    038         * Default inverse cumulative probability accuracy
    039         * @since 2.1
    040         */
    041        public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
    042
    043        /** Serializable version identifier */
    044        private static final long serialVersionUID = -8516354193418641566L;
    045
    046        /** The numerator degrees of freedom*/
    047        private double numeratorDegreesOfFreedom;
    048
    049        /** The numerator degrees of freedom*/
    050        private double denominatorDegreesOfFreedom;
    051
    052        /** Inverse cumulative probability accuracy */
    053        private final double solverAbsoluteAccuracy;
    054
    055        /**
    056         * Create a F distribution using the given degrees of freedom.
    057         * @param numeratorDegreesOfFreedom the numerator degrees of freedom.
    058         * @param denominatorDegreesOfFreedom the denominator degrees of freedom.
    059         */
    060        public FDistributionImpl(double numeratorDegreesOfFreedom,
    061                                 double denominatorDegreesOfFreedom) {
    062            this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    063        }
    064
    065        /**
    066         * Create a F distribution using the given degrees of freedom and inverse cumulative probability accuracy.
    067         * @param numeratorDegreesOfFreedom the numerator degrees of freedom.
    068         * @param denominatorDegreesOfFreedom the denominator degrees of freedom.
    069         * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates
    070         * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY})
    071         * @since 2.1
    072         */
    073        public FDistributionImpl(double numeratorDegreesOfFreedom, double denominatorDegreesOfFreedom,
    074                double inverseCumAccuracy) {
    075            super();
    076            setNumeratorDegreesOfFreedomInternal(numeratorDegreesOfFreedom);
    077            setDenominatorDegreesOfFreedomInternal(denominatorDegreesOfFreedom);
    078            solverAbsoluteAccuracy = inverseCumAccuracy;
    079        }
    080
    081        /**
    082         * Returns the probability density for a particular point.
    083         *
    084         * @param x The point at which the density should be computed.
    085         * @return The pdf at point x.
    086         * @since 2.1
    087         */
    088        @Override
    089        public double density(double x) {
    090            final double nhalf = numeratorDegreesOfFreedom / 2;
    091            final double mhalf = denominatorDegreesOfFreedom / 2;
    092            final double logx = FastMath.log(x);
    093            final double logn = FastMath.log(numeratorDegreesOfFreedom);
    094            final double logm = FastMath.log(denominatorDegreesOfFreedom);
    095            final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x + denominatorDegreesOfFreedom);
    096            return FastMath.exp(nhalf*logn + nhalf*logx - logx + mhalf*logm - nhalf*lognxm -
    097                   mhalf*lognxm - Beta.logBeta(nhalf, mhalf));
    098        }
    099
    100        /**
    101         * For this distribution, X, this method returns P(X &lt; x).
    102         *
    103         * The implementation of this method is based on:
    104         * <ul>
    105         * <li>
    106         * <a href="http://mathworld.wolfram.com/F-Distribution.html">
    107         * F-Distribution</a>, equation (4).</li>
    108         * </ul>
    109         *
    110         * @param x the value at which the CDF is evaluated.
    111         * @return CDF for this distribution.
    112         * @throws MathException if the cumulative probability can not be
    113         *            computed due to convergence or other numerical errors.
    114         */
    115        public double cumulativeProbability(double x) throws MathException {
    116            double ret;
    117            if (x <= 0.0) {
    118                ret = 0.0;
    119            } else {
    120                double n = numeratorDegreesOfFreedom;
    121                double m = denominatorDegreesOfFreedom;
    122
    123                ret = Beta.regularizedBeta((n * x) / (m + n * x),
    124                    0.5 * n,
    125                    0.5 * m);
    126            }
    127            return ret;
    128        }
    129
    130        /**
    131         * For this distribution, X, this method returns the critical point x, such
    132         * that P(X &lt; x) = <code>p</code>.
    133         * <p>
    134         * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.</p>
    135         *
    136         * @param p the desired probability
    137         * @return x, such that P(X &lt; x) = <code>p</code>
    138         * @throws MathException if the inverse cumulative probability can not be
    139         *         computed due to convergence or other numerical errors.
    140         * @throws IllegalArgumentException if <code>p</code> is not a valid
    141         *         probability.
    142         */
    143        @Override
    144        public double inverseCumulativeProbability(final double p)
    145            throws MathException {
    146            if (p == 0) {
    147                return 0d;
    148            }
    149            if (p == 1) {
    150                return Double.POSITIVE_INFINITY;
    151            }
    152            return super.inverseCumulativeProbability(p);
    153        }
    154
    155        /**
    156         * Access the domain value lower bound, based on <code>p</code>, used to
    157         * bracket a CDF root.  This method is used by
    158         * {@link #inverseCumulativeProbability(double)} to find critical values.
    159         *
    160         * @param p the desired probability for the critical value
    161         * @return domain value lower bound, i.e.
    162         *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>
    163         */
    164        @Override
    165        protected double getDomainLowerBound(double p) {
    166            return 0.0;
    167        }
    168
    169        /**
    170         * Access the domain value upper bound, based on <code>p</code>, used to
    171         * bracket a CDF root.  This method is used by
    172         * {@link #inverseCumulativeProbability(double)} to find critical values.
    173         *
    174         * @param p the desired probability for the critical value
    175         * @return domain value upper bound, i.e.
    176         *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code>
    177         */
    178        @Override
    179        protected double getDomainUpperBound(double p) {
    180            return Double.MAX_VALUE;
    181        }
    182
    183        /**
    184         * Access the initial domain value, based on <code>p</code>, used to
    185         * bracket a CDF root.  This method is used by
    186         * {@link #inverseCumulativeProbability(double)} to find critical values.
    187         *
    188         * @param p the desired probability for the critical value
    189         * @return initial domain value
    190         */
    191        @Override
    192        protected double getInitialDomain(double p) {
    193            double ret = 1.0;
    194            double d = denominatorDegreesOfFreedom;
    195            if (d > 2.0) {
    196                // use mean
    197                ret = d / (d - 2.0);
    198            }
    199            return ret;
    200        }
    201
    202        /**
    203         * Modify the numerator degrees of freedom.
    204         * @param degreesOfFreedom the new numerator degrees of freedom.
    205         * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
    206         *         positive.
    207         * @deprecated as of 2.1 (class will become immutable in 3.0)
    208         */
    209        @Deprecated
    210        public void setNumeratorDegreesOfFreedom(double degreesOfFreedom) {
    211            setNumeratorDegreesOfFreedomInternal(degreesOfFreedom);
    212        }
    213
    214        /**
    215         * Modify the numerator degrees of freedom.
    216         * @param degreesOfFreedom the new numerator degrees of freedom.
    217         * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
    218         *         positive.
    219         */
    220        private void setNumeratorDegreesOfFreedomInternal(double degreesOfFreedom) {
    221            if (degreesOfFreedom <= 0.0) {
    222                throw MathRuntimeException.createIllegalArgumentException(
    223                      LocalizedFormats.NOT_POSITIVE_DEGREES_OF_FREEDOM, degreesOfFreedom);
    224            }
    225            this.numeratorDegreesOfFreedom = degreesOfFreedom;
    226        }
    227
    228        /**
    229         * Access the numerator degrees of freedom.
    230         * @return the numerator degrees of freedom.
    231         */
    232        public double getNumeratorDegreesOfFreedom() {
    233            return numeratorDegreesOfFreedom;
    234        }
    235
    236        /**
    237         * Modify the denominator degrees of freedom.
    238         * @param degreesOfFreedom the new denominator degrees of freedom.
    239         * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
    240         *         positive.
    241         * @deprecated as of 2.1 (class will become immutable in 3.0)
    242         */
    243        @Deprecated
    244        public void setDenominatorDegreesOfFreedom(double degreesOfFreedom) {
    245            setDenominatorDegreesOfFreedomInternal(degreesOfFreedom);
    246        }
    247
    248        /**
    249         * Modify the denominator degrees of freedom.
    250         * @param degreesOfFreedom the new denominator degrees of freedom.
    251         * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
    252         *         positive.
    253         */
    254        private void setDenominatorDegreesOfFreedomInternal(double degreesOfFreedom) {
    255            if (degreesOfFreedom <= 0.0) {
    256                throw MathRuntimeException.createIllegalArgumentException(
    257                      LocalizedFormats.NOT_POSITIVE_DEGREES_OF_FREEDOM, degreesOfFreedom);
    258            }
    259            this.denominatorDegreesOfFreedom = degreesOfFreedom;
    260        }
    261
    262        /**
    263         * Access the denominator degrees of freedom.
    264         * @return the denominator degrees of freedom.
    265         */
    266        public double getDenominatorDegreesOfFreedom() {
    267            return denominatorDegreesOfFreedom;
    268        }
    269
    270        /**
    271         * Return the absolute accuracy setting of the solver used to estimate
    272         * inverse cumulative probabilities.
    273         *
    274         * @return the solver absolute accuracy
    275         * @since 2.1
    276         */
    277        @Override
    278        protected double getSolverAbsoluteAccuracy() {
    279            return solverAbsoluteAccuracy;
    280        }
    281
    282        /**
    283         * Returns the lower bound of the support for the distribution.
    284         *
    285         * The lower bound of the support is always 0, regardless of the parameters.
    286         *
    287         * @return lower bound of the support (always 0)
    288         * @since 2.2
    289         */
    290        public double getSupportLowerBound() {
    291            return 0;
    292        }
    293
    294        /**
    295         * Returns the upper bound of the support for the distribution.
    296         *
    297         * The upper bound of the support is always positive infinity,
    298         * regardless of the parameters.
    299         *
    300         * @return upper bound of the support (always Double.POSITIVE_INFINITY)
    301         * @since 2.2
    302         */
    303        public double getSupportUpperBound() {
    304            return Double.POSITIVE_INFINITY;
    305        }
    306
    307        /**
    308         * Returns the mean of the distribution.
    309         *
    310         * For denominator degrees of freedom parameter <code>b</code>,
    311         * the mean is
    312         * <ul>
    313         *  <li>if <code>b &gt; 2</code> then <code>b / (b - 2)</code></li>
    314         *  <li>else <code>undefined</code>
    315         * </ul>
    316         *
    317         * @return the mean
    318         * @since 2.2
    319         */
    320        public double getNumericalMean() {
    321            final double denominatorDF = getDenominatorDegreesOfFreedom();
    322
    323            if (denominatorDF > 2) {
    324                return denominatorDF / (denominatorDF - 2);
    325            }
    326
    327            return Double.NaN;
    328        }
    329
    330        /**
    331         * Returns the variance of the distribution.
    332         *
    333         * For numerator degrees of freedom parameter <code>a</code>
    334         * and denominator degrees of freedom parameter <code>b</code>,
    335         * the variance is
    336         * <ul>
    337         *  <li>
    338         *    if <code>b &gt; 4</code> then
    339         *    <code>[ 2 * b^2 * (a + b - 2) ] / [ a * (b - 2)^2 * (b - 4) ]</code>
    340         *  </li>
    341         *  <li>else <code>undefined</code>
    342         * </ul>
    343         *
    344         * @return the variance
    345         * @since 2.2
    346         */
    347        public double getNumericalVariance() {
    348            final double denominatorDF = getDenominatorDegreesOfFreedom();
    349
    350            if (denominatorDF > 4) {
    351                final double numeratorDF = getNumeratorDegreesOfFreedom();
    352                final double denomDFMinusTwo = denominatorDF - 2;
    353
    354                return ( 2 * (denominatorDF * denominatorDF) * (numeratorDF + denominatorDF - 2) ) /
    355                        ( (numeratorDF * (denomDFMinusTwo * denomDFMinusTwo) * (denominatorDF - 4)) );
    356            }
    357
    358            return Double.NaN;
    359        }
    360    }