1 /*
2  * Licensed to the Apache Software Foundation (ASF) under one or more
3  * contributor license agreements.  See the NOTICE file distributed with
4  * this work for additional information regarding copyright ownership.
5  * The ASF licenses this file to You under the Apache License, Version 2.0
6  * (the "License"); you may not use this file except in compliance with
7  * the License.  You may obtain a copy of the License at
8  *
9  *      http://www.apache.org/licenses/LICENSE-2.0
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  */
17 package org.apache.commons.math.distribution;
18 
19 import java.io.Serializable;
20 
21 import org.apache.commons.math.MathException;
22 import org.apache.commons.math.MathRuntimeException;
23 import org.apache.commons.math.exception.util.LocalizedFormats;
24 import org.apache.commons.math.special.Beta;
25 import org.apache.commons.math.util.FastMath;
26 
27 /**
28  * Default implementation of
29  * {@link org.apache.commons.math.distribution.FDistribution}.
30  *
31  * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $
32  */
33 public class FDistributionImpl
34     extends AbstractContinuousDistribution
35     implements FDistribution, Serializable  {
36 
37     /**
38      * Default inverse cumulative probability accuracy
39      * @since 2.1
40      */
41     public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
42 
43     /** Serializable version identifier */
44     private static final long serialVersionUID = -8516354193418641566L;
45 
46     /** The numerator degrees of freedom*/
47     private double numeratorDegreesOfFreedom;
48 
49     /** The numerator degrees of freedom*/
50     private double denominatorDegreesOfFreedom;
51 
52     /** Inverse cumulative probability accuracy */
53     private final double solverAbsoluteAccuracy;
54 
55     /**
56      * Create a F distribution using the given degrees of freedom.
57      * @param numeratorDegreesOfFreedom the numerator degrees of freedom.
58      * @param denominatorDegreesOfFreedom the denominator degrees of freedom.
59      */
FDistributionImpl(double numeratorDegreesOfFreedom, double denominatorDegreesOfFreedom)60     public FDistributionImpl(double numeratorDegreesOfFreedom,
61                              double denominatorDegreesOfFreedom) {
62         this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
63     }
64 
65     /**
66      * Create a F distribution using the given degrees of freedom and inverse cumulative probability accuracy.
67      * @param numeratorDegreesOfFreedom the numerator degrees of freedom.
68      * @param denominatorDegreesOfFreedom the denominator degrees of freedom.
69      * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates
70      * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY})
71      * @since 2.1
72      */
FDistributionImpl(double numeratorDegreesOfFreedom, double denominatorDegreesOfFreedom, double inverseCumAccuracy)73     public FDistributionImpl(double numeratorDegreesOfFreedom, double denominatorDegreesOfFreedom,
74             double inverseCumAccuracy) {
75         super();
76         setNumeratorDegreesOfFreedomInternal(numeratorDegreesOfFreedom);
77         setDenominatorDegreesOfFreedomInternal(denominatorDegreesOfFreedom);
78         solverAbsoluteAccuracy = inverseCumAccuracy;
79     }
80 
81     /**
82      * Returns the probability density for a particular point.
83      *
84      * @param x The point at which the density should be computed.
85      * @return The pdf at point x.
86      * @since 2.1
87      */
88     @Override
density(double x)89     public double density(double x) {
90         final double nhalf = numeratorDegreesOfFreedom / 2;
91         final double mhalf = denominatorDegreesOfFreedom / 2;
92         final double logx = FastMath.log(x);
93         final double logn = FastMath.log(numeratorDegreesOfFreedom);
94         final double logm = FastMath.log(denominatorDegreesOfFreedom);
95         final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x + denominatorDegreesOfFreedom);
96         return FastMath.exp(nhalf*logn + nhalf*logx - logx + mhalf*logm - nhalf*lognxm -
97                mhalf*lognxm - Beta.logBeta(nhalf, mhalf));
98     }
99 
100     /**
101      * For this distribution, X, this method returns P(X < 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      */
cumulativeProbability(double x)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
inverseCumulativeProbability(final double p)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
getDomainLowerBound(double p)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
getDomainUpperBound(double p)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
getInitialDomain(double p)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
setNumeratorDegreesOfFreedom(double degreesOfFreedom)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      */
setNumeratorDegreesOfFreedomInternal(double degreesOfFreedom)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      */
getNumeratorDegreesOfFreedom()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
setDenominatorDegreesOfFreedom(double degreesOfFreedom)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      */
setDenominatorDegreesOfFreedomInternal(double degreesOfFreedom)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      */
getDenominatorDegreesOfFreedom()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
getSolverAbsoluteAccuracy()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      */
getSupportLowerBound()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      */
getSupportUpperBound()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      */
getNumericalMean()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      */
getNumericalVariance()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 }
361