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 < 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 < 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 < <i>lower bound</i>) < <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 < <i>upper bound</i>) > <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 }