6
$\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$
  • 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 '19 at 23:00
  • $\begingroup$ stats.stackexchange.com/questions/169352/… (including some useful comments) $\endgroup$ – Glen_b -Reinstate Monica Mar 1 '19 at 1:55
6
$\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).$$

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$
  • 1
    $\begingroup$ great answer, this is exactly what I was looking for $\endgroup$ – PyRsquared Mar 1 '19 at 10:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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