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
018package org.apache.commons.math3.distribution;
019
020import org.apache.commons.math3.exception.NotStrictlyPositiveException;
021import org.apache.commons.math3.exception.util.LocalizedFormats;
022import org.apache.commons.math3.random.RandomGenerator;
023import org.apache.commons.math3.random.Well19937c;
024import org.apache.commons.math3.special.Beta;
025import org.apache.commons.math3.util.FastMath;
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 */
033public class FDistribution extends AbstractRealDistribution {
034    /**
035     * Default inverse cumulative probability accuracy.
036     * @since 2.1
037     */
038    public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
039    /** Serializable version identifier. */
040    private static final long serialVersionUID = -8516354193418641566L;
041    /** The numerator degrees of freedom. */
042    private final double numeratorDegreesOfFreedom;
043    /** The numerator degrees of freedom. */
044    private final double denominatorDegreesOfFreedom;
045    /** Inverse cumulative probability accuracy. */
046    private final double solverAbsoluteAccuracy;
047    /** Cached numerical variance */
048    private double numericalVariance = Double.NaN;
049    /** Whether or not the numerical variance has been calculated */
050    private boolean numericalVarianceIsCalculated = false;
051
052    /**
053     * Creates an F distribution using the given degrees of freedom.
054     * <p>
055     * <b>Note:</b> this constructor will implicitly create an instance of
056     * {@link Well19937c} as random generator to be used for sampling only (see
057     * {@link #sample()} and {@link #sample(int)}). In case no sampling is
058     * needed for the created distribution, it is advised to pass {@code null}
059     * as random generator via the appropriate constructors to avoid the
060     * additional initialisation overhead.
061     *
062     * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
063     * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
064     * @throws NotStrictlyPositiveException if
065     * {@code numeratorDegreesOfFreedom <= 0} or
066     * {@code denominatorDegreesOfFreedom <= 0}.
067     */
068    public FDistribution(double numeratorDegreesOfFreedom,
069                         double denominatorDegreesOfFreedom)
070        throws NotStrictlyPositiveException {
071        this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom,
072             DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
073    }
074
075    /**
076     * Creates an F distribution using the given degrees of freedom
077     * and inverse cumulative probability accuracy.
078     * <p>
079     * <b>Note:</b> this constructor will implicitly create an instance of
080     * {@link Well19937c} as random generator to be used for sampling only (see
081     * {@link #sample()} and {@link #sample(int)}). In case no sampling is
082     * needed for the created distribution, it is advised to pass {@code null}
083     * as random generator via the appropriate constructors to avoid the
084     * additional initialisation overhead.
085     *
086     * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
087     * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
088     * @param inverseCumAccuracy the maximum absolute error in inverse
089     * cumulative probability estimates.
090     * @throws NotStrictlyPositiveException if
091     * {@code numeratorDegreesOfFreedom <= 0} or
092     * {@code denominatorDegreesOfFreedom <= 0}.
093     * @since 2.1
094     */
095    public FDistribution(double numeratorDegreesOfFreedom,
096                         double denominatorDegreesOfFreedom,
097                         double inverseCumAccuracy)
098        throws NotStrictlyPositiveException {
099        this(new Well19937c(), numeratorDegreesOfFreedom,
100             denominatorDegreesOfFreedom, inverseCumAccuracy);
101    }
102
103    /**
104     * Creates an F distribution.
105     *
106     * @param rng Random number generator.
107     * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
108     * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
109     * @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0} or
110     * {@code denominatorDegreesOfFreedom <= 0}.
111     * @since 3.3
112     */
113    public FDistribution(RandomGenerator rng,
114                         double numeratorDegreesOfFreedom,
115                         double denominatorDegreesOfFreedom)
116        throws NotStrictlyPositiveException {
117        this(rng, numeratorDegreesOfFreedom, denominatorDegreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
118    }
119
120    /**
121     * Creates an F distribution.
122     *
123     * @param rng Random number generator.
124     * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
125     * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
126     * @param inverseCumAccuracy the maximum absolute error in inverse
127     * cumulative probability estimates.
128     * @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0} or
129     * {@code denominatorDegreesOfFreedom <= 0}.
130     * @since 3.1
131     */
132    public FDistribution(RandomGenerator rng,
133                         double numeratorDegreesOfFreedom,
134                         double denominatorDegreesOfFreedom,
135                         double inverseCumAccuracy)
136        throws NotStrictlyPositiveException {
137        super(rng);
138
139        if (numeratorDegreesOfFreedom <= 0) {
140            throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
141                                                   numeratorDegreesOfFreedom);
142        }
143        if (denominatorDegreesOfFreedom <= 0) {
144            throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
145                                                   denominatorDegreesOfFreedom);
146        }
147        this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom;
148        this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom;
149        solverAbsoluteAccuracy = inverseCumAccuracy;
150    }
151
152    /**
153     * {@inheritDoc}
154     *
155     * @since 2.1
156     */
157    public double density(double x) {
158        return FastMath.exp(logDensity(x));
159    }
160
161    /** {@inheritDoc} **/
162    @Override
163    public double logDensity(double x) {
164        final double nhalf = numeratorDegreesOfFreedom / 2;
165        final double mhalf = denominatorDegreesOfFreedom / 2;
166        final double logx = FastMath.log(x);
167        final double logn = FastMath.log(numeratorDegreesOfFreedom);
168        final double logm = FastMath.log(denominatorDegreesOfFreedom);
169        final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x +
170                denominatorDegreesOfFreedom);
171        return nhalf * logn + nhalf * logx - logx +
172               mhalf * logm - nhalf * lognxm - mhalf * lognxm -
173               Beta.logBeta(nhalf, mhalf);
174    }
175
176    /**
177     * {@inheritDoc}
178     *
179     * The implementation of this method is based on
180     * <ul>
181     *  <li>
182     *   <a href="http://mathworld.wolfram.com/F-Distribution.html">
183     *   F-Distribution</a>, equation (4).
184     *  </li>
185     * </ul>
186     */
187    public double cumulativeProbability(double x)  {
188        double ret;
189        if (x <= 0) {
190            ret = 0;
191        } else {
192            double n = numeratorDegreesOfFreedom;
193            double m = denominatorDegreesOfFreedom;
194
195            ret = Beta.regularizedBeta((n * x) / (m + n * x),
196                0.5 * n,
197                0.5 * m);
198        }
199        return ret;
200    }
201
202    /**
203     * Access the numerator degrees of freedom.
204     *
205     * @return the numerator degrees of freedom.
206     */
207    public double getNumeratorDegreesOfFreedom() {
208        return numeratorDegreesOfFreedom;
209    }
210
211    /**
212     * Access the denominator degrees of freedom.
213     *
214     * @return the denominator degrees of freedom.
215     */
216    public double getDenominatorDegreesOfFreedom() {
217        return denominatorDegreesOfFreedom;
218    }
219
220    /** {@inheritDoc} */
221    @Override
222    protected double getSolverAbsoluteAccuracy() {
223        return solverAbsoluteAccuracy;
224    }
225
226    /**
227     * {@inheritDoc}
228     *
229     * For denominator degrees of freedom parameter {@code b}, the mean is
230     * <ul>
231     *  <li>if {@code b > 2} then {@code b / (b - 2)},</li>
232     *  <li>else undefined ({@code Double.NaN}).
233     * </ul>
234     */
235    public double getNumericalMean() {
236        final double denominatorDF = getDenominatorDegreesOfFreedom();
237
238        if (denominatorDF > 2) {
239            return denominatorDF / (denominatorDF - 2);
240        }
241
242        return Double.NaN;
243    }
244
245    /**
246     * {@inheritDoc}
247     *
248     * For numerator degrees of freedom parameter {@code a} and denominator
249     * degrees of freedom parameter {@code b}, the variance is
250     * <ul>
251     *  <li>
252     *    if {@code b > 4} then
253     *    {@code [2 * b^2 * (a + b - 2)] / [a * (b - 2)^2 * (b - 4)]},
254     *  </li>
255     *  <li>else undefined ({@code Double.NaN}).
256     * </ul>
257     */
258    public double getNumericalVariance() {
259        if (!numericalVarianceIsCalculated) {
260            numericalVariance = calculateNumericalVariance();
261            numericalVarianceIsCalculated = true;
262        }
263        return numericalVariance;
264    }
265
266    /**
267     * used by {@link #getNumericalVariance()}
268     *
269     * @return the variance of this distribution
270     */
271    protected double calculateNumericalVariance() {
272        final double denominatorDF = getDenominatorDegreesOfFreedom();
273
274        if (denominatorDF > 4) {
275            final double numeratorDF = getNumeratorDegreesOfFreedom();
276            final double denomDFMinusTwo = denominatorDF - 2;
277
278            return ( 2 * (denominatorDF * denominatorDF) * (numeratorDF + denominatorDF - 2) ) /
279                   ( (numeratorDF * (denomDFMinusTwo * denomDFMinusTwo) * (denominatorDF - 4)) );
280        }
281
282        return Double.NaN;
283    }
284
285    /**
286     * {@inheritDoc}
287     *
288     * The lower bound of the support is always 0 no matter the parameters.
289     *
290     * @return lower bound of the support (always 0)
291     */
292    public double getSupportLowerBound() {
293        return 0;
294    }
295
296    /**
297     * {@inheritDoc}
298     *
299     * The upper bound of the support is always positive infinity
300     * no matter the parameters.
301     *
302     * @return upper bound of the support (always Double.POSITIVE_INFINITY)
303     */
304    public double getSupportUpperBound() {
305        return Double.POSITIVE_INFINITY;
306    }
307
308    /** {@inheritDoc} */
309    public boolean isSupportLowerBoundInclusive() {
310        return false;
311    }
312
313    /** {@inheritDoc} */
314    public boolean isSupportUpperBoundInclusive() {
315        return false;
316    }
317
318    /**
319     * {@inheritDoc}
320     *
321     * The support of this distribution is connected.
322     *
323     * @return {@code true}
324     */
325    public boolean isSupportConnected() {
326        return true;
327    }
328}