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
018 package org.apache.commons.math3.distribution;
019
020 import org.apache.commons.math3.exception.NotStrictlyPositiveException;
021 import org.apache.commons.math3.exception.util.LocalizedFormats;
022 import org.apache.commons.math3.special.Beta;
023 import org.apache.commons.math3.util.FastMath;
024 import org.apache.commons.math3.random.RandomGenerator;
025 import org.apache.commons.math3.random.Well19937c;
026
027 /**
028 * Implementation of the F-distribution.
029 *
030 * @see <a href="http://en.wikipedia.org/wiki/F-distribution">F-distribution (Wikipedia)</a>
031 * @see <a href="http://mathworld.wolfram.com/F-Distribution.html">F-distribution (MathWorld)</a>
032 * @version $Id: FDistribution.java 1416643 2012-12-03 19:37:14Z tn $
033 */
034 public class FDistribution extends AbstractRealDistribution {
035 /**
036 * Default inverse cumulative probability accuracy.
037 * @since 2.1
038 */
039 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
040 /** Serializable version identifier. */
041 private static final long serialVersionUID = -8516354193418641566L;
042 /** The numerator degrees of freedom. */
043 private final double numeratorDegreesOfFreedom;
044 /** The numerator degrees of freedom. */
045 private final double denominatorDegreesOfFreedom;
046 /** Inverse cumulative probability accuracy. */
047 private final double solverAbsoluteAccuracy;
048 /** Cached numerical variance */
049 private double numericalVariance = Double.NaN;
050 /** Whether or not the numerical variance has been calculated */
051 private boolean numericalVarianceIsCalculated = false;
052
053 /**
054 * Creates an F distribution using the given degrees of freedom.
055 *
056 * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
057 * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
058 * @throws NotStrictlyPositiveException if
059 * {@code numeratorDegreesOfFreedom <= 0} or
060 * {@code denominatorDegreesOfFreedom <= 0}.
061 */
062 public FDistribution(double numeratorDegreesOfFreedom,
063 double denominatorDegreesOfFreedom)
064 throws NotStrictlyPositiveException {
065 this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom,
066 DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
067 }
068
069 /**
070 * Creates an F distribution using the given degrees of freedom
071 * and inverse cumulative probability accuracy.
072 *
073 * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
074 * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
075 * @param inverseCumAccuracy the maximum absolute error in inverse
076 * cumulative probability estimates.
077 * @throws NotStrictlyPositiveException if
078 * {@code numeratorDegreesOfFreedom <= 0} or
079 * {@code denominatorDegreesOfFreedom <= 0}.
080 * @since 2.1
081 */
082 public FDistribution(double numeratorDegreesOfFreedom,
083 double denominatorDegreesOfFreedom,
084 double inverseCumAccuracy)
085 throws NotStrictlyPositiveException {
086 this(new Well19937c(), numeratorDegreesOfFreedom,
087 denominatorDegreesOfFreedom, inverseCumAccuracy);
088 }
089
090 /**
091 * Creates an F distribution.
092 *
093 * @param rng Random number generator.
094 * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
095 * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
096 * @param inverseCumAccuracy the maximum absolute error in inverse
097 * cumulative probability estimates.
098 * @throws NotStrictlyPositiveException if
099 * {@code numeratorDegreesOfFreedom <= 0} or
100 * {@code denominatorDegreesOfFreedom <= 0}.
101 * @since 3.1
102 */
103 public FDistribution(RandomGenerator rng,
104 double numeratorDegreesOfFreedom,
105 double denominatorDegreesOfFreedom,
106 double inverseCumAccuracy)
107 throws NotStrictlyPositiveException {
108 super(rng);
109
110 if (numeratorDegreesOfFreedom <= 0) {
111 throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
112 numeratorDegreesOfFreedom);
113 }
114 if (denominatorDegreesOfFreedom <= 0) {
115 throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
116 denominatorDegreesOfFreedom);
117 }
118 this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom;
119 this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom;
120 solverAbsoluteAccuracy = inverseCumAccuracy;
121 }
122
123 /**
124 * {@inheritDoc}
125 *
126 * @since 2.1
127 */
128 public double density(double x) {
129 final double nhalf = numeratorDegreesOfFreedom / 2;
130 final double mhalf = denominatorDegreesOfFreedom / 2;
131 final double logx = FastMath.log(x);
132 final double logn = FastMath.log(numeratorDegreesOfFreedom);
133 final double logm = FastMath.log(denominatorDegreesOfFreedom);
134 final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x +
135 denominatorDegreesOfFreedom);
136 return FastMath.exp(nhalf * logn + nhalf * logx - logx +
137 mhalf * logm - nhalf * lognxm - mhalf * lognxm -
138 Beta.logBeta(nhalf, mhalf));
139 }
140
141 /**
142 * {@inheritDoc}
143 *
144 * The implementation of this method is based on
145 * <ul>
146 * <li>
147 * <a href="http://mathworld.wolfram.com/F-Distribution.html">
148 * F-Distribution</a>, equation (4).
149 * </li>
150 * </ul>
151 */
152 public double cumulativeProbability(double x) {
153 double ret;
154 if (x <= 0) {
155 ret = 0;
156 } else {
157 double n = numeratorDegreesOfFreedom;
158 double m = denominatorDegreesOfFreedom;
159
160 ret = Beta.regularizedBeta((n * x) / (m + n * x),
161 0.5 * n,
162 0.5 * m);
163 }
164 return ret;
165 }
166
167 /**
168 * Access the numerator degrees of freedom.
169 *
170 * @return the numerator degrees of freedom.
171 */
172 public double getNumeratorDegreesOfFreedom() {
173 return numeratorDegreesOfFreedom;
174 }
175
176 /**
177 * Access the denominator degrees of freedom.
178 *
179 * @return the denominator degrees of freedom.
180 */
181 public double getDenominatorDegreesOfFreedom() {
182 return denominatorDegreesOfFreedom;
183 }
184
185 /** {@inheritDoc} */
186 @Override
187 protected double getSolverAbsoluteAccuracy() {
188 return solverAbsoluteAccuracy;
189 }
190
191 /**
192 * {@inheritDoc}
193 *
194 * For denominator degrees of freedom parameter {@code b}, the mean is
195 * <ul>
196 * <li>if {@code b > 2} then {@code b / (b - 2)},</li>
197 * <li>else undefined ({@code Double.NaN}).
198 * </ul>
199 */
200 public double getNumericalMean() {
201 final double denominatorDF = getDenominatorDegreesOfFreedom();
202
203 if (denominatorDF > 2) {
204 return denominatorDF / (denominatorDF - 2);
205 }
206
207 return Double.NaN;
208 }
209
210 /**
211 * {@inheritDoc}
212 *
213 * For numerator degrees of freedom parameter {@code a} and denominator
214 * degrees of freedom parameter {@code b}, the variance is
215 * <ul>
216 * <li>
217 * if {@code b > 4} then
218 * {@code [2 * b^2 * (a + b - 2)] / [a * (b - 2)^2 * (b - 4)]},
219 * </li>
220 * <li>else undefined ({@code Double.NaN}).
221 * </ul>
222 */
223 public double getNumericalVariance() {
224 if (!numericalVarianceIsCalculated) {
225 numericalVariance = calculateNumericalVariance();
226 numericalVarianceIsCalculated = true;
227 }
228 return numericalVariance;
229 }
230
231 /**
232 * used by {@link #getNumericalVariance()}
233 *
234 * @return the variance of this distribution
235 */
236 protected double calculateNumericalVariance() {
237 final double denominatorDF = getDenominatorDegreesOfFreedom();
238
239 if (denominatorDF > 4) {
240 final double numeratorDF = getNumeratorDegreesOfFreedom();
241 final double denomDFMinusTwo = denominatorDF - 2;
242
243 return ( 2 * (denominatorDF * denominatorDF) * (numeratorDF + denominatorDF - 2) ) /
244 ( (numeratorDF * (denomDFMinusTwo * denomDFMinusTwo) * (denominatorDF - 4)) );
245 }
246
247 return Double.NaN;
248 }
249
250 /**
251 * {@inheritDoc}
252 *
253 * The lower bound of the support is always 0 no matter the parameters.
254 *
255 * @return lower bound of the support (always 0)
256 */
257 public double getSupportLowerBound() {
258 return 0;
259 }
260
261 /**
262 * {@inheritDoc}
263 *
264 * The upper bound of the support is always positive infinity
265 * no matter the parameters.
266 *
267 * @return upper bound of the support (always Double.POSITIVE_INFINITY)
268 */
269 public double getSupportUpperBound() {
270 return Double.POSITIVE_INFINITY;
271 }
272
273 /** {@inheritDoc} */
274 public boolean isSupportLowerBoundInclusive() {
275 return false;
276 }
277
278 /** {@inheritDoc} */
279 public boolean isSupportUpperBoundInclusive() {
280 return false;
281 }
282
283 /**
284 * {@inheritDoc}
285 *
286 * The support of this distribution is connected.
287 *
288 * @return {@code true}
289 */
290 public boolean isSupportConnected() {
291 return true;
292 }
293 }