home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / statlib.zip / STATIST.C < prev   
C/C++ Source or Header  |  1998-02-12  |  4KB  |  147 lines
  1. /***********************************************************
  2.   Statistische Verteilungsfunktionen als Bibliotheksroutinen
  3.   Neubearbeitung meines Omikron-Basic-Merge-Files SIGNIF.BAS
  4.   optimiert für C und Coprozessor
  5.   (c) 1998 by Heinz Repp
  6.   zuletzt bearbeitet am 12.02.98
  7. ***********************************************************/
  8.  
  9. #include <math.h>
  10.  
  11. #if ! defined PI
  12. #define PI 3.14159265358979323846
  13. #endif
  14.  
  15. #if (defined CHIQUADRAT || defined STANDARDNORMAL) && ! defined EPSILON
  16. #define EPSILON 5E-9
  17. #endif
  18.  
  19. #if defined FTEST || defined TSTUDENT
  20.  
  21. static double sum_sub (register double f, double a,
  22.                        register double b, double d)
  23. { register double s;
  24.   if (b < a) return 0.0;
  25.   s = 1.0;
  26.   for (b -= 2.0; b >= a; b -= 2.0)
  27.   { s *= (b + d) / b * f;
  28.     s += 1.0;
  29.   } /* endfor */
  30.   return s;
  31. } /* sum_sub */
  32.  
  33. #endif
  34.  
  35. #if defined FTEST
  36.  
  37. static double pow_i (register double f, register unsigned int i)
  38. { register double p;
  39.   p = f;    /* i >= 1 required */
  40.   while (--i) p *= f;
  41.   return p;
  42. } /* pow_i */
  43.  
  44. double distribution_of_F (double f, unsigned int f1, unsigned int f2)
  45. { double i, p;
  46.   if (! (f1 && f2))
  47.     return -1;
  48.   f *= (double) f1 / (double) f2;
  49.   if (f1 & f2 & 1)      /* f1 and f2 odd */
  50.   { p = - sum_sub (f / (f + 1.0), 3.0, (double) f1, (double) f2 - 2.0);
  51.     for (i = (double) f2 - 2.0; i >= 1.0; i -= 2.0)
  52.       p *= (i + 1.0) / i;
  53.     p *= sqrt (f) / pow_i (f + 1.0, (f2 + 1) >> 1);
  54.     p += sum_sub (1.0 / (f + 1.0), 3.0, (double) f2, -1.0)
  55.          * sqrt (f) / (f + 1.0);
  56.     p += atan (sqrt (f));
  57.     p = (p + p) / PI;
  58.   } else
  59.   { if (f1 & 1 || ! (f2 & 1) && f1 > f2)
  60.       p = sum_sub (1.0 / (f + 1.0), 2.0, (double) f2, (double) f1 - 2.0)
  61.           * pow_i (sqrt (f / (f + 1.0)), f1);
  62.     else
  63.       p = 1.0 - sum_sub (f / (f + 1.0), 2.0, (double) f1, (double) f2 - 2.0)
  64.                 / pow_i (sqrt (f + 1.0), f2);
  65.   } /* endif */
  66.   return p;
  67. } /* distribution_of_F */
  68.  
  69. #endif
  70.  
  71. #if defined TSTUDENT
  72.  
  73. double distribution_of_t (double t, unsigned int f)
  74. { double p;
  75.   if (! f)
  76.     return -1;
  77.   t /= sqrt ((double) f);
  78.   if (f & 1)      /* df odd */
  79.     p = (sum_sub (1.0 / (t * t + 1.0), 3.0, (double) f, -1.0)
  80.          * t / (t * t + 1.0) + atan (t)) / PI;
  81.   else            /* df even */
  82.     p = sum_sub (1.0 / (t * t + 1.0), 2.0, (double) f, -1.0)
  83.         * t / sqrt (t * t + 1.0) / 2.0;
  84.   return 0.5 + p;
  85. } /* distribution_of_t */
  86.  
  87. #endif
  88.  
  89. #if defined CHIQUADRAT
  90.  
  91. double distribution_of_Chi2 (double x, unsigned int f)
  92. { double i, j, p, g;
  93.   if (! f)
  94.     return -1;
  95.   if (x <= 0.0)
  96.     return 0;
  97.   x /= 2.0;
  98.   if (f & 1)                  /* df odd */
  99.   { p = exp (- x) / sqrt (PI * x);
  100.     if (p == 0.0) return 1.0;          /* underflow */
  101.     j = (double) f / 2;
  102.     for (i = 0.5; i <= j; i += 1.0)    /* 1/2 to f/2 */
  103.       p *= x / i;
  104.     for (j = p; i <= x; i += 1.0)      /* f/2 to x/2 */
  105.     { j *= x / i;
  106.       p += j;
  107.     } /* endfor2 */
  108.     g = EPSILON / x;
  109.     for ( ; j > (i - x) * g; i += 1.0) /* x/2 to accuracy limit */
  110.     { j *= x / i;
  111.       p += j;
  112.     } /* endfor3 */
  113.     return p;
  114.   } else                      /* df even */
  115.   { p = 1.0;
  116.     for (f >>= 1; --f; )
  117.       p = p * x / (double) f + 1.0;
  118.     return 1.0 - p * exp (- x);
  119.   } /* endif */
  120. } /* distribution_of_Chi2 */
  121.  
  122. #endif
  123.  
  124. #if defined STANDARDNORMAL
  125.  
  126. double distribution_of_Z (double z)
  127. { double g, i, j, p;
  128.   if (z == 0) return 0.5;
  129.   p = z;
  130.   g = z < 0 ? - EPSILON : EPSILON;
  131.   z = z * z / 2.0;
  132.   p *= exp (- z) / sqrt ((PI + PI));
  133.   if (p == 0.0) return g < 0 ? 0.0 : 1.0; /* underflow */
  134.   for (i = p, j = 1.5; j <= z; j += 1)
  135.   { i *= z / j;
  136.     p += i;
  137.   } /* endfor1 */
  138.   g /= z;
  139.   for ( ; i / g > j - z; j += 1)
  140.   { i *= z / j;
  141.     p += i;
  142.   } /* endfor2 */
  143.   return 0.5 + p;
  144. } /* distribution_of_Z */
  145.  
  146. #endif
  147.