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
































































