Trailing-Edge
-
PDP-10 Archives
-
-
There are no other files named in the archive.
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