math.stats.distribution.StudentT Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of finwhale Show documentation
Statistical distributions library (in statu nascendi)
The newest version!
/*
Copyright ? 1999 CERN - European Organization for Nuclear Research.
Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
is hereby granted without fee, provided that the above copyright notice appear in all copies and
that both that copyright notice and this permission notice appear in supporting documentation.
CERN makes no representations about the suitability of this software for any purpose.
It is provided "as is" without expressed or implied warranty.
 */
package math.stats.distribution;

import math.FastGamma;
import math.FastMath;
import math.rng.DefaultRng;
import math.rng.PseudoRandom;
import math.stats.ProbabilityFuncs;

/**
 * StudentT distribution (a.k.a "T distribution").
 * 
 * p(x) = const  *  (1 + x^2/ν) ^ -(ν+1)/2 where
 * const = Γ((ν+1)/2) / (√(Π*ν) * Γ(ν/2))
 * and Γ(a) being the Gamma function and ν being the
 * degrees of freedom.
 * 
 * Implementation: This is a port of RandStudentT used in CLHEP 1.4.0 (C++). CLHEP's
 * implementation, in turn, is based on tpol.c from the C-RAND /
 * WIN-RAND library. C-RAND's implementation, in turn, is based upon
 * 
 * R.W. Bailey (1994): Polar generation of random variates with the
 * t-distribution, Mathematics of Computation 62, 779-781.
 * 
 * @author [email protected]
 * @version 1.0, 09/24/99
 */
public class StudentT extends AbstractContinuousDistribution {

    private final double df;
    private final double pdfConst;

    public StudentT(double df) {
        this(DefaultRng.newPseudoRandom(), df);
    }

    public StudentT(PseudoRandom prng, double df) {
        super(prng);
        if (df <= 0.0) {
            throw new IllegalArgumentException("df <= 0.0 : " + df);
        }
        final double tmp = FastGamma.logGamma((df + 1.0) / 2.0)
                - FastGamma.logGamma(df / 2.0);
        this.pdfConst = FastMath.exp(tmp) / Math.sqrt(Math.PI * df);
        this.df = df;
    }

    @Override
    public double pdf(double x) {
        return pdfConst * FastMath.pow((1.0 + x * x / df), -(df + 1.0) * 0.5);
    }

    @Override
    public double cdf(double x) {
        return ProbabilityFuncs.studentT(df, x);
    }

    @Override
    public double sample() {
        /*
         * Marsaglia's formulation of the Box/Muller polar method for generating
         * Normal variates is adapted to the Student-t distribution. The two
         * generated variates are not independent and the expected number of
         * uniforms per variate is 2.5464.
         * 
         * Reference:
         * 
         * R.W. Bailey (1994): Polar generation of random variates with the
         * t-distribution, Mathematics of Computation 62, 779-781.
         */
        double u1, u2, q;
        do {
            u1 = 2.0 * prng.nextDouble() - 1.0; // between -1 and 1
            u2 = 2.0 * prng.nextDouble() - 1.0; // between -1 and 1
            q = u1 * u1 + u2 * u2;
        } while (q > 1.0);
        return u1
                * Math.sqrt(df * (FastMath.exp(-2.0 / df * Math.log(q)) - 1.0)
                        / q);
    }

    @Override
    public double mean() {
        if (df <= 1.0) {
            return Double.NaN;
        }
        return 0.0;
    }

    @Override
    public double variance() {
        if (df > 2.0) {
            return df / ((double) df - 2.0);
        }
        if (df == 2.0) {
            return Double.POSITIVE_INFINITY;
        }
        return Double.NaN;
    }

    public double getDegreesOfFreedom() {
        return df;
    }

    @Override
    public String toString() {
        return getSimpleName(df);
    }
}