lib/libm/bdtrf.c

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