PDP-10 Archive: decus/20-0137/oneaov/oneaov.for from decus_20tap5

Trailing-Edge - PDP-10 Archives - decus_20tap5_198111 - 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