home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / cprob / stdtr.c < prev   
Encoding:
C/C++ Source or Header  |  1992-11-17  |  2.9 KB  |  153 lines
  1. /*                            stdtr.c
  2.  *
  3.  *    Student's t distribution
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double t, stdtr();
  10.  * short k;
  11.  *
  12.  * y = stdtr( k, t );
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Computes the integral from minus infinity to t of the Student
  18.  * t distribution with integer k > 0 degrees of freedom:
  19.  *
  20.  *                                      t
  21.  *                                      -
  22.  *                                     | |
  23.  *              -                      |         2   -(k+1)/2
  24.  *             | ( (k+1)/2 )           |  (     x   )
  25.  *       ----------------------        |  ( 1 + --- )        dx
  26.  *                     -               |  (      k  )
  27.  *       sqrt( k pi ) | ( k/2 )        |
  28.  *                                   | |
  29.  *                                    -
  30.  *                                   -inf.
  31.  * 
  32.  * Relation to incomplete beta integral:
  33.  *
  34.  *        1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z )
  35.  * where
  36.  *        z = k/(k + t**2).
  37.  *
  38.  * For t < -1, this is the method of computation.  For higher t,
  39.  * a direct method is derived from integration by parts.
  40.  * Since the function is symmetric about t=0, the area under the
  41.  * right tail of the density is found by calling the function
  42.  * with -t instead of t.
  43.  * 
  44.  * ACCURACY:
  45.  *
  46.  * Tested at random 1 <= k <= 25.  The "range" refers to t:
  47.  *                      Relative error:
  48.  * arithmetic   domain     # trials      peak         rms
  49.  *    DEC       0,24        12000       4.7e-17     8.9e-18
  50.  *    DEC       -24,0       11000       2.3e-15     2.7e-16
  51.  *    IEEE      0,24        30000       4.5e-16     8.0e-17
  52.  *    IEEE      -24,0       30000       1.9e-14     2.3e-15
  53.  */
  54.  
  55.  
  56. /*
  57. Cephes Math Library Release 2.0:  April, 1987
  58. Copyright 1984, 1987 by Stephen L. Moshier
  59. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  60. */
  61.  
  62. #include "mconf.h"
  63.  
  64. extern double PI, MACHEP;
  65.  
  66. double stdtr( k, t )
  67. int k;
  68. double t;
  69. {
  70. double x, rk, z, f, tz, p, xsqk;
  71. double sqrt(), atan(), incbet();
  72. int j;
  73.  
  74. if( k <= 0 )
  75.     {
  76.     mtherr( "stdtr", DOMAIN );
  77.     return(0.0);
  78.     }
  79.  
  80. if( t == 0 )
  81.     return( 0.5 );
  82.  
  83. if( t < -1.0 )
  84.     {
  85.     rk = k;
  86.     z = rk / (rk + t * t);
  87.     p = 0.5 * incbet( 0.5*rk, 0.5, z );
  88.     return( p );
  89.     }
  90.  
  91. /*    compute integral from -t to + t */
  92.  
  93. if( t < 0 )
  94.     x = -t;
  95. else
  96.     x = t;
  97.  
  98. rk = k;    /* degrees of freedom */
  99. z = 1.0 + ( x * x )/rk;
  100.  
  101. /* test if k is odd or even */
  102. if( (k & 1) != 0)
  103.     {
  104.  
  105.     /*    computation for odd k    */
  106.  
  107.     xsqk = x/sqrt(rk);
  108.     p = atan( xsqk );
  109.     if( k > 1 )
  110.         {
  111.         f = 1.0;
  112.         tz = 1.0;
  113.         j = 3;
  114.         while(  (j<=(k-2)) && ( (tz/f) > MACHEP )  )
  115.             {
  116.             tz *= (j-1)/( z * j );
  117.             f += tz;
  118.             j += 2;
  119.             }
  120.         p += f * xsqk/z;
  121.         }
  122.     p *= 2.0/PI;
  123.     }
  124.  
  125.  
  126. else
  127.     {
  128.  
  129.     /*    computation for even k    */
  130.  
  131.     f = 1.0;
  132.     tz = 1.0;
  133.     j = 2;
  134.  
  135.     while(  ( j <= (k-2) ) && ( (tz/f) > MACHEP )  )
  136.         {
  137.         tz *= (j - 1)/( z * j );
  138.         f += tz;
  139.         j += 2;
  140.         }
  141.     p = f * x/sqrt(z*rk);
  142.     }
  143.  
  144. /*    common exit    */
  145.  
  146.  
  147. if( t < 0 )
  148.     p = -p;    /* note destruction of relative accuracy */
  149.  
  150.     p = 0.5 + 0.5 * p;
  151. return(p);
  152. }
  153.