home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / gmt_stat.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-04  |  5.8 KB  |  248 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)gmt_stat.c    2.3  2/4/95
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /* stat_functions.c
  8.    some stuff for significance tests.
  9. */
  10.  
  11. #include "gmt.h"
  12.  
  13. int    inc_beta(a,b,x,ibeta)
  14. double    a, b, x, *ibeta;
  15. {
  16.     double    bt, gama, gamb, gamab, cf_beta ();
  17.  
  18.     if (a <= 0.0) {
  19.         fprintf(stderr,"Inc Beta:  Bad a (a <= 0).\n");
  20.         return(-1);
  21.     }
  22.     if (b <= 0.0) {
  23.         fprintf(stderr,"Inc Beta:  Bad b (b <= 0).\n");
  24.         return(-1);
  25.     }
  26.     if (x > 0.0 && x < 1.0) {
  27.         ln_gamma_r(a, &gama);
  28.         ln_gamma_r(b, &gamb);
  29.         ln_gamma_r( (a+b), &gamab);
  30.         bt = exp(gamab - gama - gamb
  31.             + a * d_log(x) + b * d_log(1.0 - x) );
  32.  
  33.         /* Here there is disagreement on the range of x which
  34.             converges efficiently.  Abramowitz and Stegun
  35.             say to use x < (a - 1) / (a + b - 2).  Editions
  36.             of Numerical Recipes thru mid 1987 say
  37.             x < ( (a + 1) / (a + b + 1), but the code has
  38.             x < ( (a + 1) / (a + b + 2).  Editions printed
  39.             late 1987 and after say x < ( (a + 1) / (a + b + 2)
  40.             in text as well as code.  What to do ? */
  41.  
  42.         if (x < ( (a + 1) / (a + b + 2) ) ) {
  43.             *ibeta = bt * cf_beta(a, b, x) / a;
  44.         }
  45.         else {
  46.             *ibeta = 1.0 - bt * cf_beta(b, a, (1.0 - x) ) / b;
  47.         }
  48.         return(0);
  49.     }
  50.     else if (x == 0.0) {
  51.         *ibeta = 0.0;
  52.         return(0);
  53.     }
  54.     else if (x == 1.0) {
  55.         *ibeta = 1.0;
  56.         return(0);
  57.     }
  58.     else if (x < 0.0) {
  59.         fprintf(stderr,"Inc Beta:  Bad x (x < 0).\n");
  60.         *ibeta = 0.0;
  61.     }
  62.     else if (x > 1.0) {
  63.         fprintf(stderr,"Inc Beta:  Bad x (x > 1).\n");
  64.         *ibeta = 1.0;
  65.     }
  66.     return(-1);
  67. }
  68.  
  69. int    ln_gamma_r(x, lngam)
  70. double    x, *lngam;
  71. {
  72.     /* Get natural logrithm of Gamma(x), x > 0.
  73.         To maintain full accuracy, this
  74.         routine uses Gamma(1 + x) / x when
  75.         x < 1.  This routine in turn calls
  76.         ln_gamma(x), which computes the
  77.         actual function value.  ln_gamma
  78.         assumes it is being called in a
  79.         smart way, and does not check the
  80.         range of x.  */
  81.  
  82.     if (x > 1.0) {
  83.         *lngam = ln_gamma(x);
  84.         return(0);
  85.     }
  86.     if (x > 0.0 && x < 1.0) {
  87.         *lngam = ln_gamma(1.0 + x) - d_log(x);
  88.         return(0);
  89.     }
  90.     if (x == 1.0) {
  91.         *lngam = 0.0;
  92.         return(0);
  93.     }
  94.     fprintf(stderr,"Ln Gamma:  Bad x (x <= 0).\n");
  95.     return(-1);
  96. }
  97.  
  98. double    ln_gamma(xx)
  99. double    xx;
  100. {
  101.     /* Routine to compute natural log of Gamma(x)
  102.         by Lanczos approximation.  Most accurate
  103.         for x > 1; fails for x <= 0.  No error
  104.         checking is done here; it is assumed
  105.         that this is called by ln_gamma_r()  */
  106.  
  107.     static double    cof[6] = {
  108.         76.18009173,
  109.         -86.50532033,
  110.         24.01409822,
  111.         -1.231739516,
  112.         0.120858003e-2,
  113.         -0.536382e-5
  114.     };
  115.  
  116.     static double    stp = 2.50662827465, half = 0.5, one = 1.0, fpf = 5.5;
  117.     double    x, tmp, ser;
  118.  
  119.     int    i;
  120.  
  121.     x = xx - one;
  122.     tmp = x + fpf;
  123.     tmp = (x + half) * d_log(tmp) - tmp;
  124.     ser = one;
  125.     for (i = 0; i < 6; i++) {
  126.         x += one;
  127.         ser += (cof[i]/x);
  128.     }
  129.     return(tmp + d_log(stp*ser) );
  130. }
  131.  
  132. double    cf_beta(a,b,x)
  133. double    a,b,x;
  134. {
  135.     /* Continued fraction method called by inc_beta.  */
  136.  
  137.     static int    itmax = 100;
  138.     static double    eps = 3.0e-7;
  139.  
  140.     double    am = 1.0, bm = 1.0, az = 1.0;
  141.     double    qab, qap, qam, bz, em, tem, d;
  142.     double    ap, bp, app, bpp, aold;
  143.  
  144.     int    m = 0;
  145.  
  146.     qab = a + b;
  147.     qap = a + 1.0;
  148.     qam = a - 1.0;
  149.     bz = 1.0 - qab * x / qap;
  150.  
  151.     do {
  152.         m++;
  153.         em = m;
  154.         tem = em + em;
  155.         d = em*(b-m)*x/((qam+tem)*(a+tem));
  156.         ap = az+d*am;
  157.         bp = bz+d*bm;
  158.         d = -(a+m)*(qab+em)*x/((a+tem)*(qap+tem));
  159.         app = ap+d*az;
  160.         bpp = bp+d*bz;
  161.         aold = az;
  162.         am = ap/bpp;
  163.         bm = bp/bpp;
  164.         az = app/bpp;
  165.         bz = 1.0;
  166.     } while ( ( (fabs(az-aold) ) >= (eps * fabs(az) ) ) && (m < itmax) );
  167.  
  168.     if (m == itmax)
  169.         fprintf(stderr,"BetaCF:  A or B too big, or ITMAX too small.\n");
  170.  
  171.     return(az);
  172. }
  173.  
  174. int    f_test(chisq1, nu1, chisq2, nu2, prob)
  175. double    chisq1, chisq2, *prob;
  176. int    nu1, nu2;
  177. {
  178.     /* Routine to compute the probability that
  179.         two variances are the same.
  180.         chisq1 is distributed as chisq with
  181.         nu1 degrees of freedom; ditto for
  182.         chisq2 and nu2.  If these are independent
  183.         and we form the ratio
  184.         F = max(chisq1,chisq2)/min(chisq1,chisq2)
  185.         then we can ask what is the probability 
  186.         that an F greater than this would occur
  187.         by chance.  It is this probability that
  188.         is returned in prob.  When prob is small,
  189.         it is likely that the two chisq represent
  190.         two different populations; the confidence
  191.         that the two do not represent the same pop
  192.         is 1.0 - prob.  This is a two-sided test.
  193.     This follows some ideas in Numerical Recipes, CRC Handbook,
  194.     and Abramowitz and Stegun.  */
  195.  
  196.     double    f, df1, df2, p1, p2;
  197.  
  198.     if ( chisq1 <= 0.0) {
  199.         fprintf(stderr,"f_test:  Chi-Square One <= 0.0\n");
  200.         return(-1);
  201.     }
  202.     if ( chisq2 <= 0.0) {
  203.         fprintf(stderr,"f_test:  Chi-Square Two <= 0.0\n");
  204.         return(-1);
  205.     }
  206.     if (chisq1 > chisq2) {
  207.         f = chisq1/chisq2;
  208.         df1 = nu1;
  209.         df2 = nu2;
  210.     }
  211.     else {
  212.         f = chisq2/chisq1;
  213.         df1 = nu2;
  214.         df2 = nu1;
  215.     }
  216.     if (inc_beta(0.5*df2, 0.5*df1, df2/(df2+df1*f), &p1) ) {
  217.         fprintf(stderr,"f_test:  Trouble on 1st inc_beta call.\n");
  218.         return(-1);
  219.     }
  220.     if (inc_beta(0.5*df1, 0.5*df2, df1/(df1+df2/f), &p2) ) {
  221.         fprintf(stderr,"f_test:  Trouble on 2nd inc_beta call.\n");
  222.         return(-1);
  223.     }
  224.     *prob = p1 + (1.0 - p2);
  225.     return(0);
  226. }
  227.  
  228. int    sig_f(chi1, n1, chi2, n2, level, prob)
  229. double    chi1, chi2, level, *prob;
  230. int    n1, n2;
  231. {
  232.     /* Returns TRUE if chi1 significantly less than chi2
  233.         at the level level.  Returns FALSE if:
  234.             chi1 > chi2;
  235.             error occurs in f_test();
  236.             chi1 < chi2 but not significant at level.  */
  237.  
  238.     int    trouble;
  239.     double    p;
  240.  
  241.     trouble = f_test(chi1, n1, chi2, n2, &p);
  242.     *prob = (1.0 - 0.5*p);    /* Make it one-sided  */
  243.     if (trouble) return(0);
  244.     if (chi1 > chi2) return (0);
  245.     return((*prob) >= level);
  246. }
  247.  
  248.