C ANALYSIS OF VARIANCE FOR ONE-WAY DESIGN JUNE 15, 1966 C THIS IS A SIFTED VERSION OF BMD01V ORIGINALLY WRITTEN IN C FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE C AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION. DOUBLE PRECISION PROBLM,CODE,FINISH,SAMSIZ,SPECTG,PROBNO DIMENSION DATE(2), X(9000), NC(5000), FMT(180), FT(12), 1 SDEV(12),XBAR(12),NUSE(20),LCODE(72),CONST(72),L1(8),D2(8) COMMON X , NCI , NCUE , FMT , NTAPE , LUMP COMMON LCODE , CONST , NH , FM , NTG , CASE C 200 FORMAT('1BMD01V - ANALYSIS OF VARIANCE FOR ONE-WAY DESIGN - REVISE 1D JUNE 24, 1969'/ 341H HEALTH SCIENCES COMPUTING FACILITY, UCLA//) C DATA YES,XNO,PROBLM,FINISH,SAMSIZ,SPECTG/3HYES,2HNO,6HPROBLM, 16HFINISH,6HSAMSIZ,6HSPECTG/ CALL USAGEB('BMD01V') DO 666 I = 1, 20 666 NUSE(I) = 0 NEVER = 0 WRITE (6,200) 1111 READ (5,300)CODE,PROBNO,NP,OUT,NTG,BIN,WIND,NTAPE,IVF IF(CODE.EQ.FINISH)GO TO 1914 IF(CODE.EQ.PROBLM) GO TO 2 608 WRITE(6,6080) 6080 FORMAT(87H0PROBLM OR FINISH IS PUNCHED WRONG ON THE PROPER CARD OR 1 EITHER CARD IS OUT OF SEQUENCE) STOP 2 IF((NP-1)*(NP-5001)) 6001,6002,6002 6002 WRITE (6,6102) 6102 FORMAT(72H0ERROR ON PROBLM CARD, CHECK PROGRAM LIMITS ON NUMBER OF 1 TREATMENT GROUP) STOP 6005 WRITE(6,6015) 6015 FORMAT(58H0 WRONG PROBLEM CARD, THE NUMBER OF SPECTG CARD IS TOO B 1IG) STOP 609 WRITE(6,6019) 6019 FORMAT(55H00 AND 5 CAN NOT BE SPECIFIED AS INPUT BINARY TAPE UNIT) STOP 610 WRITE(6,6110) 6110 FORMAT(33H0TOTAL NUMBER OF CASES IS TOO BIG) STOP 406 WRITE(6,4016)J,NI 4016 FORMAT(11H0THE NUMBER,I3,15H SPECTG CODE ON,I3,23H SPECTG CARD IS 1ILLEGAL) STOP 605 WRITE(6,6050) 6050 FORMAT(31H0ILLEGAL NUMBER OF SPECTG CARDS) STOP 613 WRITE(6,6013)I 6013 FORMAT(' THE NUMBER OF CASES IN GROUP',I5,' IS OUTSIDE THE RANGE O *F BETWEEN 1 AND 9000, PROBLEM WILL BE IGNORED') STOP 6001 IF(NTG-9) 6003,6003,6005 6003 NCUE=1 IF(BTWEEN(NP,2,500).LE.0) GO TO 6002 IF(BIN.NE.YES)GO TO 1113 1112 NCUE=2 1113 IF (BTWEEN(NP, 2, 5000)*BTWEEN(IVF, 0, 10)*BTWEEN(NCUE, 1, 2)* 1 (BTWEEN(NTAPE, 0, 5) + BTWEEN(NTAPE, 7, 15))) 607, 607, 3 3 NUTS=NTAPE*(NTAPE-5) GO TO (4, 5), NCUE 5 IF (NUTS) 7, 609, 7 4 IF(IVF.GT.0.AND.IVF.LE.10)GO TO 6 IVF=1 WRITE(6,4000) 6 IF (NUTS) 7, 8, 7 8 NTAPE = 5 GO TO 9 7 NUSE(NTAPE) = 1 IF(WIND.EQ.XNO)GO TO 9 667 REWIND NTAPE 9 IF (NEVER) 10, 10, 11 11 WRITE (6,200) 10 NEVER = NEVER + 1 WRITE (6,400)PROBNO, NP NCARDS=(NP+12)/13 I2=0 DO 17 J=1,NCARDS I1=I2+1 I2=I2+13 IF(NP-I2)15,16,16 15 I2=NP 16 READ (5,700)CODE,(NC(I),I=I1,I2) IF(CODE.EQ.SAMSIZ)GO TO 17 WRITE(6,7000) 7000 FORMAT(55H0SAMSIZ IS SPELLED WRONG OR SAMSIZ CARD OUT OF SEQUENCE) CALL EXIT 165 NCARDS=-NCARDS 17 CONTINUE IF(-NCARDS)18,18,6002 18 CASE=0.0 DO 19 I=1,NP IF(BTWEEN(NC(I),1,9000))613,613,181 181 CASE1=NC(I) CASE=CASE+CASE1 19 CONTINUE IF(CASE-(10.0**8))195,195,610 195 IF(NTG)605,499,401 401 NH=0 DO 445 NI=1,NTG READ(5,425)CODE,M,(L1(I),D2(I),I=1,M) IF(CODE.EQ.SPECTG)GO TO 404 403 WRITE (6,405) STOP 404 DO 6004 J=1,M IF(L1(J)*(L1(J)-11)) 6004,406,406 6004 CONTINUE DO 415 NG=1,M NG1=NH+NG LCODE(NG1)=L1(NG) 415 CONST(NG1)=D2(NG) 445 NH=NH+M 499 GO TO (12,13),NCUE 12 WRITE (6,500)IVF, NTAPE WRITE(6,750) 750 FORMAT(20H0THE VARIABLE FORMAT) IVF = IVF*18 READ (5,100)(FMT(I), I = 1, IVF) WRITE(6,101)(FMT(I),I=1,IVF) 101 FORMAT(1X,18A4) GO TO 14 13 READ (5,601)LUMP WRITE (6,600)LUMP, NTAPE 14 IF(NTG)605,439,426 426 WRITE (6,427) DO 428 NG=1,NH 428 WRITE (6,431)LCODE(NG),CONST(NG) 439 SPSQ=0.0 SMSQ = 0.0 GT = 0.0 NP1 = NP - 1 NCI = NC(1) TOBS = NCI CALL READIN DO 20 I = 1, NCI 20 GT = GT + X(I) Y = GT/TOBS DO 21 J = 1, NCI 21 SPSQ = SPSQ + (X(J) - Y)**2 NGO = 1 IF(OUT.EQ.YES)GO TO 31 30 WRITE (6,801) GO TO 40 31 SDV = SQRT( SPSQ/(TOBS - 1.0) ) CALL DECFND(Y, LY) CALL DECFND(SDV, LSDV) IF (NP - 12) 32, 32, 33 33 WRITE (6,900) LY= MAX0(0, 7 - LY) LSDV= MAX0(0, 7 - LSDV) WRITE (6,901)NGO, NCI, Y, SDV NGO = 3 GO TO 40 32 LY= MAX0(LY, LSDV) LSDV= MAX0(0, 5 - LY) XBAR(1) = Y SDEV(1) = SDV NGO = 2 40 DO 50 I = 2, NP NCI = NC(I) XNI = NCI TOBS = TOBS + XNI CALL READIN AVE = 0.0 DO 51 J = 1, NCI 51 AVE = AVE + X(J) GT = GT + AVE AVE = AVE/XNI SDV = 0.0 DO 52 J = 1, NCI 52 SDV = (X(J) - AVE)**2 + SDV SPSQ = SPSQ + SDV SMSQ = SMSQ + XNI*(AVE - Y)**2 SDV = SQRT( SDV/(XNI - 1.0) ) GO TO (50, 61, 62), NGO 61 XBAR(I) = AVE SDEV(I) = SDV GO TO 50 62 WRITE (6,901)I, NCI, AVE, SDV 50 CONTINUE GO TO (70, 71, 70), NGO 71 WRITE (6,1000)(I, I = 1, NP) 1000 FORMAT(////22H TREATMENT GROUP , 12(I7, 2X) ) WRITE (6,1001)(NC(I), I = 1, NP) WRITE (6,1002)(XBAR(I), I = 1, NP) WRITE (6,1003)(SDEV(I), I = 1, NP) 70 WRITE (6,1100) GT = GT/TOBS SMSQ = SMSQ - TOBS*(GT - Y)**2 GT = SMSQ + SPSQ TOBS = TOBS - 1.0 TLEVEL = NP1 DF = TOBS - TLEVEL AVE = SMSQ/TLEVEL SDV = SPSQ/DF FRATIO = AVE/SDV WRITE (6,1200)SMSQ, NP1, AVE, FRATIO INTDF=DF +.0000001 INTOBS=TOBS+.0000001 WRITE (6,1500)SPSQ, INTDF , SDV WRITE (6,1600)GT, INTOBS GO TO 1111 607 WRITE (6,6070) STOP 1914 DO 73 I = 1, 20 IF (NUSE(I)) 72, 73, 72 72 IF(I.NE.5)REWIND I 73 CONTINUE WRITE (6,1900) 1915 STOP 100 FORMAT(18A4) 300 FORMAT(2A6,I4,A3,I1,43X,A3,A2,2I2) 400 FORMAT(///// 15H PROBLEM CODE , A6, / 20H NUMBER OF TREATMENT , 1 9H GROUPS , I4 ) 405 FORMAT(51H0ERROR ON TRANS-GENERATION CARD, SPECTG IS MISSPELT) 425 FORMAT(A6,I1,8(I2,F6.0)) 427 FORMAT(1H023HSPECIAL TRANSGENERATION/1H 3X4HCODE4X8HCONSTANT) 431 FORMAT(1H 5XI2,1XF11.5) 500 FORMAT(34H NUMBER OF VARIABLE FORMAT CARDS , I2 / 1 18H DATA INPUT TAPE , I2 ) 600 FORMAT(23H BINARY RECORD LENGTH , I3 / 18H DATA INPUT TAPE , I2) 601 FORMAT(I3) 700 FORMAT(A6, 1X, 13I5) 801 FORMAT(////41H LISTING OF TREATMENT MEANS NOT REQUESTED ) 900 FORMAT(//// 49H TREATMENT SAMPLE SIZE MEAN STANDARD 1 10H DEVIATION //) 901 FORMAT (1X,I6,8X,I6,F16.5, F17.5) 1001 FORMAT(/ 22H SAMPLE SIZE , 12(I7, 2X) ) 1002 FORMAT (/22H MEAN , 12F9.4) 1003 FORMAT (/22H STANDARD DEVIATION , 12F9.4) 1100 FORMAT(///// 33X, 21H ANALYSIS OF VARIANCE /// 1 19X, 15H SUM OF SQUARES, 6X, 2HDF, 5X, 12H MEAN SQUARE , 5X, 2 8H F RATIO / ) 1200 FORMAT(15H0BETWEEN GROUPS , F17.4, I10, 2F15.4 ) 1300 FORMAT(2F11.0) 1400 FORMAT(2A5, 1X, 2A5) 1500 FORMAT(15H0WITHIN GROUPS , F17.4, I10, F15.4) 1600 FORMAT(// 1H , 5X, 5HTOTAL, F21.4, I10) 1900 FORMAT(40H1COMPUTATION HALTED BY FINISH CARD. ) 4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF 1IED, ASSUMED TO BE 1.) 6070 FORMAT(123H0PROBLEM CARD IS SET UP WRONG, EITHER INPUT TAPE UNIT I 1S SPECIFIED WRONG OR THE NUMBER OF F-TYPE VARIABLE FORMAT IS TOO B 2IG) END C FUNCTION BTWEEN FOR BMD01V JUNE 15, 1966 C BTWEEN = 1.0 IF N IS BETWEN I AND J FUNCTION BTWEEN(N, I, J) BTWEEN=0.0 IF(N.GE.I.AND.N.LE.J)BTWEEN=1.0 RETURN END C SUBROUTINE DECFND FOR BMD01V JUNE 15, 1966 C X HAS A MAGNITUDE BETWEEN 10**(LEFT-1) AND 10**(LEFT) SUBROUTINE DECFND(X, LEFT) A=10. DO 1 LEFT=1,19 IF(A.GT.X) GO TO 2 1 A=A*10. 2 RETURN END C SUBROUTINE DIVALG FOR BMD01V JUNE 15, 1966 C DIVISION ALGORITHM SUBROUTINE DIVALG(NA, NB, NQ, NR, NCODE) NN = NA IF (NCODE) 1, 1, 2 2 NN = NN - 1 1 NQ = NN/NB NR = NA - NB*NQ RETURN END C SUBROUTINE READIN FOR BMD01V JUNE 15, 1966 SUBROUTINE READIN DIMENSION X(9000),LCODE(72),CONST(72),FMT(180) COMMON X , NCI , NCUE , FMT , NTAPE , LUMP COMMON LCODE , CONST , NH , FM , NTG GO TO (1, 2), NCUE 1 READ (NTAPE,FMT)(X(I), I = 1, NCI) 6 IF(NTG)7,8,7 7 CALL TRNGEN 8 RETURN 2 CALL DIVALG(NCI, LUMP, NDO, LEFT, 1) N2 = 0 IF (NDO) 3, 3, 4 4 DO 5 J = 1, NDO N1 = N2 + 1 N2 = N2 + LUMP 5 READ(NTAPE,9876)(X(I),I=N1,N2) 3 N1 = N2 + 1 READ(NTAPE,9876)(X(I),I=N1,NCI) GO TO 6 9876 FORMAT(20A4) END C SUBROUTINE TRNGEN FOR BMD01V JUNE 15, 1966 SUBROUTINE TRNGEN DIMENSION X(9000),LCODE(72),CONST(72),FMT(180) COMMON X , NCI , NCUE , FMT , NTAPE , LUMP COMMON LCODE , CONST , NH , FM , NTG , SAMP ASN(XX)=ATAN(XX/SQRT(1.0-XX**2)) DO 210 I=1,NH FM=CONST(I) JUMP=LCODE(I) DO 85 J=1,NCI D=X(J) GO TO(10,20,30,40,50,60,70,80,90,110),JUMP 10 IF(D)99,11,12 11 D=0.0 GO TO 8 12 D=SQRT(D) GO TO 8 20 IF(D)99,21,22 21 D=1.0 GO TO 8 22 D=SQRT(D)+SQRT(D+1.0) GO TO 8 30 IF(D)99,99,31 31 D=ALOG10(D) GO TO 8 40 D=EXP(D) GO TO 8 50 IF(-D)52,11,99 52 IF(D-1.0)53,54,99 54 D=3.14159265/2.0 GO TO 8 53 D=ASN(SQRT(D)) GO TO 8 60 A=D/(SAMP+1.0) B=A+1.0/(SAMP+1.0) IF(A)99,61,62 61 IF(-B)64,11,99 64 D=ASN(SQRT(B)) GO TO 8 62 IF(B)99,65,66 65 D=ASN(SQRT(A)) GO TO 8 66 D=ASN(SQRT(A))+ASN(SQRT(B)) GO TO 8 70 IF(D)71,99,71 71 D=1.0/D GO TO 8 80 D=D+FM GO TO 8 90 D=D*FM GO TO 8 110 IF(D)111,11,111 111 D=D**FM 8 X(J)=D 85 CONTINUE 210 CONTINUE GO TO 1000 99 WRITE(6,105)I,D STOP 105 FORMAT(23H0TRANS-GENERATION ERROR//10H PASS NO.=I3,9H X VALUE=F10. 15) 1000 RETURN END