C DISCRIMINANT ANALYSIS - TWO GROUPS JUNE 9, 1966 C THIS IS A SIFTED VERSION OF BMD04M ORIGINALLY WRITTEN IN C FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE C AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION. DIMENSION X(2,20,250),SUM(2,25,25),TOTAL(2,25),D(25),B(25),Z(2,300 1), MMM(25), N(2),LL(25),DIV(2),C(25,25),DD(25),ZBAR(2),VARZ(2), 2LLL(25),IRANK(2,300),LX(25),LY(25),FMT(180) C 343 FORMAT(55H1BMD04M - DISCRIMINANT ANALYSIS-TWO GROUPS - VERSION OF, 118H JUNE 9, 1966 ,/ 241H HEALTH SCIENCES COMPUTING FACILITY, UCLA) C DOUBLE PRECISION A1,A2,A3,FINISH DATA A2,FINISH,A3/6HPROBLM,6HFINISH,6HSELECT/ MTAPE=5 CALL USAGEB('BMD04M') 6 READ (5,310)A1,PROB,K,N(1),N(2),NVG,NADD,NTAPE,NUM,KVR IF(A1.EQ.A2)GO TO 110 IF(A1.EQ.FINISH)GO TO 100 107 WRITE (6,345) A1 GO TO 100 110 CALL TPWD(NTAPE,MTAPE) 116 WRITE (6,343) IF((K-1)*(K-26))112,348,348 112 K1=K+NADD IF((K1-1)*(K1-26))114,350,350 114 WRITE (6,339)PROB,K1 DO 117 I=1,2 IF((N(I)-K1+1)*(N(I)-301))117,352,352 117 CONTINUE NTRAN=0 IERROR=0 TYPE 9876,KVR IF(KVR.GT.0.AND.KVR.LE.10)GO TO 118 KVR=1 WRITE(6,4000) 118 KVR=KVR*18 TYPE 9876,KVR 9876 FORMAT(I7) READ (5,344)(FMT(I),I=1,KVR) WRITE (6,347) (FMT(I),I=1,KVR) DO 410 I=1,2 NN=N(I) DO 410 J=1,NN 410 READ (MTAPE,FMT)(X(I,L,J), L=1,K) IF(NVG)21,21,500 500 CALL TRNGEN(X,K,N,IERROR,2,NVG) IF(IERROR) 10, 21, 21 10 DO 15 I=1,NUM 15 READ (5,340)A1 GO TO 6 21 K=K1 NSEL=0 DO 35 M=1,2 DIV(M)=N(M) DO 30 I=1,K LL(I)=I TOTAL(M,I)=0.0 DO 30 L=I,K C(I,L)=0.0 SUM(M,I,L)=0.0 NN=N(M) DO 25 J=1,NN 25 SUM(M,I,L)=SUM(M,I,L)+X(M,I,J)*X(M,L,J) 30 SUM(M,L,I)=SUM(M,I,L) DO 35 I=1,K NN=N(M) DO 35 J=1,NN 35 TOTAL(M,I)=TOTAL(M,I)+X(M,I,J) DO 37 I=1,K 37 D(I)=TOTAL(1,I)/DIV(1)-TOTAL(2,I)/DIV(2) DO 38 I=1,K 38 DD(I)=D(I) WRITE (6,313) WRITE (6,321) DO 36 I=1,2 DO 36 J=1,K 36 Z(I,J)=TOTAL(I,J)/DIV(I) DO 39 I=1,K 39 WRITE (6,314)I,Z(1,I),Z(2,I),D(I) DO 40 I=1,K DO 40 L=1,K C(I,L)=SUM(1,I,L)+SUM(2,I,L)-TOTAL(1,I)*TOTAL(1,L)/DIV(1)-TOTAL(2, 1I)*TOTAL(2,L)/DIV(2) 40 SUM(1,I,L)=C(I,L) WRITE (6,301) DO 42 I=1,K 42 WRITE (6,302)(SUM(1,I,L),L=1,K) 43 CALL INVERT (C,K,DET,LX,LY) WRITE (6,323) DO 44 I=1,K 44 WRITE (6,302)(C(I,J),J=1,K) DO 45 I=1,K B(I)=0.0 DO 45 L=1,K 45 B(I)=B(I)+C(I,L)*D(L) WRITE (6,303) WRITE (6,302)(B(I),I=1,K) DSQ=0.0 DO 48I=1,K DO 48 J=1,K 48 DSQ=DSQ+D(I)*D(J)*C(I,J) DSQ=DSQ*(DIV(1)+DIV(2)-2.0) WRITE (6,320)DSQ IDF=N(1)+N(2)-K-1 DK=IDF FK=K VAL=DIV(1)*DIV(2)*DK/(FK*(DIV(1)+DIV(2))*(DIV(1)+DIV(2)-2.0)) VAL=VAL*DSQ WRITE (6,341)K,IDF,VAL DO 50 M=1,2 NN=N(M) DO 50 J=1,NN Z(M,J)=0.0 DO 50 I=1,K LI=LL(I) 50 Z(M,J)=Z(M,J)+B(I)*X(M,LI,J) WRITE (6,322) DO 52 M=1,2 NN=N(M) SUMZ=0.0 SUMZSQ=0.0 ZBAR(M)=0.0 DIVN=N(M) VARZ(M)=0.0 DO 51 I=1,NN SUMZ=SUMZ+Z(M,I) 51 SUMZSQ=SUMZSQ+Z(M,I)**2 ZBAR(M)=SUMZ/DIVN VARZ(M)=(SUMZSQ-SUMZ**2/DIVN)/(DIVN-1.0) STDVZ=SQRT(VARZ(M)) 52 WRITE (6,319)M,NN,ZBAR(M),VARZ(M),STDVZ WRITE (6,304) WRITE (6,324) DO 53 I=1,2 DO 53 J=1,300 53 IRANK (I,J)=0 NTOTAL=N(1)+N(2) DO 70 I=1,NTOTAL HOLD=-(10.0**35.0) DO 59 M=1,2 NN=N(M) DO 59 J=1,NN IF(Z(M,J)-HOLD)59,59,54 54 IF(IRANK(M,J))55,55,59 55 MM=M JJ=J HOLD=Z(M,J) 59 CONTINUE IRANK (MM,JJ)=999 IF(MM-1)60,60,61 60 WRITE (6,305)I,HOLD,JJ GO TO 70 61 WRITE (6,325)I,HOLD,JJ 70 CONTINUE 71 IF(NUM) 6, 6, 72 72 NSEL=NSEL+1 WRITE (6,317)NSEL READ (5,340)A1,K,(LL(I), I=1,25) IF(A1.EQ.A3)GO TO 74 73 WRITE (6,346)NSEL WRITE (6,354) A1 NUM=NUM-1 GO TO 71 74 WRITE (6,309) WRITE (6,307)(LL(I),I=1,K) DO 76 I=1,K LLI=LL(I) DO 75 L=1,K LLLL=LL(L) 75 C(I,L)=SUM(1,LLI,LLLL) 76 D(I)=DD(LLI) NUM=NUM-1 GO TO 43 301 FORMAT(36H0 SUM OF PRODUCTS OF DEV. FROM MEANS) 302 FORMAT(7F16.5) 303 FORMAT(36H0 DISCRIMINANT FUNCTION COEFFICIENTS) 304 FORMAT(69H0 FIRST GROUP SECOND GROUP FIRST GROUP 1 SECOND GROUP) 305 FORMAT(1H I4,F17.5,17X,I12) 307 FORMAT(10I5) 309 FORMAT(28H VARIABLES USED IN FUNCTION) 310 FORMAT(A6,A2,I2,2I3,4I2,46X,I2) 313 FORMAT(49H0 VARIABLE MEANS BY GROUP AND DIFFERENCE IN MEANS) 314 FORMAT(I5,2X,3F16.5) 317 FORMAT(1H0//15H SELECTION NO.I4) 319 FORMAT(I7,I14,F17.5,2F20.5) 320 FORMAT(22H0 MAHALANOBIS DSQUARE=F16.5) 321 FORMAT(57H VARIABLE MEAN 1 MEAN 2 DIFFERENC 1E) 322 FORMAT(79H0 POP. NO. SAMPLE SIZE MEAN Z VARIANC 1E Z STD. DEV. Z) 323 FORMAT(47H0 INVERSE OF SUM OF PRODUCTS OF DEV. FROM MEANS) 324 FORMAT(66H RANK VALUES VALUES ITEM NO. 1 ITEM NO.) 325 FORMAT(1H I4,17X,F17.5,12X,I13) 339 FORMAT(15H0 PROBLEM NO. A2/21H NUMBER OF VARIABLESI4) 340 FORMAT(A6,26I2) 341 FORMAT(4H0 F(I2,1H,I3,2H)=F16.5) 344 FORMAT(18A4) 345 FORMAT (' PROBLM OR FINISH CARD EXPECTED,BUT READ ',A6) 346 FORMAT(40H0 ERROR ON SELECTION CARD OR DECK SET-UPI4) 347 FORMAT (' VARIABLE FORMAT IS '/(5X,18A4)) 348 WRITE (6,349) K 349 FORMAT (1X,I3,' ORIGINAL VARIABLES IS ILLEGAL') GO TO 100 350 WRITE (6,351) K1 351 FORMAT (1X,I4,' IS AN ILLEGAL NUMBER OF VARIABLES AFTER TRANSFORMA 1TIONS') GO TO 100 352 WRITE (6,353) N(I),I 353 FORMAT (1X,I4,' SAMPLE SIZE FOR GROUP ',I2,' IS ILLEGAL') GO TO 100 354 FORMAT (' CARD READ WAS ',A6) 4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF 1IED, ASSUMED TO BE 1.) 100 IF(MTAPE-5)102,102,101 101 REWIND MTAPE 102 STOP END C SUBROUTINE INVERT FOR BMD04M JUNE 9, 1966 SUBROUTINE INVERT (A,N,D,L,M) C PROGRAM FOR FINDING THE INVERSE OF A NXN MATRIX DIMENSION A(25,25),L(25),M(25) C SEARCH FOR LARGEST ELEMENT D=1.0 DO80 K=1,N L(K)=K M(K)=K BIGA=A(K,K) DO20 I=K,N DO20 J=K,N IF(ABS(BIGA)-ABS(A(I,J))) 10,20,20 10 BIGA=A(I,J) L(K)=I M(K)=J 20 CONTINUE C INTERCHANGE ROWS J=L(K) IF(L(K)-K) 35,35,25 25 DO30 I=1,N HOLD=-A(K,I) A(K,I)=A(J,I) 30 A(J,I)=HOLD C INTERCHANGE COLUMNS 35 I=M(K) IF(M(K)-K) 45,45,37 37 DO40 J=1,N HOLD=-A(J,K) A(J,K)=A(J,I) 40 A(J,I)=HOLD C DIVIDE COLUMN BY MINUS PIVOT 45 DO55 I=1,N 46 IF(I-K)50,55,50 50 A(I,K)=A(I,K)/(-A(K,K)) 55 CONTINUE C REDUCE MATRIX DO65 I=1,N DO65 J=1,N 56 IF(I-K) 57,65,57 57 IF(J-K) 60,65,60 60 A(I,J)=A(I,K)*A(K,J)+A(I,J) 65 CONTINUE C DIVIDE ROW BY PIVOT DO75 J=1,N 68 IF(J-K)70,75,70 70 A(K,J)=A(K,J)/A(K,K) 75 CONTINUE C CONTINUED PRODUCT OF PIVOTS D=D*A(K,K) C REPLACE PIVOT BY RECIPROCAL A(K,K)=1.0/A(K,K) 80 CONTINUE C FINAL ROW AND COLUMN INTERCHANGE K=N 100 K=(K-1) IF(K) 150,150,103 103 I=L(K) IF(I-K) 120,120,105 105 DO110 J=1,N HOLD=A(J,K) A(J,K)=-A(J,I) 110 A(J,I)=HOLD 120 J=M(K) IF(J-K) 100,100,125 125 DO130 I=1,N HOLD=A(K,I) A(K,I)=-A(J,I) 130 A(J,I)=HOLD GO TO 100 150 RETURN END C SUBROUTINE TPWD FOR BMD04M JUNE 9, 1966 SUBROUTINE TPWD(NT1,NT2) IF(NT1)40,10,12 10 NT1=5 12 IF(NT1-NT2)14,19,14 14 IF(NT2.EQ.5)GO TO 18 17 REWIND NT2 19 IF(NT1-5)18,24,18 18 IF(NT1-6)22,40,22 22 REWIND NT1 24 NT2=NT1 28 RETURN 40 WRITE (6,49) STOP 49 FORMAT(25H ERROR ON TAPE ASSIGNMENT) END C SUBROUTINE TRNGEN FOR BMD04M JUNE 9, 1966 SUBROUTINE TRNGEN (DATA,NVAR,NSAM,IERROR,NGRP,NVG) DIMENSION DATA(2,20,250),NSAM(2) ASN(XX)=ATAN(XX/SQRT(1.0-XX**2)) DOUBLE PRECISION A1,A2 DATA A2/6HTRNGEN/ MARY=0 WRITE (6,1403) WRITE (6,1400) DO1000J=1,NVG READ (5,1100)A1,NEWA,LCODE,LVA,BNEW IF(A1.EQ.A2)GO TO 250 200 WRITE (6,1401)J IERROR=-999 250 WRITE (6,1402)J,NEWA,LCODE,LVA,BNEW IF(-IERROR)251,251,1000 251 IF(LCODE*(LCODE-15))255,500,500 255 IF(LCODE-10)4,5,5 5 NEWB=BNEW 4 DO3000JJ=1,NGRP NN=NSAM(JJ) DO 300 I=1,NN D=DATA(JJ,LVA,I) GOTO(10,20,30,40,50,60,70,80,90,100,110,120,130,140),LCODE 10 IF(D)99,7,8 7 D2=0.0 GOTO3 8 D2=SQRT(D) GOTO3 20 IF(D)99,11,12 11 D2=1.0 GOTO3 12 D2=SQRT(D)+SQRT(D+1.0) GOTO3 30 IF(-D)14,99,99 14 D2=ALOG10(D) GO TO 3 40 D2=EXP(D) GOTO3 50 IF(-D)17,7,99 17 IF(D-1.0)18,19,99 19 D2=3.14159/2.0 GOTO3 18 D2=ASN(SQRT(D)) GOTO3 60 FN=NN A=D/(FN+1.0) B=A+1.0/(FN+1.0) IF(A)99,23,24 23 IF(-B)27,7,99 27 D2=ASN(SQRT(B)) GOTO3 24 IF(B)99,28,29 28 D2=ASN(SQRT(A)) GOTO3 29 A=SQRT(A) B=SQRT(B) D2=ASN(A)+ASN(B) GOTO3 70 IF(D)31,99,31 31 D2=1.0/D GOTO3 80 D2=D+BNEW GOTO3 90 D2=D*BNEW GOTO3 100 IF(-D)33,7,99 33 D2=D**NEWB GOTO3 110 D2=D+DATA(JJ,NEWB,I) GOTO3 120 D2=D-DATA(JJ,NEWB,I) GOTO3 130 D2=D*DATA(JJ,NEWB,I) GOTO3 140 IF(DATA(JJ,NEWB,I))34,99,34 34 D2=D/DATA(JJ,NEWB,I) GOTO3 99 IF(MARY)43,44,44 44 MARY=-999 IERROR=-999 WRITE (6,1404)J 43 WRITE (6,1405)I GO TO 300 3 DATA(JJ,NEWA,I)=D2 300 CONTINUE IF(IERROR)45,3000,3000 45 WRITE (6,1406)JJ 3000 CONTINUE IF(IERROR)42,1000,1000 500 WRITE (6,1408)NEWA 1000 CONTINUE IF(IERROR) 42,1111,1111 42 WRITE (6,1407) 1100 FORMAT(A6,I3,I2,I3,F6.0) 1400 FORMAT(46H0CARD NEW TRANS ORIG. ORIG. VAR(B)/45H NO. 1VARIABLE CODE VAR(A) OR CONSTANT) 1401 FORMAT(30H0ERROR ON TRANSGENERATION CARDI4) 1402 FORMAT(2H I2,I8,2I9,F15.5) 1403 FORMAT(1H06X,21HTRANSGENERATION CARDS) 1404 FORMAT(30H0THE INSTRUCTIONS INDICATED ON 25H TRANSGENERATION CARD 1NO.I2,4H RE-/29H SULTED IN THE VIOLATION OF A 31H RESTRICTION FOR 2THIS TRANSFOR-/31H MATION. THE VIOLATION OCCURED 29H FOR THE ITEM 3S LISTED BELOW./) 1405 FORMAT(10H ITEM NO. I3) 1406 FORMAT(32H0ITEMS ABOVE ARE FROM GROUP NO. I2//) 1407 FORMAT(41H0PROGRAM CANNOT CONTINUE FOR THIS PROBLEM) 1408 FORMAT(1H0,35X,95HILLEGAL TRANSGENERATION CODE SPECIFIED ON PRECED 1ING CARD. PROGRAM WILL PROCEED LEAVING VARIABLE,I3,11H UNCHANGED.) 1111 RETURN END