/** This implements the statistical "Student's T" distribution which is tricky
    to implement, especially for arbitrary-precision and calculating the
    cumulative distribution function (CDF).

    TODO:  Integrate this with statistics.frink ? That intentionally has no
    dependencies on arbitrary-precision libraries nor numerical integration.

    TODO:  Implement an efficient inverse CDF function?  See the paper below.

    See:
    https://www.homepages.ucl.ac.uk/~ucahwts/lgsnotes/JCF_Student.pdf
    especially eq. (19)
*/

use pi2.frink
use NumericalIntegrationArbitrary.frink
use ArbitraryPrecision.frink

/** This gives the cumulative distribution function (CDF) of the Student-t
    distribution with the specified normalized t parameter (called x here) and n
    degrees of freedom. */
studentTProbability[x, n, digits=getPrecision[]] :=
{
   digits = max[14, digits + 5]
   origPrec = getPrecision[]
   try
   {
      setPrecision[digits]
      left = 1/sqrt[Pi.getPi[digits] * n, digits]
      left = left * gamma[(n+1)/2, digits] / gamma[n/2, digits]
      f = {|t,n| arbitraryPow[1/(1 + t^2/n), (n+1)/2] }
      res = left * NIntegrateData[f, "-inf", x, n, digits]
      return res
   }
   finally
      setPrecision[origPrec]
}

/** This gives the probability density function (PDF) of the Student-t
    distribution with the specified normalized t parameter (called x) and n
    degrees of freedom. */
studentTDensity[x, n, digits=getPrecision[]] :=
{
   digits = max[14, digits + 5]
   origPrec = getPrecision[]
   try
   {
      setPrecision[digits]
      left = 1/sqrt[Pi.getPi[digits] * n, digits]
      left = left * gamma[(n+1)/2, digits] / gamma[n/2, digits]
      right = arbitraryPow[1/(1 + x^2/n), (n+1)/2]
      res = left * right
      return res
   }
   finally
      setPrecision[origPrec]
}