Trailing-Edge
-
PDP-10 Archives
-
decuslib10-10
-
43,50520/stp16.stp
There is 1 other file named stp16.stp in the archive. Click here to see a list.
C ***** STAT PACK *****
C SUBROUTINE FOR PERFORMING DISCRIMINATE ANALYSIS
C CALLING SEQUENCE: CALL DISCR(NV,NC,MV,MC,DATA,AV,JV,NAMES)
C WHERE: NV - NUMBER OF VARIABLES
C NC - NUMBER OF OBSERVATIONS
C MV - MAXIMUM NUMBER OF VARIABLES POSSIBLE
C MC - MAXIMUM NUMBER OF OBSERVATIONS POSSIBLE
C DATA - MATRIX CONTAINING DATA
C AV - EXTRA VECTOR AT LEAST NC LONG
C JV - EXTRA VECTOR AT LEAST NC LONG
C NAMES - VECTOR CONTAINING NAMES OF VARIABLES
C
C DISCRIMINATE ANALYSIS WAS IN ORIGINAL REQUEST FOR ROUTINES
C WHICH PROMPTED THE CREATION OF STAT PACK, SUBMITTED
C BY AL LEADER (MANAGEMENT). MORE RECENTLY SAM ANEMA
C (COMPUTER CENTER FACULTY,STAFF COORDINATOR) EXPRESSED
C AN OPINION THAT STAT PACK SHOULD CONTAIN A DISCRIMANNTE
C ANALYSIS. IN DISCUSSIONS WITH ULDIS SMMIDCHENS THE NEED
C WAS ALSO INDICATED. THIS ROUTINE IS ADOPTED FROM A PROGRAM
C ORIGINALLY WRITTEN BY SAM ANEMA.
C
SUBROUTINE DISCR(NV,NC,MV,MC,DATA,AV,JV,NAMES)
DIMENSION XBAR(25,25),D(25,25),C(25,25),IP(25,25)
DIMENSION OPT(8),IVAR(25),LOPT(8),RANGE(3,25),INPUT(50)
DIMENSION IVALUE(15),SPACE(3),IV(25),XC(25),CON(25),XVAL(25)
DIMENSION JV(1),AV(1),NAMES(1),DATA(MC,MV),NGRP(25)
EQUIVALENCE (XBAR,IP),(INPUT,IV)
COMMON /DEV/ ICC,IDATA,IOUT,IDLG,IDSK
COMMON /PRNT/ LINPP,ICOPS,RUNPRG
COMMON /EXTRA/ HEDR(70),NSZ
DATA LOPT/'DISCR','RANGE','MEANS','DISPM','INVDM','EVALF',
1'XPROD','CLASS'/
1 IF(ICC.NE.2) WRITE(IDLG,2)
2 FORMAT('0ENTER OPTIONS SEPARATED BY COMMAS'/)
READ(ICC,3,END=1)(IVAR(I),I=1,10)
3 FORMAT(10(A5,1X))
IF(IVAR(1).EQ.'!') RETURN
DO 13 I=1,8
13 OPT(I)=0
DO 12 I=1,10
IF(IVAR(I).EQ.' ') GO TO 20
IF(IVAR(I).NE.'HELP,') GO TO 6
4 WRITE(IDLG,5)
5 FORMAT('0THE DISCRIMINATE ANALYSIS MAY BE USED FOR UP TO'/
1' 25 VARIABLES WITH UP TO 25 GROUPS. PLACEMENT OF OBSERVATIONS'/
2' INTO GROUPS IS BASED ON THE VALUE OF A BREAKDOWN VARIABLE.'/
3' RANGES FOR GROUPING OBSERVATIONS ARE ASSUMED TO BE ENTERED BY'/
4' THE USER, HOWEVER, OPTIONS ARE AVAILABLE TO FORM RANGES'/
5' AUTOMATICALLY. IF THE BREAKDOWN VARIABLE HAS 25 OR FEWER'/
6' DISCRETE VALUES, A SEPARATE RANGE WILL BE CREATED FOR EACH'/
7' INDIVIDUAL VALUE, OTHERWISE THE RANGE OF VALUES FOR THE '/
8' BREAKDOWN VARIABLE WILL BE BROKEN INTO 25 SEPARATE AND EQUAL'/
9' RANGES.'/
1'0OPTIONS AVAILABLE ARE:'/
2' "DISCR" - AUTOMATICALLY FORM GROUPINGS'/
3' "RANGE" - LIST RANGES USED IN BREAKDOWNS (AVAILABLE ONLY IF'/
4' "DISCR" IS USED.)'/
5' "AUTO" - SAME AS "DISCR" EXCEPT SPECIFIED WHEN RANGES ARE'/
6' REQUESTED.'/
7' "MEANS" - PRINT A MEANS TABLE FOR VARIABLE GROUP COMBINATIONS'/
8' "XPROD" - PRINT THE CROSS PRODUCTS OF VARIANCE FROM MEAN FOR'/
9' EACH GROUP')
WRITE(IDLG,16)
16 FORMAT(' "DISPM" - PRINT THE DISPERSION MATRIX'/
1' "INVDM" - PRINT THE INVERSE OF THE DISPERSION MATRIX'/
2' "EVALF" - EVALUATE THE DISCRIMINATE FUNCTION FOR EACH OBSER.'/
4' "CLASS" - PRINT THE CLASSIFICATION MATIRX'/
5' "OUTPT" - "MEANS", "DISPM", "XPROD", "INVDM", "EVALF", AND'/
6' "CLASS"'/
7' IF NO OPTIONS ARE DESIRED TYPE A CARRIAGE RETURN')
GO TO 1
6 IF(IVAR(I).EQ.'AUTO') GO TO 7
IF(IVAR(I).NE.'AUTO,') GO TO 9
7 WRITE(IDLG,8)
8 FORMAT(' "AUTO" IS ONLY SPECIFIED WHEN RANGES ARE REQUESTED')
GO TO 1
9 IOPT=0
DO 10 J=1,8
IF(IVAR(I).NE.LOPT(J)) GO TO 10
OPT(J)=1
GO TO 12
10 CONTINUE
IF(IVAR(I).NE.'OUTPT') GO TO 15
DO 14 J=3,8
14 OPT(J)=1
GO TO 12
15 WRITE(IDLG,11) IVAR(I)
11 FORMAT(' OPTION "',A5,'" DOES NOT EXIST')
GO TO 1
12 CONTINUE
C
C OPTIONS TAKEN CRE OF NOW FIND VARIABLES TO BE ANALYSED
C
20 IF(ICC.NE.2) WRITE(IDLG,21)
21 FORMAT(' WHICH VARIABLES ARE TO BE ANALYSED?'/)
CALL ALPHA(IVAR,25,NN,IRET,IHELP,IERR,NAMES,NV)
IF(IRET.EQ.1) RETURN
IF(IERR.EQ.1) GO TO 20
IF(IHELP.EQ.1) GO TO 20
IF(NN.LT.NV) GO TO 23
WRITE(IDLG,22)
22 FORMAT(' MORE VARIABLES SPECIFIED THAN POSSIBLE FOR THIS DATA')
GO TO 20
23 DO 24 I=1,NN-1
IF(IVAR(I).LT.0) GO TO 24
DO 924 J=I,NN
IF(IVAR(I).NE.IVAR(N)) GO TO 924
WRITE(IDLG,25)NAMES(IVAR(I))
25 FORMAT(' VARIABLE "',A5,'" SPECIFIED TWICE')
GO TO 20
924 CONTINUE
24 CONTINUE
30 IF(ICC.NE.2) WRITE(IDLG,31)
31 FORMAT(' ENTER THE BREAKDOWN VARIABLE? ',$)
CALL ALPHA(IBREAK,1,I,IRET,IHELP,IERR,NAMES,NV)
IF(IERR.EQ.1) GO TO 30
IF(IRET.EQ.1) RETURN
IF(IHELP.EQ.1) GO TO 30
IF(IBREAK.GT.0) GO TO 33
WRITE(IDLG,32)
32 FORMAT(' "*" AND "ALL" MAY NOT BE USED FOR A BREAKDOWN')
GO TO 30
33 DO 34 I=1,NN
IF(IBREAK.NE.IVAR(I)) GO TO 34
WRITE(IDLG,35)
35 FORMAT(' THE BREAKDOWN VARIABLE WAS ALSO SPECIFIED AS ONE OF'/
1' THE VARIALBES TO BE ANALYSED')
GO TO 30
34 CONTINUE
C
C NOW WORK ON BREAKDOWN VARIABLE
C
DO 40 J=1,NC
JV(J)=J
40 AV(J)=DATA(J,IBREAK)
CALL SORT(AV,JV,NC)
IF(OPT(1).EQ.1) GO TO 100
C
C RANGES ENTERED BY USER
C
70 IF(ICC.NE.2) WRITE(IDLG,71)
71 FORMAT(' ENTER RANGES FOR BREAKDOWNS'/)
I=1
72 IF(ICC.NE.2) WRITE(IDLG,73)
73 FORMAT('+? ',$)
READ(ICC,74,END=60) INPUT
74 FORMAT(50A1)
IF(INPUT(1).EQ.'!') RETURN
IF(INPUT(1).EQ.' ') GO TO 60
IF((INPUT(1).EQ.'A').AND.(INPUT(2).EQ.'U').AND.(INPUT(3).EQ.'T')
1.AND.(INPUT(4).EQ.'O')) GO TO 100
IF((INPUT(1).NE.'H').OR.(INPUT(2).NE.'E').OR.(INPUT(3).NE.'P')
1.OR.(INPUT(4).NE.'P')) GO TO 65
WRITE(IDLG,75)
75 FORMAT('0ENTER RANGES 1 PER LINE IN RESPONSE TO A QUESTION MARK.'/
1' EACH RANGE IS COMPRISED OF A MINIMUM, A MAXIMUM, AND AN'/
2' OPTIONAL RANGE NAME (UP TO 5 CHARACTERS). EACH ENTRY IN'/
3' THE RANGE IS SEPARATED FROM THE NEXT ENTRY BY A COMMA.'/
4' FOR EXAMPLE: 1,7,LOW; THE MINIMUM IS 1, THE MAXIMUM FOR THE'/
5' RANGE IS 7, AND THE NAME OF THAT GROUPING IS LOW. WHEN THE'/
6' LAST RANGE HAS BEEN ENTERED TYPE A CONTROL Z (^Z) OR A BLANK'/
7' LINE.'//)
GO TO 72
C
C MINIMUM OF RANGE
C
65 DO 76 J=1,15
76 IVALUE(J)=' '
K=1
J=1
77 IF(INPUT(J).EQ.' ') GO TO 80
IF(INPUT(J).EQ.',') GO TO 82
IF(INPUT(J).EQ.'E') GO TO 79
IF(INPUT(J).EQ.'-') GO TO 79
IF(INPUT(J).EQ.'.') GO TO 79
IF((INPUT(J).LE.'9').AND.(INPUT(J).GE.'0')) GO TO 79
WRITE(IDLG,78) INPUT(J)
78 FORMAT(' ILLEGAL CHARACTER "',A1,'" IN MINIMUM OF RANGE'/)
GO TO 72
79 IF(K.GT.15) GO TO 80
IVALUE(K)=INPUT(J)
K=K+1
J=J+1
GO TO 77
80 WRITE(IDLG,81)
81 FORMAT(' NO COMMA BETWEEN MINIMUM AND MAXIMUM OR EXTRA SPACES'/)
GO TO 72
82 ENCODE(15,74,SPACE) IVALUE
DECODE(15,83,SPACE) RANGE(1,I)
83 FORMAT(G)
DO 84 K=1,15
84 IVALUE(K)=' '
K=1
J=J+1
C
C MAXIMUM OF RANGE
C
85 IF(INPUT(J).EQ.' ') GO TO 90
IF(INPUT(J).EQ.',') GO TO 90
IF(INPUT(J).EQ.'E') GO TO 87
IF(INPUT(J).EQ.'-') GO TO 87
IF(INPUT(J).EQ.'.') GO TO 87
IF((INPUT(J).LE.'9').AND.(INPUT(J).GE.'0')) GO TO 87
WRITE(IDLG,86) INPUT(J)
86 FORMAT(' ILLEGAL CHARACTER "',A1,'" IN MAXIMUM OF RANGE'/)
GO TO 72
87 IF(K.GT.15) GO TO 88
IVALUE(K)=INPUT(J)
K=K+1
J=J+1
GO TO 85
88 WRITE(IDLG,89)
89 FORMAT(' NO COMMA BETWEEN MAXIMUM AND NAME OR TOO MANY CHAR.'/)
GO TO 72
90 ENCODE(15,74,SPACE) IVALUE
DECODE(15,83,SPACE) RANGE(2,I)
IF(INPUT(J).NE.',') GO TO 91
IF(INPUT(J+1).EQ.' ') GO TO 91
ENCODE(5,74,RANGE(3,I))(INPUT(K),K=J+1,J+5)
GO TO 93
91 ENCODE(5,92,RANGE(3,I)) I
92 FORMAT('GRP',I2)
93 I=I+1
IF(I.LE.25) GO TO 72
WRITE(IDLG,94)
94 FORMAT(' MAXIMUM OF 25 RANGES')
GO TO 60
60 NBRK=I-1
GO TO 108
C
C DISCR OR AUTO USED
C
100 I=1
J=2
RANGE(1,I)=AV(1)
RANGE(2,I)=AV(1)
ENCODE(5,92,RANGE(3,I)) I
101 IF(RANGE(1,I).EQ.AV(J)) GO TO 102
I=I+1
IF(I.GT.25) GO TO 103
RANGE(1,I)=AV(J)
RANGE(2,I)=AV(J)
ENCODE(5,92,RANGE(3,I)) I
102 J=J+1
IF(J.LE.NC) GO TO 101
NBRK=I
GO TO 105
103 XNC=(AV(NC)-AV(1))/25.
RANGE(1,I)=AV(1)
DO 104 I=1,24
RANGE(2,I)=RANGE(1,I)+XNC
RANGE(1,I+1)=RANGE(2,I)
104 CONTINUE
RANGE(2,25)=AV(NC)
NBRK=25
C
C OK RANGES ESTABLISHED SEE IF USER WANTS LIST
C
105 IF(OPT(2).NE.1) GO TO 108
WRITE(IDLG,95)
95 FORMAT(' RANGES USED:')
DO 106 I=1,NBRK
106 WRITE(IDLG,107) RANGE(1,I),RANGE(2,I),RANGE(3,I)
107 FORMAT(1X,G,',',G,',',A5)
C
C NOW ARRANGE GROUPS
C
108 DO 109 I=1,NBRK
109 NGRP(I)=0
DO 110 I=1,NC
DO 111 J=1,NBRK
IF(AV(I).LT.RANGE(1,J)) GO TO 111
IF(AV(I).GT.RANGE(2,J)) GO TO 111
AV(I)=J
NGRP(J)=NGRP(J)+1
GO TO 110
111 CONTINUE
115 AV(I)=99
110 CONTINUE
CALL SORT(AV,JV,NC)
I=1
112 IF(NGRP(I).NE.0) GO TO 117
WRITE(IDLG,113) RANGE(3,I)
113 FORMAT(' RANGE "',A5,'" HAD NO OBSERVATIONS - DISCARDED')
IF(I.EQ.NBRK) GO TO 116
DO 114 J=I+1,NBRK
NGRP(J-1)=NGRP(J)
DO 114 K=1,3
114 RANGE(K,J-1)=RANGE(K,J)
116 NBRK=NBRK-1
IF(I.LE.NBRK) GO TO 112
GO TO 118
117 I=I+1
IF(I.LE.NBRK) GO TO 112
118 IF(NBRK.GT.1) GO TO 120
WRITE(IDLG,119)
119 FORMAT(' ONLY 1 GROUP - NO DISCRIMINATE ANALYSIS PERFORMED')
RETURN
C
C OK NOW GET VARIABLES TO BE USED.
C
120 K=1
DO 121 I=1,NN
IV(I)=IVAR(I)
IF(IVAR(I).GT.0) GO TO 121
IV(I)=K
K=K+1
121 CONTINUE
GO TO 135
130 J=NN
131 IF(IVAR(J).GT.0) GO TO 132
IV(J)=IV(J)+1
IF(IV(J).LE.NV) GO TO 133
132 J=J-1
IF(J.GE.1) GO TO 131
RETURN
133 K=IV(J)
IF(J.EQ.NN) GO TO 135
DO 134 I=J+1,NN
IF(IVAR(I).GT.0) GO TO 134
K=K+1
IF(K.GT.NV) GO TO 132
IV(I)=K
134 CONTINUE
135 DO 136 I=1,NN-1
DO 136 K=I+1,NN
IF(IV(I).EQ.IV(K)) GO TO 130
136 CONTINUE
C
C BEGIN ANALYSIS
C
IF(IOUT.NE.21) WRITE(IOUT,5566) (HEDR(J),J=1,NSZ)
5566 FORMAT('1',70A1)
IF(IOUT.EQ.21) CALL PRNTHD
WRITE(IOUT,137) NAMES(IBREAK),(NAMES(IV(I)),I=1,NN)
137 FORMAT('0',15X,'***** DISCRIMINATE ANALYSIS *****'/
1'0VARIABLE ',A5,' IS USED TO DETERMINE GROUPINGS. THE '/
2' VARIABLES BEING ANALYSED ARE:',7(1X,A5)/11(1X,A5))
LINES=9
LINPO=(NN+7)/8+1
IF(NN.EQ.8) LINPO=3
LINPOB=(NBRK+7)/8+1
IF(NBRK.EQ.8) LINPOB=3
DO 140 J=1,NN
XC(J)=0
DO 141 K=1,NN
141 D(J,K)=0.
DO 142 K=1,NBRK
142 XBAR(J,K)=0.
140 CONTINUE
N=0
NTOT=0
DO 150 I=1,NBRK
NTOT=NTOT+NGRP(I)
DO 151 J=1,NN
DO 151 K=1,NN
151 C(J,K)=0.
XN=NGRP(I)
DO 152 J=N+1,N+NGRP(I)
JJ=JV(J)
DO 153 K=1,NN
X=DATA(JJ,IV(K))
XBAR(K,I)=XBAR(K,I)+X
XC(K)=XC(K)+X
DO 154 L=K,NN
154 C(K,L)=C(K,L)+X*DATA(JJ,IV(L))
153 CONTINUE
152 CONTINUE
DO 155 K=1,NN
155 XBAR(K,I)=XBAR(K,I)/NGRP(I)
DO 156 K=1,NN
DO 156 L=K,NN
C(K,L)=C(K,L)-XN*XBAR(K,I)*XBAR(L,I)
156 C(L,K)=C(K,L)
IF(OPT(7).EQ.0) GO TO 165
IF(IOUT.NE.21) GO TO 400
LINES=LINES+4+LINPO
IF(LINES.LE.(LINPP-LINPO)) GO TO 400
CALL PRNTHD
LINES=6+LINPO
400 WRITE(IOUT,157) RANGE(3,I)
157 FORMAT(//'0MATRIX OF CROSS PRODUCTS OF DEVIATIONS FROM THE MEANS'
1,' FOR GROUP: ',A5)
IF(IOUT.NE.21) GO TO 401
WRITE(IOUT,160) (NAMES(IV(K)),K=1,NN)
160 FORMAT('0',7X,8(2X,A5,8X)/(8X,8(2X,A5,8X)))
GO TO 402
401 WRITE(IOUT,161) (NAMES(IV(K)),K=1,NN)
161 FORMAT('0',7X,5(2X,A5,6X)/(8X,5(2X,A5,6X)))
402 DO 158 K=1,NN
IF(IOUT.NE.21) GO TO 404
LINES=LINES+LINPO
IF(LINES.LE.LINPP) GO TO 403
CALL PRNTHD
WRITE(IOUT,160)(NAMES(IV(KLP)),KLP=1,NN)
LINES=2+2*LINPO
403 WRITE(IOUT,159) NAMES(IV(K)),(C(K,L),L=1,NN)
GO TO 158
159 FORMAT('0',A5,2X,8G15.7/(8X,8G15.7))
404 WRITE(IOUT,162) NAMES(IV(K)),(C(K,L),L=1,NN)
162 FORMAT('0',A5,1X,5G13.7/(7X,5G13.7))
158 CONTINUE
165 DO 166 K=1,NN
DO 166 L=1,NN
166 D(K,L)=D(K,L)+C(K,L)
N=N+NGRP(I)
150 CONTINUE
XN=N-NBRK
IF(OPT(3).NE.1) GO TO 172
IF(IOUT.NE.21) GO TO 405
LINES=LINES+4+LINPOB
IF(LINES.LE.(LINPP-LINPOB)) GO TO 405
CALL PRNTHD
LINES=6+LINPOB
405 WRITE(IOUT,171)
171 FORMAT(//'0','MEANS')
IF(IOUT.NE.21) GO TO 407
406 WRITE(IOUT,160)(RANGE(3,K),K=1,NBRK)
GO TO 408
407 WRITE(IOUT,161)(RANGE(3,K),K=1,NBRK)
408 DO 170 I=1,NN
IF(IOUT.NE.21) GO TO 410
LINES=LINES+LINPOB
IF(LINES.LE.LINPP) GO TO 409
CALL PRNTHD
WRITE(IOUT,160)(RANGE(3,K),K=1,NBRK)
LINES=2+2*LINPOB
409 WRITE(IOUT,159) NAMES(IV(I)),(XBAR(I,J),J=1,NBRK)
GO TO 170
410 WRITE(IOUT,162) NAMES(IV(I)),(XBAR(I,J),J=1,NBRK)
170 CONTINUE
172 DO 175 K=1,NN
DO 175 L=1,NN
175 C(K,L)=D(K,L)/XN
IF(OPT(4).NE.1) GO TO 176
IF(IOUT.NE.21) GO TO 411
LINES=LINES+4+LINPO
IF(LINES.LE.(LINPP-LINPO)) GO TO 411
CALL PRNTHD
LINES=6+LINPO
411 WRITE(IOUT,173)
173 FORMAT(//'0',5X,'DISPERSION MATRIX')
IF(IOUT.NE.21) GO TO 413
412 WRITE(IOUT,160)(NAMES(IV(K)),K=1,NN)
GO TO 414
413 WRITE(IOUT,161)(NAMES(IV(K)),K=1,NN)
414 DO 174 I=1,NN
IF(IOUT.NE.21) GO TO 416
LINES=LINES+LINPO
IF(LINES.LE.LINPP) GO TO 415
CALL PRNTHD
WRITE(IOUT,160)(NAMES(IV(K)),K=1,NN)
LINES=2+2*LINPO
415 WRITE(IOUT,159) NAMES(IV(I)),(C(I,J),J=1,NN)
GO TO 174
416 WRITE(IOUT,162) NAMES(IV(I)),(C(I,J),J=1,NN)
174 CONTINUE
C
C
C INVERSE OF DISPERSION MATRIX
C
176 DO 202 I=1,NN
DO 201 J=1,NN
201 D(I,J)=0.
202 D(I,I)=1.0
DO 210 I=1,NN
IF(((C(I,I)+100.)-100.).NE.0.0) GO TO 220
IF(I.EQ.N) GO TO 300
DO 211 J=I+1,NN
IF(((C(J,I)+100.)-100.).NE.0.0) GO TO 212
211 CONTINUE
GO TO 300
212 DO 213 K=1,NN
D(I,K)=D(I,K)+D(J,K)
213 C(I,K)=C(I,K)+C(J,K)
220 G=C(I,I)
DO 221 J=1,NN
D(I,J)=D(I,J)/G
221 C(I,J)=C(I,J)/G
DO 230 L=1,NN
IF(L.EQ.I) GO TO 230
G=C(L,I)
DO 231 J=1,NN
D(L,J)=D(L,J)-G*D(I,J)
231 C(L,J)=C(L,J)-G*C(I,J)
230 CONTINUE
210 CONTINUE
C
C
IF(OPT(5).NE.1) GO TO 242
IF(IOUT.NE.21) GO TO 417
LINES=LINES+4+LINPO
IF(LINES.LE.(LINPP-LINPO)) GO TO 417
CALL PRNTHD
LINES=6+LINPO
417 WRITE(IOUT,240)
240 FORMAT(//'0','INVERSE OF DISPERSION MATRIX')
IF(IOUT.NE.21) GO TO 419
418 WRITE(IOUT,160)(NAMES(IV(K)),K=1,NN)
GO TO 420
419 WRITE(IOUT,161)(NAMES(IV(K)),K=1,NN)
420 DO 241 I=1,NN
IF(IOUT.NE.21) GO TO 422
LINES=LINES+LINPO
IF(LINES.LE.LINPP) GO TO 421
CALL PRNTHD
WRITE(IOUT,160)(NAMES(IV(K)),K=1,NN)
LINES=2+2*LINPO
421 WRITE(IOUT,159) NAMES(IV(I)),(D(I,J),J=1,NN)
GO TO 241
422 WRITE(IOUT,162) NAMES(IV(I)),(D(I,J),J=1,NN)
241 CONTINUE
242 DO 243 I=1,NN
243 XC(I)=XC(I)/NTOT
DSQ=0.
DO 244 I=1,NN
DO 244 J=1,NN
DS=0
DO 245 K=1,NBRK
245 DS=DS+NGRP(K)*(XBAR(I,K)-XC(I))*(XBAR(J,K)-XC(J))
244 DSQ=DSQ+DS*D(I,J)
IDF=NN*(NBRK-1)
PROB=FISHER(IDF,1000,DSQ)
IF(IOUT.NE.21) GO TO 423
LINES=LINES+8
IF(LINES.LE.LINPP) GO TO 423
CALL PRNTHD
LINES=10
423 WRITE(IOUT,246) DSQ,DSQ,IDF,NBRK,NN,PROB
246 FORMAT(//'0GENERALIZED MAHALANOBIS D- SQUARE =',G/
1'0USING THE VALUE ',G,' AS A CHI-SQUARE WITH ',I4/
1' DEGREES OF FREEDOM TO TEST THE HYPOTHESIS THAT THE MEAN'/
1' VALUES ARE THE SAME IN ALL ',I3,' GROUPS FOR THESE ',I3/
4' VARIABLES, THE PROBABILITY IS ',F5.2)
DO 250 I=1,NBRK
DSS=0
DO 251 J=1,NN
DS=0
DO 252 K=1,NN
DS=DS+D(J,K)*XBAR(K,I)
252 DSS=DSS+D(J,K)*XBAR(J,I)*XBAR(K,I)
251 C(J,I)=DS
250 CON(I)=DSS*-.5
C** 250+1
IF(IOUT.NE.21) GO TO 424
LINES=LINES+4+LINPOB
IF(LINES.LE.(LINPP-LINPOB)) GO TO 424
CALL PRNTHD
LINES=6+LINPOB
424 WRITE(IOUT,260)
260 FORMAT(//'0FUNCTION COEFFICIENTS')
IF(IOUT.NE.21) GO TO 425
WRITE(IOUT,160)(RANGE(3,K),K=1,NBRK)
GO TO 426
425 WRITE(IOUT,161)(RANGE(3,K),K=1,NBRK)
426 DO 261 I=1,NN
IF(IOUT.NE.21) GO TO 428
LINES=LINES+LINPOB
IF(LINES.LE.LINPP) GO TO 427
CALL PRNTHD
WRITE(IOUT,160)(RANGE(3,K),K=1,NBRK)
LINES=2+2*LINPOB
427 WRITE(IOUT,159) NAMES(IV(I)),(C(I,J),J=1,NBRK)
GO TO 261
428 WRITE(IOUT,162) NAMES(IV(I)),(C(I,J),J=1,NBRK)
261 CONTINUE
LINES=LINES+1+LINPOB
IF(LINES.LE.LINPP) GO TO 429
CALL PRNTHD
WRITE(IOUT,160)(RANGE(3,K),K=1,NBRK)
LINES=2+2*LINPOB
429 WRITE(IOUT,262)
262 FORMAT(/)
CONST='CONST'
IF(IOUT.EQ.21) WRITE(IOUT,159) CONST,(CON(I),I=1,NBRK)
IF(IOUT.NE.21) WRITE(IOUT,162) CONST,(CON(I),I=1,NBRK)
IF((OPT(6).NE.1).AND.(OPT(8).NE.1)) GO TO 130
IF(OPT(6).NE.1) GO TO 430
LINES=LINES+9
LINPOB=(NBRK+9)/10
IF(NBRK.EQ.10) LINPOB=2
IF(LINES.LE.(LINPP-2-LINPOB)) GO TO 431
CALL PRNTHD
LINES=11
431 WRITE(IOUT,272)
272 FORMAT(//'0','EVALUATION OF CLASSIFICATION FUNCTIONS',
1' FOR EACH CASE')
WRITE(IOUT,285)
285 FORMAT('0',9X,'FUNCTION'/12X,'WITH'/10X,'LARGEST',4X,'LARGEST'/
32X,'OBS.',6X,'PROB.',6X,'PROB.')
430 DO 270 I=1,NBRK
DO 270 J=1,NBRK
270 IP(I,J)=0
N=0
DO 271 I=1,NBRK
IF(IOUT.NE.21) GO TO 432
LINES=LINES+2
IF(LINES.LE.(LINPP-LINPOB)) GO TO 432
CALL PRNTHD
WRITE(IOUT,285)
LINES=9
432 WRITE(IOUT,273) RANGE(3,I)
273 FORMAT('0GROUP: ',A5)
DO 274 J=1+N,NGRP(I)+N
M=JV(J)
XX=-1.E30
DO 275 K=1,NBRK
E=CON(K)
DO 276 L=1,NN
276 E=E+C(L,K)*DATA(M,IV(L))
IF(E.LE.XX) GO TO 275
XX=E
IX=K
275 XVAL(K)=E
DS=0
DO 277 JJ=1,NBRK
XVAL(JJ)=EXP(XVAL(JJ)-XX)
277 DS=DS+XVAL(JJ)
XX=1./DS
DO 278 JJ=1,NBRK
278 XVAL(JJ)=XVAL(JJ)/DS
IF(OPT(6).NE.1) GO TO 274
IF(IOUT.NE.21) GO TO 434
LINES=LINES+LINPOB
IF(LINES.LE.LINPP) GO TO 433
CALL PRNTHD
WRITE(IOUT,285)
LINES=7+LINPOB
433 IF(IOUT.EQ.21) WRITE(IOUT,279) JV(J),IX,XX,(XVAL(JJ),JJ=1,NBRK)
279 FORMAT(1X,I5,7X,I2,4X,F8.5,4X,10(2X,F8.5)/(31X,10(2X,F8.5)))
GO TO 274
434 WRITE(IOUT,280) JV(J),IX,XX,(XVAL(JJ),JJ=1,NBRK)
280 FORMAT(1X,I5,7X,I2,4X,F8.5,4X,5(1X,F7.4)/(31X,5(F7.4)))
274 IP(I,IX)=IP(I,IX)+1
N=N+NGRP(I)
271 CONTINUE
IF(OPT(8).NE.1) GO TO 130
LINPOB=(NBRK+19)/20
IF(NBRK.EQ.20) LINPOB=2
IF(IOUT.NE.21) GO TO 435
LINES=LINES+7+LINPOB
IF(LINES.LE.(LINPP-LINPOB-1)) GO TO 435
CALL PRNTHD
LINES=9+LINPOB
435 WRITE(IOUT,284)
284 FORMAT(//'0CLASSIFICATION MATRIX'/8X,'GROUP PLACED IN AS A',
1' RESULT OF EVALUATING THE FUNCTION')
IF(IOUT.EQ.21) WRITE(IOUT,281)(RANGE(3,I),I=1,NBRK)
281 FORMAT('0GROUP ',20(1X,A5)/(8X,20(1X,A5)))
IF(IOUT.NE.21) WRITE(IOUT,282) (RANGE(3,I),I=1,NBRK)
282 FORMAT('0GROUP ',10(1X,A5)/(8X,10(1X,A5)))
WRITE(IOUT,283)
283 FORMAT(1X)
DO 290 I=1,NBRK
IF(IOUT.NE.21) GO TO 437
LINES=LINES+LINPOB
IF(LINES.LE.LINPP) GO TO 436
CALL PRNTHD
WRITE(IOUT,281)(RANGE(3,J),J=1,NBRK)
WRITE(IOUT,283)
LINES=4+LINPOB*2
436 WRITE(IOUT,291) RANGE(3,I),(IP(I,J),J=1,NBRK)
291 FORMAT(1X,A5,2X,20(1X,I5)/(8X,20(1X,I5)))
GO TO 290
437 WRITE(IOUT,292) RANGE(3,I),(IP(I,J),J=1,NBRK)
292 FORMAT(1X,A5,2X,10(1X,I5)/(8X,10(1X,I5)))
290 CONTINUE
GO TO 130
300 WRITE(IDLG,301)
IF(IOUT.EQ.21) WRITE(IOUT,301)
301 FORMAT(' INVERSE OF DISPERSION MATRIX DOES NOT EXIST')
GO TO 130
END
SUBROUTINE SORT(AV,SP,NC)
DIMENSION SP(1),AV(1),IL(20),IU(20)
M=1
II=1
J=NC
71 IF(II.GE.J) GO TO 78
72 K=II
IJ=(J+II)/2
T=AV(IJ)
S=SP(IJ)
IF(AV(II).GT.T) GO TO 83
IF(AV(II).LT.T) GO TO 73
IF(SP(II).LE.S) GO TO 73
83 AV(IJ)=AV(II)
AV(II)=T
T=AV(IJ)
SP(IJ)=SP(II)
SP(II)=S
S=SP(IJ)
73 LL=J
IF(AV(J).LT.T) GO TO 84
IF(AV(J).GT.T) GO TO 75
IF(SP(J).GE.S) GO TO 75
84 AV(IJ)=AV(J)
AV(J)=T
T=AV(IJ)
SP(IJ)=SP(J)
SP(J)=S
S=SP(IJ)
IF(AV(II).GT.T) GO TO 85
IF(AV(II).LT.T) GO TO 75
IF(SP(II).LE.S) GO TO 75
85 AV(IJ)=AV(II)
AV(II)=T
T=AV(IJ)
SP(IJ)=SP(II)
SP(II)=S
S=SP(IJ)
GO TO 75
74 AV(LL)=AV(K)
AV(K)=TT
SP(LL)=SP(K)
SP(K)=SS
75 LL=LL-1
IF(AV(LL).GT.T) GO TO 75
IF((AV(LL).EQ.T).AND.(SP(LL).GT.S)) GO TO 75
TT=AV(LL)
SS=SP(LL)
76 K=K+1
IF(AV(K).LT.T) GO TO 76
IF((AV(K).EQ.T).AND.(SP(K).LT.S)) GO TO 76
IF(K.LE.LL) GO TO74
IF((LL-II).LE.(J-K)) GO TO 77
IL(M)=II
IU(M)=LL
II=K
M=M+1
GO TO 79
77 IL(M)=K
IU(M)=J
J=LL
M=M+1
GO TO 79
78 M=M-1
IF(M.EQ.0) RETURN
II=IL(M)
J=IU(M)
79 IF((J-II).GE.11) GO TO 72
IF(II.EQ.1) GO TO 71
II=II-1
80 II=II+1
IF(II.EQ.J) GO TO 78
S=SP(II+1)
T=AV(II+1)
IF(AV(II).LT.T) GO TO 80
IF((AV(II).EQ.T).AND.(SP(II).LE.S)) GO TO 80
K=II
81 AV(K+1)=AV(K)
SP(K+1)=SP(K)
K=K-1
IF(T.LT.AV(K)) GO TO 81
IF((T.EQ.AV(K)).AND.(S.LT.SP(K))) GO TO 81
AV(K+1)=T
SP(K+1)=S
GO TO 80
END