Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/oneaov/oneaov.for
There is 1 other file named oneaov.for in the archive. Click here to see a list.
C WESTERN MICHIGAN UNIVERSITY
C ONEAOV.F4 (FILENAME ON LIBRARY DECTAPE)
C ONEAOV, 1.9.5 (CALLING NAME, SUBLST. NO.)
C ONE WAY ANALYSIS OF VARIANCE (UNBALANCED)
C ONEAOV WAS PROGRAMMED BY R.R. BARR III, LATER MODIFIED BY M.T.
C O'BRYAN
C CONSULTANT: DR. MIKE STOLINE
C LIBRARY DECTAPE PROGRAMS USED: USAGE.MAC
C APLIB PROGS. USED: IOB, GETFOR, FISHER
C FORWMU PROGS. USED: TTYPTY, ALLCOR, DEVCH, EXISTS, PRINTS
C INTERNAL SUBR. USED: CUNO, CUCS, MNL
C METHOD FROM B. J. WINER, "STATISTICAL PRINCIPLES IN
C EXPERIEMENTAL METHODS"
C ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C---------------ID(16) LIMITS USER SPEC. HEADING FOR OUTPUT TO 80 CH.,
C--------------- IFMT(48) LIMITS OBJ. TIME FORMAT TO 240 CH.
DIMENSION ID(16),IFMT(48)
C---------------XBR(20) (SEE ST. 312 AND 312+1)
DIMENSION XBR(20)
C---------------S(1) IS USED IN CALL MNL
DIMENSION S(1)
C---------------SUBR. IOB IN APLIB SHARES COMMON /IOBLK/
C--------------- AND COMMON /IOBLKA/
COMMON/IOBLK/IDLG,IRSP,IIN,IOUT,IDVI,IDVO,ICODE,IBNK,NAMI(2)
COMMON/IOBLKA/NAMO(2),IPJ,IPG,NCOPYS,ITYCH
IIN=4
IOUT=6
ITYCH=7
CALL ERRSET(1000)
IDLG=-1
IRSP=-4
WRITE(IDLG,100)
100 FORMAT(///,' ---WMU ONE-WAY ANALYSIS OF VARIANCE---',/)
C CALL USAGE('ONEAOV')
C
C PARAMETER INPUT SECTION
C
C RELATIVE LINES .+1,3,4,5 PLUS 208+2,210+3 AND 321+0,326+1
C CONTAIN MODIFICATIONS TO ALLOW SAME OPTION FOR TTY:
C---------------TTYPTY RETURNS ZERO - TTY, MINUS 0NE - BATCH JOB
CALL TTYPTY(ICODE)
C---------------1 MEANS OUTPUT? PRINTS, O MEANS INTPUT? PRINTS.
C--------------- IDLG, IRSP, IIN, IOUT, IDVI, IDVO, ICODE, ARE
C--------------- INPUT THRU COMMON /IOBLK/. ITYCH IS RETURNED
C---------------THRU COMMON /IOBLKA/
CALL IOB(1)
102 CALL IOB(0)
NBK=0
C---------------ITYCH RETURNED THRU COMMON /IOBLKA/ BY SUBR. IOB(APLIB)
IF(IMTH.NE.1.OR.IIN.NE.ITYCH)GO TO 103
WRITE(IDLG,105)
105 FORMAT(' ?CANNOT USE "SAME" OPTION FOR METHOD ONE AND',
1 ' TTY: DATA',/)
GO TO 102
103 WRITE(IDLG,104)
104 FORMAT(' ENTER OUTPUT IDENTIFICATION, IF DESIRED',/)
READ(IRSP,106)ID
106 FORMAT(16A5)
IF(IIN.EQ.ITYCH)GO TO 300
C---------------IFMT, ISTD RETURNED. OTHER ARGS. ARE INPUT.
C--------------- 48=NO. OF OBJ. TIME FORMAT WORDS (3 LINES).
C--------------- 2 MEANS F - TYPE FORMAT ONLY.
CALL GETFOR(-1,-4,IFMT,ISTD,48,2)
IF(ISTD.EQ.1)IFMT(1)='(10F)'
107 WRITE(IDLG,108)
108 FORMAT(' HOW MANY VARIABLES? ',$)
READ(IRSP,110)NVAR
NVARA=NVAR
110 FORMAT(10I)
IF(NVAR.GT.0)GO TO 112
WRITE(IDLG,114)
IF(ICODE.NE.0)CALL EXIT
GO TO 107
112 WRITE(IDLG,124)
124 FORMAT(' ENTER METHOD OF INPUT(1 OR 2) ',$)
READ(IRSP,110)IMTH
IF(IMTH.NE.1.OR.IIN.NE.ITYCH)GO TO 126
WRITE(IDLG,105)
GO TO 112
126 IF(IMTH.NE.2.OR.NVARA.NE.1)GO TO 127
WRITE(IDLG,114)
GO TO 102
127 IF(IMTH.LE.2)IF(IMTH-1)128,200,300
128 WRITE(IDLG,114)
114 FORMAT(' ?RESPONSE OUT OF BOUNDS',/)
IF(ICODE.NE.0)CALL EXIT
GO TO 112
130 WRITE(IDLG,132)MAX
132 FORMAT(' ?PROBLEM TOO LARGE(SIZE=',I5,')',/)
GO TO 107
C
C METHOD ONE FOR INPUT
C
200 WRITE(IDLG,202)
202 FORMAT(' HOW MANY GROUPS? ',$)
READ(IRSP,110)NGR
IF(NGR.GT.0)GO TO 320
WRITE(IDLG,114)
IF(ICODE.NE.0)CALL EXIT
GO TO 200
C
C METHOD TWO FOR INPUT
C
300 IF(IIN.EQ.ITYCH)NVAR=NVARA
WRITE(IDLG,302)
302 FORMAT(' WHICH VARIABLE IS THE BREAKDOWN VARIABLE? ',$)
READ(IRSP,110)NBK
IF(NBK.GE.1.AND.NBK.LE.NVAR)GO TO 310
WRITE(IDLG,114)
IF(ICODE.NE.0)CALL EXIT
GO TO 300
310 WRITE(IDLG,312)
312 FORMAT(' ENTER BREAKDOWN LIMITS(MAX. 20)',/)
READ(IRSP,121)(XBR(I),I=1,20)
120 FORMAT(10F)
121 FORMAT(20F)
DO 314 NGR=20,1,-1
IF(XBR(NGR))316,314,320
314 CONTINUE
316 WRITE(IDLG,318)
318 FORMAT(' ?BREAKDOWN LIMIT ERROR')
GO TO 310
320 CONTINUE
C
C CORE ALLOCATION SECTION
C
MAX=3*NGR*NVAR+2*NGR+3*NVAR
CALL ALLCOR(MAX,IERR,ISTT,S)
IF(IERR.LT.0)GO TO 130
IJMAX=NGR*NVAR
ISTX2=ISTT+IJMAX
ISTN=ISTX2+IJMAX
ISTNEL=ISTN+IJMAX
ISTAMN=ISTNEL+NGR
ISTX=ISTAMN+NGR
ISTSYM=ISTX+NVAR
ISTVAR=ISTSYM+NVAR
C FOR EXPANSION
ISTXXX=ISTVAR+NGR
CALL MNL(S(ISTT),S(ISTX2),S(ISTN),S(ISTNEL),S(ISTAMN),
1 S(ISTX),S(ISTSYM),S(ISTVAR),NVAR,NGR,IMTH,IFMT,IDVI,IDVO,
2 ID,NBK,XBR,IIN,ITYCH)
GO TO 102
END
C
C
C
C---------------NVAR, NGR, IMTH, IFMT, IDVI, IDVO, ID, NBK,
C--------------- XBR, IIN, ITYCH ARE INPUT. THE OTHER ARGS.
C--------------- ARE SPACES RESERVED BY DYN. ALLOC.
SUBROUTINE MNL(T,X2,N,NEL,AMN,X,SYM,VAR,NVAR,NGR,IMTH,IFMT,
1 IDVI,IDVO,ID,NBK,XBR,IIN,ITYCH)
DIMENSION T(1),X2(1),N(1),NEL(1),AMN(1),X(1),SYM(1),VAR(1),
1 IFMT(48),ID(16),XBR(1)
IDLG=-1
IRSP=-4
WRITE(IDLG,116)
116 FORMAT(' ARE ANY SYMBOLS TO BE USED FOR MISSING',
1 ' DATA?(YES OR NO) ',$)
ISYM=0
READ(IRSP,114)ANS
114 FORMAT(A3)
IF(ANS.NE.'YES')GO TO 122
ISYM=1
WRITE(IDLG,118)
118 FORMAT(' TYPE MISSING DATA SYMBOL FOR EACH VARIABLE(10 PER',
1 ' LINE)',/)
READ(IRSP,120)(SYM(I),I=1,NVAR)
120 FORMAT(10F)
122 IF(IMTH-1)203,203,321
C
C METHOD ONE FOR INPUT
C
203 WRITE(IDLG,204)
204 FORMAT(' ENTER NUMBER OF ELEMENTS PER GROUP(10 PER LINE)',/)
READ(IRSP,205)(NEL(I),I=1,NGR)
205 FORMAT(10I)
IF(IDVI.NE.'TTY')WRITE(IDLG,206)IDVI
206 FORMAT(' DATA IS BEING READ FROM ',A5,/)
DO 212 I=1,NGR
DO 208 J=1,NVAR
IJ=(I-1)*NVAR+J
T(IJ)=0
X2(IJ)=0
N(IJ)=0
208 CONTINUE
IF(NEL(I).LE.0)GO TO 212
IF(IIN.NE.ITYCH.AND.IDVI.EQ.'TTY')WRITE(IDLG,210)I
210 FORMAT(' ENTER DATA FOR GROUP 'I2,/)
DO 212 K=1,NEL(I)
IF(IIN.EQ.ITYCH)GO TO 213
READ(IIN,IFMT,END=250)(X(J),J=1,NVAR)
IF(IDVI.EQ.'TTY')WRITE(ITYCH)(X(J),J=1,NVAR)
GO TO 214
213 READ(IIN,END=250)(X(J),J=1,NVAR)
214 DO 212 J=1,NVAR
IJ=(I-1)*NVAR+J
IF(ISYM.EQ.1.AND.X(J).EQ.SYM(J))GO TO 211
T(IJ)=T(IJ)+X(J)
X2(IJ)=X2(IJ)+X(J)*X(J)
N(IJ)=N(IJ)+1
211 CONTINUE
212 CONTINUE
GO TO 336
250 WRITE(IDLG,252)
252 FORMAT(' ?NOT ENOUGH DATA',/)
RETURN
C
C METHOD TWO FOR INPUT
C
321 IF(IIN.NE.ITYCH.AND.IDVI.EQ.'TTY')WRITE(IDLG,324)
IF(IDVI.NE.'TTY')WRITE(IDLG,206)IDVI
320 DO 322 I=1,NGR
DO 322 J=1,NVAR
IJ=(I-1)*NVAR+J
T(IJ)=0
X2(IJ)=0
N(IJ)=0
322 CONTINUE
IRJ=0
324 FORMAT(' ENTER DATA',/)
326 IF(IIN.EQ.ITYCH)GO TO 327
READ(IIN,IFMT,END=332)(X(I),I=1,NVAR)
IF(IDVI.EQ.'TTY')WRITE(ITYCH)(X(I),I=1,NVAR)
GO TO 329
327 READ(IIN,END=332)(X(I),I=1,NVAR)
329 DO 330 I=1,NGR
IF(ISYM.EQ.1.AND.X(NBK).EQ.SYM(NBK))GO TO 326
IF(X(NBK).GT.XBR(I))GO TO 330
DO 328 JX=1,NVAR
IF(JX.EQ.NBK)GO TO 328
IF(ISYM.EQ.1.AND.X(JX).EQ.SYM(JX))GO TO 328
J=JX
IF(JX.GT.NBK)J=JX-1
IJ=(I-1)*(NVAR-1)+J
T(IJ)=T(IJ)+X(JX)
X2(IJ)=X2(IJ)+X(JX)*X(JX)
N(IJ)=N(IJ)+1
328 CONTINUE
GO TO 326
330 CONTINUE
IRJ=IRJ+1
GO TO 326
332 NVAR=NVAR-1
WRITE(IDLG,334)IRJ
334 FORMAT(' NUMBER OF REJECTED SAMPLES=',I4,/)
336 WRITE(IDLG,338)
338 FORMAT(' WOULD YOU LIKE TWO-SAMPLE T''S?(YES OR NO) ',$)
READ(IRSP,114)SCHANS
IF(SCHANS.NE.'YES')GO TO 346
IBB=1
340 WRITE(IDLG,342)NGR
342 FORMAT(' ENTER A 1 TO USE THE POOLED MEAN SQUARE ',
1 'FOR JUST THE TWO GROUPS,',/,' ENTER A 2 TO USE',
2 ' THE POOLED MEAN SQUARE FOR ALL ',I4,' GROUPS',/)
READ(IRSP,205)IBB
IF(IBB.LT.1.OR.IBB.GT.2)GO TO 340
346 CONTINUE
C
C FIRST WRITE TO OUTPUT FILE
C
WRITE(6,350)ID
350 FORMAT('1',16A5)
DIMENSION DAY(2)
CALL DATE(DAY)
CALL TIME(HM)
WRITE(6,351)HM,DAY
IF(IDVO.NE.'TTY')WRITE(IDLG,351)HM,DAY
351 FORMAT(/,1X,A5,1X,2A5)
DO 590 J=1,NVAR
C
C MEANS SECTION
C
IHEAD=0
V=0
VI=0
BART=0
NOBART=0
SPSQ=0
NGRA=0
JZ=J
IF(J.GE.NBK.AND.NBK.NE.0)JZ=J+1
DO 414 I=1,NGR
IJ=(I-1)*NVAR+J
IF(N(IJ).LT.1)GO TO 414
NGRA=NGRA+1
IF(IHEAD.EQ.0)WRITE(6,352)JZ
352 FORMAT(///,1X,T21,'*** VARIABLE ',I3,' ***'//,T9,'GROUP',T22,
1 'SIZE',T30,'MEANS',T41,'VARIANCE',T54,'STD DEV',/)
IHEAD=1
IF(N(IJ).GT.1)GO TO 410
AMN(I)=T(IJ)
NGRA=NGRA-1
NOBART=1
WRITE(6,409)I,N(IJ),AMN(I)
409 FORMAT(1X,T10,I3,T20,I5,T26,F10.3,T40,10('-'),T52,10('-'),/)
GO TO 414
410 AMN(I)=T(IJ)/N(IJ)
VAR(I)=(X2(IJ)-T(IJ)*T(IJ)/N(IJ))/(N(IJ)-1)
C**********************************************************
C MTO: 26-JUL-76. ROUNDOFF ERRORS MADE VAR NEGATIVE.
C
IF (VAR(I).LT.0.0) VAR(I)=0.0
C
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
STD=SQRT(VAR(I))
IF (VAR(I).GT.1E-9) GOTO 356
NOBART=1
NGRA=NGRA-1
GO TO 413
356 V=V+N(IJ)-1
VI=VI+1./(N(IJ)-1)
BART=BART+(N(IJ)-1)*ALOG(VAR(I))
SPSQ=SPSQ+VAR(I)*(N(IJ)-1)
413 WRITE(6,412)I,N(IJ),AMN(I),VAR(I),STD
412 FORMAT(1X,T10,I3,T20,I5,T26,F10.3,T40,F10.4,T52,F10.4,/)
414 CONTINUE
IF(NGRA.GT.0.OR.NOBART.EQ.1)GO TO 440
WRITE(6,416)JZ
416 FORMAT(///,' ALL GROUPS EMPTY FOR VARIABLE ',I2,/)
GO TO 590
C
C BARTLETT TEST SECTION
C
440 IF(NOBART.NE.1)GO TO 444
WRITE(6,442)
442 FORMAT(5X,'BARTLETT''S TEST STATISTIC IS NOT POSSIBLE',/,
1 5X,'SINCE AT LEAST ONE SAMPLE SIZE IS 1 OR VARIANCE IS 0.',/)
GO TO 480
444 SPSQ=SPSQ/V
IDF=NGRA-1
IF(IDF.GT.0)GO TO 448
WRITE(6,446)
446 FORMAT(5X,'NO BARTLETT TEST IF NUMBER OF GROUPS IS LESS',
1 ' THAN 2',/)
GO TO 480
448 BART=(V*ALOG(SPSQ)-BART)/(1.+(VI-1./V)/(3.*IDF))
WRITE(6,450)NGRA,IDF,BART
450 FORMAT(//,5X,'NUMBER OF GROUPS=',T36,I4,T47,'DF= ',I4,/,
1 5X,'BARTLETT''S STATISTIC=',T35,F9.3)
CALL CUCS(IDF,BART)
480 CONTINUE
C
C ANOVA TABLE SECTION
C
KDF=NGRA-1
IF(KDF.GT.0)GO TO 486
WRITE(6,482)KDF
482 FORMAT(/,' ANALYSIS OF VARIANCE TABLE DELETED BECAUSE',/,
1 ' DEGREES OF FREEDOM FOR GROUPS IS ',I2,/)
GO TO 590
486 G=0
NT=0
THREE=0
TWO=0
DO 490 I=1,NGR
IJ=(I-1)*NVAR+J
G=G+T(IJ)
IF(N(IJ).LT.1)GO TO 490
NT=NT+N(IJ)
TWO=TWO+X2(IJ)
IF(N(IJ).EQ.0)GO TO 490
THREE=THREE+T(IJ)*T(IJ)/N(IJ)
490 CONTINUE
GBAR=G/NT
ONE=G*G/NT
SSTR=THREE-ONE
SSE=TWO-THREE
SST=TWO-ONE
KDFER=NT-NGRA
KDFT=NT-1
TRMS=SSTR/KDF
ERMS=SSE/KDFER
TRF=TRMS/ERMS
PROB=FISHER(KDF,KDFER,TRF)
WRITE(6,500)
500 FORMAT(///,17X,'ANALYSIS OF VARIANCE TABLE',//,
1 ' SOURCE OF VAR.',T18,'D.F.',T27,'SUM OF SQ.',T40,'MEAN SQ.',
2 T58,'F',T66,'PROB',/)
WRITE(6,502)KDF,SSTR,TRMS,TRF,PROB
502 FORMAT(' GROUPS',T18,I5,T26,F11.3,T39,F9.3,T51,
1 F9.3,T65,F6.3/)
WRITE(6,504)KDFER,SSE,ERMS
504 FORMAT(' ERROR',T18,I5,T26,F11.3,T39,F9.3,/)
WRITE(6,506)KDFT,SST
506 FORMAT(' TOTAL',T18,I5,T26,F11.3,/)
C
C TWO-SAMPLE T ANALYSIS
C
IF(SCHANS.NE.'YES')GO TO 590
ALPHA=.95
PCT=ALPHA*100.
ALPHA=1.-ALPHA
WRITE(6,520)PCT
520 FORMAT(///,21X,'TWO-SAMPLE T ANALYSIS'///,T17,'MEAN',T30,'T',
1 T59,F3.0,'% C. I.',/,1X,'GRP-GRP',T13,'DIFFERANCE',
2 T28,'VALUE',T39,'PROB',T47,'DF',T51,'LOWER LIMIT,',
3 'UPPER LIMIT',/)
CALL RELEAS(-1)
IFLG=0
DO 562 I=1,NGR-1
DO 562 K=I+1,NGR
IJ=(I-1)*NVAR+J
KJ=(K-1)*NVAR+J
FLAG=' '
IF(N(IJ)*N(KJ).EQ.0)GO TO 562
IF((N(IJ).EQ.1).OR.(N(KJ).EQ.1))FLAG='#'
IF(FLAG.EQ.'#')IFLG=1
DMN=AMN(I)-AMN(K)
RKPI=(N(IJ)+N(KJ))/FLOAT(N(IJ)*N(KJ))
IF(IBB.EQ.2)GO TO 521
KP=N(IJ)+N(KJ)-2
C OMMIT GRP-GRP IF KP=0
IF(KP.LE.0)GO TO 562
XTMP=RKPI*((N(IJ)-1)*VAR(I)+(N(KJ)-1)*VAR(K))
IF(XTMP.EQ.0)GO TO 562
XTMP=XTMP/KP
SSH=SQRT(XTMP)
GO TO 522
521 SSH=SQRT(ERMS*RKPI)
KP=KDFER
522 TSH=TPCT(ALPHA,KP)
AL=DMN-TSH*SSH
U=DMN+TSH*SSH
TII=DMN/SSH
TPR=FISHER(1,KP,TII*TII)
WRITE(6,524)I,K,FLAG,DMN,TII,TPR,KP,AL,U
524 FORMAT(1X,I3,'-',I3,T10,A1,T13,F8.2,T25,
1 F8.3,T37,F6.3,T45,I4,T52,'(',F9.2,',',F9.2,')',/)
562 CONTINUE
IF(IFLG.EQ.1)WRITE(6,560)
560 FORMAT(' # - INDICATES AT LEAST ONE GROUP OF SIZE - 1',/)
IF(IBB.EQ.1)WRITE(6,564)
564 FORMAT(' USING POOLED MEAN SQUARE FOR JUST TWO GROUPS',/)
IF(IBB.EQ.2)WRITE(6,566)NGR
566 FORMAT(' USING POOLED MEAN SQUARE FOR ALL ',I4,' GROUPS',/)
590 CONTINUE
WRITE(IDLG,600)
600 FORMAT(//)
RETURN
END
C
C CUMULATIVE CHI SQUARE SUBROUTINE
C
C SOURCE-CHARLES A. NAGY
C
C---------------BOTH PARAM. INPUT. THIS IS SAME AS CHISQR IN APLIB.
SUBROUTINE CUCS(K,X)
IF(K.LE.0.OR.K.GT.100)GO TO 14
IF(X.GT.141)GO TO 14
Y=0
IF(X.LE.0)GO TO 13
IF(X.LE.0)GO TO 13
IF(K.GE.4)GO TO 4
GO TO (1,2,3),K
1 P=SQRT(X)
CALL CUNO(P,S)
IF(P.LT.0)RETURN
Y=2.*S-1.
GO TO 13
2 Y=1.-(1./EXP(X/2.))
GO TO 13
3 P=SQRT(X)
CALL CUNO(P,S)
IF(P.LT.0)RETURN
Y=(2.*S-1.)-P/(1.25331414*EXP(X/2.))
GO TO 13
4 M=K/2
IF(K-2*M)5,6,5
5 P=SQRT(X)
CALL CUNO(P,S)
IF(P.LT.0)RETURN
Y=2.*S-1.
S=X/2.
C=1./(.62665707*P*EXP(S))
P=S
T=.5
GO TO 7
6 C=1./EXP(X/2.)
Y=1.-C
S=0
P=1
T=0
7 DO 8 I=1,M-1
T=T+1.
P=P*(X/(T*2.))
8 S=S+P
Y=Y-C*S
13 Y=1.-Y
WRITE(6,11)Y
11 FORMAT(5X,'WITH CHI-SQUARE PROBABILITY=',T39,F5.3)
RETURN
14 WRITE(6,15)
15 FORMAT(5X,'NO PROBABILITY CALC. DF OR CHI SQ. OUTSIDE LIMITS')
X=-1
RETURN
END
C
C NORMAL ROUTINE FOR CUCS
C
C SOURCE CHARLES A. NAGY
C
C---------------X IS INPUT. Y IS RETURNED. THIS IS SAME AS
C--------------- CUNO IN APLIB.
SUBROUTINE CUNO(X,Y)
DIMENSION F(25)
DATA (F(I),I=1,25)/.5,.598706326,.691462461,.773372648,
1 .841344746,.894350226,.933192799,.959940843,.977249868,
2 .987775527,.993790335,.997020237,.998650102,.999422975,
3 .999767371,.999911583,.999968329,.999989311,.999996602,
4 .999998983,.999999713,.999999924,.999999981,.999999996,
5 .999999999/
W=ABS(X)
IF(W-6.125)2,2,1
1 Z=1
GO TO 7
2 K=INT(4.*W)+1
A=.25*FLOAT(K-1)
IF(W-A)10,3,4
3 Z=F(K)
GO TO 7
4 IF(W-(A+.125))6,6,5
5 K=K+1
A=A+.25
6 H=W-A
C1=((-A*A+10.)*A*A-15.)*A
C2=(6.*A*A-36.)*A*A+18.
C3=(-30.*A*A+90.)*A
C4=120.*(A*A-1.)
C5=-360.*A
C6=(((((C1*H+C2)*H+C3)*H+C4)*H+C5)*H+720.)*H
Z=F(K)+C6/(720.*SQRT(6.28318531*EXP(A*A)))
7 IF(X)8,9,9
8 Y=1.-Z
RETURN
9 Y=Z
RETURN
10 WRITE(30,11)
11 FORMAT(' ?ERROR IN SUBROUTINE CUNO',/)
X=-1
RETURN
END
C
C FISHER PROBABILITY FUNCTION
C
C
C T-VALUE FOR T-TEST ANALYSIS
C
C
FUNCTION TPCT(ALPHA,KDF)
Z=ALPHA*ALPHA
Q=1./Z
T=SQRT(ALOG(Q))
U=((.010328*T+.802853)*T+2.515517)
V=(((.001380*T+.189269)*T+1.432788)*T+1.)
XP=T-U/V
X2=XP*XP
A=XP*(X2+1.)/4.
B=((5.*X2+16.)*X2+3.)*XP/96.
C=(((3.*X2+19.)*X2+17.)*X2-15.)*XP/384.
D=((((79.*X2+776.)*X2+1482.)*X2-1920.)*X2-945.)*XP/92160.
E=(((((27.*X2+339.)*X2+930.)*X2-1782.)*X2-765.)*X2+17955.)*XP/3686
140.
V=1./KDF
TPCT=XP+(A+(B+(C+(D+E*V)*V)*V)*V)*V
RETURN
END