3 #define TWOVRPI 0.636619772367581343
4 #define HALF_PI 1.5707963268
5 #define EPSILON .000001
7 extern double ppnd(), ppbeta();
9 /* CACM Algorithm 395, by G. W. Hill */
11 studentbase(x
, df
, cdf
)
14 double t
, y
, b
, a
, z
, j
, n
;
22 if (n
> floor(n
) || (n
>= 20 && t
< n
) || (n
> 20)) {
23 if (n
< 2.0 && n
!= 1.0) {
24 /* beta integral aproximation for small df */
25 double da
= 0.5, db
= 0.5 * n
, dx
, dp
;
26 int ia
= 0, ib
= floor(db
);
28 dx
= db
/ (db
+ da
* t
);
29 betabase(&dx
, &db
, &da
, &ib
, &ia
, &dp
);
30 *cdf
= (*x
>= 0) ? 1.0 - .5 * dp
: .5 * dp
;
33 /* asymptotic series for large or non-integer df */
34 if(y
> EPSILON
) y
= log(b
);
38 y
= (((((-0.4 * y
- 3.3) * y
- 24.0) * y
- 85.5)
39 / (0.8 * y
* y
+ 100.0 + b
) + y
+ 3.0) / b
+ 1.0) * sqrt(y
);
43 if (*x
> 0.0) *cdf
= 1.0 - *cdf
;
47 /* nested summation of cosine series */
48 if (n
< 20 && t
< 4.0) {
55 for(j
= 2; fabs(a
- z
) > EPSILON
; j
+= 2.0) {
57 y
= (y
* (j
- 1)) / (b
* j
);
64 for(n
= n
- 2.0; n
> 1.0; n
-= 2.0) a
= ((n
- 1.0) / (b
* n
)) * a
+ y
;
65 a
= (fabs(n
) < EPSILON
) ? a
/sqrt(b
) : TWOVRPI
* (atan(y
) + a
/ b
);
67 if(*x
> 0.0) *cdf
= 1.0 - 0.5 * *cdf
;
68 else *cdf
= 0.5 * *cdf
;
72 /* CACM Algorithm 396, by G. W. Hill */
74 double ppstudent(pp
, n
, ifault
)
78 double sq
, p
, a
, b
, c
, d
, x
, y
;
80 /* convert to double upper tailed probability */
81 p
= (pp
< 0.5) ? 2.0 * pp
: 2.0 * (1.0 - pp
);
84 if (n
== 1) sq
= tan(HALF_PI
* (1.0 - p
));
85 else if (n
== 2.0) sq
= sqrt(2.0 / (p
* (2.0 - p
)) - 2.0);
87 sq
= ppbeta(p
, 0.5 * n
, 0.5, ifault
);
88 if (sq
!= 0.0) sq
= sqrt((n
/ sq
) - n
);
94 c
= ((20700.0 * a
/ b
- 98.0) * a
- 16) * a
+ 96.36;
95 d
= ((94.5 / (b
+ c
) - 3.0) / b
+ 1.0) * sqrt(a
* HALF_PI
) * n
;
99 /* asymptotic inverse expansion about normal */
100 x
= ppnd(0.5 * p
, ifault
);
102 if (n
< 5) c
= c
+ 0.3 * (n
- 4.5) * (x
+ 0.6);
103 c
= (((0.05 * d
* x
- 5.0) * x
- 7.0) * x
- 2.0) * x
+ b
+ c
;
104 y
= (((((0.4 * y
+ 6.3) * y
+ 36.0) * y
+ 94.5) / c
- y
- 3.0) / b
107 y
= (y
> .002) ? exp(y
) - 1.0 : 0.5 * y
* y
+ y
;
110 y
= ((1.0 / (((n
+ 6.0) / (n
* y
) - 0.089 * d
- 0.822)
111 * (n
+ 2.0) * 3.0) + 0.5 / (n
+ 4.0)) * y
- 1.0)
112 * (n
+ 1.0) / (n
+ 2.0) + 1.0 / y
;
117 /* decode based on tail */
118 if (pp
< 0.5) sq
= -sq
;