TDistribution.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math3.distribution;

  19. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  20. import org.apache.commons.math3.exception.util.LocalizedFormats;
  21. import org.apache.commons.math3.random.RandomGenerator;
  22. import org.apache.commons.math3.random.Well19937c;
  23. import org.apache.commons.math3.special.Beta;
  24. import org.apache.commons.math3.special.Gamma;
  25. import org.apache.commons.math3.util.FastMath;

  27. /**
  28.  * Implementation of Student's t-distribution.
  29.  *
  30.  * @see "<a href='http://en.wikipedia.org/wiki/Student&apos;s_t-distribution'>Student's t-distribution (Wikipedia)</a>"
  31.  * @see "<a href='http://mathworld.wolfram.com/Studentst-Distribution.html'>Student's t-distribution (MathWorld)</a>"
  32.  */
  33. public class TDistribution extends AbstractRealDistribution {
  34.     /**
  35.      * Default inverse cumulative probability accuracy.
  36.      * @since 2.1
  37.      */
  38.     public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
  39.     /** Serializable version identifier */
  40.     private static final long serialVersionUID = -5852615386664158222L;
  41.     /** The degrees of freedom. */
  42.     private final double degreesOfFreedom;
  43.     /** Inverse cumulative probability accuracy. */
  44.     private final double solverAbsoluteAccuracy;
  45.     /** Static computation factor based on degreesOfFreedom. */
  46.     private final double factor;

  48.     /**
  49.      * Create a t distribution using the given degrees of freedom.
  50.      * <p>
  51.      * <b>Note:</b> this constructor will implicitly create an instance of
  52.      * {@link Well19937c} as random generator to be used for sampling only (see
  53.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  54.      * needed for the created distribution, it is advised to pass {@code null}
  55.      * as random generator via the appropriate constructors to avoid the
  56.      * additional initialisation overhead.
  57.      *
  58.      * @param degreesOfFreedom Degrees of freedom.
  59.      * @throws NotStrictlyPositiveException if {@code degreesOfFreedom <= 0}
  60.      */
  61.     public TDistribution(double degreesOfFreedom)
  62.         throws NotStrictlyPositiveException {
  63.         this(degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
  64.     }

  66.     /**
  67.      * Create a t distribution using the given degrees of freedom and the
  68.      * specified inverse cumulative probability absolute accuracy.
  69.      * <p>
  70.      * <b>Note:</b> this constructor will implicitly create an instance of
  71.      * {@link Well19937c} as random generator to be used for sampling only (see
  72.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  73.      * needed for the created distribution, it is advised to pass {@code null}
  74.      * as random generator via the appropriate constructors to avoid the
  75.      * additional initialisation overhead.
  76.      *
  77.      * @param degreesOfFreedom Degrees of freedom.
  78.      * @param inverseCumAccuracy the maximum absolute error in inverse
  79.      * cumulative probability estimates
  80.      * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
  81.      * @throws NotStrictlyPositiveException if {@code degreesOfFreedom <= 0}
  82.      * @since 2.1
  83.      */
  84.     public TDistribution(double degreesOfFreedom, double inverseCumAccuracy)
  85.         throws NotStrictlyPositiveException {
  86.         this(new Well19937c(), degreesOfFreedom, inverseCumAccuracy);
  87.     }

  89.     /**
  90.      * Creates a t distribution.
  91.      *
  92.      * @param rng Random number generator.
  93.      * @param degreesOfFreedom Degrees of freedom.
  94.      * @throws NotStrictlyPositiveException if {@code degreesOfFreedom <= 0}
  95.      * @since 3.3
  96.      */
  97.     public TDistribution(RandomGenerator rng, double degreesOfFreedom)
  98.         throws NotStrictlyPositiveException {
  99.         this(rng, degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
  100.     }

  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);

  119.         if (degreesOfFreedom <= 0) {
  120.             throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
  121.                                                    degreesOfFreedom);
  122.         }
  123.         this.degreesOfFreedom = degreesOfFreedom;
  124.         solverAbsoluteAccuracy = inverseCumAccuracy;

  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.     }

  133.     /**
  134.      * Access the degrees of freedom.
  135.      *
  136.      * @return the degrees of freedom.
  137.      */
  138.     public double getDegreesOfFreedom() {
  139.         return degreesOfFreedom;
  140.     }

  142.     /** {@inheritDoc} */
  143.     public double density(double x) {
  144.         return FastMath.exp(logDensity(x));
  145.     }

  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.     }

  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.         }

  173.         return ret;
  174.     }

  176.     /** {@inheritDoc} */
  177.     @Override
  178.     protected double getSolverAbsoluteAccuracy() {
  179.         return solverAbsoluteAccuracy;
  180.     }

  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();

  194.         if (df > 1) {
  195.             return 0;
  196.         }

  198.         return Double.NaN;
  199.     }

  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();

  215.         if (df > 2) {
  216.             return df / (df - 2);
  217.         }

  219.         if (df > 1 && df <= 2) {
  220.             return Double.POSITIVE_INFINITY;
  221.         }

  223.         return Double.NaN;
  224.     }

  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.     }

  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.     }

  252.     /** {@inheritDoc} */
  253.     public boolean isSupportLowerBoundInclusive() {
  254.         return false;
  255.     }

  257.     /** {@inheritDoc} */
  258.     public boolean isSupportUpperBoundInclusive() {
  259.         return false;
  260.     }

  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. }
Created with JaCoCo 0.7.2.201409121644