summaryrefslogtreecommitdiff
path: root/libm/ldouble/bdtrl.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/ldouble/bdtrl.c')
-rw-r--r--libm/ldouble/bdtrl.c260
1 files changed, 260 insertions, 0 deletions
diff --git a/libm/ldouble/bdtrl.c b/libm/ldouble/bdtrl.c
new file mode 100644
index 000000000..aca9577d1
--- /dev/null
+++ b/libm/ldouble/bdtrl.c
@@ -0,0 +1,260 @@
+/* bdtrl.c
+ *
+ * Binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * long double p, y, bdtrl();
+ *
+ * y = bdtrl( k, n, p );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the sum of the terms 0 through k of the Binomial
+ * probability density:
+ *
+ * k
+ * -- ( n ) j n-j
+ * > ( ) p (1-p)
+ * -- ( j )
+ * j=0
+ *
+ * The terms are not summed directly; instead the incomplete
+ * beta integral is employed, according to the formula
+ *
+ * y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ).
+ *
+ * The arguments must be positive, with p ranging from 0 to 1.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Tested at random points (k,n,p) with a and b between 0
+ * and 10000 and p between 0 and 1.
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0,10000 3000 1.6e-14 2.2e-15
+ *
+ * ERROR MESSAGES:
+ *
+ * message condition value returned
+ * bdtrl domain k < 0 0.0
+ * n < k
+ * x < 0, x > 1
+ *
+ */
+ /* bdtrcl()
+ *
+ * Complemented binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * long double p, y, bdtrcl();
+ *
+ * y = bdtrcl( k, n, p );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the sum of the terms k+1 through n of the Binomial
+ * probability density:
+ *
+ * n
+ * -- ( n ) j n-j
+ * > ( ) p (1-p)
+ * -- ( j )
+ * j=k+1
+ *
+ * The terms are not summed directly; instead the incomplete
+ * beta integral is employed, according to the formula
+ *
+ * y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ).
+ *
+ * The arguments must be positive, with p ranging from 0 to 1.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * See incbet.c.
+ *
+ * ERROR MESSAGES:
+ *
+ * message condition value returned
+ * bdtrcl domain x<0, x>1, n<k 0.0
+ */
+ /* bdtril()
+ *
+ * Inverse binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * long double p, y, bdtril();
+ *
+ * p = bdtril( k, n, y );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Finds the event probability p such that the sum of the
+ * terms 0 through k of the Binomial probability density
+ * is equal to the given cumulative probability y.
+ *
+ * This is accomplished using the inverse beta integral
+ * function and the relation
+ *
+ * 1 - p = incbi( n-k, k+1, y ).
+ *
+ * ACCURACY:
+ *
+ * See incbi.c.
+ * Tested at random k, n between 1 and 10000. The "domain" refers to p:
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0,1 3500 2.0e-15 8.2e-17
+ *
+ * ERROR MESSAGES:
+ *
+ * message condition value returned
+ * bdtril domain k < 0, n <= k 0.0
+ * x < 0, x > 1
+ */
+
+/* bdtr() */
+
+
+/*
+Cephes Math Library Release 2.3: March, 1995
+Copyright 1984, 1995 by Stephen L. Moshier
+*/
+
+#include <math.h>
+#ifdef ANSIPROT
+extern long double incbetl ( long double, long double, long double );
+extern long double incbil ( long double, long double, long double );
+extern long double powl ( long double, long double );
+extern long double expm1l ( long double );
+extern long double log1pl ( long double );
+#else
+long double incbetl(), incbil(), powl(), expm1l(), log1pl();
+#endif
+
+long double bdtrcl( k, n, p )
+int k, n;
+long double p;
+{
+long double dk, dn;
+
+if( (p < 0.0L) || (p > 1.0L) )
+ goto domerr;
+if( k < 0 )
+ return( 1.0L );
+
+if( n < k )
+ {
+domerr:
+ mtherr( "bdtrcl", DOMAIN );
+ return( 0.0L );
+ }
+
+if( k == n )
+ return( 0.0L );
+dn = n - k;
+if( k == 0 )
+ {
+ if( p < .01L )
+ dk = -expm1l( dn * log1pl(-p) );
+ else
+ dk = 1.0L - powl( 1.0L-p, dn );
+ }
+else
+ {
+ dk = k + 1;
+ dk = incbetl( dk, dn, p );
+ }
+return( dk );
+}
+
+
+
+long double bdtrl( k, n, p )
+int k, n;
+long double p;
+{
+long double dk, dn, q;
+
+if( (p < 0.0L) || (p > 1.0L) )
+ goto domerr;
+if( (k < 0) || (n < k) )
+ {
+domerr:
+ mtherr( "bdtrl", DOMAIN );
+ return( 0.0L );
+ }
+
+if( k == n )
+ return( 1.0L );
+
+q = 1.0L - p;
+dn = n - k;
+if( k == 0 )
+ {
+ dk = powl( q, dn );
+ }
+else
+ {
+ dk = k + 1;
+ dk = incbetl( dn, dk, q );
+ }
+return( dk );
+}
+
+
+long double bdtril( k, n, y )
+int k, n;
+long double y;
+{
+long double dk, dn, p;
+
+if( (y < 0.0L) || (y > 1.0L) )
+ goto domerr;
+if( (k < 0) || (n <= k) )
+ {
+domerr:
+ mtherr( "bdtril", DOMAIN );
+ return( 0.0L );
+ }
+
+dn = n - k;
+if( k == 0 )
+ {
+ if( y > 0.8L )
+ p = -expm1l( log1pl(y-1.0L) / dn );
+ else
+ p = 1.0L - powl( y, 1.0L/dn );
+ }
+else
+ {
+ dk = k + 1;
+ p = incbetl( dn, dk, y );
+ if( p > 0.5 )
+ p = incbil( dk, dn, 1.0L-y );
+ else
+ p = 1.0 - incbil( dn, dk, y );
+ }
+return( p );
+}