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