← Back to branch summary

~ubuntu-branches/ubuntu/quantal/commons-math/quantal

« back to all changes in this revision

Viewing changes to src/main/java/org/apache/commons/math/distribution/FDistributionImpl.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2010-04-05 23:33:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20100405233302-gpqlceked76nw28a
Tags: 2.1-1
* New upstream release.
* Bump Standards-Version to 3.8.4: no changes needed
* Bump debhelper to >= 7
* Switch to 3.0 (quilt) source format:
  - Remove B-D on quilt
  - Add d/source/format
  - Remove d/README.source

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/main/java/org/apache/commons/math/distribution/FDistributionImpl.java
26
26
 * Default implementation of
27
27
 * {@link org.apache.commons.math.distribution.FDistribution}.
28
28
 *
29
 
 * @version $Revision: 772119 $ $Date: 2009-05-06 05:43:28 -0400 (Wed, 06 May 2009) $
 
29
 * @version $Revision: 925897 $ $Date: 2010-03-21 17:06:46 -0400 (Sun, 21 Mar 2010) $
30
30
 */
31
31
public class FDistributionImpl
32
32
    extends AbstractContinuousDistribution
33
33
    implements FDistribution, Serializable  {
34
34
 
 
35
    /**
 
36
     * Default inverse cumulative probability accuracy
 
37
     * @since 2.1
 
38
     */
 
39
    public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
 
40
 
 
41
    /** Message for non positive degrees of freddom. */
 
42
    private static final String NON_POSITIVE_DEGREES_OF_FREEDOM_MESSAGE =
 
43
        "degrees of freedom must be positive ({0})";
 
44
 
35
45
    /** Serializable version identifier */
36
46
    private static final long serialVersionUID = -8516354193418641566L;
37
47
 
40
50
 
41
51
    /** The numerator degrees of freedom*/
42
52
    private double denominatorDegreesOfFreedom;
43
 
    
 
53
 
 
54
    /** Inverse cumulative probability accuracy */
 
55
    private final double solverAbsoluteAccuracy;
 
56
 
44
57
    /**
45
58
     * Create a F distribution using the given degrees of freedom.
46
59
     * @param numeratorDegreesOfFreedom the numerator degrees of freedom.
47
60
     * @param denominatorDegreesOfFreedom the denominator degrees of freedom.
48
61
     */
49
62
    public FDistributionImpl(double numeratorDegreesOfFreedom,
50
 
            double denominatorDegreesOfFreedom) {
 
63
                             double denominatorDegreesOfFreedom) {
 
64
        this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
 
65
    }
 
66
 
 
67
    /**
 
68
     * Create a F distribution using the given degrees of freedom and inverse cumulative probability accuracy.
 
69
     * @param numeratorDegreesOfFreedom the numerator degrees of freedom.
 
70
     * @param denominatorDegreesOfFreedom the denominator degrees of freedom.
 
71
     * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates
 
72
     * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY})
 
73
     * @since 2.1
 
74
     */
 
75
    public FDistributionImpl(double numeratorDegreesOfFreedom, double denominatorDegreesOfFreedom,
 
76
            double inverseCumAccuracy) {
51
77
        super();
52
 
        setNumeratorDegreesOfFreedom(numeratorDegreesOfFreedom);
53
 
        setDenominatorDegreesOfFreedom(denominatorDegreesOfFreedom);
54
 
    }
55
 
    
 
78
        setNumeratorDegreesOfFreedomInternal(numeratorDegreesOfFreedom);
 
79
        setDenominatorDegreesOfFreedomInternal(denominatorDegreesOfFreedom);
 
80
        solverAbsoluteAccuracy = inverseCumAccuracy;
 
81
    }
 
82
 
 
83
    /**
 
84
     * Returns the probability density for a particular point.
 
85
     *
 
86
     * @param x The point at which the density should be computed.
 
87
     * @return The pdf at point x.
 
88
     * @since 2.1
 
89
     */
 
90
    @Override
 
91
    public double density(double x) {
 
92
        final double nhalf = numeratorDegreesOfFreedom / 2;
 
93
        final double mhalf = denominatorDegreesOfFreedom / 2;
 
94
        final double logx = Math.log(x);
 
95
        final double logn = Math.log(numeratorDegreesOfFreedom);
 
96
        final double logm = Math.log(denominatorDegreesOfFreedom);
 
97
        final double lognxm = Math.log(numeratorDegreesOfFreedom * x + denominatorDegreesOfFreedom);
 
98
        return Math.exp(nhalf*logn + nhalf*logx - logx + mhalf*logm - nhalf*lognxm -
 
99
               mhalf*lognxm - Beta.logBeta(nhalf, mhalf));
 
100
    }
 
101
 
56
102
    /**
57
103
     * For this distribution, X, this method returns P(X < x).
58
 
     * 
 
104
     *
59
105
     * The implementation of this method is based on:
60
106
     * <ul>
61
107
     * <li>
62
108
     * <a href="http://mathworld.wolfram.com/F-Distribution.html">
63
109
     * F-Distribution</a>, equation (4).</li>
64
110
     * </ul>
65
 
     * 
 
111
     *
66
112
     * @param x the value at which the CDF is evaluated.
67
 
     * @return CDF for this distribution. 
 
113
     * @return CDF for this distribution.
68
114
     * @throws MathException if the cumulative probability can not be
69
115
     *            computed due to convergence or other numerical errors.
70
116
     */
73
119
        if (x <= 0.0) {
74
120
            ret = 0.0;
75
121
        } else {
76
 
            double n = getNumeratorDegreesOfFreedom();
77
 
            double m = getDenominatorDegreesOfFreedom();
78
 
            
 
122
            double n = numeratorDegreesOfFreedom;
 
123
            double m = denominatorDegreesOfFreedom;
 
124
 
79
125
            ret = Beta.regularizedBeta((n * x) / (m + n * x),
80
126
                0.5 * n,
81
127
                0.5 * m);
82
128
        }
83
129
        return ret;
84
130
    }
85
 
    
 
131
 
86
132
    /**
87
133
     * For this distribution, X, this method returns the critical point x, such
88
134
     * that P(X &lt; x) = <code>p</code>.
97
143
     *         probability.
98
144
     */
99
145
    @Override
100
 
    public double inverseCumulativeProbability(final double p) 
 
146
    public double inverseCumulativeProbability(final double p)
101
147
        throws MathException {
102
148
        if (p == 0) {
103
149
            return 0d;
107
153
        }
108
154
        return super.inverseCumulativeProbability(p);
109
155
    }
110
 
        
 
156
 
111
157
    /**
112
158
     * Access the domain value lower bound, based on <code>p</code>, used to
113
159
     * bracket a CDF root.  This method is used by
114
160
     * {@link #inverseCumulativeProbability(double)} to find critical values.
115
 
     * 
 
161
     *
116
162
     * @param p the desired probability for the critical value
117
163
     * @return domain value lower bound, i.e.
118
 
     *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code> 
 
164
     *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>
119
165
     */
120
166
    @Override
121
167
    protected double getDomainLowerBound(double p) {
126
172
     * Access the domain value upper bound, based on <code>p</code>, used to
127
173
     * bracket a CDF root.  This method is used by
128
174
     * {@link #inverseCumulativeProbability(double)} to find critical values.
129
 
     * 
 
175
     *
130
176
     * @param p the desired probability for the critical value
131
177
     * @return domain value upper bound, i.e.
132
 
     *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code> 
 
178
     *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code>
133
179
     */
134
180
    @Override
135
181
    protected double getDomainUpperBound(double p) {
140
186
     * Access the initial domain value, based on <code>p</code>, used to
141
187
     * bracket a CDF root.  This method is used by
142
188
     * {@link #inverseCumulativeProbability(double)} to find critical values.
143
 
     * 
 
189
     *
144
190
     * @param p the desired probability for the critical value
145
191
     * @return initial domain value
146
192
     */
147
193
    @Override
148
194
    protected double getInitialDomain(double p) {
149
195
        double ret = 1.0;
150
 
        double d = getDenominatorDegreesOfFreedom();
 
196
        double d = denominatorDegreesOfFreedom;
151
197
        if (d > 2.0) {
152
198
            // use mean
153
199
            ret = d / (d - 2.0);
154
200
        }
155
201
        return ret;
156
202
    }
157
 
    
 
203
 
158
204
    /**
159
205
     * Modify the numerator degrees of freedom.
160
206
     * @param degreesOfFreedom the new numerator degrees of freedom.
161
207
     * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
162
208
     *         positive.
 
209
     * @deprecated as of 2.1 (class will become immutable in 3.0)
163
210
     */
 
211
    @Deprecated
164
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) {
165
223
        if (degreesOfFreedom <= 0.0) {
166
224
            throw MathRuntimeException.createIllegalArgumentException(
167
 
                  "degrees of freedom must be positive ({0})",
168
 
                  degreesOfFreedom);
 
225
                  NON_POSITIVE_DEGREES_OF_FREEDOM_MESSAGE, degreesOfFreedom);
169
226
        }
170
227
        this.numeratorDegreesOfFreedom = degreesOfFreedom;
171
228
    }
172
 
    
 
229
 
173
230
    /**
174
231
     * Access the numerator degrees of freedom.
175
232
     * @return the numerator degrees of freedom.
177
234
    public double getNumeratorDegreesOfFreedom() {
178
235
        return numeratorDegreesOfFreedom;
179
236
    }
180
 
    
 
237
 
181
238
    /**
182
239
     * Modify the denominator degrees of freedom.
183
240
     * @param degreesOfFreedom the new denominator degrees of freedom.
184
241
     * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
185
242
     *         positive.
 
243
     * @deprecated as of 2.1 (class will become immutable in 3.0)
186
244
     */
 
245
    @Deprecated
187
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) {
188
257
        if (degreesOfFreedom <= 0.0) {
189
258
            throw MathRuntimeException.createIllegalArgumentException(
190
 
                  "degrees of freedom must be positive ({0})",
191
 
                  degreesOfFreedom);
 
259
                  NON_POSITIVE_DEGREES_OF_FREEDOM_MESSAGE, degreesOfFreedom);
192
260
        }
193
261
        this.denominatorDegreesOfFreedom = degreesOfFreedom;
194
262
    }
195
 
    
 
263
 
196
264
    /**
197
265
     * Access the denominator degrees of freedom.
198
266
     * @return the denominator degrees of freedom.
200
268
    public double getDenominatorDegreesOfFreedom() {
201
269
        return denominatorDegreesOfFreedom;
202
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
    }
203
283
}