How to approximate the student-t CDF at a point without the hypergeometric function?

Question

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.

Answer 1

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).$$

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 */