libm/ldouble/bdtrl.c

   1 /*                                                      bdtrl.c
   2  *
   3  *      Binomial distribution
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * int k, n;
  10  * long double p, y, bdtrl();
  11  *
  12  * y = bdtrl( 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  * Tested at random points (k,n,p) with a and b between 0
  39  * and 10000 and p between 0 and 1.
  40  *    Relative error:
  41  * arithmetic   domain     # trials      peak         rms
  42  *    IEEE      0,10000      3000       1.6e-14     2.2e-15
  43  *
  44  * ERROR MESSAGES:
  45  *
  46  *   message         condition      value returned
  47  * bdtrl domain        k < 0            0.0
  48  *                     n < k
  49  *                     x < 0, x > 1
  50  *
  51  */
  52 \f/*                                                     bdtrcl()
  53  *
  54  *      Complemented binomial distribution
  55  *
  56  *
  57  *
  58  * SYNOPSIS:
  59  *
  60  * int k, n;
  61  * long double p, y, bdtrcl();
  62  *
  63  * y = bdtrcl( k, n, p );
  64  *
  65  *
  66  *
  67  * DESCRIPTION:
  68  *
  69  * Returns the sum of the terms k+1 through n of the Binomial
  70  * probability density:
  71  *
  72  *   n
  73  *   --  ( n )   j      n-j
  74  *   >   (   )  p  (1-p)
  75  *   --  ( j )
  76  *  j=k+1
  77  *
  78  * The terms are not summed directly; instead the incomplete
  79  * beta integral is employed, according to the formula
  80  *
  81  * y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ).
  82  *
  83  * The arguments must be positive, with p ranging from 0 to 1.
  84  *
  85  *
  86  *
  87  * ACCURACY:
  88  *
  89  * See incbet.c.
  90  *
  91  * ERROR MESSAGES:
  92  *
  93  *   message         condition      value returned
  94  * bdtrcl domain     x<0, x>1, n<k       0.0
  95  */
  96 \f/*                                                     bdtril()
  97  *
  98  *      Inverse binomial distribution
  99  *
 100  *
 101  *
 102  * SYNOPSIS:
 103  *
 104  * int k, n;
 105  * long double p, y, bdtril();
 106  *
 107  * p = bdtril( 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  * ACCURACY:
 123  *
 124  * See incbi.c.
 125  * Tested at random k, n between 1 and 10000.  The "domain" refers to p:
 126  *                      Relative error:
 127  * arithmetic   domain     # trials      peak         rms
 128  *    IEEE       0,1        3500       2.0e-15     8.2e-17
 129  *
 130  * ERROR MESSAGES:
 131  *
 132  *   message         condition      value returned
 133  * bdtril domain     k < 0, n <= k         0.0
 134  *                  x < 0, x > 1
 135  */
 136 \f
 137 /*                                                              bdtr() */
 138
 139
 140 /*
 141 Cephes Math Library Release 2.3:  March, 1995
 142 Copyright 1984, 1995 by Stephen L. Moshier
 143 */
 144
 145 #include <math.h>
 146 #ifdef ANSIPROT
 147 extern long double incbetl ( long double, long double, long double );
 148 extern long double incbil ( long double, long double, long double );
 149 extern long double powl ( long double, long double );
 150 extern long double expm1l ( long double );
 151 extern long double log1pl ( long double );
 152 #else
 153 long double incbetl(), incbil(), powl(), expm1l(), log1pl();
 154 #endif
 155
 156 long double bdtrcl( k, n, p )
 157 int k, n;
 158 long double p;
 159 {
 160 long double dk, dn;
 161
 162 if( (p < 0.0L) || (p > 1.0L) )
 163         goto domerr;
 164 if( k < 0 )
 165         return( 1.0L );
 166
 167 if( n < k )
 168         {
 169 domerr:
 170         mtherr( "bdtrcl", DOMAIN );
 171         return( 0.0L );
 172         }
 173
 174 if( k == n )
 175         return( 0.0L );
 176 dn = n - k;
 177 if( k == 0 )
 178         {
 179         if( p < .01L )
 180                 dk = -expm1l( dn * log1pl(-p) );
 181         else
 182                 dk = 1.0L - powl( 1.0L-p, dn );
 183         }
 184 else
 185         {
 186         dk = k + 1;
 187         dk = incbetl( dk, dn, p );
 188         }
 189 return( dk );
 190 }
 191
 192
 193
 194 long double bdtrl( k, n, p )
 195 int k, n;
 196 long double p;
 197 {
 198 long double dk, dn, q;
 199
 200 if( (p < 0.0L) || (p > 1.0L) )
 201         goto domerr;
 202 if( (k < 0) || (n < k) )
 203         {
 204 domerr:
 205         mtherr( "bdtrl", DOMAIN );
 206         return( 0.0L );
 207         }
 208
 209 if( k == n )
 210         return( 1.0L );
 211
 212 q = 1.0L - p;
 213 dn = n - k;
 214 if( k == 0 )
 215         {
 216         dk = powl( q, dn );
 217         }
 218 else
 219         {
 220         dk = k + 1;
 221         dk = incbetl( dn, dk, q );
 222         }
 223 return( dk );
 224 }
 225
 226
 227 long double bdtril( k, n, y )
 228 int k, n;
 229 long double y;
 230 {
 231 long double dk, dn, p;
 232
 233 if( (y < 0.0L) || (y > 1.0L) )
 234         goto domerr;
 235 if( (k < 0) || (n <= k) )
 236         {
 237 domerr:
 238         mtherr( "bdtril", DOMAIN );
 239         return( 0.0L );
 240         }
 241
 242 dn = n - k;
 243 if( k == 0 )
 244         {
 245         if( y > 0.8L )
 246                 p = -expm1l( log1pl(y-1.0L) / dn );
 247         else
 248                 p = 1.0L - powl( y, 1.0L/dn );
 249         }
 250 else
 251         {
 252         dk = k + 1;
 253         p = incbetl( dn, dk, y );
 254         if( p > 0.5 )
 255                 p = incbil( dk, dn, 1.0L-y );
 256         else
 257                 p = 1.0 - incbil( dn, dk, y );
 258         }
 259 return( p );
 260 }