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 *      http://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 */
17package org.apache.commons.math.distribution;
18
19import java.io.Serializable;
20
21import org.apache.commons.math.MathException;
22import org.apache.commons.math.MathRuntimeException;
23import org.apache.commons.math.exception.util.LocalizedFormats;
24import org.apache.commons.math.special.Beta;
25import org.apache.commons.math.special.Gamma;
26import org.apache.commons.math.util.FastMath;
27
28/**
29 * Default implementation of
30 * {@link org.apache.commons.math.distribution.TDistribution}.
31 *
32 * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $
33 */
34public class TDistributionImpl
35    extends AbstractContinuousDistribution
36    implements TDistribution, Serializable  {
37
38    /**
39     * Default inverse cumulative probability accuracy
40     * @since 2.1
41    */
42    public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
43
44    /** Serializable version identifier */
45    private static final long serialVersionUID = -5852615386664158222L;
46
47    /** The degrees of freedom*/
48    private double degreesOfFreedom;
49
50    /** Inverse cumulative probability accuracy */
51    private final double solverAbsoluteAccuracy;
52
53    /**
54     * Create a t distribution using the given degrees of freedom and the
55     * specified inverse cumulative probability absolute accuracy.
56     *
57     * @param degreesOfFreedom the degrees of freedom.
58     * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates
59     * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY})
60     * @since 2.1
61     */
62    public TDistributionImpl(double degreesOfFreedom, double inverseCumAccuracy) {
63        super();
64        setDegreesOfFreedomInternal(degreesOfFreedom);
65        solverAbsoluteAccuracy = inverseCumAccuracy;
66    }
67
68    /**
69     * Create a t distribution using the given degrees of freedom.
70     * @param degreesOfFreedom the degrees of freedom.
71     */
72    public TDistributionImpl(double degreesOfFreedom) {
73        this(degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
74    }
75
76    /**
77     * Modify the degrees of freedom.
78     * @param degreesOfFreedom the new degrees of freedom.
79     * @deprecated as of 2.1 (class will become immutable in 3.0)
80     */
81    @Deprecated
82    public void setDegreesOfFreedom(double degreesOfFreedom) {
83        setDegreesOfFreedomInternal(degreesOfFreedom);
84    }
85
86    /**
87     * Modify the degrees of freedom.
88     * @param newDegreesOfFreedom the new degrees of freedom.
89     */
90    private void setDegreesOfFreedomInternal(double newDegreesOfFreedom) {
91        if (newDegreesOfFreedom <= 0.0) {
92            throw MathRuntimeException.createIllegalArgumentException(
93                  LocalizedFormats.NOT_POSITIVE_DEGREES_OF_FREEDOM,
94                  newDegreesOfFreedom);
95        }
96        this.degreesOfFreedom = newDegreesOfFreedom;
97    }
98
99    /**
100     * Access the degrees of freedom.
101     * @return the degrees of freedom.
102     */
103    public double getDegreesOfFreedom() {
104        return degreesOfFreedom;
105    }
106
107    /**
108     * Returns the probability density for a particular point.
109     *
110     * @param x The point at which the density should be computed.
111     * @return The pdf at point x.
112     * @since 2.1
113     */
114    @Override
115    public double density(double x) {
116        final double n = degreesOfFreedom;
117        final double nPlus1Over2 = (n + 1) / 2;
118        return FastMath.exp(Gamma.logGamma(nPlus1Over2) - 0.5 * (FastMath.log(FastMath.PI) + FastMath.log(n)) -
119                Gamma.logGamma(n/2) - nPlus1Over2 * FastMath.log(1 + x * x /n));
120    }
121
122    /**
123     * For this distribution, X, this method returns P(X &lt; <code>x</code>).
124     * @param x the value at which the CDF is evaluated.
125     * @return CDF evaluated at <code>x</code>.
126     * @throws MathException if the cumulative probability can not be
127     *            computed due to convergence or other numerical errors.
128     */
129    public double cumulativeProbability(double x) throws MathException{
130        double ret;
131        if (x == 0.0) {
132            ret = 0.5;
133        } else {
134            double t =
135                Beta.regularizedBeta(
136                    degreesOfFreedom / (degreesOfFreedom + (x * x)),
137                    0.5 * degreesOfFreedom,
138                    0.5);
139            if (x < 0.0) {
140                ret = 0.5 * t;
141            } else {
142                ret = 1.0 - 0.5 * t;
143            }
144        }
145
146        return ret;
147    }
148
149    /**
150     * For this distribution, X, this method returns the critical point x, such
151     * that P(X &lt; x) = <code>p</code>.
152     * <p>
153     * Returns <code>Double.NEGATIVE_INFINITY</code> for p=0 and
154     * <code>Double.POSITIVE_INFINITY</code> for p=1.</p>
155     *
156     * @param p the desired probability
157     * @return x, such that P(X &lt; x) = <code>p</code>
158     * @throws MathException if the inverse cumulative probability can not be
159     *         computed due to convergence or other numerical errors.
160     * @throws IllegalArgumentException if <code>p</code> is not a valid
161     *         probability.
162     */
163    @Override
164    public double inverseCumulativeProbability(final double p)
165    throws MathException {
166        if (p == 0) {
167            return Double.NEGATIVE_INFINITY;
168        }
169        if (p == 1) {
170            return Double.POSITIVE_INFINITY;
171        }
172        return super.inverseCumulativeProbability(p);
173    }
174
175    /**
176     * Access the domain value lower bound, based on <code>p</code>, used to
177     * bracket a CDF root.  This method is used by
178     * {@link #inverseCumulativeProbability(double)} to find critical values.
179     *
180     * @param p the desired probability for the critical value
181     * @return domain value lower bound, i.e.
182     *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>
183     */
184    @Override
185    protected double getDomainLowerBound(double p) {
186        return -Double.MAX_VALUE;
187    }
188
189    /**
190     * Access the domain value upper bound, based on <code>p</code>, used to
191     * bracket a CDF root.  This method is used by
192     * {@link #inverseCumulativeProbability(double)} to find critical values.
193     *
194     * @param p the desired probability for the critical value
195     * @return domain value upper bound, i.e.
196     *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code>
197     */
198    @Override
199    protected double getDomainUpperBound(double p) {
200        return Double.MAX_VALUE;
201    }
202
203    /**
204     * Access the initial domain value, based on <code>p</code>, used to
205     * bracket a CDF root.  This method is used by
206     * {@link #inverseCumulativeProbability(double)} to find critical values.
207     *
208     * @param p the desired probability for the critical value
209     * @return initial domain value
210     */
211    @Override
212    protected double getInitialDomain(double p) {
213        return 0.0;
214    }
215
216    /**
217     * Return the absolute accuracy setting of the solver used to estimate
218     * inverse cumulative probabilities.
219     *
220     * @return the solver absolute accuracy
221     * @since 2.1
222     */
223    @Override
224    protected double getSolverAbsoluteAccuracy() {
225        return solverAbsoluteAccuracy;
226    }
227
228    /**
229     * Returns the lower bound of the support for the distribution.
230     *
231     * The lower bound of the support is always negative infinity
232     * no matter the parameters.
233     *
234     * @return lower bound of the support (always Double.NEGATIVE_INFINITY)
235     * @since 2.2
236     */
237    public double getSupportLowerBound() {
238        return Double.NEGATIVE_INFINITY;
239    }
240
241    /**
242     * Returns the upper bound of the support for the distribution.
243     *
244     * The upper bound of the support is always positive infinity
245     * no matter the parameters.
246     *
247     * @return upper bound of the support (always Double.POSITIVE_INFINITY)
248     * @since 2.2
249     */
250    public double getSupportUpperBound() {
251        return Double.POSITIVE_INFINITY;
252    }
253
254    /**
255     * Returns the mean.
256     *
257     * For degrees of freedom parameter df, the mean is
258     * <ul>
259     *  <li>if <code>df &gt; 1</code> then <code>0</code></li>
260     * <li>else <code>undefined</code></li>
261     * </ul>
262     *
263     * @return the mean
264     * @since 2.2
265     */
266    public double getNumericalMean() {
267        final double df = getDegreesOfFreedom();
268
269        if (df > 1) {
270            return 0;
271        }
272
273        return Double.NaN;
274    }
275
276    /**
277     * Returns the variance.
278     *
279     * For degrees of freedom parameter df, the variance is
280     * <ul>
281     *  <li>if <code>df &gt; 2</code> then <code>df / (df - 2)</code> </li>
282     *  <li>if <code>1 &lt; df &lt;= 2</code> then <code>positive infinity</code></li>
283     *  <li>else <code>undefined</code></li>
284     * </ul>
285     *
286     * @return the variance
287     * @since 2.2
288     */
289    public double getNumericalVariance() {
290        final double df = getDegreesOfFreedom();
291
292        if (df > 2) {
293            return df / (df - 2);
294        }
295
296        if (df > 1 && df <= 2) {
297            return Double.POSITIVE_INFINITY;
298        }
299
300        return Double.NaN;
301    }
302
303}
304