Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : STDTR.C

/* stdtr.c
*
* Student's t distribution
*
*
*
* SYNOPSIS:
*
* double t, stdtr();
* short k;
*
* y = stdtr( k, t );
*
*
* DESCRIPTION:
*
* Computes the integral from minus infinity to t of the Student
* t distribution with integer k > 0 degrees of freedom:
*
* t
* -
* | |
* - | 2 -(k+1)/2
* | ( (k+1)/2 ) | ( x )
* ---------------------- | ( 1 + --- ) dx
* - | ( k )
* sqrt( k pi ) | ( k/2 ) |
* | |
* -
* -inf.
*
* Relation to incomplete beta integral:
*
* 1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z )
* where
* z = k/(k + t**2).
*
* For t < -1, this is the method of computation. For higher t,
* a direct method is derived from integration by parts.
* Since the function is symmetric about t=0, the area under the
* right tail of the density is found by calling the function
* with -t instead of t.
*
* ACCURACY:
*
* Tested at random 1 <= k <= 25. The "range" refers to t:
* Relative error:
* arithmetic domain # trials peak rms
* DEC 0,24 12000 4.7e-17 8.9e-18
* DEC -24,0 11000 2.3e-15 2.7e-16
* IEEE 0,24 30000 4.5e-16 8.0e-17
* IEEE -24,0 30000 1.9e-14 2.3e-15
*/


/*
Cephes Math Library Release 2.0: April, 1987
Copyright 1984, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

#include "mconf.h"

extern double PI, MACHEP;

double stdtr( k, t )
int k;
double t;
{
double x, rk, z, f, tz, p, xsqk;
double sqrt(), atan(), incbet();
int j;

if( k <= 0 )
{
mtherr( "stdtr", DOMAIN );
return(0.0);
}

if( t == 0 )
return( 0.5 );

if( t < -1.0 )
{
rk = k;
z = rk / (rk + t * t);
p = 0.5 * incbet( 0.5*rk, 0.5, z );
return( p );
}

/* compute integral from -t to + t */

if( t < 0 )
x = -t;
else
x = t;

rk = k; /* degrees of freedom */
z = 1.0 + ( x * x )/rk;

/* test if k is odd or even */
if( (k & 1) != 0)
{

/* computation for odd k */

xsqk = x/sqrt(rk);
p = atan( xsqk );
if( k > 1 )
{
f = 1.0;
tz = 1.0;
j = 3;
while( (j<=(k-2)) && ( (tz/f) > MACHEP ) )
{
tz *= (j-1)/( z * j );
f += tz;
j += 2;
}
p += f * xsqk/z;
}
p *= 2.0/PI;
}


else
{

/* computation for even k */

f = 1.0;
tz = 1.0;
j = 2;

while( ( j <= (k-2) ) && ( (tz/f) > MACHEP ) )
{
tz *= (j - 1)/( z * j );
f += tz;
j += 2;
}
p = f * x/sqrt(z*rk);
}

/* common exit */


if( t < 0 )
p = -p; /* note destruction of relative accuracy */

p = 0.5 + 0.5 * p;
return(p);
}