3 #define TWOVRPI 0.636619772367581343
4 #define HALF_PI 1.5707963268
5 #define EPSILON .000001
7 extern void betabase(double *, double *, double *,int *, int *,double *);
8 extern void normbase(double *,double *);
10 extern double ppnd(), ppbeta();
12 /* CACM Algorithm 395, by G. W. Hill */
15 studentbase(double *x
, double *df
, double *cdf
)
17 double t
, y
, b
, a
, z
, j
, n
;
25 if (n
> floor(n
) || (n
>= 20 && t
< n
) || (n
> 20)) {
26 if (n
< 2.0 && n
!= 1.0) {
27 /* beta integral aproximation for small df */
28 double da
= 0.5, db
= 0.5 * n
, dx
, dp
;
29 int ia
= 0, ib
= floor(db
);
31 dx
= db
/ (db
+ da
* t
);
32 betabase(&dx
, &db
, &da
, &ib
, &ia
, &dp
);
33 *cdf
= (*x
>= 0) ? 1.0 - .5 * dp
: .5 * dp
;
35 /* asymptotic series for large or non-integer df */
42 y
= (((((-0.4 * y
- 3.3) * y
- 24.0) * y
- 85.5)
43 / (0.8 * y
* y
+ 100.0 + b
) + y
+ 3.0) / b
+ 1.0) * sqrt(y
);
47 if (*x
> 0.0) *cdf
= 1.0 - *cdf
;
50 /* nested summation of cosine series */
51 if (n
< 20 && t
< 4.0) {
57 for(j
= 2; fabs(a
- z
) > EPSILON
; j
+= 2.0) {
59 y
= (y
* (j
- 1)) / (b
* j
);
66 for(n
= n
- 2.0; n
> 1.0; n
-= 2.0) a
= ((n
- 1.0) / (b
* n
)) * a
+ y
;
67 a
= (fabs(n
) < EPSILON
) ? a
/sqrt(b
) : TWOVRPI
* (atan(y
) + a
/ b
);
70 *cdf
= 1.0 - 0.5 * *cdf
;
77 /* CACM Algorithm 396, by G. W. Hill */
79 double ppstudent(pp
, n
, ifault
)
83 double sq
, p
, a
, b
, c
, d
, x
, y
;
85 /* convert to double upper tailed probability */
86 p
= (pp
< 0.5) ? 2.0 * pp
: 2.0 * (1.0 - pp
);
89 if (n
== 1) sq
= tan(HALF_PI
* (1.0 - p
));
90 else if (n
== 2.0) sq
= sqrt(2.0 / (p
* (2.0 - p
)) - 2.0);
92 sq
= ppbeta(p
, 0.5 * n
, 0.5, ifault
);
93 if (sq
!= 0.0) sq
= sqrt((n
/ sq
) - n
);
99 c
= ((20700.0 * a
/ b
- 98.0) * a
- 16) * a
+ 96.36;
100 d
= ((94.5 / (b
+ c
) - 3.0) / b
+ 1.0) * sqrt(a
* HALF_PI
) * n
;
104 /* asymptotic inverse expansion about normal */
105 x
= ppnd(0.5 * p
, ifault
);
107 if (n
< 5) c
= c
+ 0.3 * (n
- 4.5) * (x
+ 0.6);
108 c
= (((0.05 * d
* x
- 5.0) * x
- 7.0) * x
- 2.0) * x
+ b
+ c
;
109 y
= (((((0.4 * y
+ 6.3) * y
+ 36.0) * y
+ 94.5) / c
- y
- 3.0) / b
112 y
= (y
> .002) ? exp(y
) - 1.0 : 0.5 * y
* y
+ y
;
115 y
= ((1.0 / (((n
+ 6.0) / (n
* y
) - 0.089 * d
- 0.822)
116 * (n
+ 2.0) * 3.0) + 0.5 / (n
+ 4.0)) * y
- 1.0)
117 * (n
+ 1.0) / (n
+ 2.0) + 1.0 / y
;
122 /* decode based on tail */
123 if (pp
< 0.5) sq
= -sq
;