home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / math / ols / f.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-07-28  |  2.9 KB  |  143 lines
  1.  
  2. #include <stdio.h>
  3. #include <math.h>
  4. #include "distribs.h"
  5.  
  6. /* The F distribution. */
  7.  
  8. #ifdef TESTING
  9. int 
  10. main ()
  11. {
  12.   double x, p, invp;
  13.   int df1, df2;
  14.  
  15.   for (df1 = 1; df1 < 10; df1 += 4)
  16.     for (df2 = 1; df2 < 10; df2 += 4)
  17.       {
  18.     printf ("\nF(%d, %d)\n", df1, df2);
  19.     for (x = 0.5; x < 5.0; x += 2.0)
  20.       {
  21.         p = pF (x, df1, df2);
  22.         invp = critF (p, df1, df2);
  23.         printf ("x = %2.1f, Prob(>F) = %10.8f, InvF = %10.8f\n",
  24.             x, p, invp);
  25.       }
  26.       }
  27.  
  28.   return 0;
  29. }
  30.  
  31. #endif
  32.  
  33. /* Rest of this file essentially cleaned up version of Gary
  34. Perlman's source on StatLib. */
  35.  
  36. #define I_PI 0.3183098861837906715377675    /* 1/pi */
  37. #define F_EPSILON 1.0e-9    /* accuracy of critF approximation */
  38. #define F_MAX 9999.0        /* maximum F ratio */
  39.  
  40. /*FUNCTION pF: probability of F.
  41.   ALGORITHM Compute probability of F ratio.
  42.     Adapted from Collected Algorithms of the CACM
  43.     Algorithm 322
  44.     Egon Dorrer
  45. */
  46.  
  47. double 
  48. pF (double F, int df1, int df2)
  49. {
  50.   int i, j;
  51.   int a, b;
  52.   double w, y, z, d, p;
  53.  
  54.   if (F < F_EPSILON || df1 < 1 || df2 < 1)
  55.     return (1.0);
  56.  
  57.   a = df1 % 2 ? 1 : 2;
  58.   b = df2 % 2 ? 1 : 2;
  59.   w = (F * df1) / df2;
  60.   z = 1.0 / (1.0 + w);
  61.   if (a == 1)
  62.     if (b == 1)
  63.       {
  64.     p = sqrt (w);
  65.     y = I_PI;
  66.     d = y * z / p;
  67.     p = 2.0 * y * atan (p);
  68.       }
  69.     else
  70.       {
  71.     p = sqrt (w * z);
  72.     d = 0.5 * p * z / w;
  73.       }
  74.   else if (b == 1)
  75.     {
  76.       p = sqrt (z);
  77.       d = 0.5 * z * p;
  78.       p = 1.0 - p;
  79.     }
  80.   else
  81.     {
  82.       d = z * z;
  83.       p = w * z;
  84.     }
  85.   y = 2.0 * w / z;
  86.  
  87.   for (j = b + 2; j <= df2; j += 2)
  88.     {
  89.       d *= (1.0 + a / (j - 2.0)) * z;
  90.       p = (a == 1 ? p + d * y / (j - 1.0) : (p + w) * z);
  91.     }
  92.   y = w * z;
  93.   z = 2.0 / z;
  94.   b = df2 - 2;
  95.   for (i = a + 2; i <= df1; i += 2)
  96.     {
  97.       j = i + b;
  98.       d *= y * j / (i - 2.0);
  99.       p -= z * d / j;
  100.     }
  101.   /* correction for approximation errors suggested in certification */
  102.   if (p < 0.0)
  103.     p = 0.0;
  104.   else if (p > 1.0)
  105.     p = 1.0;
  106.  
  107.   return (1.0 - p);
  108. }
  109.  
  110. /* critF: compute critical F value t produce given probability.
  111.    ALGORITHM
  112.     Begin with upper and lower limits for F values (maxf and minf)
  113.     set to extremes.  Choose an f value (fval) between the extremes.
  114.     Compute the probability of the f value.  Set minf or maxf, based
  115.     on whether the probability is less than or greater than the
  116.     desired p.  Continue adjusting the extremes until they are
  117.     within F_EPSILON of each other.
  118. */
  119.  
  120. double 
  121. critF (double p, int df1, int df2)
  122. {
  123.   double fval;
  124.   double maxf = F_MAX;        /* maximum F ratio */
  125.   double minf = 0.0;        /* minimum F ratio */
  126.  
  127.   if (p <= 0.0 || p >= 1.0)
  128.     return (0.0);
  129.  
  130.   fval = 1.0 / p;        /* the smaller the p, the larger the F */
  131.  
  132.   while (fabs (maxf - minf) > F_EPSILON)
  133.     {
  134.       if (pF (fval, df1, df2) < p)    /* F too large */
  135.     maxf = fval;
  136.       else            /* F too small */
  137.     minf = fval;
  138.       fval = (maxf + minf) * 0.5;
  139.     }
  140.  
  141.   return (fval);
  142. }
  143.