123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222 |
- /* nbdtr.c
- *
- * Negative binomial distribution
- *
- *
- *
- * SYNOPSIS:
- *
- * int k, n;
- * double p, y, nbdtr();
- *
- * y = nbdtr( k, n, p );
- *
- * DESCRIPTION:
- *
- * Returns the sum of the terms 0 through k of the negative
- * binomial distribution:
- *
- * k
- * -- ( n+j-1 ) n j
- * > ( ) p (1-p)
- * -- ( j )
- * j=0
- *
- * In a sequence of Bernoulli trials, this is the probability
- * that k or fewer failures precede the nth success.
- *
- * The terms are not computed individually; instead the incomplete
- * beta integral is employed, according to the formula
- *
- * y = nbdtr( k, n, p ) = incbet( n, k+1, p ).
- *
- * The arguments must be positive, with p ranging from 0 to 1.
- *
- * ACCURACY:
- *
- * Tested at random points (a,b,p), with p between 0 and 1.
- *
- * a,b Relative error:
- * arithmetic domain # trials peak rms
- * IEEE 0,100 100000 1.7e-13 8.8e-15
- * See also incbet.c.
- *
- */
- /* nbdtrc.c
- *
- * Complemented negative binomial distribution
- *
- *
- *
- * SYNOPSIS:
- *
- * int k, n;
- * double p, y, nbdtrc();
- *
- * y = nbdtrc( k, n, p );
- *
- * DESCRIPTION:
- *
- * Returns the sum of the terms k+1 to infinity of the negative
- * binomial distribution:
- *
- * inf
- * -- ( n+j-1 ) n j
- * > ( ) p (1-p)
- * -- ( j )
- * j=k+1
- *
- * The terms are not computed individually; instead the incomplete
- * beta integral is employed, according to the formula
- *
- * y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ).
- *
- * The arguments must be positive, with p ranging from 0 to 1.
- *
- * ACCURACY:
- *
- * Tested at random points (a,b,p), with p between 0 and 1.
- *
- * a,b Relative error:
- * arithmetic domain # trials peak rms
- * IEEE 0,100 100000 1.7e-13 8.8e-15
- * See also incbet.c.
- */
- /* nbdtrc
- *
- * Complemented negative binomial distribution
- *
- *
- *
- * SYNOPSIS:
- *
- * int k, n;
- * double p, y, nbdtrc();
- *
- * y = nbdtrc( k, n, p );
- *
- * DESCRIPTION:
- *
- * Returns the sum of the terms k+1 to infinity of the negative
- * binomial distribution:
- *
- * inf
- * -- ( n+j-1 ) n j
- * > ( ) p (1-p)
- * -- ( j )
- * j=k+1
- *
- * The terms are not computed individually; instead the incomplete
- * beta integral is employed, according to the formula
- *
- * y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ).
- *
- * The arguments must be positive, with p ranging from 0 to 1.
- *
- * ACCURACY:
- *
- * See incbet.c.
- */
- /* nbdtri
- *
- * Functional inverse of negative binomial distribution
- *
- *
- *
- * SYNOPSIS:
- *
- * int k, n;
- * double p, y, nbdtri();
- *
- * p = nbdtri( k, n, y );
- *
- * DESCRIPTION:
- *
- * Finds the argument p such that nbdtr(k,n,p) is equal to y.
- *
- * ACCURACY:
- *
- * Tested at random points (a,b,y), with y between 0 and 1.
- *
- * a,b Relative error:
- * arithmetic domain # trials peak rms
- * IEEE 0,100 100000 1.5e-14 8.5e-16
- * See also incbi.c.
- */
- /*
- Cephes Math Library Release 2.8: June, 2000
- Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
- */
- #include <math.h>
- #ifdef ANSIPROT
- extern double incbet ( double, double, double );
- extern double incbi ( double, double, double );
- #else
- double incbet(), incbi();
- #endif
- double nbdtrc( k, n, p )
- int k, n;
- double p;
- {
- double dk, dn;
- if( (p < 0.0) || (p > 1.0) )
- goto domerr;
- if( k < 0 )
- {
- domerr:
- mtherr( "nbdtr", DOMAIN );
- return( 0.0 );
- }
- dk = k+1;
- dn = n;
- return( incbet( dk, dn, 1.0 - p ) );
- }
- double nbdtr( k, n, p )
- int k, n;
- double p;
- {
- double dk, dn;
- if( (p < 0.0) || (p > 1.0) )
- goto domerr;
- if( k < 0 )
- {
- domerr:
- mtherr( "nbdtr", DOMAIN );
- return( 0.0 );
- }
- dk = k+1;
- dn = n;
- return( incbet( dn, dk, p ) );
- }
- double nbdtri( k, n, p )
- int k, n;
- double p;
- {
- double dk, dn, w;
- if( (p < 0.0) || (p > 1.0) )
- goto domerr;
- if( k < 0 )
- {
- domerr:
- mtherr( "nbdtri", DOMAIN );
- return( 0.0 );
- }
- dk = k+1;
- dn = n;
- w = incbi( dn, dk, p );
- return( w );
- }
|