    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    }