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    }