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 < 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    }