3 * Binomial distribution
10 * long double p, y, bdtrl();
12 * y = bdtrl( k, n, p );
18 * Returns the sum of the terms 0 through k of the Binomial
19 * probability density:
27 * The terms are not summed directly; instead the incomplete
28 * beta integral is employed, according to the formula
30 * y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ).
32 * The arguments must be positive, with p ranging from 0 to 1.
38 * Tested at random points (k,n,p) with a and b between 0
39 * and 10000 and p between 0 and 1.
41 * arithmetic domain # trials peak rms
42 * IEEE 0,10000 3000 1.6e-14 2.2e-15
46 * message condition value returned
47 * bdtrl domain k < 0 0.0
54 * Complemented binomial distribution
61 * long double p, y, bdtrcl();
63 * y = bdtrcl( k, n, p );
69 * Returns the sum of the terms k+1 through n of the Binomial
70 * probability density:
78 * The terms are not summed directly; instead the incomplete
79 * beta integral is employed, according to the formula
81 * y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ).
83 * The arguments must be positive, with p ranging from 0 to 1.
93 * message condition value returned
94 * bdtrcl domain x<0, x>1, n<k 0.0
98 * Inverse binomial distribution
105 * long double p, y, bdtril();
107 * p = bdtril( k, n, y );
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.
117 * This is accomplished using the inverse beta integral
118 * function and the relation
120 * 1 - p = incbi( n-k, k+1, y ).
125 * Tested at random k, n between 1 and 10000. The "domain" refers to p:
127 * arithmetic domain # trials peak rms
128 * IEEE 0,1 3500 2.0e-15 8.2e-17
132 * message condition value returned
133 * bdtril domain k < 0, n <= k 0.0
141 Cephes Math Library Release 2.3: March, 1995
142 Copyright 1984, 1995 by Stephen L. Moshier
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 );
153 long double incbetl(), incbil(), powl(), expm1l(), log1pl();
156 long double bdtrcl( k, n, p )
162 if( (p < 0.0L) || (p > 1.0L) )
170 mtherr( "bdtrcl", DOMAIN );
180 dk = -expm1l( dn * log1pl(-p) );
182 dk = 1.0L - powl( 1.0L-p, dn );
187 dk = incbetl( dk, dn, p );
194 long double bdtrl( k, n, p )
198 long double dk, dn, q;
200 if( (p < 0.0L) || (p > 1.0L) )
202 if( (k < 0) || (n < k) )
205 mtherr( "bdtrl", DOMAIN );
221 dk = incbetl( dn, dk, q );
227 long double bdtril( k, n, y )
231 long double dk, dn, p;
233 if( (y < 0.0L) || (y > 1.0L) )
235 if( (k < 0) || (n <= k) )
238 mtherr( "bdtril", DOMAIN );
246 p = -expm1l( log1pl(y-1.0L) / dn );
248 p = 1.0L - powl( y, 1.0L/dn );
253 p = incbetl( dn, dk, y );
255 p = incbil( dk, dn, 1.0L-y );
257 p = 1.0 - incbil( dn, dk, y );