com.opengamma.analytics.math.statistics.distribution.StudentTDistribution.java Source code

Introduction

Here is the source code for com.opengamma.analytics.math.statistics.distribution.StudentTDistribution.java

Source

/**
 * Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
 * 
 * Please see distribution for license.
 */
package com.opengamma.analytics.math.statistics.distribution;

import java.util.Date;

import org.apache.commons.lang.Validate;

import cern.jet.random.StudentT;
import cern.jet.random.engine.MersenneTwister64;
import cern.jet.random.engine.RandomEngine;

import com.opengamma.analytics.math.function.Function1D;
import com.opengamma.analytics.math.function.special.InverseIncompleteBetaFunction;

/**
 * Student's T-distribution is a continuous probability distribution with probability density function
 * $$
 * \begin{align*}
 * f(x) = \frac{\Gamma\left(\frac{\nu + 1}{2}\right)}{\sqrt{\nu\pi}\Gamma(\left(\frac{\nu}{2}\right)}\left(1 + \frac{x^2}{\nu}\right)^{-\frac{1}{2}(\nu + 1)}
 * \end{align*}
 * $$
 * where $\nu$ is the number of degrees of freedom and $\Gamma$ is the Gamma function ({@link com.opengamma.analytics.math.function.special.GammaFunction}).
 * <p>
 * This class is a wrapper for the <a href="http://acs.lbl.gov/software/colt/api/cern/jet/random/StudentT.html">Colt</a> implementation of Student's T-distribution for cdf and pdf calculations
 * and Student's T-distributed random number generation.
 */
public class StudentTDistribution implements ProbabilityDistribution<Double> {
    // TODO need a better seed
    private final double _degFreedom;
    private final StudentT _dist;
    private final Function1D<Double, Double> _beta;

    /**
     * @param degFreedom The number of degrees of freedom, not negative or zero
     */
    public StudentTDistribution(final double degFreedom) {
        this(degFreedom, new MersenneTwister64(new Date()));
    }

    /**
     * @param degFreedom The number of degrees of freedom, not negative or zero
     * @param engine A generator of uniform random numbers, not null
     */
    public StudentTDistribution(final double degFreedom, final RandomEngine engine) {
        Validate.isTrue(degFreedom > 0, "degrees of freedom");
        Validate.notNull(engine);
        _degFreedom = degFreedom;
        _dist = new StudentT(degFreedom, engine);
        _beta = new InverseIncompleteBetaFunction(degFreedom / 2., 0.5);
    }

    /**
     * {@inheritDoc}
     */
    @Override
    public double getCDF(final Double x) {
        Validate.notNull(x);
        return _dist.cdf(x);
    }

    /**
     * {@inheritDoc}
     */
    @Override
    public double getPDF(final Double x) {
        Validate.notNull(x);
        return _dist.pdf(x);
    }

    /**
     * {@inheritDoc}
     */
    @Override
    public double nextRandom() {
        return _dist.nextDouble();
    }

    /**
     * {@inheritDoc}
     * The inverse cdf is given by:
     * $$
     * \begin{align*}
     * F(P) &= \mathrm{sign}(p - \frac{1}{2})\sqrt{\frac{\nu}{x - 1}}\\
     * x &= B(2 \min(p, 1-p)) 
     * \end{align*}
     * $$
     * where $B$ is the inverse incomplete Beta function ({@link com.opengamma.analytics.math.function.special.InverseIncompleteBetaFunction}).
     */
    @Override
    public double getInverseCDF(final Double p) {
        Validate.notNull(p);
        Validate.isTrue(p >= 0 && p <= 1, "Probability must be >= 0 and <= 1");
        final double x = _beta.evaluate(2 * Math.min(p, 1 - p));
        return Math.signum(p - 0.5) * Math.sqrt(_degFreedom * (1. / x - 1));
    }

    /**
     * @return The number of degrees of freedom
     */
    public double getDegreesOfFreedom() {
        return _degFreedom;
    }

    @Override
    public int hashCode() {
        final int prime = 31;
        int result = 1;
        long temp;
        temp = Double.doubleToLongBits(_degFreedom);
        result = prime * result + (int) (temp ^ (temp >>> 32));
        return result;
    }

    @Override
    public boolean equals(final Object obj) {
        if (this == obj) {
            return true;
        }
        if (obj == null) {
            return false;
        }
        if (getClass() != obj.getClass()) {
            return false;
        }
        final StudentTDistribution other = (StudentTDistribution) obj;
        return Double.doubleToLongBits(_degFreedom) == Double.doubleToLongBits(other._degFreedom);
    }
}