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    package org.apache.commons.math.distribution;
018    
019    import java.io.Serializable;
020    
021    import org.apache.commons.math.MathException;
022    import org.apache.commons.math.MathRuntimeException;
023    import org.apache.commons.math.exception.util.LocalizedFormats;
024    import org.apache.commons.math.special.Beta;
025    import org.apache.commons.math.special.Gamma;
026    import org.apache.commons.math.util.FastMath;
027    
028    /**
029     * Default implementation of
030     * {@link org.apache.commons.math.distribution.TDistribution}.
031     *
032     * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $
033     */
034    public class TDistributionImpl
035        extends AbstractContinuousDistribution
036        implements TDistribution, Serializable  {
037    
038        /**
039         * Default inverse cumulative probability accuracy
040         * @since 2.1
041        */
042        public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
043    
044        /** Serializable version identifier */
045        private static final long serialVersionUID = -5852615386664158222L;
046    
047        /** The degrees of freedom*/
048        private double degreesOfFreedom;
049    
050        /** Inverse cumulative probability accuracy */
051        private final double solverAbsoluteAccuracy;
052    
053        /**
054         * Create a t distribution using the given degrees of freedom and the
055         * specified inverse cumulative probability absolute accuracy.
056         *
057         * @param degreesOfFreedom the degrees of freedom.
058         * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates
059         * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY})
060         * @since 2.1
061         */
062        public TDistributionImpl(double degreesOfFreedom, double inverseCumAccuracy) {
063            super();
064            setDegreesOfFreedomInternal(degreesOfFreedom);
065            solverAbsoluteAccuracy = inverseCumAccuracy;
066        }
067    
068        /**
069         * Create a t distribution using the given degrees of freedom.
070         * @param degreesOfFreedom the degrees of freedom.
071         */
072        public TDistributionImpl(double degreesOfFreedom) {
073            this(degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
074        }
075    
076        /**
077         * Modify the degrees of freedom.
078         * @param degreesOfFreedom the new degrees of freedom.
079         * @deprecated as of 2.1 (class will become immutable in 3.0)
080         */
081        @Deprecated
082        public void setDegreesOfFreedom(double degreesOfFreedom) {
083            setDegreesOfFreedomInternal(degreesOfFreedom);
084        }
085    
086        /**
087         * Modify the degrees of freedom.
088         * @param newDegreesOfFreedom the new degrees of freedom.
089         */
090        private void setDegreesOfFreedomInternal(double newDegreesOfFreedom) {
091            if (newDegreesOfFreedom <= 0.0) {
092                throw MathRuntimeException.createIllegalArgumentException(
093                      LocalizedFormats.NOT_POSITIVE_DEGREES_OF_FREEDOM,
094                      newDegreesOfFreedom);
095            }
096            this.degreesOfFreedom = newDegreesOfFreedom;
097        }
098    
099        /**
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    }