1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22
23 package org.hipparchus.distribution.continuous;
24
25 import org.hipparchus.exception.LocalizedCoreFormats;
26 import org.hipparchus.exception.MathIllegalArgumentException;
27 import org.hipparchus.special.Beta;
28 import org.hipparchus.util.FastMath;
29
30 /**
31 * Implementation of the F-distribution.
32 *
33 * @see <a href="http://en.wikipedia.org/wiki/F-distribution">F-distribution (Wikipedia)</a>
34 * @see <a href="http://mathworld.wolfram.com/F-Distribution.html">F-distribution (MathWorld)</a>
35 */
36 public class FDistribution extends AbstractRealDistribution {
37 /** Serializable version identifier. */
38 private static final long serialVersionUID = 20160320L;
39 /** The numerator degrees of freedom. */
40 private final double numeratorDegreesOfFreedom;
41 /** The numerator degrees of freedom. */
42 private final double denominatorDegreesOfFreedom;
43 /** Cached numerical variance */
44 private final double numericalVariance;
45
46 /**
47 * Creates an F distribution using the given degrees of freedom.
48 *
49 * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
50 * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
51 * @throws MathIllegalArgumentException if
52 * {@code numeratorDegreesOfFreedom <= 0} or
53 * {@code denominatorDegreesOfFreedom <= 0}.
54 */
55 public FDistribution(double numeratorDegreesOfFreedom,
56 double denominatorDegreesOfFreedom)
57 throws MathIllegalArgumentException {
58 this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom,
59 DEFAULT_SOLVER_ABSOLUTE_ACCURACY);
60 }
61
62
63 /**
64 * Creates an F distribution.
65 *
66 * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
67 * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
68 * @param inverseCumAccuracy the maximum absolute error in inverse
69 * cumulative probability estimates.
70 * @throws MathIllegalArgumentException if {@code numeratorDegreesOfFreedom <= 0} or
71 * {@code denominatorDegreesOfFreedom <= 0}.
72 */
73 public FDistribution(double numeratorDegreesOfFreedom,
74 double denominatorDegreesOfFreedom,
75 double inverseCumAccuracy)
76 throws MathIllegalArgumentException {
77 super(inverseCumAccuracy);
78
79 if (numeratorDegreesOfFreedom <= 0) {
80 throw new MathIllegalArgumentException(LocalizedCoreFormats.DEGREES_OF_FREEDOM,
81 numeratorDegreesOfFreedom);
82 }
83 if (denominatorDegreesOfFreedom <= 0) {
84 throw new MathIllegalArgumentException(LocalizedCoreFormats.DEGREES_OF_FREEDOM,
85 denominatorDegreesOfFreedom);
86 }
87 this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom;
88 this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom;
89 this.numericalVariance = calculateNumericalVariance();
90 }
91
92 /**
93 * {@inheritDoc}
94 */
95 @Override
96 public double density(double x) {
97 return FastMath.exp(logDensity(x));
98 }
99
100 /** {@inheritDoc} **/
101 @Override
102 public double logDensity(double x) {
103 final double nhalf = numeratorDegreesOfFreedom / 2;
104 final double mhalf = denominatorDegreesOfFreedom / 2;
105 final double logx = FastMath.log(x);
106 final double logn = FastMath.log(numeratorDegreesOfFreedom);
107 final double logm = FastMath.log(denominatorDegreesOfFreedom);
108 final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x +
109 denominatorDegreesOfFreedom);
110 return nhalf * logn + nhalf * logx - logx +
111 mhalf * logm - nhalf * lognxm - mhalf * lognxm -
112 Beta.logBeta(nhalf, mhalf);
113 }
114
115 /**
116 * {@inheritDoc}
117 *
118 * The implementation of this method is based on
119 * <ul>
120 * <li>
121 * <a href="http://mathworld.wolfram.com/F-Distribution.html">
122 * F-Distribution</a>, equation (4).
123 * </li>
124 * </ul>
125 */
126 @Override
127 public double cumulativeProbability(double x) {
128 double ret;
129 if (x <= 0) {
130 ret = 0;
131 } else {
132 double n = numeratorDegreesOfFreedom;
133 double m = denominatorDegreesOfFreedom;
134
135 ret = Beta.regularizedBeta((n * x) / (m + n * x),
136 0.5 * n,
137 0.5 * m);
138 }
139 return ret;
140 }
141
142 /**
143 * Access the numerator degrees of freedom.
144 *
145 * @return the numerator degrees of freedom.
146 */
147 public double getNumeratorDegreesOfFreedom() {
148 return numeratorDegreesOfFreedom;
149 }
150
151 /**
152 * Access the denominator degrees of freedom.
153 *
154 * @return the denominator degrees of freedom.
155 */
156 public double getDenominatorDegreesOfFreedom() {
157 return denominatorDegreesOfFreedom;
158 }
159
160 /**
161 * {@inheritDoc}
162 *
163 * For denominator degrees of freedom parameter {@code b}, the mean is
164 * <ul>
165 * <li>if {@code b > 2} then {@code b / (b - 2)},</li>
166 * <li>else undefined ({@code Double.NaN}).
167 * </ul>
168 */
169 @Override
170 public double getNumericalMean() {
171 final double denominatorDF = getDenominatorDegreesOfFreedom();
172
173 if (denominatorDF > 2) {
174 return denominatorDF / (denominatorDF - 2);
175 }
176
177 return Double.NaN;
178 }
179
180 /**
181 * {@inheritDoc}
182 *
183 * For numerator degrees of freedom parameter {@code a} and denominator
184 * degrees of freedom parameter {@code b}, the variance is
185 * <ul>
186 * <li>
187 * if {@code b > 4} then
188 * {@code [2 * b^2 * (a + b - 2)] / [a * (b - 2)^2 * (b - 4)]},
189 * </li>
190 * <li>else undefined ({@code Double.NaN}).
191 * </ul>
192 */
193 @Override
194 public double getNumericalVariance() {
195 return numericalVariance;
196 }
197
198 /**
199 * Calculates the numerical variance.
200 *
201 * @return the variance of this distribution
202 */
203 private double calculateNumericalVariance() {
204 final double denominatorDF = getDenominatorDegreesOfFreedom();
205
206 if (denominatorDF > 4) {
207 final double numeratorDF = getNumeratorDegreesOfFreedom();
208 final double denomDFMinusTwo = denominatorDF - 2;
209
210 return ( 2 * (denominatorDF * denominatorDF) * (numeratorDF + denominatorDF - 2) ) /
211 ( (numeratorDF * (denomDFMinusTwo * denomDFMinusTwo) * (denominatorDF - 4)) );
212 }
213
214 return Double.NaN;
215 }
216
217 /**
218 * {@inheritDoc}
219 *
220 * The lower bound of the support is always 0 no matter the parameters.
221 *
222 * @return lower bound of the support (always 0)
223 */
224 @Override
225 public double getSupportLowerBound() {
226 return 0;
227 }
228
229 /**
230 * {@inheritDoc}
231 *
232 * The upper bound of the support is always positive infinity
233 * no matter the parameters.
234 *
235 * @return upper bound of the support (always Double.POSITIVE_INFINITY)
236 */
237 @Override
238 public double getSupportUpperBound() {
239 return Double.POSITIVE_INFINITY;
240 }
241
242 /**
243 * {@inheritDoc}
244 *
245 * The support of this distribution is connected.
246 *
247 * @return {@code true}
248 */
249 @Override
250 public boolean isSupportConnected() {
251 return true;
252 }
253 }