/[gsl]/gsl/cdf/f.c
ViewVC logotype

Diff of /gsl/cdf/f.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.2.7 by bjg, Mon Apr 7 09:51:02 2003 UTC revision 1.1.2.8 by bjg, Sun Jun 22 21:20:36 2003 UTC
# Line 19  Line 19 
19    
20  #include <config.h>  #include <config.h>
21  #include <gsl/gsl_cdf.h>  #include <gsl/gsl_cdf.h>
22  #include <gsl/gsl_sf.h>  #include <gsl/gsl_sf_gamma.h>
 #include <gsl/gsl_sf_result.h>  
23  #include <gsl/gsl_errno.h>  #include <gsl/gsl_errno.h>
24    
25  /*  /*
26   * Lower tail.   * Lower tail.
27   */   */
28  double  double
29  gsl_cdf_F_ratio_P(double x, double nu1, double nu2)  gsl_cdf_F_ratio_P (double x, double nu1, double nu2)
30  {  {
31        double P;
32    int rc;    double r = nu2 / nu1;
33    double F, nu1_2, nu2_2;  
34    gsl_sf_result result;    if (x < r)
35          {
36    F = nu2 / (nu2 + nu1 * x);        double u = x / (r + x);
37    nu1_2 = nu1 / 2.0;  
38    nu2_2 = nu2 / 2.0;        P = gsl_sf_beta_inc (nu1 / 2.0, nu2 / 2.0, u);
39          }
40    rc = gsl_sf_beta_inc_e(nu2_2, nu1_2, F, &result);    else
41        {
42    /*        double u = r / (r + x);
43     * I don't like subtracting from 1.0 here. It could cause  
44     * too much cancellation error. A direct approximation for        P = 1 - gsl_sf_beta_inc (nu2 / 2.0, nu1 / 2.0, u);
    * the lower tail would be better.  
    */  
   if(rc == GSL_SUCCESS)  
     {  
       result.val = 1.0 - result.val;  
45      }      }
46        
47    return result.val;    return P;
48  }  }
49    
50  /*  /*
51   * Upper tail.   * Upper tail.
52   */   */
53  double  double
54  gsl_cdf_F_ratio_Q(double x, double nu1, double nu2)  gsl_cdf_F_ratio_Q (double x, double nu1, double nu2)
55  {  {
56    int rc;    double P;
57    double F, nu1_2, nu2_2;    double r = nu2 / nu1;
   gsl_sf_result result;  
     
   F = nu2 / (nu2 + nu1 * x);  
   nu1_2 = nu1 / 2.0;  
   nu2_2 = nu2 / 2.0;  
     
   rc = gsl_sf_beta_inc_e(nu2_2, nu1_2, F, &result);  
58    
59    return result.val;    if (x < r)
60        {
61          double u = x / (r + x);
62    
63  }        P = 1 - gsl_sf_beta_inc (nu1 / 2.0, nu2 / 2.0, u);
64        }
65      else
66        {
67          double u = r / (r + x);
68    
69  /*        P = gsl_sf_beta_inc (nu2 / 2.0, nu1 / 2.0, u);
70   * Compute the F distribution function by transforming      }
  * to the incomplete beta function.  
  */  
 /* int gsl_cdf_F_e ( double x, double nu1, double nu2,  */  
 /*                gsl_cdf_result * result, gsl_cdf_tail_t tail) */  
 /* { */  
 /*   int rc; */  
     
 /*   if( tail == GSL_CDF_LOWER ) */  
 /*     { */  
 /*       if( x <= 0.0 ) */  
 /*      { */  
 /*        result->val = 0.0; */  
 /*        result->err = 0.0; */  
 /*        return GSL_SUCCESS; */  
 /*      } */  
 /*       else */  
 /*      { */  
 /*        rc = gsl_F_ratio_cdf_P( x, nu1, nu2, result); */  
 /*      } */  
 /*     } */  
 /*   else if (tail == GSL_CDF_UPPER ) */  
 /*     { */  
 /*       if( x <= 0.0 ) */  
 /*      { */  
 /*        result->val = 1.0; */  
 /*        result->err = 0.0; */  
 /*        return GSL_SUCCESS; */  
 /*      } */  
 /*       else */  
 /*      { */  
 /*        rc = gsl_F_ratio_cdf_Q (x, nu1, nu2, result); */  
 /*      } */  
 /*     } */  
 /*   else */  
 /*     { */  
 /*       return GSL_EINVAL; */  
 /*     } */  
71    
72  /*   return rc; */    return P;
73  /* } */  }

Legend:
Removed from v.1.1.2.7  
changed lines
  Added in v.1.1.2.8

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26