    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.math3.distribution;
    018
    019    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    020    import org.apache.commons.math3.exception.util.LocalizedFormats;
    021    import org.apache.commons.math3.special.Beta;
    022    import org.apache.commons.math3.special.Gamma;
    023    import org.apache.commons.math3.util.FastMath;
    024    import org.apache.commons.math3.random.RandomGenerator;
    025    import org.apache.commons.math3.random.Well19937c;
    026
    027    /**
    028     * Implementation of Student's t-distribution.
    029     *
    030     * @see "<a href='http://en.wikipedia.org/wiki/Student&apos;s_t-distribution'>Student's t-distribution (Wikipedia)</a>"
    031     * @see "<a href='http://mathworld.wolfram.com/Studentst-Distribution.html'>Student's t-distribution (MathWorld)</a>"
    032     * @version $Id: TDistribution.java 1416643 2012-12-03 19:37:14Z tn $
    033     */
    034    public class TDistribution extends AbstractRealDistribution {
    035        /**
    036         * Default inverse cumulative probability accuracy.
    037         * @since 2.1
    038         */
    039        public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
    040        /** Serializable version identifier */
    041        private static final long serialVersionUID = -5852615386664158222L;
    042        /** The degrees of freedom. */
    043        private final double degreesOfFreedom;
    044        /** Inverse cumulative probability accuracy. */
    045        private final double solverAbsoluteAccuracy;
    046
    047        /**
    048         * Create a t distribution using the given degrees of freedom.
    049         *
    050         * @param degreesOfFreedom Degrees of freedom.
    051         * @throws NotStrictlyPositiveException if {@code degreesOfFreedom <= 0}
    052         */
    053        public TDistribution(double degreesOfFreedom)
    054            throws NotStrictlyPositiveException {
    055            this(degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    056        }
    057
    058        /**
    059         * Create a t distribution using the given degrees of freedom and the
    060         * specified inverse cumulative probability absolute accuracy.
    061         *
    062         * @param degreesOfFreedom Degrees of freedom.
    063         * @param inverseCumAccuracy the maximum absolute error in inverse
    064         * cumulative probability estimates
    065         * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
    066         * @throws NotStrictlyPositiveException if {@code degreesOfFreedom <= 0}
    067         * @since 2.1
    068         */
    069        public TDistribution(double degreesOfFreedom, double inverseCumAccuracy)
    070            throws NotStrictlyPositiveException {
    071            this(new Well19937c(), degreesOfFreedom, inverseCumAccuracy);
    072        }
    073
    074        /**
    075         * Creates a t distribution.
    076         *
    077         * @param rng Random number generator.
    078         * @param degreesOfFreedom Degrees of freedom.
    079         * @param inverseCumAccuracy the maximum absolute error in inverse
    080         * cumulative probability estimates
    081         * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
    082         * @throws NotStrictlyPositiveException if {@code degreesOfFreedom <= 0}
    083         * @since 3.1
    084         */
    085        public TDistribution(RandomGenerator rng,
    086                             double degreesOfFreedom,
    087                             double inverseCumAccuracy)
    088            throws NotStrictlyPositiveException {
    089            super(rng);
    090
    091            if (degreesOfFreedom <= 0) {
    092                throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
    093                                                       degreesOfFreedom);
    094            }
    095            this.degreesOfFreedom = degreesOfFreedom;
    096            solverAbsoluteAccuracy = inverseCumAccuracy;
    097        }
    098
    099        /**
    100         * Access the degrees of freedom.
    101         *
    102         * @return the degrees of freedom.
    103         */
    104        public double getDegreesOfFreedom() {
    105            return degreesOfFreedom;
    106        }
    107
    108        /** {@inheritDoc} */
    109        public double density(double x) {
    110            final double n = degreesOfFreedom;
    111            final double nPlus1Over2 = (n + 1) / 2;
    112            return FastMath.exp(Gamma.logGamma(nPlus1Over2) -
    113                                0.5 * (FastMath.log(FastMath.PI) +
    114                                       FastMath.log(n)) -
    115                                Gamma.logGamma(n / 2) -
    116                                nPlus1Over2 * FastMath.log(1 + x * x / n));
    117        }
    118
    119        /** {@inheritDoc} */
    120        public double cumulativeProbability(double x) {
    121            double ret;
    122            if (x == 0) {
    123                ret = 0.5;
    124            } else {
    125                double t =
    126                    Beta.regularizedBeta(
    127                        degreesOfFreedom / (degreesOfFreedom + (x * x)),
    128                        0.5 * degreesOfFreedom,
    129                        0.5);
    130                if (x < 0.0) {
    131                    ret = 0.5 * t;
    132                } else {
    133                    ret = 1.0 - 0.5 * t;
    134                }
    135            }
    136
    137            return ret;
    138        }
    139
    140        /** {@inheritDoc} */
    141        @Override
    142        protected double getSolverAbsoluteAccuracy() {
    143            return solverAbsoluteAccuracy;
    144        }
    145
    146        /**
    147         * {@inheritDoc}
    148         *
    149         * For degrees of freedom parameter {@code df}, the mean is
    150         * <ul>
    151         *  <li>if {@code df > 1} then {@code 0},</li>
    152         * <li>else undefined ({@code Double.NaN}).</li>
    153         * </ul>
    154         */
    155        public double getNumericalMean() {
    156            final double df = getDegreesOfFreedom();
    157
    158            if (df > 1) {
    159                return 0;
    160            }
    161
    162            return Double.NaN;
    163        }
    164
    165        /**
    166         * {@inheritDoc}
    167         *
    168         * For degrees of freedom parameter {@code df}, the variance is
    169         * <ul>
    170         *  <li>if {@code df > 2} then {@code df / (df - 2)},</li>
    171         *  <li>if {@code 1 < df <= 2} then positive infinity
    172         *  ({@code Double.POSITIVE_INFINITY}),</li>
    173         *  <li>else undefined ({@code Double.NaN}).</li>
    174         * </ul>
    175         */
    176        public double getNumericalVariance() {
    177            final double df = getDegreesOfFreedom();
    178
    179            if (df > 2) {
    180                return df / (df - 2);
    181            }
    182
    183            if (df > 1 && df <= 2) {
    184                return Double.POSITIVE_INFINITY;
    185            }
    186
    187            return Double.NaN;
    188        }
    189
    190        /**
    191         * {@inheritDoc}
    192         *
    193         * The lower bound of the support is always negative infinity no matter the
    194         * parameters.
    195         *
    196         * @return lower bound of the support (always
    197         * {@code Double.NEGATIVE_INFINITY})
    198         */
    199        public double getSupportLowerBound() {
    200            return Double.NEGATIVE_INFINITY;
    201        }
    202
    203        /**
    204         * {@inheritDoc}
    205         *
    206         * The upper bound of the support is always positive infinity no matter the
    207         * parameters.
    208         *
    209         * @return upper bound of the support (always
    210         * {@code Double.POSITIVE_INFINITY})
    211         */
    212        public double getSupportUpperBound() {
    213            return Double.POSITIVE_INFINITY;
    214        }
    215
    216        /** {@inheritDoc} */
    217        public boolean isSupportLowerBoundInclusive() {
    218            return false;
    219        }
    220
    221        /** {@inheritDoc} */
    222        public boolean isSupportUpperBoundInclusive() {
    223            return false;
    224        }
    225
    226        /**
    227         * {@inheritDoc}
    228         *
    229         * The support of this distribution is connected.
    230         *
    231         * @return {@code true}
    232         */
    233        public boolean isSupportConnected() {
    234            return true;
    235        }
    236    }