summaryrefslogtreecommitdiff
path: root/libm/float/bdtrf.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/float/bdtrf.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/float/bdtrf.c')
-rw-r--r--libm/float/bdtrf.c247
1 files changed, 247 insertions, 0 deletions
diff --git a/libm/float/bdtrf.c b/libm/float/bdtrf.c
new file mode 100644
index 000000000..e063f1c77
--- /dev/null
+++ b/libm/float/bdtrf.c
@@ -0,0 +1,247 @@
+/* bdtrf.c
+ *
+ * Binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * float p, y, bdtrf();
+ *
+ * y = bdtrf( 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:
+ *
+ * Relative error (p varies from 0 to 1):
+ * arithmetic domain # trials peak rms
+ * IEEE 0,100 2000 6.9e-5 1.1e-5
+ *
+ * ERROR MESSAGES:
+ *
+ * message condition value returned
+ * bdtrf domain k < 0 0.0
+ * n < k
+ * x < 0, x > 1
+ *
+ */
+ /* bdtrcf()
+ *
+ * Complemented binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * float p, y, bdtrcf();
+ *
+ * y = bdtrcf( 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:
+ *
+ * Relative error (p varies from 0 to 1):
+ * arithmetic domain # trials peak rms
+ * IEEE 0,100 2000 6.0e-5 1.2e-5
+ *
+ * ERROR MESSAGES:
+ *
+ * message condition value returned
+ * bdtrcf domain x<0, x>1, n<k 0.0
+ */
+ /* bdtrif()
+ *
+ * Inverse binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * float p, y, bdtrif();
+ *
+ * p = bdtrf( 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:
+ *
+ * Relative error (p varies from 0 to 1):
+ * arithmetic domain # trials peak rms
+ * IEEE 0,100 2000 3.5e-5 3.3e-6
+ *
+ * ERROR MESSAGES:
+ *
+ * message condition value returned
+ * bdtrif domain k < 0, n <= k 0.0
+ * x < 0, x > 1
+ *
+ */
+
+/* bdtr() */
+
+
+/*
+Cephes Math Library Release 2.2: July, 1992
+Copyright 1984, 1987, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+#include <math.h>
+
+#ifdef ANSIC
+float incbetf(float, float, float), powf(float, float);
+float incbif( float, float, float );
+#else
+float incbetf(), powf(), incbif();
+#endif
+
+float bdtrcf( int k, int n, float pp )
+{
+float p, dk, dn;
+
+p = pp;
+if( (p < 0.0) || (p > 1.0) )
+ goto domerr;
+if( k < 0 )
+ return( 1.0 );
+
+if( n < k )
+ {
+domerr:
+ mtherr( "bdtrcf", DOMAIN );
+ return( 0.0 );
+ }
+
+if( k == n )
+ return( 0.0 );
+dn = n - k;
+if( k == 0 )
+ {
+ dk = 1.0 - powf( 1.0-p, dn );
+ }
+else
+ {
+ dk = k + 1;
+ dk = incbetf( dk, dn, p );
+ }
+return( dk );
+}
+
+
+
+float bdtrf( int k, int n, float pp )
+{
+float p, dk, dn;
+
+p = pp;
+if( (p < 0.0) || (p > 1.0) )
+ goto domerr;
+if( (k < 0) || (n < k) )
+ {
+domerr:
+ mtherr( "bdtrf", DOMAIN );
+ return( 0.0 );
+ }
+
+if( k == n )
+ return( 1.0 );
+
+dn = n - k;
+if( k == 0 )
+ {
+ dk = powf( 1.0-p, dn );
+ }
+else
+ {
+ dk = k + 1;
+ dk = incbetf( dn, dk, 1.0 - p );
+ }
+return( dk );
+}
+
+
+float bdtrif( int k, int n, float yy )
+{
+float y, dk, dn, p;
+
+y = yy;
+if( (y < 0.0) || (y > 1.0) )
+ goto domerr;
+if( (k < 0) || (n <= k) )
+ {
+domerr:
+ mtherr( "bdtrif", DOMAIN );
+ return( 0.0 );
+ }
+
+dn = n - k;
+if( k == 0 )
+ {
+ p = 1.0 - powf( y, 1.0/dn );
+ }
+else
+ {
+ dk = k + 1;
+ p = 1.0 - incbif( dn, dk, y );
+ }
+return( p );
+}