lib/libm/bdtr.c

   1 /*                                                      bdtr.c
   2  *
   3  *      Binomial distribution
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * int k, n;
  10  * double p, y, bdtr();
  11  *
  12  * y = bdtr( k, n, p );
  13  *
  14  * DESCRIPTION:
  15  *
  16  * Returns the sum of the terms 0 through k of the Binomial
  17  * probability density:
  18  *
  19  *   k
  20  *   --  ( n )   j      n-j
  21  *   >   (   )  p  (1-p)
  22  *   --  ( j )
  23  *  j=0
  24  *
  25  * The terms are not summed directly; instead the incomplete
  26  * beta integral is employed, according to the formula
  27  *
  28  * y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ).
  29  *
  30  * The arguments must be positive, with p ranging from 0 to 1.
  31  *
  32  * ACCURACY:
  33  *
  34  * Tested at random points (a,b,p), with p between 0 and 1.
  35  *
  36  *               a,b                     Relative error:
  37  * arithmetic  domain     # trials      peak         rms
  38  *  For p between 0.001 and 1:
  39  *    IEEE     0,100       100000      4.3e-15     2.6e-16
  40  * See also incbet.c.
  41  *
  42  * ERROR MESSAGES:
  43  *
  44  *   message         condition      value returned
  45  * bdtr domain         k < 0            0.0
  46  *                     n < k
  47  *                     x < 0, x > 1
  48  */
  49 \f/*                                                     bdtrc()
  50  *
  51  *      Complemented binomial distribution
  52  *
  53  *
  54  *
  55  * SYNOPSIS:
  56  *
  57  * int k, n;
  58  * double p, y, bdtrc();
  59  *
  60  * y = bdtrc( k, n, p );
  61  *
  62  * DESCRIPTION:
  63  *
  64  * Returns the sum of the terms k+1 through n of the Binomial
  65  * probability density:
  66  *
  67  *   n
  68  *   --  ( n )   j      n-j
  69  *   >   (   )  p  (1-p)
  70  *   --  ( j )
  71  *  j=k+1
  72  *
  73  * The terms are not summed directly; instead the incomplete
  74  * beta integral is employed, according to the formula
  75  *
  76  * y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ).
  77  *
  78  * The arguments must be positive, with p ranging from 0 to 1.
  79  *
  80  * ACCURACY:
  81  *
  82  * Tested at random points (a,b,p).
  83  *
  84  *               a,b                     Relative error:
  85  * arithmetic  domain     # trials      peak         rms
  86  *  For p between 0.001 and 1:
  87  *    IEEE     0,100       100000      6.7e-15     8.2e-16
  88  *  For p between 0 and .001:
  89  *    IEEE     0,100       100000      1.5e-13     2.7e-15
  90  *
  91  * ERROR MESSAGES:
  92  *
  93  *   message         condition      value returned
  94  * bdtrc domain      x<0, x>1, n<k       0.0
  95  */
  96 \f/*                                                     bdtri()
  97  *
  98  *      Inverse binomial distribution
  99  *
 100  *
 101  *
 102  * SYNOPSIS:
 103  *
 104  * int k, n;
 105  * double p, y, bdtri();
 106  *
 107  * p = bdtr( k, n, y );
 108  *
 109  * DESCRIPTION:
 110  *
 111  * Finds the event probability p such that the sum of the
 112  * terms 0 through k of the Binomial probability density
 113  * is equal to the given cumulative probability y.
 114  *
 115  * This is accomplished using the inverse beta integral
 116  * function and the relation
 117  *
 118  * 1 - p = incbi( n-k, k+1, y ).
 119  *
 120  * ACCURACY:
 121  *
 122  * Tested at random points (a,b,p).
 123  *
 124  *               a,b                     Relative error:
 125  * arithmetic  domain     # trials      peak         rms
 126  *  For p between 0.001 and 1:
 127  *    IEEE     0,100       100000      2.3e-14     6.4e-16
 128  *    IEEE     0,10000     100000      6.6e-12     1.2e-13
 129  *  For p between 10^-6 and 0.001:
 130  *    IEEE     0,100       100000      2.0e-12     1.3e-14
 131  *    IEEE     0,10000     100000      1.5e-12     3.2e-14
 132  * See also incbi.c.
 133  *
 134  * ERROR MESSAGES:
 135  *
 136  *   message         condition      value returned
 137  * bdtri domain     k < 0, n <= k         0.0
 138  *                  x < 0, x > 1
 139  */
 140 \f
 141 /*                                                              bdtr() */
 142
 143
 144 /*
 145 Cephes Math Library Release 2.8:  June, 2000
 146 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
 147 */
 148
 149 #include "mconf.h"
 150 #ifdef ANSIPROT
 151 extern double incbet ( double, double, double );
 152 extern double incbi ( double, double, double );
 153 extern double pow ( double, double );
 154 extern double log1p ( double );
 155 extern double expm1 ( double );
 156 #else
 157 double incbet(), incbi(), pow(), log1p(), expm1();
 158 #endif
 159
 160 double bdtrc( k, n, p )
 161 int k, n;
 162 double p;
 163 {
 164 double dk, dn;
 165
 166 if( (p < 0.0) || (p > 1.0) )
 167         goto domerr;
 168 if( k < 0 )
 169         return( 1.0 );
 170
 171 if( n < k )
 172         {
 173 domerr:
 174         mtherr( "bdtrc", DOMAIN );
 175         return( 0.0 );
 176         }
 177
 178 if( k == n )
 179         return( 0.0 );
 180 dn = n - k;
 181 if( k == 0 )
 182         {
 183         if( p < .01 )
 184                 dk = -expm1( dn * log1p(-p) );
 185         else
 186                 dk = 1.0 - pow( 1.0-p, dn );
 187         }
 188 else
 189         {
 190         dk = k + 1;
 191         dk = incbet( dk, dn, p );
 192         }
 193 return( dk );
 194 }
 195
 196
 197
 198 double bdtr( k, n, p )
 199 int k, n;
 200 double p;
 201 {
 202 double dk, dn;
 203
 204 if( (p < 0.0) || (p > 1.0) )
 205         goto domerr;
 206 if( (k < 0) || (n < k) )
 207         {
 208 domerr:
 209         mtherr( "bdtr", DOMAIN );
 210         return( 0.0 );
 211         }
 212
 213 if( k == n )
 214         return( 1.0 );
 215
 216 dn = n - k;
 217 if( k == 0 )
 218         {
 219         dk = pow( 1.0-p, dn );
 220         }
 221 else
 222         {
 223         dk = k + 1;
 224         dk = incbet( dn, dk, 1.0 - p );
 225         }
 226 return( dk );
 227 }
 228
 229
 230 double bdtri( k, n, y )
 231 int k, n;
 232 double y;
 233 {
 234 double dk, dn, p;
 235
 236 if( (y < 0.0) || (y > 1.0) )
 237         goto domerr;
 238 if( (k < 0) || (n <= k) )
 239         {
 240 domerr:
 241         mtherr( "bdtri", DOMAIN );
 242         return( 0.0 );
 243         }
 244
 245 dn = n - k;
 246 if( k == 0 )
 247         {
 248         if( y > 0.8 )
 249                 p = -expm1( log1p(y-1.0) / dn );
 250         else
 251                 p = 1.0 - pow( y, 1.0/dn );
 252         }
 253 else
 254         {
 255         dk = k + 1;
 256         p = incbet( dn, dk, 0.5 );
 257         if( p > 0.5 )
 258                 p = incbi( dk, dn, 1.0-y );
 259         else
 260                 p = 1.0 - incbi( dn, dk, y );
 261         }
 262 return( p );
 263 }