libm/ldouble/nbdtrl.c

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