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