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 */
017package org.apache.commons.math3.distribution;
018
019import org.apache.commons.math3.exception.NotStrictlyPositiveException;
020import org.apache.commons.math3.exception.util.LocalizedFormats;
021import org.apache.commons.math3.random.RandomGenerator;
022import org.apache.commons.math3.random.Well19937c;
023import org.apache.commons.math3.special.Beta;
024import org.apache.commons.math3.special.Gamma;
025import org.apache.commons.math3.util.FastMath;
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 */
033public class TDistribution extends AbstractRealDistribution {
034    /**
035     * Default inverse cumulative probability accuracy.
036     * @since 2.1
037     */
038    public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
039    /** Serializable version identifier */
040    private static final long serialVersionUID = -5852615386664158222L;
041    /** The degrees of freedom. */
042    private final double degreesOfFreedom;
043    /** Inverse cumulative probability accuracy. */
044    private final double solverAbsoluteAccuracy;
045    /** Static computation factor based on degreesOfFreedom. */
046    private final double factor;
047
048    /**
049     * Create a t distribution using the given degrees of freedom.
050     * <p>
051     * <b>Note:</b> this constructor will implicitly create an instance of
052     * {@link Well19937c} as random generator to be used for sampling only (see
053     * {@link #sample()} and {@link #sample(int)}). In case no sampling is
054     * needed for the created distribution, it is advised to pass {@code null}
055     * as random generator via the appropriate constructors to avoid the
056     * additional initialisation overhead.
057     *
058     * @param degreesOfFreedom Degrees of freedom.
059     * @throws NotStrictlyPositiveException if {@code degreesOfFreedom <= 0}
060     */
061    public TDistribution(double degreesOfFreedom)
062        throws NotStrictlyPositiveException {
063        this(degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
064    }
065
066    /**
067     * Create a t distribution using the given degrees of freedom and the
068     * specified inverse cumulative probability absolute accuracy.
069     * <p>
070     * <b>Note:</b> this constructor will implicitly create an instance of
071     * {@link Well19937c} as random generator to be used for sampling only (see
072     * {@link #sample()} and {@link #sample(int)}). In case no sampling is
073     * needed for the created distribution, it is advised to pass {@code null}
074     * as random generator via the appropriate constructors to avoid the
075     * additional initialisation overhead.
076     *
077     * @param degreesOfFreedom Degrees of freedom.
078     * @param inverseCumAccuracy the maximum absolute error in inverse
079     * cumulative probability estimates
080     * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
081     * @throws NotStrictlyPositiveException if {@code degreesOfFreedom <= 0}
082     * @since 2.1
083     */
084    public TDistribution(double degreesOfFreedom, double inverseCumAccuracy)
085        throws NotStrictlyPositiveException {
086        this(new Well19937c(), degreesOfFreedom, inverseCumAccuracy);
087    }
088
089    /**
090     * Creates a t distribution.
091     *
092     * @param rng Random number generator.
093     * @param degreesOfFreedom Degrees of freedom.
094     * @throws NotStrictlyPositiveException if {@code degreesOfFreedom <= 0}
095     * @since 3.3
096     */
097    public TDistribution(RandomGenerator rng, double degreesOfFreedom)
098        throws NotStrictlyPositiveException {
099        this(rng, degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
100    }
101
102    /**
103     * Creates a t distribution.
104     *
105     * @param rng Random number generator.
106     * @param degreesOfFreedom Degrees of freedom.
107     * @param inverseCumAccuracy the maximum absolute error in inverse
108     * cumulative probability estimates
109     * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
110     * @throws NotStrictlyPositiveException if {@code degreesOfFreedom <= 0}
111     * @since 3.1
112     */
113    public TDistribution(RandomGenerator rng,
114                         double degreesOfFreedom,
115                         double inverseCumAccuracy)
116        throws NotStrictlyPositiveException {
117        super(rng);
118
119        if (degreesOfFreedom <= 0) {
120            throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
121                                                   degreesOfFreedom);
122        }
123        this.degreesOfFreedom = degreesOfFreedom;
124        solverAbsoluteAccuracy = inverseCumAccuracy;
125
126        final double n = degreesOfFreedom;
127        final double nPlus1Over2 = (n + 1) / 2;
128        factor = Gamma.logGamma(nPlus1Over2) -
129                 0.5 * (FastMath.log(FastMath.PI) + FastMath.log(n)) -
130                 Gamma.logGamma(n / 2);
131    }
132
133    /**
134     * Access the degrees of freedom.
135     *
136     * @return the degrees of freedom.
137     */
138    public double getDegreesOfFreedom() {
139        return degreesOfFreedom;
140    }
141
142    /** {@inheritDoc} */
143    public double density(double x) {
144        return FastMath.exp(logDensity(x));
145    }
146
147    /** {@inheritDoc} */
148    @Override
149    public double logDensity(double x) {
150        final double n = degreesOfFreedom;
151        final double nPlus1Over2 = (n + 1) / 2;
152        return factor - nPlus1Over2 * FastMath.log(1 + x * x / n);
153    }
154
155    /** {@inheritDoc} */
156    public double cumulativeProbability(double x) {
157        double ret;
158        if (x == 0) {
159            ret = 0.5;
160        } else {
161            double t =
162                Beta.regularizedBeta(
163                    degreesOfFreedom / (degreesOfFreedom + (x * x)),
164                    0.5 * degreesOfFreedom,
165                    0.5);
166            if (x < 0.0) {
167                ret = 0.5 * t;
168            } else {
169                ret = 1.0 - 0.5 * t;
170            }
171        }
172
173        return ret;
174    }
175
176    /** {@inheritDoc} */
177    @Override
178    protected double getSolverAbsoluteAccuracy() {
179        return solverAbsoluteAccuracy;
180    }
181
182    /**
183     * {@inheritDoc}
184     *
185     * For degrees of freedom parameter {@code df}, the mean is
186     * <ul>
187     *  <li>if {@code df > 1} then {@code 0},</li>
188     * <li>else undefined ({@code Double.NaN}).</li>
189     * </ul>
190     */
191    public double getNumericalMean() {
192        final double df = getDegreesOfFreedom();
193
194        if (df > 1) {
195            return 0;
196        }
197
198        return Double.NaN;
199    }
200
201    /**
202     * {@inheritDoc}
203     *
204     * For degrees of freedom parameter {@code df}, the variance is
205     * <ul>
206     *  <li>if {@code df > 2} then {@code df / (df - 2)},</li>
207     *  <li>if {@code 1 < df <= 2} then positive infinity
208     *  ({@code Double.POSITIVE_INFINITY}),</li>
209     *  <li>else undefined ({@code Double.NaN}).</li>
210     * </ul>
211     */
212    public double getNumericalVariance() {
213        final double df = getDegreesOfFreedom();
214
215        if (df > 2) {
216            return df / (df - 2);
217        }
218
219        if (df > 1 && df <= 2) {
220            return Double.POSITIVE_INFINITY;
221        }
222
223        return Double.NaN;
224    }
225
226    /**
227     * {@inheritDoc}
228     *
229     * The lower bound of the support is always negative infinity no matter the
230     * parameters.
231     *
232     * @return lower bound of the support (always
233     * {@code Double.NEGATIVE_INFINITY})
234     */
235    public double getSupportLowerBound() {
236        return Double.NEGATIVE_INFINITY;
237    }
238
239    /**
240     * {@inheritDoc}
241     *
242     * The upper bound of the support is always positive infinity no matter the
243     * parameters.
244     *
245     * @return upper bound of the support (always
246     * {@code Double.POSITIVE_INFINITY})
247     */
248    public double getSupportUpperBound() {
249        return Double.POSITIVE_INFINITY;
250    }
251
252    /** {@inheritDoc} */
253    public boolean isSupportLowerBoundInclusive() {
254        return false;
255    }
256
257    /** {@inheritDoc} */
258    public boolean isSupportUpperBoundInclusive() {
259        return false;
260    }
261
262    /**
263     * {@inheritDoc}
264     *
265     * The support of this distribution is connected.
266     *
267     * @return {@code true}
268     */
269    public boolean isSupportConnected() {
270        return true;
271    }
272}