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: 772119 $ $Date: 2009-05-06 05:43:28 -0400 (Wed, 06 May 2009) $
030     */
031    public class FDistributionImpl
032        extends AbstractContinuousDistribution
033        implements FDistribution, Serializable  {
034    
035        /** Serializable version identifier */
036        private static final long serialVersionUID = -8516354193418641566L;
037    
038        /** The numerator degrees of freedom*/
039        private double numeratorDegreesOfFreedom;
040    
041        /** The numerator degrees of freedom*/
042        private double denominatorDegreesOfFreedom;
043        
044        /**
045         * Create a F distribution using the given degrees of freedom.
046         * @param numeratorDegreesOfFreedom the numerator degrees of freedom.
047         * @param denominatorDegreesOfFreedom the denominator degrees of freedom.
048         */
049        public FDistributionImpl(double numeratorDegreesOfFreedom,
050                double denominatorDegreesOfFreedom) {
051            super();
052            setNumeratorDegreesOfFreedom(numeratorDegreesOfFreedom);
053            setDenominatorDegreesOfFreedom(denominatorDegreesOfFreedom);
054        }
055        
056        /**
057         * For this distribution, X, this method returns P(X < x).
058         * 
059         * The implementation of this method is based on:
060         * <ul>
061         * <li>
062         * <a href="http://mathworld.wolfram.com/F-Distribution.html">
063         * F-Distribution</a>, equation (4).</li>
064         * </ul>
065         * 
066         * @param x the value at which the CDF is evaluated.
067         * @return CDF for this distribution. 
068         * @throws MathException if the cumulative probability can not be
069         *            computed due to convergence or other numerical errors.
070         */
071        public double cumulativeProbability(double x) throws MathException {
072            double ret;
073            if (x <= 0.0) {
074                ret = 0.0;
075            } else {
076                double n = getNumeratorDegreesOfFreedom();
077                double m = getDenominatorDegreesOfFreedom();
078                
079                ret = Beta.regularizedBeta((n * x) / (m + n * x),
080                    0.5 * n,
081                    0.5 * m);
082            }
083            return ret;
084        }
085        
086        /**
087         * For this distribution, X, this method returns the critical point x, such
088         * that P(X &lt; x) = <code>p</code>.
089         * <p>
090         * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.</p>
091         *
092         * @param p the desired probability
093         * @return x, such that P(X &lt; x) = <code>p</code>
094         * @throws MathException if the inverse cumulative probability can not be
095         *         computed due to convergence or other numerical errors.
096         * @throws IllegalArgumentException if <code>p</code> is not a valid
097         *         probability.
098         */
099        @Override
100        public double inverseCumulativeProbability(final double p) 
101            throws MathException {
102            if (p == 0) {
103                return 0d;
104            }
105            if (p == 1) {
106                return Double.POSITIVE_INFINITY;
107            }
108            return super.inverseCumulativeProbability(p);
109        }
110            
111        /**
112         * Access the domain value lower bound, based on <code>p</code>, used to
113         * bracket a CDF root.  This method is used by
114         * {@link #inverseCumulativeProbability(double)} to find critical values.
115         * 
116         * @param p the desired probability for the critical value
117         * @return domain value lower bound, i.e.
118         *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code> 
119         */
120        @Override
121        protected double getDomainLowerBound(double p) {
122            return 0.0;
123        }
124    
125        /**
126         * Access the domain value upper bound, based on <code>p</code>, used to
127         * bracket a CDF root.  This method is used by
128         * {@link #inverseCumulativeProbability(double)} to find critical values.
129         * 
130         * @param p the desired probability for the critical value
131         * @return domain value upper bound, i.e.
132         *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code> 
133         */
134        @Override
135        protected double getDomainUpperBound(double p) {
136            return Double.MAX_VALUE;
137        }
138    
139        /**
140         * Access the initial domain value, based on <code>p</code>, used to
141         * bracket a CDF root.  This method is used by
142         * {@link #inverseCumulativeProbability(double)} to find critical values.
143         * 
144         * @param p the desired probability for the critical value
145         * @return initial domain value
146         */
147        @Override
148        protected double getInitialDomain(double p) {
149            double ret = 1.0;
150            double d = getDenominatorDegreesOfFreedom();
151            if (d > 2.0) {
152                // use mean
153                ret = d / (d - 2.0);
154            }
155            return ret;
156        }
157        
158        /**
159         * Modify the numerator degrees of freedom.
160         * @param degreesOfFreedom the new numerator degrees of freedom.
161         * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
162         *         positive.
163         */
164        public void setNumeratorDegreesOfFreedom(double degreesOfFreedom) {
165            if (degreesOfFreedom <= 0.0) {
166                throw MathRuntimeException.createIllegalArgumentException(
167                      "degrees of freedom must be positive ({0})",
168                      degreesOfFreedom);
169            }
170            this.numeratorDegreesOfFreedom = degreesOfFreedom;
171        }
172        
173        /**
174         * Access the numerator degrees of freedom.
175         * @return the numerator degrees of freedom.
176         */
177        public double getNumeratorDegreesOfFreedom() {
178            return numeratorDegreesOfFreedom;
179        }
180        
181        /**
182         * Modify the denominator degrees of freedom.
183         * @param degreesOfFreedom the new denominator degrees of freedom.
184         * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
185         *         positive.
186         */
187        public void setDenominatorDegreesOfFreedom(double degreesOfFreedom) {
188            if (degreesOfFreedom <= 0.0) {
189                throw MathRuntimeException.createIllegalArgumentException(
190                      "degrees of freedom must be positive ({0})",
191                      degreesOfFreedom);
192            }
193            this.denominatorDegreesOfFreedom = degreesOfFreedom;
194        }
195        
196        /**
197         * Access the denominator degrees of freedom.
198         * @return the denominator degrees of freedom.
199         */
200        public double getDenominatorDegreesOfFreedom() {
201            return denominatorDegreesOfFreedom;
202        }
203    }