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.special.Beta;
024    
025    /**
026     * Default implementation of
027     * {@link org.apache.commons.math.distribution.FDistribution}.
028     *
029     * @version $Revision: 925897 $ $Date: 2010-03-21 17:06:46 -0400 (Sun, 21 Mar 2010) $
030     */
031    public class FDistributionImpl
032        extends AbstractContinuousDistribution
033        implements FDistribution, Serializable  {
034    
035        /**
036         * Default inverse cumulative probability accuracy
037         * @since 2.1
038         */
039        public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
040    
041        /** Message for non positive degrees of freddom. */
042        private static final String NON_POSITIVE_DEGREES_OF_FREEDOM_MESSAGE =
043            "degrees of freedom must be positive ({0})";
044    
045        /** Serializable version identifier */
046        private static final long serialVersionUID = -8516354193418641566L;
047    
048        /** The numerator degrees of freedom*/
049        private double numeratorDegreesOfFreedom;
050    
051        /** The numerator degrees of freedom*/
052        private double denominatorDegreesOfFreedom;
053    
054        /** Inverse cumulative probability accuracy */
055        private final double solverAbsoluteAccuracy;
056    
057        /**
058         * Create a F distribution using the given degrees of freedom.
059         * @param numeratorDegreesOfFreedom the numerator degrees of freedom.
060         * @param denominatorDegreesOfFreedom the denominator degrees of freedom.
061         */
062        public FDistributionImpl(double numeratorDegreesOfFreedom,
063                                 double denominatorDegreesOfFreedom) {
064            this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
065        }
066    
067        /**
068         * Create a F distribution using the given degrees of freedom and inverse cumulative probability accuracy.
069         * @param numeratorDegreesOfFreedom the numerator degrees of freedom.
070         * @param denominatorDegreesOfFreedom the denominator degrees of freedom.
071         * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates
072         * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY})
073         * @since 2.1
074         */
075        public FDistributionImpl(double numeratorDegreesOfFreedom, double denominatorDegreesOfFreedom,
076                double inverseCumAccuracy) {
077            super();
078            setNumeratorDegreesOfFreedomInternal(numeratorDegreesOfFreedom);
079            setDenominatorDegreesOfFreedomInternal(denominatorDegreesOfFreedom);
080            solverAbsoluteAccuracy = inverseCumAccuracy;
081        }
082    
083        /**
084         * Returns the probability density for a particular point.
085         *
086         * @param x The point at which the density should be computed.
087         * @return The pdf at point x.
088         * @since 2.1
089         */
090        @Override
091        public double density(double x) {
092            final double nhalf = numeratorDegreesOfFreedom / 2;
093            final double mhalf = denominatorDegreesOfFreedom / 2;
094            final double logx = Math.log(x);
095            final double logn = Math.log(numeratorDegreesOfFreedom);
096            final double logm = Math.log(denominatorDegreesOfFreedom);
097            final double lognxm = Math.log(numeratorDegreesOfFreedom * x + denominatorDegreesOfFreedom);
098            return Math.exp(nhalf*logn + nhalf*logx - logx + mhalf*logm - nhalf*lognxm -
099                   mhalf*lognxm - Beta.logBeta(nhalf, mhalf));
100        }
101    
102        /**
103         * For this distribution, X, this method returns P(X < x).
104         *
105         * The implementation of this method is based on:
106         * <ul>
107         * <li>
108         * <a href="http://mathworld.wolfram.com/F-Distribution.html">
109         * F-Distribution</a>, equation (4).</li>
110         * </ul>
111         *
112         * @param x the value at which the CDF is evaluated.
113         * @return CDF for this distribution.
114         * @throws MathException if the cumulative probability can not be
115         *            computed due to convergence or other numerical errors.
116         */
117        public double cumulativeProbability(double x) throws MathException {
118            double ret;
119            if (x <= 0.0) {
120                ret = 0.0;
121            } else {
122                double n = numeratorDegreesOfFreedom;
123                double m = denominatorDegreesOfFreedom;
124    
125                ret = Beta.regularizedBeta((n * x) / (m + n * x),
126                    0.5 * n,
127                    0.5 * m);
128            }
129            return ret;
130        }
131    
132        /**
133         * For this distribution, X, this method returns the critical point x, such
134         * that P(X &lt; x) = <code>p</code>.
135         * <p>
136         * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.</p>
137         *
138         * @param p the desired probability
139         * @return x, such that P(X &lt; x) = <code>p</code>
140         * @throws MathException if the inverse cumulative probability can not be
141         *         computed due to convergence or other numerical errors.
142         * @throws IllegalArgumentException if <code>p</code> is not a valid
143         *         probability.
144         */
145        @Override
146        public double inverseCumulativeProbability(final double p)
147            throws MathException {
148            if (p == 0) {
149                return 0d;
150            }
151            if (p == 1) {
152                return Double.POSITIVE_INFINITY;
153            }
154            return super.inverseCumulativeProbability(p);
155        }
156    
157        /**
158         * Access the domain value lower bound, based on <code>p</code>, used to
159         * bracket a CDF root.  This method is used by
160         * {@link #inverseCumulativeProbability(double)} to find critical values.
161         *
162         * @param p the desired probability for the critical value
163         * @return domain value lower bound, i.e.
164         *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>
165         */
166        @Override
167        protected double getDomainLowerBound(double p) {
168            return 0.0;
169        }
170    
171        /**
172         * Access the domain value upper bound, based on <code>p</code>, used to
173         * bracket a CDF root.  This method is used by
174         * {@link #inverseCumulativeProbability(double)} to find critical values.
175         *
176         * @param p the desired probability for the critical value
177         * @return domain value upper bound, i.e.
178         *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code>
179         */
180        @Override
181        protected double getDomainUpperBound(double p) {
182            return Double.MAX_VALUE;
183        }
184    
185        /**
186         * Access the initial domain value, based on <code>p</code>, used to
187         * bracket a CDF root.  This method is used by
188         * {@link #inverseCumulativeProbability(double)} to find critical values.
189         *
190         * @param p the desired probability for the critical value
191         * @return initial domain value
192         */
193        @Override
194        protected double getInitialDomain(double p) {
195            double ret = 1.0;
196            double d = denominatorDegreesOfFreedom;
197            if (d > 2.0) {
198                // use mean
199                ret = d / (d - 2.0);
200            }
201            return ret;
202        }
203    
204        /**
205         * Modify the numerator degrees of freedom.
206         * @param degreesOfFreedom the new numerator degrees of freedom.
207         * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
208         *         positive.
209         * @deprecated as of 2.1 (class will become immutable in 3.0)
210         */
211        @Deprecated
212        public void setNumeratorDegreesOfFreedom(double degreesOfFreedom) {
213            setNumeratorDegreesOfFreedomInternal(degreesOfFreedom);
214        }
215    
216        /**
217         * Modify the numerator degrees of freedom.
218         * @param degreesOfFreedom the new numerator degrees of freedom.
219         * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
220         *         positive.
221         */
222        private void setNumeratorDegreesOfFreedomInternal(double degreesOfFreedom) {
223            if (degreesOfFreedom <= 0.0) {
224                throw MathRuntimeException.createIllegalArgumentException(
225                      NON_POSITIVE_DEGREES_OF_FREEDOM_MESSAGE, degreesOfFreedom);
226            }
227            this.numeratorDegreesOfFreedom = degreesOfFreedom;
228        }
229    
230        /**
231         * Access the numerator degrees of freedom.
232         * @return the numerator degrees of freedom.
233         */
234        public double getNumeratorDegreesOfFreedom() {
235            return numeratorDegreesOfFreedom;
236        }
237    
238        /**
239         * Modify the denominator degrees of freedom.
240         * @param degreesOfFreedom the new denominator degrees of freedom.
241         * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
242         *         positive.
243         * @deprecated as of 2.1 (class will become immutable in 3.0)
244         */
245        @Deprecated
246        public void setDenominatorDegreesOfFreedom(double degreesOfFreedom) {
247            setDenominatorDegreesOfFreedomInternal(degreesOfFreedom);
248        }
249    
250        /**
251         * Modify the denominator degrees of freedom.
252         * @param degreesOfFreedom the new denominator degrees of freedom.
253         * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
254         *         positive.
255         */
256        private void setDenominatorDegreesOfFreedomInternal(double degreesOfFreedom) {
257            if (degreesOfFreedom <= 0.0) {
258                throw MathRuntimeException.createIllegalArgumentException(
259                      NON_POSITIVE_DEGREES_OF_FREEDOM_MESSAGE, degreesOfFreedom);
260            }
261            this.denominatorDegreesOfFreedom = degreesOfFreedom;
262        }
263    
264        /**
265         * Access the denominator degrees of freedom.
266         * @return the denominator degrees of freedom.
267         */
268        public double getDenominatorDegreesOfFreedom() {
269            return denominatorDegreesOfFreedom;
270        }
271    
272        /**
273         * Return the absolute accuracy setting of the solver used to estimate
274         * inverse cumulative probabilities.
275         *
276         * @return the solver absolute accuracy
277         * @since 2.1
278         */
279        @Override
280        protected double getSolverAbsoluteAccuracy() {
281            return solverAbsoluteAccuracy;
282        }
283    }