]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:02:42 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | FUNCTION STUDIN(Q,N) | |
11 | C THIS ROUTINE COMPUTES THE INVERSE OF THE DISTRIBUTION | |
12 | C FUNCTION OF THE STUDENT DISTRIBUTION WITH N DEGREES OF FREEDOM. | |
13 | C Q=PROBABILITY AT WHICH THE FUNCTION IS INVERTED, | |
14 | C N=INTEGER GREATER THAN ZERO. | |
15 | C THE ROUTINE WAS WRITTEN BY G.W.HILL IN ALGOL | |
16 | C C.A.C.M. ALGORITHM 396 | |
17 | C AT LEAST 6 SIGNIFICANT FIGURES ARE CORRECT. | |
18 | DATA HP/1.5707963268/ | |
19 | IF(N. LT. 1) GO TO 10 | |
20 | RL=1. | |
21 | IF(Q .GE. 0.5) GO TO 1 | |
22 | RL=-1. | |
23 | P=2.*Q | |
24 | GO TO 2 | |
25 | 1 P=2.*(1.-Q) | |
26 | 2 IF(P .LE. 0. .OR. P .GT. 1.) GO TO 20 | |
27 | IF(N .GT. 1) GO TO 3 | |
28 | PP=COS(HP*P) | |
29 | STUDIN=PP/SQRT(1.-PP*PP)*RL | |
30 | RETURN | |
31 | 3 IF(N .GT. 2) GO TO 4 | |
32 | STUDIN=SQRT(2./(P*(2.-P))-2.)*RL | |
33 | RETURN | |
34 | 4 RN=N | |
35 | A=1./(RN-0.5) | |
36 | B=48./(A*A) | |
37 | C=((20700.*A/B-98.)*A-16.)*A+96.36 | |
38 | D=((94.5/(B+C)-3.)/B+1.)*SQRT(A*HP)*RN | |
39 | X=D*P | |
40 | Y=X**(2./RN) | |
41 | IF(Y .LE. 0.05+A) GO TO 5 | |
42 | PP=0.5*P | |
43 | X=GAUSIN(PP) | |
44 | Y=X*X | |
45 | IF(N .GE. 5) C=C+0.3*(RN-4.5)*(X+0.6) | |
46 | C=(((0.05*D*X-5.)*X-7.)*X-2.)*X+B+C | |
47 | Y=(((((0.4*Y+6.3)*Y+36.)*Y+94.5)/C-Y-3.)/B+1.)*X | |
48 | Y=A*Y*Y | |
49 | IF(Y .LE. 0.002) Y=0.5*Y*Y+Y | |
50 | IF(Y .GT. 0.002) Y=EXP(Y)-1. | |
51 | GO TO 6 | |
52 | 5 Y=((1./(((RN+6.)/(RN*Y)-0.089*D-0.822)*(RN+2.)*3.)+0.5/(RN+4.))*Y- | |
53 | 11.)*(RN+1.)/(RN+2.)+1./Y | |
54 | 6 STUDIN=SQRT(RN*Y)*RL | |
55 | RETURN | |
56 | 10 WRITE(6,7) N | |
57 | 20 WRITE(6,8) Q | |
58 | STOP | |
59 | 7 FORMAT(/10X,'DEGREE OF FREEDOM N=',I5,' IN STUDIN ILLEGAL'/) | |
60 | 8 FORMAT(/10X,'ARGUMENT Q=',E15.5,' IN STUDIN ILLEGAL'/) | |
61 | END |