lib/libm/fdtr.c

   1 /*                                                      fdtr.c
   2  *
   3  *      F distribution
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * int df1, df2;
  10  * double x, y, fdtr();
  11  *
  12  * y = fdtr( df1, df2, x );
  13  *
  14  * DESCRIPTION:
  15  *
  16  * Returns the area from zero to x under the F density
  17  * function (also known as Snedcor's density or the
  18  * variance ratio density).  This is the density
  19  * of x = (u1/df1)/(u2/df2), where u1 and u2 are random
  20  * variables having Chi square distributions with df1
  21  * and df2 degrees of freedom, respectively.
  22  *
  23  * The incomplete beta integral is used, according to the
  24  * formula
  25  *
  26  *      P(x) = incbet( df1/2, df2/2, (df1*x/(df2 + df1*x) ).
  27  *
  28  *
  29  * The arguments a and b are greater than zero, and x is
  30  * nonnegative.
  31  *
  32  * ACCURACY:
  33  *
  34  * Tested at random points (a,b,x).
  35  *
  36  *                x     a,b                     Relative error:
  37  * arithmetic  domain  domain     # trials      peak         rms
  38  *    IEEE      0,1    0,100       100000      9.8e-15     1.7e-15
  39  *    IEEE      1,5    0,100       100000      6.5e-15     3.5e-16
  40  *    IEEE      0,1    1,10000     100000      2.2e-11     3.3e-12
  41  *    IEEE      1,5    1,10000     100000      1.1e-11     1.7e-13
  42  * See also incbet.c.
  43  *
  44  *
  45  * ERROR MESSAGES:
  46  *
  47  *   message         condition      value returned
  48  * fdtr domain     a<0, b<0, x<0         0.0
  49  *
  50  */
  51 \f/*                                                     fdtrc()
  52  *
  53  *      Complemented F distribution
  54  *
  55  *
  56  *
  57  * SYNOPSIS:
  58  *
  59  * int df1, df2;
  60  * double x, y, fdtrc();
  61  *
  62  * y = fdtrc( df1, df2, x );
  63  *
  64  * DESCRIPTION:
  65  *
  66  * Returns the area from x to infinity under the F density
  67  * function (also known as Snedcor's density or the
  68  * variance ratio density).
  69  *
  70  *
  71  *                      inf.
  72  *                       -
  73  *              1       | |  a-1      b-1
  74  * 1-P(x)  =  ------    |   t    (1-t)    dt
  75  *            B(a,b)  | |
  76  *                     -
  77  *                      x
  78  *
  79  *
  80  * The incomplete beta integral is used, according to the
  81  * formula
  82  *
  83  *      P(x) = incbet( df2/2, df1/2, (df2/(df2 + df1*x) ).
  84  *
  85  *
  86  * ACCURACY:
  87  *
  88  * Tested at random points (a,b,x) in the indicated intervals.
  89  *                x     a,b                     Relative error:
  90  * arithmetic  domain  domain     # trials      peak         rms
  91  *    IEEE      0,1    1,100       100000      3.7e-14     5.9e-16
  92  *    IEEE      1,5    1,100       100000      8.0e-15     1.6e-15
  93  *    IEEE      0,1    1,10000     100000      1.8e-11     3.5e-13
  94  *    IEEE      1,5    1,10000     100000      2.0e-11     3.0e-12
  95  * See also incbet.c.
  96  *
  97  * ERROR MESSAGES:
  98  *
  99  *   message         condition      value returned
 100  * fdtrc domain    a<0, b<0, x<0         0.0
 101  *
 102  */
 103 \f/*                                                     fdtri()
 104  *
 105  *      Inverse of complemented F distribution
 106  *
 107  *
 108  *
 109  * SYNOPSIS:
 110  *
 111  * int df1, df2;
 112  * double x, p, fdtri();
 113  *
 114  * x = fdtri( df1, df2, p );
 115  *
 116  * DESCRIPTION:
 117  *
 118  * Finds the F density argument x such that the integral
 119  * from x to infinity of the F density is equal to the
 120  * given probability p.
 121  *
 122  * This is accomplished using the inverse beta integral
 123  * function and the relations
 124  *
 125  *      z = incbi( df2/2, df1/2, p )
 126  *      x = df2 (1-z) / (df1 z).
 127  *
 128  * Note: the following relations hold for the inverse of
 129  * the uncomplemented F distribution:
 130  *
 131  *      z = incbi( df1/2, df2/2, p )
 132  *      x = df2 z / (df1 (1-z)).
 133  *
 134  * ACCURACY:
 135  *
 136  * Tested at random points (a,b,p).
 137  *
 138  *              a,b                     Relative error:
 139  * arithmetic  domain     # trials      peak         rms
 140  *  For p between .001 and 1:
 141  *    IEEE     1,100       100000      8.3e-15     4.7e-16
 142  *    IEEE     1,10000     100000      2.1e-11     1.4e-13
 143  *  For p between 10^-6 and 10^-3:
 144  *    IEEE     1,100        50000      1.3e-12     8.4e-15
 145  *    IEEE     1,10000      50000      3.0e-12     4.8e-14
 146  * See also fdtrc.c.
 147  *
 148  * ERROR MESSAGES:
 149  *
 150  *   message         condition      value returned
 151  * fdtri domain   p <= 0 or p > 1       0.0
 152  *                     v < 1
 153  *
 154  */
 155 \f
 156
 157 /*
 158 Cephes Math Library Release 2.8:  June, 2000
 159 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
 160 */
 161
 162
 163 #include "mconf.h"
 164 #ifdef ANSIPROT
 165 extern double incbet ( double, double, double );
 166 extern double incbi ( double, double, double );
 167 #else
 168 double incbet(), incbi();
 169 #endif
 170
 171 double fdtrc( ia, ib, x )
 172 int ia, ib;
 173 double x;
 174 {
 175 double a, b, w;
 176
 177 if( (ia < 1) || (ib < 1) || (x < 0.0) )
 178         {
 179         mtherr( "fdtrc", DOMAIN );
 180         return( 0.0 );
 181         }
 182 a = ia;
 183 b = ib;
 184 w = b / (b + a * x);
 185 return( incbet( 0.5*b, 0.5*a, w ) );
 186 }
 187
 188
 189
 190 double fdtr( ia, ib, x )
 191 int ia, ib;
 192 double x;
 193 {
 194 double a, b, w;
 195
 196 if( (ia < 1) || (ib < 1) || (x < 0.0) )
 197         {
 198         mtherr( "fdtr", DOMAIN );
 199         return( 0.0 );
 200         }
 201 a = ia;
 202 b = ib;
 203 w = a * x;
 204 w = w / (b + w);
 205 return( incbet(0.5*a, 0.5*b, w) );
 206 }
 207
 208
 209 double fdtri( ia, ib, y )
 210 int ia, ib;
 211 double y;
 212 {
 213 double a, b, w, x;
 214
 215 if( (ia < 1) || (ib < 1) || (y <= 0.0) || (y > 1.0) )
 216         {
 217         mtherr( "fdtri", DOMAIN );
 218         return( 0.0 );
 219         }
 220 a = ia;
 221 b = ib;
 222 /* Compute probability for x = 0.5.  */
 223 w = incbet( 0.5*b, 0.5*a, 0.5 );
 224 /* If that is greater than y, then the solution w < .5.
 225    Otherwise, solve at 1-y to remove cancellation in (b - b*w).  */
 226 if( w > y || y < 0.001)
 227         {
 228         w = incbi( 0.5*b, 0.5*a, y );
 229         x = (b - b*w)/(a*w);
 230         }
 231 else
 232         {
 233         w = incbi( 0.5*a, 0.5*b, 1.0-y );
 234         x = b*w/(a*(1.0-w));
 235         }
 236 return(x);
 237 }