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.statistics.distribution;
019
020import org.apache.commons.numbers.gamma.LogBeta;
021import org.apache.commons.numbers.gamma.RegularizedBeta;
022
023/**
024 * Implementation of the F-distribution.
025 *
026 * @see <a href="http://en.wikipedia.org/wiki/F-distribution">F-distribution (Wikipedia)</a>
027 * @see <a href="http://mathworld.wolfram.com/F-Distribution.html">F-distribution (MathWorld)</a>
028 */
029public class FDistribution extends AbstractContinuousDistribution {
030    /** Support lower bound. */
031    private static final double SUPPORT_LO = 0;
032    /** Support upper bound. */
033    private static final double SUPPORT_HI = Double.POSITIVE_INFINITY;
034    /** The minimum degrees of freedom for the denominator when computing the mean. */
035    private static final double MIN_DENOMINATOR_DF_FOR_MEAN = 2.0;
036    /** The minimum degrees of freedom for the denominator when computing the variance. */
037    private static final double MIN_DENOMINATOR_DF_FOR_VARIANCE = 4.0;
038
039    /** The numerator degrees of freedom. */
040    private final double numeratorDegreesOfFreedom;
041    /** The numerator degrees of freedom. */
042    private final double denominatorDegreesOfFreedom;
043
044    /**
045     * Creates a distribution.
046     *
047     * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
048     * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
049     * @throws IllegalArgumentException if {@code numeratorDegreesOfFreedom <= 0} or
050     * {@code denominatorDegreesOfFreedom <= 0}.
051     */
052    public FDistribution(double numeratorDegreesOfFreedom,
053                         double denominatorDegreesOfFreedom) {
054        if (numeratorDegreesOfFreedom <= 0) {
055            throw new DistributionException(DistributionException.NEGATIVE,
056                                            numeratorDegreesOfFreedom);
057        }
058        if (denominatorDegreesOfFreedom <= 0) {
059            throw new DistributionException(DistributionException.NEGATIVE,
060                                            denominatorDegreesOfFreedom);
061        }
062        this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom;
063        this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom;
064    }
065
066    /**
067     * {@inheritDoc}
068     */
069    @Override
070    public double density(double x) {
071        return Math.exp(logDensity(x));
072    }
073
074    /** {@inheritDoc} **/
075    @Override
076    public double logDensity(double x) {
077        if (x <= SUPPORT_LO ||
078            x >= SUPPORT_HI) {
079            return Double.NEGATIVE_INFINITY;
080        }
081
082        final double nhalf = numeratorDegreesOfFreedom / 2;
083        final double mhalf = denominatorDegreesOfFreedom / 2;
084        final double logx = Math.log(x);
085        final double logn = Math.log(numeratorDegreesOfFreedom);
086        final double logm = Math.log(denominatorDegreesOfFreedom);
087        final double lognxm = Math.log(numeratorDegreesOfFreedom * x +
088                denominatorDegreesOfFreedom);
089        return nhalf * logn + nhalf * logx - logx +
090               mhalf * logm - nhalf * lognxm - mhalf * lognxm -
091               LogBeta.value(nhalf, mhalf);
092    }
093
094    /**
095     * {@inheritDoc}
096     *
097     * The implementation of this method is based on
098     * <ul>
099     *  <li>
100     *   <a href="http://mathworld.wolfram.com/F-Distribution.html">
101     *   F-Distribution</a>, equation (4).
102     *  </li>
103     * </ul>
104     */
105    @Override
106    public double cumulativeProbability(double x)  {
107        if (x <= SUPPORT_LO) {
108            return 0;
109        } else if (x >= SUPPORT_HI) {
110            return 1;
111        }
112
113        final double n = numeratorDegreesOfFreedom;
114        final double m = denominatorDegreesOfFreedom;
115
116        return RegularizedBeta.value((n * x) / (m + n * x),
117                                     0.5 * n,
118                                     0.5 * m);
119    }
120
121    /**
122     * Access the numerator degrees of freedom.
123     *
124     * @return the numerator degrees of freedom.
125     */
126    public double getNumeratorDegreesOfFreedom() {
127        return numeratorDegreesOfFreedom;
128    }
129
130    /**
131     * Access the denominator degrees of freedom.
132     *
133     * @return the denominator degrees of freedom.
134     */
135    public double getDenominatorDegreesOfFreedom() {
136        return denominatorDegreesOfFreedom;
137    }
138
139    /**
140     * {@inheritDoc}
141     *
142     * For denominator degrees of freedom parameter {@code b}, the mean is
143     * <ul>
144     *  <li>if {@code b > 2} then {@code b / (b - 2)},</li>
145     *  <li>else undefined ({@code Double.NaN}).
146     * </ul>
147     */
148    @Override
149    public double getMean() {
150        final double denominatorDF = getDenominatorDegreesOfFreedom();
151
152        if (denominatorDF > MIN_DENOMINATOR_DF_FOR_MEAN) {
153            return denominatorDF / (denominatorDF - 2);
154        }
155
156        return Double.NaN;
157    }
158
159    /**
160     * {@inheritDoc}
161     *
162     * For numerator degrees of freedom parameter {@code a} and denominator
163     * degrees of freedom parameter {@code b}, the variance is
164     * <ul>
165     *  <li>
166     *    if {@code b > 4} then
167     *    {@code [2 * b^2 * (a + b - 2)] / [a * (b - 2)^2 * (b - 4)]},
168     *  </li>
169     *  <li>else undefined ({@code Double.NaN}).
170     * </ul>
171     */
172    @Override
173    public double getVariance() {
174        final double denominatorDF = getDenominatorDegreesOfFreedom();
175
176        if (denominatorDF > MIN_DENOMINATOR_DF_FOR_VARIANCE) {
177            final double numeratorDF = getNumeratorDegreesOfFreedom();
178            final double denomDFMinusTwo = denominatorDF - 2;
179
180            return (2 * (denominatorDF * denominatorDF) * (numeratorDF + denominatorDF - 2)) /
181                   (numeratorDF * (denomDFMinusTwo * denomDFMinusTwo) * (denominatorDF - 4));
182        }
183
184        return Double.NaN;
185    }
186
187    /**
188     * {@inheritDoc}
189     *
190     * The lower bound of the support is always 0 no matter the parameters.
191     *
192     * @return lower bound of the support (always 0)
193     */
194    @Override
195    public double getSupportLowerBound() {
196        return SUPPORT_LO;
197    }
198
199    /**
200     * {@inheritDoc}
201     *
202     * The upper bound of the support is always positive infinity
203     * no matter the parameters.
204     *
205     * @return upper bound of the support (always Double.POSITIVE_INFINITY)
206     */
207    @Override
208    public double getSupportUpperBound() {
209        return SUPPORT_HI;
210    }
211
212    /**
213     * {@inheritDoc}
214     *
215     * The support of this distribution is connected.
216     *
217     * @return {@code true}
218     */
219    @Override
220    public boolean isSupportConnected() {
221        return true;
222    }
223}