lib/libm/nbdtr.c

   1 /*                                                      nbdtr.c
   2  *
   3  *      Negative binomial distribution
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * int k, n;
  10  * double p, y, nbdtr();
  11  *
  12  * y = nbdtr( k, n, p );
  13  *
  14  * DESCRIPTION:
  15  *
  16  * Returns the sum of the terms 0 through k of the negative
  17  * binomial distribution:
  18  *
  19  *   k
  20  *   --  ( n+j-1 )   n      j
  21  *   >   (       )  p  (1-p)
  22  *   --  (   j   )
  23  *  j=0
  24  *
  25  * In a sequence of Bernoulli trials, this is the probability
  26  * that k or fewer failures precede the nth success.
  27  *
  28  * The terms are not computed individually; instead the incomplete
  29  * beta integral is employed, according to the formula
  30  *
  31  * y = nbdtr( k, n, p ) = incbet( n, k+1, p ).
  32  *
  33  * The arguments must be positive, with p ranging from 0 to 1.
  34  *
  35  * ACCURACY:
  36  *
  37  * Tested at random points (a,b,p), with p between 0 and 1.
  38  *
  39  *               a,b                     Relative error:
  40  * arithmetic  domain     # trials      peak         rms
  41  *    IEEE     0,100       100000      1.7e-13     8.8e-15
  42  * See also incbet.c.
  43  *
  44  */
  45 \f/*                                                     nbdtrc.c
  46  *
  47  *      Complemented negative binomial distribution
  48  *
  49  *
  50  *
  51  * SYNOPSIS:
  52  *
  53  * int k, n;
  54  * double p, y, nbdtrc();
  55  *
  56  * y = nbdtrc( k, n, p );
  57  *
  58  * DESCRIPTION:
  59  *
  60  * Returns the sum of the terms k+1 to infinity of the negative
  61  * binomial distribution:
  62  *
  63  *   inf
  64  *   --  ( n+j-1 )   n      j
  65  *   >   (       )  p  (1-p)
  66  *   --  (   j   )
  67  *  j=k+1
  68  *
  69  * The terms are not computed individually; instead the incomplete
  70  * beta integral is employed, according to the formula
  71  *
  72  * y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ).
  73  *
  74  * The arguments must be positive, with p ranging from 0 to 1.
  75  *
  76  * ACCURACY:
  77  *
  78  * Tested at random points (a,b,p), with p between 0 and 1.
  79  *
  80  *               a,b                     Relative error:
  81  * arithmetic  domain     # trials      peak         rms
  82  *    IEEE     0,100       100000      1.7e-13     8.8e-15
  83  * See also incbet.c.
  84  */
  85 \f
  86 /*                                                      nbdtrc
  87  *
  88  *      Complemented negative binomial distribution
  89  *
  90  *
  91  *
  92  * SYNOPSIS:
  93  *
  94  * int k, n;
  95  * double p, y, nbdtrc();
  96  *
  97  * y = nbdtrc( k, n, p );
  98  *
  99  * DESCRIPTION:
 100  *
 101  * Returns the sum of the terms k+1 to infinity of the negative
 102  * binomial distribution:
 103  *
 104  *   inf
 105  *   --  ( n+j-1 )   n      j
 106  *   >   (       )  p  (1-p)
 107  *   --  (   j   )
 108  *  j=k+1
 109  *
 110  * The terms are not computed individually; instead the incomplete
 111  * beta integral is employed, according to the formula
 112  *
 113  * y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ).
 114  *
 115  * The arguments must be positive, with p ranging from 0 to 1.
 116  *
 117  * ACCURACY:
 118  *
 119  * See incbet.c.
 120  */
 121 \f/*                                                     nbdtri
 122  *
 123  *      Functional inverse of negative binomial distribution
 124  *
 125  *
 126  *
 127  * SYNOPSIS:
 128  *
 129  * int k, n;
 130  * double p, y, nbdtri();
 131  *
 132  * p = nbdtri( k, n, y );
 133  *
 134  * DESCRIPTION:
 135  *
 136  * Finds the argument p such that nbdtr(k,n,p) is equal to y.
 137  *
 138  * ACCURACY:
 139  *
 140  * Tested at random points (a,b,y), with y between 0 and 1.
 141  *
 142  *               a,b                     Relative error:
 143  * arithmetic  domain     # trials      peak         rms
 144  *    IEEE     0,100       100000      1.5e-14     8.5e-16
 145  * See also incbi.c.
 146  */
 147 \f
 148 /*
 149 Cephes Math Library Release 2.8:  June, 2000
 150 Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
 151 */
 152
 153 #include "mconf.h"
 154 #ifdef ANSIPROT
 155 extern double incbet ( double, double, double );
 156 extern double incbi ( double, double, double );
 157 #else
 158 double incbet(), incbi();
 159 #endif
 160
 161 double nbdtrc( k, n, p )
 162 int k, n;
 163 double p;
 164 {
 165 double dk, dn;
 166
 167 if( (p < 0.0) || (p > 1.0) )
 168         goto domerr;
 169 if( k < 0 )
 170         {
 171 domerr:
 172         mtherr( "nbdtr", DOMAIN );
 173         return( 0.0 );
 174         }
 175
 176 dk = k+1;
 177 dn = n;
 178 return( incbet( dk, dn, 1.0 - p ) );
 179 }
 180
 181
 182
 183 double nbdtr( k, n, p )
 184 int k, n;
 185 double p;
 186 {
 187 double dk, dn;
 188
 189 if( (p < 0.0) || (p > 1.0) )
 190         goto domerr;
 191 if( k < 0 )
 192         {
 193 domerr:
 194         mtherr( "nbdtr", DOMAIN );
 195         return( 0.0 );
 196         }
 197 dk = k+1;
 198 dn = n;
 199 return( incbet( dn, dk, p ) );
 200 }
 201
 202
 203
 204 double nbdtri( k, n, p )
 205 int k, n;
 206 double p;
 207 {
 208 double dk, dn, w;
 209
 210 if( (p < 0.0) || (p > 1.0) )
 211         goto domerr;
 212 if( k < 0 )
 213         {
 214 domerr:
 215         mtherr( "nbdtri", DOMAIN );
 216         return( 0.0 );
 217         }
 218 dk = k+1;
 219 dn = n;
 220 w = incbi( dn, dk, p );
 221 return( w );
 222 }