summaryrefslogtreecommitdiff
path: root/libm/double/stdtr.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/double/stdtr.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/double/stdtr.c')
-rw-r--r--libm/double/stdtr.c225
1 files changed, 225 insertions, 0 deletions
diff --git a/libm/double/stdtr.c b/libm/double/stdtr.c
new file mode 100644
index 000000000..743e01704
--- /dev/null
+++ b/libm/double/stdtr.c
@@ -0,0 +1,225 @@
+/* stdtr.c
+ *
+ * Student's t distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double t, stdtr();
+ * short k;
+ *
+ * y = stdtr( k, t );
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Computes the integral from minus infinity to t of the Student
+ * t distribution with integer k > 0 degrees of freedom:
+ *
+ * t
+ * -
+ * | |
+ * - | 2 -(k+1)/2
+ * | ( (k+1)/2 ) | ( x )
+ * ---------------------- | ( 1 + --- ) dx
+ * - | ( k )
+ * sqrt( k pi ) | ( k/2 ) |
+ * | |
+ * -
+ * -inf.
+ *
+ * Relation to incomplete beta integral:
+ *
+ * 1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z )
+ * where
+ * z = k/(k + t**2).
+ *
+ * For t < -2, this is the method of computation. For higher t,
+ * a direct method is derived from integration by parts.
+ * Since the function is symmetric about t=0, the area under the
+ * right tail of the density is found by calling the function
+ * with -t instead of t.
+ *
+ * ACCURACY:
+ *
+ * Tested at random 1 <= k <= 25. The "domain" refers to t.
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE -100,-2 50000 5.9e-15 1.4e-15
+ * IEEE -2,100 500000 2.7e-15 4.9e-17
+ */
+
+/* stdtri.c
+ *
+ * Functional inverse of Student's t distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double p, t, stdtri();
+ * int k;
+ *
+ * t = stdtri( k, p );
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Given probability p, finds the argument t such that stdtr(k,t)
+ * is equal to p.
+ *
+ * ACCURACY:
+ *
+ * Tested at random 1 <= k <= 100. The "domain" refers to p:
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE .001,.999 25000 5.7e-15 8.0e-16
+ * IEEE 10^-6,.001 25000 2.0e-12 2.9e-14
+ */
+
+
+/*
+Cephes Math Library Release 2.8: June, 2000
+Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
+*/
+
+#include <math.h>
+
+extern double PI, MACHEP, MAXNUM;
+#ifdef ANSIPROT
+extern double sqrt ( double );
+extern double atan ( double );
+extern double incbet ( double, double, double );
+extern double incbi ( double, double, double );
+extern double fabs ( double );
+#else
+double sqrt(), atan(), incbet(), incbi(), fabs();
+#endif
+
+double stdtr( k, t )
+int k;
+double t;
+{
+double x, rk, z, f, tz, p, xsqk;
+int j;
+
+if( k <= 0 )
+ {
+ mtherr( "stdtr", DOMAIN );
+ return(0.0);
+ }
+
+if( t == 0 )
+ return( 0.5 );
+
+if( t < -2.0 )
+ {
+ rk = k;
+ z = rk / (rk + t * t);
+ p = 0.5 * incbet( 0.5*rk, 0.5, z );
+ return( p );
+ }
+
+/* compute integral from -t to + t */
+
+if( t < 0 )
+ x = -t;
+else
+ x = t;
+
+rk = k; /* degrees of freedom */
+z = 1.0 + ( x * x )/rk;
+
+/* test if k is odd or even */
+if( (k & 1) != 0)
+ {
+
+ /* computation for odd k */
+
+ xsqk = x/sqrt(rk);
+ p = atan( xsqk );
+ if( k > 1 )
+ {
+ f = 1.0;
+ tz = 1.0;
+ j = 3;
+ while( (j<=(k-2)) && ( (tz/f) > MACHEP ) )
+ {
+ tz *= (j-1)/( z * j );
+ f += tz;
+ j += 2;
+ }
+ p += f * xsqk/z;
+ }
+ p *= 2.0/PI;
+ }
+
+
+else
+ {
+
+ /* computation for even k */
+
+ f = 1.0;
+ tz = 1.0;
+ j = 2;
+
+ while( ( j <= (k-2) ) && ( (tz/f) > MACHEP ) )
+ {
+ tz *= (j - 1)/( z * j );
+ f += tz;
+ j += 2;
+ }
+ p = f * x/sqrt(z*rk);
+ }
+
+/* common exit */
+
+
+if( t < 0 )
+ p = -p; /* note destruction of relative accuracy */
+
+ p = 0.5 + 0.5 * p;
+return(p);
+}
+
+double stdtri( k, p )
+int k;
+double p;
+{
+double t, rk, z;
+int rflg;
+
+if( k <= 0 || p <= 0.0 || p >= 1.0 )
+ {
+ mtherr( "stdtri", DOMAIN );
+ return(0.0);
+ }
+
+rk = k;
+
+if( p > 0.25 && p < 0.75 )
+ {
+ if( p == 0.5 )
+ return( 0.0 );
+ z = 1.0 - 2.0 * p;
+ z = incbi( 0.5, 0.5*rk, fabs(z) );
+ t = sqrt( rk*z/(1.0-z) );
+ if( p < 0.5 )
+ t = -t;
+ return( t );
+ }
+rflg = -1;
+if( p >= 0.5)
+ {
+ p = 1.0 - p;
+ rflg = 1;
+ }
+z = incbi( 0.5*rk, 0.5, 2.0*p );
+
+if( MAXNUM * z < rk )
+ return(rflg* MAXNUM);
+t = sqrt( rk/z - rk );
+return( rflg * t );
+}