lib/libm/nbdtrf.c

   1 /*                                                      nbdtrf.c
   2  *
   3  *      Negative binomial distribution
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * int k, n;
  10  * float p, y, nbdtrf();
  11  *
  12  * y = nbdtrf( 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  *        Relative error:
  42  * arithmetic   domain     # trials      peak         rms
  43  *    IEEE       0,100       5000       1.5e-4      1.9e-5
  44  *
  45  */
  46 \f/*                                                     nbdtrcf.c
  47  *
  48  *      Complemented negative binomial distribution
  49  *
  50  *
  51  *
  52  * SYNOPSIS:
  53  *
  54  * int k, n;
  55  * float p, y, nbdtrcf();
  56  *
  57  * y = nbdtrcf( k, n, p );
  58  *
  59  *
  60  *
  61  * DESCRIPTION:
  62  *
  63  * Returns the sum of the terms k+1 to infinity of the negative
  64  * binomial distribution:
  65  *
  66  *   inf
  67  *   --  ( n+j-1 )   n      j
  68  *   >   (       )  p  (1-p)
  69  *   --  (   j   )
  70  *  j=k+1
  71  *
  72  * The terms are not computed individually; instead the incomplete
  73  * beta integral is employed, according to the formula
  74  *
  75  * y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ).
  76  *
  77  * The arguments must be positive, with p ranging from 0 to 1.
  78  *
  79  *
  80  *
  81  * ACCURACY:
  82  *
  83  *        Relative error:
  84  * arithmetic   domain     # trials      peak         rms
  85  *    IEEE       0,100       5000       1.4e-4      2.0e-5
  86  *
  87  */
  88 \f
  89 /*
  90 Cephes Math Library Release 2.2:  July, 1992
  91 Copyright 1984, 1987 by Stephen L. Moshier
  92 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  93 */
  94
  95 #include "mconf.h"
  96
  97 #ifdef ANSIC
  98 float incbetf(float, float, float);
  99 #else
 100 float incbetf();
 101 #endif
 102
 103
 104 #ifdef ANSIC
 105 float nbdtrcf( int k, int n, float pp )
 106 #else
 107 float nbdtrcf( k, n, pp )
 108 int k, n;
 109 double pp;
 110 #endif
 111 {
 112 float dk, dn, p;
 113
 114 p = pp;
 115 if( (p < 0.0) || (p > 1.0) )
 116         goto domerr;
 117 if( k < 0 )
 118         {
 119 domerr:
 120         mtherr( "nbdtrf", DOMAIN );
 121         return( 0.0 );
 122         }
 123
 124 dk = k+1;
 125 dn = n;
 126 return( incbetf( dk, dn, 1.0 - p ) );
 127 }
 128
 129
 130
 131 #ifdef ANSIC
 132 float nbdtrf( int k, int n, float pp )
 133 #else
 134 float nbdtrf( k, n, pp )
 135 int k, n;
 136 double pp;
 137 #endif
 138 {
 139 float dk, dn, p;
 140
 141 p = pp;
 142 if( (p < 0.0) || (p > 1.0) )
 143         goto domerr;
 144 if( k < 0 )
 145         {
 146 domerr:
 147         mtherr( "nbdtrf", DOMAIN );
 148         return( 0.0 );
 149         }
 150 dk = k+1;
 151 dn = n;
 152 return( incbetf( dn, dk, p ) );
 153 }