1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16 package org.apache.commons.math.distribution;
17
18 import java.io.Serializable;
19
20 import org.apache.commons.math.MathException;
21 import org.apache.commons.math.special.Beta;
22
23 /**
24 * Default implementation of
25 * {@link org.apache.commons.math.distribution.FDistribution}.
26 *
27 * @version $Revision: 348519 $ $Date: 2005-11-23 12:12:18 -0700 (Wed, 23 Nov 2005) $
28 */
29 public class FDistributionImpl
30 extends AbstractContinuousDistribution
31 implements FDistribution, Serializable {
32
33 /** Serializable version identifier */
34 private static final long serialVersionUID = -8516354193418641566L;
35
36 /** The numerator degrees of freedom*/
37 private double numeratorDegreesOfFreedom;
38
39 /** The numerator degrees of freedom*/
40 private double denominatorDegreesOfFreedom;
41
42 /**
43 * Create a F distribution using the given degrees of freedom.
44 * @param numeratorDegreesOfFreedom the numerator degrees of freedom.
45 * @param denominatorDegreesOfFreedom the denominator degrees of freedom.
46 */
47 public FDistributionImpl(double numeratorDegreesOfFreedom,
48 double denominatorDegreesOfFreedom) {
49 super();
50 setNumeratorDegreesOfFreedom(numeratorDegreesOfFreedom);
51 setDenominatorDegreesOfFreedom(denominatorDegreesOfFreedom);
52 }
53
54 /**
55 * For this disbution, X, this method returns P(X < x).
56 *
57 * The implementation of this method is based on:
58 * <ul>
59 * <li>
60 * <a href="http://mathworld.wolfram.com/F-Distribution.html">
61 * F-Distribution</a>, equation (4).</li>
62 * </ul>
63 *
64 * @param x the value at which the CDF is evaluated.
65 * @return CDF for this distribution.
66 * @throws MathException if the cumulative probability can not be
67 * computed due to convergence or other numerical errors.
68 */
69 public double cumulativeProbability(double x) throws MathException {
70 double ret;
71 if (x <= 0.0) {
72 ret = 0.0;
73 } else {
74 double n = getNumeratorDegreesOfFreedom();
75 double m = getDenominatorDegreesOfFreedom();
76
77 ret = Beta.regularizedBeta((n * x) / (m + n * x),
78 0.5 * n,
79 0.5 * m);
80 }
81 return ret;
82 }
83
84 /**
85 * For this distribution, X, this method returns the critical point x, such
86 * that P(X < x) = <code>p</code>.
87 * <p>
88 * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.
89 *
90 * @param p the desired probability
91 * @return x, such that P(X < x) = <code>p</code>
92 * @throws MathException if the inverse cumulative probability can not be
93 * computed due to convergence or other numerical errors.
94 * @throws IllegalArgumentException if <code>p</code> is not a valid
95 * probability.
96 */
97 public double inverseCumulativeProbability(final double p)
98 throws MathException {
99 if (p == 0) {
100 return 0d;
101 }
102 if (p == 1) {
103 return Double.POSITIVE_INFINITY;
104 }
105 return super.inverseCumulativeProbability(p);
106 }
107
108 /**
109 * Access the domain value lower bound, based on <code>p</code>, used to
110 * bracket a CDF root. This method is used by
111 * {@link #inverseCumulativeProbability(double)} to find critical values.
112 *
113 * @param p the desired probability for the critical value
114 * @return domain value lower bound, i.e.
115 * P(X < <i>lower bound</i>) < <code>p</code>
116 */
117 protected double getDomainLowerBound(double p) {
118 return 0.0;
119 }
120
121 /**
122 * Access the domain value upper bound, based on <code>p</code>, used to
123 * bracket a CDF root. This method is used by
124 * {@link #inverseCumulativeProbability(double)} to find critical values.
125 *
126 * @param p the desired probability for the critical value
127 * @return domain value upper bound, i.e.
128 * P(X < <i>upper bound</i>) > <code>p</code>
129 */
130 protected double getDomainUpperBound(double p) {
131 return Double.MAX_VALUE;
132 }
133
134 /**
135 * Access the initial domain value, based on <code>p</code>, used to
136 * bracket a CDF root. This method is used by
137 * {@link #inverseCumulativeProbability(double)} to find critical values.
138 *
139 * @param p the desired probability for the critical value
140 * @return initial domain value
141 */
142 protected double getInitialDomain(double p) {
143 return getDenominatorDegreesOfFreedom() /
144 (getDenominatorDegreesOfFreedom() - 2.0);
145 }
146
147 /**
148 * Modify the numerator degrees of freedom.
149 * @param degreesOfFreedom the new numerator degrees of freedom.
150 * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
151 * positive.
152 */
153 public void setNumeratorDegreesOfFreedom(double degreesOfFreedom) {
154 if (degreesOfFreedom <= 0.0) {
155 throw new IllegalArgumentException(
156 "degrees of freedom must be positive.");
157 }
158 this.numeratorDegreesOfFreedom = degreesOfFreedom;
159 }
160
161 /**
162 * Access the numerator degrees of freedom.
163 * @return the numerator degrees of freedom.
164 */
165 public double getNumeratorDegreesOfFreedom() {
166 return numeratorDegreesOfFreedom;
167 }
168
169 /**
170 * Modify the denominator degrees of freedom.
171 * @param degreesOfFreedom the new denominator degrees of freedom.
172 * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
173 * positive.
174 */
175 public void setDenominatorDegreesOfFreedom(double degreesOfFreedom) {
176 if (degreesOfFreedom <= 0.0) {
177 throw new IllegalArgumentException(
178 "degrees of freedom must be positive.");
179 }
180 this.denominatorDegreesOfFreedom = degreesOfFreedom;
181 }
182
183 /**
184 * Access the denominator degrees of freedom.
185 * @return the denominator degrees of freedom.
186 */
187 public double getDenominatorDegreesOfFreedom() {
188 return denominatorDegreesOfFreedom;
189 }
190 }