home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / math / ols / t.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-07-28  |  3.5 KB  |  153 lines
  1. /* Algorithm AS 3 from
  2.    Applied Statistics (1968) Vol. 17, p. 189.
  3.    Student's t distribution probability (lower tail).
  4.  
  5.    Fortran source from statlib,
  6.    translated by f2c,
  7.    cleaned up A LOT by hand by Ajay Shah, Sun May  5 03:48:17 PDT 1991
  8. */
  9.  
  10. #include <math.h>
  11. #include "distribs.h"
  12.  
  13. double critt (double p, int dof);
  14.  
  15. #ifdef TESTING
  16. int 
  17. main ()
  18. {
  19.   int dof;
  20.   double p, t, probability;
  21.   void printf ();
  22.  
  23.   printf ("\nt Probabilities:\n\n");
  24.   printf ("dof     Pr(t<0.5)   Pr(t<1.5)   Pr(t<2.5)   Pr(t<3.5)   Pr(t<4.5)\n\n");
  25.   for (dof = 1; dof < 35; dof += 7)
  26.     {
  27.       printf ("%2d     ", dof);
  28.       for (t = 0.5; t <= 5.0; t += 1.0)
  29.     if (0 == tprob (t, dof, &probability))
  30.       printf ("  %f  ", probability);
  31.     else
  32.       {
  33.         printf ("Ouch!\n");
  34.         return 1;
  35.       }
  36.       printf ("\n");
  37.     }
  38.   printf ("\nNormal:");
  39.   for (t = 0.5; t <= 5.0; t += 1.0)
  40.     printf ("  %f  ", pnorm1 (t));
  41.   printf ("\n\n\n");
  42.  
  43.   printf ("Inverse probabilities:\n\n");
  44.   printf ("       %10s%10s%10s%10s\n\n",
  45.       "Prob=0.2", "Prob=0.4", "Prob=0.6", "Prob=0.8");
  46.   for (dof = 1; dof < 35; dof += 6)
  47.     {
  48.       printf ("%2d     ", dof);
  49.       for (p = 0.2; p <= 0.8; p += 0.2)
  50.     printf ("%10.6f", critt (p, dof));
  51.       printf ("\n");
  52.     }
  53.   printf ("\nNormal:");
  54.   for (p = 0.2; p <= 0.8; p += 0.2)
  55.     printf ("%10.6f", critz (p));
  56.   printf ("\n");
  57.  
  58.   return 0;
  59. }
  60.  
  61. #endif
  62.  
  63. int 
  64. tprob (double t, int dof, double *probability)
  65.  
  66.      /* Arguments:
  67.     t : wanna evaluate integral from t to +infinity.
  68.     dof: degrees of freedom of this t distribution.
  69.     probability : computed value might be returned.
  70.  
  71.     function returns error code 1 if something is wrong,
  72.     0 if all is well.
  73. */
  74.  
  75. #define ONEBYPI ((double) 0.3183098861837906715377675)
  76.  
  77. {
  78.   double d_dof, s, c, f, a, b;
  79.   int fk, ks, im2, ioe, k;
  80.  
  81.   if (dof < 1)
  82.     return 1;
  83.   d_dof = (double) dof;        /* d_dof is F of fortran code. */
  84.  
  85.   a = t / sqrt (d_dof);
  86.   b = d_dof / (d_dof + (t * t));
  87.   im2 = dof - 2;
  88.   ioe = dof % 2;
  89.   s = c = f = 1.0;
  90.   fk = ks = 2 + ioe;
  91.   if (im2 > 2)
  92.     for (k = ks; k <= im2; k += 2)
  93.       {
  94.     c = c * b * (fk - 1.0) / fk;
  95.     s += c;
  96.     if (s == f)
  97.       break;
  98.     f = s;
  99.     fk += 2.0;
  100.       }
  101.   if (ioe != 1)
  102.     {                /* Label 20 of fortran code. */
  103.       *probability = 0.5 + (0.5 * a * sqrt (b) * s);
  104.       return 0;
  105.     }
  106.   else
  107.     {                /* Label 30 of fortran code. */
  108.       if (dof == 1)
  109.     s = 0.0;
  110.       *probability = 0.5 + ((a * b * s + atan (a)) * ONEBYPI);
  111.       return 0;
  112.     }
  113. }
  114.  
  115. #undef ONEBYPI
  116.  
  117. #define T_MAX 35.0
  118. #define T_EPSILON 0.0001
  119.  
  120. /* critt: compute critical t value to produce given probability.
  121.    ALGORITHM
  122.     Begin with upper and lower limits for t values (maxt and mint)
  123.     set to extremes.  Choose a t value (tval) between the extremes.
  124.     Compute the probability of the t value.  Set mint or maxt, based
  125.     on whether the probability is less than or greater than the
  126.     desired p.  Continue adjusting the extremes until they are
  127.     within T_EPSILON of each other.
  128. */
  129.  
  130. double 
  131. critt (double p, int dof)
  132. {
  133.   double pt, tval;
  134.   double maxt = T_MAX;        /* maximum t ratio */
  135.   double mint = -T_MAX;        /* minimum t ratio */
  136.  
  137.   if (p <= 0.0 || p >= 1.0)
  138.     return (0.0);
  139.  
  140.   tval = 0.0;
  141.  
  142.   while (fabs (maxt - mint) > T_EPSILON)
  143.     {
  144.       tprob (tval, dof, &pt);
  145.       if (pt > p)        /* t too large */
  146.     maxt = tval;
  147.       else            /* t too small */
  148.     mint = tval;
  149.       tval = (maxt + mint) / 2.0;
  150.     }
  151.   return (tval);
  152. }
  153.