8
$\begingroup$

Is there a way to closely approximate the CDF of a student-t distribution at a point $x$ without involving the hypergeometric function? For example, by using a series expansion, or expressing the CDF in terms of simpler functions?

The t-distribution can be expressed in terms of a normal and chi-squared distributions:

$$t = \frac{X\sqrt{n}}{Y}$$

where $X$ is normally distributed with mean 0 and variance $\sigma^2$; $\frac{Y^2}{\sigma^2}$ has a chi-squared distribution. Maybe this expression can help simplify the t CDF?

I am trying to write an algorithm to compute the student-t CDF for any real number, so using the hypergeometric function will be extremely difficult. I am also trying to avoid integrating the PDF to get the CDF.

$\endgroup$
2
  • 1
    $\begingroup$ There are standard, efficient algorithms to compute this CDF. They were worked out in the 60's through the 80's and published in TOMS (Transactions on Mathematical Software), typically in Fortran. It's usually straightforward to port them to any numerically-oriented computing platform. Series expansions (except around infinity for the tails) won't work for CDFs in general. Numerical integration often does well, depending on how much accuracy you need, for what values you typically need high accuracy, and how you measure that accuracy. $\endgroup$
    – whuber
    Feb 28, 2019 at 23:00
  • $\begingroup$ stats.stackexchange.com/questions/169352/… (including some useful comments) $\endgroup$
    – Glen_b
    Mar 1, 2019 at 1:55

1 Answer 1

8
$\begingroup$

The Wikipedia article on the $t$ distribution helpfully informs us that

If $X$ has a Student's t-distribution with degree of freedom $\nu$ then one can obtain a Beta distribution: $$\frac {\nu }{\nu +X^{2}}\sim \mathrm {B} \left({\frac {1}{2}},{\frac {\nu }{2}}\right).$$

Equivalently, $1-\nu/(\nu + X^2) = X^2/(\nu + X^2)$ has a $B\left(\frac{\nu}{2},\frac{1}{2}\right)$ distribution. This is referred to as a "symmetry relation" below.

Here is an implementation in a C-like language.

extern FLOAT tcum(FLOAT t, DEGREES m) {
    if (m == INFINITY) return zcum(t);
    assert(m > 0);

    if (t == 0.0) return 0.5;
    if (t < 0.0)  return 0.5 * (1.0 - inbet(t*t / (m + t*t), 1, m));
    else          return 0.5 * (1.0 + inbet(t*t / (m + t*t), 1, m));
} 

(It falls back to the standard Normal CDF zcum for arbitrarily large $\nu$.)

A good implementation of the Beta cumulative distribution inbet is given in Numerical Recipes, where inbet(x,a,b) is the incomplete Beta function $I_x(a,b).$

Long ago I ported their Fortran version of the function to C (cleaning up a few problems along the way). It wasn't too painful. It's a continued fraction expansion. Numerical Recipes remarks

This continued fraction converges rapidly for $x\lt (a+1)/(a+b+2),$ taking in the worst case $O\left(\sqrt{\max(a,b)}\right)$ iterations. But for $x \gt (a+1)/(a+b+2)$ we can just use the symmetry relation $I_x(a,b)=1-I_{1-x}(b,a)$ to obtain an equivalent computation where the continued fraction will also converge rapidly.

You can see this strategy in the extensively tested version below. You don't need the details of the header (.h) files to port this. You do need enough experience writing scientific software to know that you must test your port thoroughly and to know how to do the testing and debugging. (Reproducing extensive tables of the function to perfect accuracy is one approach. It helps to graph the results, too: that catches all kinds of erratic errors that can afflict numerical programs.)

/* 
 *  INBET.C
 *
 * The incomplete beta function with parameters a/2 and b/2, i.e.,
 * integral, 0 to x, of 
 *     pow(t,(FLOAT)a/2-1)*pow(1-t,(FLOAT)b/2-1),
 * divided by beta(a/2,b/2).
 * x must be >= 0, <= 1, a, b must be > 0.  
 */                       
#include <math.h>
#include <float.h>
#include "stats.h"

static FLOAT betacf(FLOAT a, FLOAT b, FLOAT x);
/*--------------------------------------------------------------------*/

extern FLOAT inbet(FLOAT x, DEGREES a, DEGREES b) {
    FLOAT bt;

    assert((FLOAT)0.0<=x && x<=(FLOAT)1.0);
    if (x==0.0 || x==1.0) {
        bt = 0.0;
    } else {
        bt = exp( lngamma2(a+b)-lngamma2(a)-lngamma2(b) +
                  0.5 * (a*log(x) + b*log(1.0-x)) );
    }
    if (x < (a+2.0) / (a+b+4.0))
        return 2 * bt * betacf((FLOAT)a/2.0, (FLOAT)b/2.0, x) / a;
    else
        return 1.0 - 2 * bt * betacf((FLOAT)b/2.0, (FLOAT)a/2.0, 1.0-x) / b;
} /* inbet() */
/*--------------------------------------------------------------------*/

FLOAT betacf(FLOAT a, FLOAT b, FLOAT x) {
    FLOAT am = 1.0,
          bm = 1.0,
          az = 1.0,
          qab = a+b,
          qap = a+1.0,
          qam = a-1.0,
          bz  = 1.0 - qab * x / qap;
    FLOAT em, tem, d, ap, bp, app, bpp, aold;
    long m=1;

    do {
        em = m; tem = em+em;
        d = em * (b-m) * x / ((qam+tem) * (a+tem));
        ap = az + d*am; bp = bz + d*bm;
        d = -(a+em) * (qab+em) * x / ((a+tem) * (qap+tem));
        app = ap + d*az; bpp = bp + d*bz;
        aold = az;
        am = ap/bpp; bm = bp/bpp; az = app/bpp;
        bz = 1.0;
    } while (fabs(az-aold) >= stat_eps * fabs(az) && m++ <= _stat_iter);

    if (m > _stat_iter)
        _stat_err = _ERR_STAT_EPS;

    return az;
} /* betacf() */
/* eof inbet.c */
$\endgroup$
7
  • 2
    $\begingroup$ great answer, this is exactly what I was looking for $\endgroup$
    – PyRsquared
    Mar 1, 2019 at 10:06
  • $\begingroup$ @whuber I feel like you skip a few steps going from the relationship of the t-distribution to beta and then to the C implementation. For instance, why is it $t*t / (m + t*t)$ instead of $m / (m + t * t)$? $\endgroup$
    – John
    Oct 21, 2022 at 14:22
  • 2
    $\begingroup$ @John I see what you mean. It comes down to a simple fact which is clear from inspecting the Beta density function: when $X$ has a Beta$(a,b)$ distribution, $1-X$ has a Beta$(b,a)$ distribution. $\endgroup$
    – whuber
    Oct 21, 2022 at 18:41
  • $\begingroup$ @PyRsquared I think what this answer brushes over and what your original post didn't elaborate on is the need for "any Real value". Such language makes me think you are trying to some kind of automatic differentiation, in which case this solution might not get you what you need. $\endgroup$
    – bdeonovic
    Oct 21, 2022 at 20:00
  • 1
    $\begingroup$ @mwater That's the log gamma function: it's used to compute the Normalizing constant. $\endgroup$
    – whuber
    Nov 25, 2022 at 15:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.