root/trunk/mathextra/student.d

/**
 * Macros:
 *  GAMMA =  Γ
 *  INTEGRATE = $(BIG &#8747;<sub>$(SMALL $1)</sub><sup>$2</sup>)
 *  POWER = $1<sup>$2</sup>
 */
module mathextra.student;
private import std.math;
private import mathextra.betadist;
private import mathextra.normaldist;
debug private import std.stdio;

unittest {
    debug writefln("--- UnitTest: " __FILE__ " ---");
} 


/**************************************************
 * Computes the integral from minus infinity to t of the Student
 * t distribution with integer nu > 0 degrees of freedom:
 *
 *   $(GAMMA)( (nu+1)/2) / ( sqrt(nu &pi;) $(GAMMA)(nu/2) ) *
 * $(INTEGRATE -&infin;, t) $(POWER (1+$(POWER x, 2)/nu), -(nu+1)/2) dx
 *   
 * It is related to the incomplete beta integral:
 *        1 - studentsDistribution(nu,t) = 0.5 * betaDistribution( nu/2, 1/2, z )
 * where
 *        z = nu/(nu + t<sup>2</sup>).
 *
 * For t < -1.6, 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.
 */
real studentsDistribution(int nu, real t)
{
  // Author: Don Clugston. Public domain.
  /* Based on code from Cephes Math Library Release 2.3:  January, 1995
     Copyright 1984, 1995 by Stephen L. Moshier
 */


    if ( nu <= 0 ) return real.nan; // domain error -- or should it return 0?
    if ( t == 0.0 )  return 0.5;

    real rk, z, p;
    
    if ( t < -1.6 ) {
        rk = nu;
        z = rk / (rk + t * t);
        return 0.5L * betaDistribution( 0.5L*rk, 0.5L, z );
    }

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

    rk = nu;    /* degrees of freedom */
    
    real x;
    if (t < 0) x = -t; else x = t;
    z = 1.0L + ( x * x )/rk;

    real f, tz;
    int j;

    if ( nu & 1)    {
        /*  computation for odd nu  */
        real xsqk = x/sqrt(rk);
        p = atan( xsqk );
        if ( nu > 1 )   {
            f = 1.0L;
            tz = 1.0L;
            j = 3;
            while(  (j<=(nu-2)) && ( (tz/f) > real.epsilon )  ) {
                tz *= (j-1)/( z * j );
                f += tz;
                j += 2;
            }
            p += f * xsqk/z;
            }
        p *= 2.0L/PI;
    } else {
        /*  computation for even nu */
        f = 1.0L;
        tz = 1.0L;
        j = 2;
    
        while ( ( j <= (nu-2) ) && ( (tz/f) > real.epsilon )  ) {
            tz *= (j - 1)/( z * j );
            f += tz;
            j += 2;
        }
        p = f * x/sqrt(z*rk);
    }
    if ( t < 0.0L )
        p = -p; /* note destruction of relative accuracy */

    p = 0.5L + 0.5L * p;
    return p;
}

/******************************************
*   Inverse of Student's t distribution
*
* Given probability p and degrees of freedom nu,
* finds the argument t such that the one-sided
* studentsDistribution(nu,t) is equal to p.
* Used to test whether two distributions have the same
* standard deviation.
*
* Params:
* nu = degrees of freedom. Must be >1
* p  = probability. 0 < p < 1
*/
real studentsDistributionInv(int nu, real p )
// Author: Don Clugston. Public domain.
in {
   assert(nu>0);
   assert(p>=0.0L && p<=1.0L);
}
body
{
    if (p==0) return -real.infinity;
    if (p==1) return real.infinity;
    
    real rk, z;
    rk =  nu;
    
    if ( p > 0.25L && p < 0.75L ) {
        if ( p == 0.5L ) return 0;
        z = 1.0L - 2.0L * p;
        z = betaDistributionInv( 0.5L, 0.5L*rk, fabs(z) );
        real t = sqrt( rk*z/(1.0L-z) );
        if( p < 0.5L )
            t = -t;
        return t;
    }
    int rflg = -1; // sign of the result
    if (p >= 0.5L) {
        p = 1.0L - p;
        rflg = 1;
    }
    z = betaDistributionInv( 0.5L*rk, 0.5L, 2.0L*p );
    
    if (z<0) return rflg * real.infinity;
    return rflg * sqrt( rk/z - rk );
}

unittest {

// There are simple forms for nu = 1 and nu = 2.

// if (nu == 1), tDistribution(x) = 0.5 + atan(x)/PI
//              so tDistributionInv(p) = tan( PI * (p-0.5) );
// nu==2: tDistribution(x) = 0.5 * (1 + x/ sqrt(2+x*x) )

assert( studentsDistribution(1, -0.4)== 0.5 + atan(-0.4)/PI);
assert(studentsDistribution(2, 0.9) == 0.5L * (1 + 0.9L/sqrt(2.0L + 0.9*0.9)) );
assert(studentsDistribution(2, -5.4) == 0.5L * (1 - 5.4L/sqrt(2.0L + 5.4*5.4)) );

// Check a few spot values with statsoft.com (Mathworld values are wrong!!)
// According to statsoft.com, studentsDistributionInv(10, 0.995)= 3.16927.

// The remaining values listed here are from Excel, and are unlikely to be accurate
// in the last decimal places. However, they are helpful as a sanity check.

//  Microsoft Excel 2003 gives TINV(2*(1-0.995), 10) == 3.16927267160917
assert(isfeqabs(studentsDistributionInv(10, 0.995), 3.169_272_67L, 0.000_000_005L));

 
assert(isfeqabs(studentsDistributionInv(8, 0.6), 0.261_921_096_769_043L,0.000_000_000_05L));
// -TINV(2*0.4, 18) ==  -0.257123042655869

assert(isfeqabs(studentsDistributionInv(18, 0.4), -0.257_123_042_655_869L, 0.000_000_000_05L));
assert( feqrel(studentsDistribution(18, studentsDistributionInv(18, 0.4L)),0.4L)
 > real.mant_dig-2 );
assert( feqrel(studentsDistribution(11, studentsDistributionInv(11, 0.9L)),0.9L)
  > real.mant_dig-2);
//writefln("%.20f", studentsDistribution(18, -0.257_123_042_655_869));

}

// return true if a==b to given number of places.
bool isfeqabs(real a, real b, real diff)
{
  return fabs(a-b) < diff;
}