1 SUBROUTINE MHIST(MM, M, N, A, CLAB, RLAB, TITLE, XMIN, XMAX, NV, 2 * NDIAGS, DMWORK, WORK, DMIWRK, IWORK, IWORK2, XLL, 3 * OUNIT) 4C 5C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> 6C 7C PURPOSE 8C ------- 9C 10C PRODUCES MULTIVARIATE HISTOGRAMS 11C 12C DESCRIPTION 13C ----------- 14C 15C 1. THE DIAGRAMS ARE COMPUTED SEQUENTIALLY. THE FIRST DIAGRAM 16C BEGINS WITH THE SMALLEST RANGE FOR EACH VARIABLE WHICH CONTAINS 17C EVERY OBSERVATION. THEN THE RANGES FOR EACH VARIABLE ARE 18C VARIED AND THE NEW RANGES ARE ACCEPTED IF THEY INCREASE THE 19C LOG-LIKELIHOOD OF THE OBSERVATIONS WITHIN THOSE RANGES. EACH 20C RANGE IS VARIED UNTIL NO FURTHER CHANGE INCREASES THE 21C LIKELIHOOD. THE FIRST DIAGRAM IS COMPLETED AND PRINTED AND THE 22C OBSERVATIONS WITHIN THE FINAL RANGES ARE DELETED AND THE VOLUME 23C OF THE BLOCKS ARE IGNORED IN LATER DIAGRAM CALCULATIONS. A 24C SECOND DIAGRAM, AGAIN STARTING WITH THE SMALLEST RANGE FOR EACH 25C VARIABLE WHICH CONTAINS EVERY REMAINING OBSERVATION, IS ADDED 26C AND VARIED TO MAXIMIZE THE LIKELIHOOD OF BOTH BLOCKS. THE 27C OTHER DIAGRAMS ARE COMPUTED SIMILARLY. 28C 29C 2. THE OUTPUT IS ON FORTRAN UNIT OUNIT AND CONSISTS OF THE SET OF 30C DIAGRAMS. FOR EACH DIAGRAM, THE RANGES FOR EACH VARIABLE ARE 31C ENCLOSED IN A RECTANGULAR BLOCK. TO THE RIGHT OF EACH DIAGRAM 32C IS A RECTANGLE WHOSE LENGTH IN THE HORIZONTAL DIRECTION IS THE 33C NUMBER OF CASES WITHIN THE SET OF RANGES FOR THAT PARTICULAR 34C DIAGRAM. THE SUM OF THE HORIZONTAL DISTANCES FOR ALL DIAGRAMS 35C SHOULD ADD UP TO THE TOTAL NUMBER OF CASES. 36C 37C INPUT PARAMETERS 38C ---------------- 39C 40C MM INTEGER SCALAR (UNCHANGED ON OUTPUT). 41C THE FIRST DIMENSION OF THE MATRIX A. MUST BE AT LEAST M. 42C 43C M INTEGER SCALAR (UNCHANGED ON OUTPUT). 44C THE NUMBER OF CASES. 45C 46C N INTEGER SCALAR (UNCHANGED ON OUTPUT). 47C THE NUMBER OF VARIABLES. 48C 49C A REAL MATRIX WHOSE FIRST DIMENSION MUST BE MM AND WHOSE SECOND 50C DIMENSION MUST BE AT LEAST N (UNCHANGED ON OUTPUT). 51C THE MATRIX OF DATA VALUES. 52C 53C A(I,J) IS THE VALUE FOR THE J-TH VARIABLE FOR THE I-TH CASE. 54C 55C CLAB VECTOR OF 4-CHARACTER VARIABLES DIMENSIONED AT LEAST N. 56C (UNCHANGED ON OUTPUT). 57C THE LABELS OF THE VARIABLES. 58C 59C RLAB VECTOR OF 4-CHARACTER VARIABLES DIMENSIONED AT LEAST M. 60C (UNCHANGED ON OUTPUT). 61C THE LABELS OF THE CASES. 62C 63C TITLE 10-CHARACTER VARIABLE (UNCHANGED ON OUTPUT). 64C TITLE OF THE DATA SET. 65C 66C XMIN INTEGER VECTOR DIMENSIONED AT LEAST N (UNCHANGED ON OUTPUT). 67C XMIN(I) IS THE MINIMUM VALUE TO BE PLOTTED FOR VARIABLE I. 68C 69C XMAX INTEGER VECTOR DIMENSIONED AT LEAST N (UNCHANGED ON OUTPUT). 70C XMAX(I) IS THE MAXIMUM VALUE TO BE PLOTTED FOR VARIABLE I. 71C 72C IF XMIN(I) .GE. XMAX(I) ON INPUT, THEIR VALUES WILL BE 73C DETERMINED BY THE ROUTINE. 74C 75C NV INTEGER VECTOR DIMENSIONED AT LEAST N (UNCHANGED ON OUTPUT). 76C VECTOR DEFINING THE ORDER OF THE VARIABLES. 77C 78C NDIAGS INTEGER SCALAR (UNCHANGED ON OUTPUT). 79C THE NUMBER OF DIAGRAMS. 80C 81C DMWORK INTEGER SCALAR (UNCHANGED ON OUTPUT). 82C THE FIRST DIMENSION OF THE MATRIX WORK. MUST BE AT LEAST 83C 2*NDIAGS. 84C 85C WORK REAL MATRIX WHOSE FIRST DIMENSION MUST BE DMWORK AND WHOSE 86C SECOND DIMENSION MUST BE AT LEAST N+1. 87C WORK MATRIX. 88C 89C DMIWRK INTEGER SCALAR (UNCHANGED ON OUTPUT). 90C THE FIRST DIMENSION OF THE MATRIX IWORK. MUST BE AT LEAST 91C 2*NDIAGS. 92C 93C IWORK INTEGER MATRIX WHOSE FIRST DIMENSION MUST BE DMIWRK AND WHOSE 94C SECOND DIMENSION MUST BE AT LEAST N. 95C WORK MATRIX. 96C 97C IWORK2 INTEGER VECTOR DIMENSIONED AT LEAST M. 98C WORK VECTOR. 99C 100C OUNIT INTEGER SCALAR (UNCHANGED ON OUTPUT). 101C UNIT NUMBER FOR OUTPUT. 102C 103C OUTPUT PARAMETER 104C ---------------- 105C 106C XLL REAL SCALAR (CHANGED ON OUTPUT). 107C THE LOG LIKELIHOOD OF THE FINAL SET OF BLOCKS. 108C 109C REFERENCES 110C ---------- 111C 112C HARTIGAN, J. A. (1975). CLUSTERING ALGORITHMS, JOHN WILEY & 113C SONS, INC., NEW YORK. PAGES 40, 54-55. 114C 115C HARTIGAN, J. A. (1975) PRINTER GRAPHICS FOR CLUSTERING. JOURNAL OF 116C STATISTICAL COMPUTATION AND SIMULATION. VOLUME 4,PAGES 187-213. 117C 118C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> 119C 120 INTEGER DMWORK, DMIWRK, OUNIT 121 DIMENSION A(MM,*), XMIN(*), XMAX(*), NV(*), WORK(DMWORK,*), 122 * IWORK(DMIWRK,*), IWORK2(*) 123 CHARACTER*4 RLAB(*), CLAB(*) 124 CHARACTER*10 TITLE 125 CHARACTER*1 BL, X, H, P, AA(112) 126 DATA BL,X,H,P/' ','!','_','.'/ 127C 128 DO 10 I = 1 , N 129 IF (XMIN(I) .GE. XMAX(I))CALL RANGE(M, A(1,I), XMIN(I), XMAX(I)) 130 DO 10 L = 1 , NDIAGS 131 WORK(L,I) = XMIN(I) 132 WORK(NDIAGS+L,I) = XMAX(I) 133 10 CONTINUE 134C 135C COMPUTE OPTIMAL BLOCKS 136C 137 CALL TRY1(MM,M,N,A,NDIAGS,NV,DMWORK,WORK,XLL,WORK(1,N+1)) 138C 139C OUTPUT BLOCKS 140C 141 IF (OUNIT .GT. 0) 142 * WRITE(OUNIT,1) TITLE,(CLAB(J),XMIN(J),XMAX(J),J=1,N) 143 1 FORMAT('1 MODAL HISTOGRAM OF ',A10/ 144 *' INTERVALS ...! !.... FOR EACH VARIABLE DEFINE MODAL BLOCKS '/ 145 *' RECTANGLES TO RIGHT HAVE AREA PROPORTIONAL TO NUMBER IN MODE '// 146 *' VARIABLE RANGES '/(1X,A4,F9.3,2X,16('.'),1X,F10.3)) 147 DO 20 J=1,N 148 DO 20 K=1,NDIAGS 149 IWORK(K,J)=INT(15*(WORK(K,J)-XMIN(J))/(XMAX(J)-XMIN(J)))+1 150 IWORK(NDIAGS+K,J)=INT(15*(WORK(NDIAGS+K,J)-XMIN(J))/(XMAX(J) 151 * -XMIN(J)))+1 152 20 CONTINUE 153 DO 30 I=1,M 154 30 IWORK2(I)=0 155 DO 130 K=1,NDIAGS 156C 157C COUNT NUMBER OF CASES WHOSE VALUES ARE COMPLETELY WITHIN BLOCK 158C THAT HAVEN'T BEEN COUNTED BEFORE 159C 160 NK=0 161 DO 50 I=1,M 162 IF(IWORK2(I).EQ.0) THEN 163 DO 40 J=1,N 164 40 IF(A(I,J).GT.WORK(NDIAGS+K,J).OR.A(I,J).LT.WORK(K,J)) 165 * GOTO 50 166 IWORK2(I)=1 167 NK=NK+1 168 ENDIF 169 50 CONTINUE 170C 171C GENERATE HEIGHT 172C 173 NK=NK+19 174 DO 60 I=1,112 175 60 AA(I)=BL 176 DO 70 I=3,NK 177 70 IF(I.NE.17.AND.I.NE.18.AND.I.NE.NK) AA(I)=H 178 IF (OUNIT .GT. 0) WRITE(OUNIT,2) (AA(I),I=1,112) 179 2 FORMAT(15X,112A1) 180 DO 80 I=19,NK 181 80 AA(I)=BL 182C 183C WORK ON PROFILE 184C 185 DO 110 JJ=1,N 186 J=NV(JJ) 187 DO 90 I=1,18 188 90 AA(I)=BL 189 AA(19)=X 190 AA(NK)=X 191 K1=IWORK(K,J)+1 192 K2=IWORK(NDIAGS+K,J)+1 193 AA(K1)=X 194 AA(K2)=X 195 J1=JJ+1 196 IF(JJ.EQ.N) J1=JJ 197 J1=NV(J1) 198 K3=IWORK(K,J1)+1 199 K4=IWORK(NDIAGS+K,J1)+1 200 DO 100 KK=2,17 201 IF(KK.NE.K1.AND.KK.NE.K2) THEN 202 IF(KK.LT.K1) AA(KK)=P 203 IF(KK.GT.K2) AA(KK)=P 204 IF((KK-K3)*(KK-K1).LT.0) AA(KK)=H 205 IF((KK-K4)*(KK-K2).LT.0) AA(KK)=H 206 ENDIF 207 100 CONTINUE 208C 209C PRINT A LINE 210C 211 IF (OUNIT .GT. 0) WRITE(OUNIT,3) CLAB(J),(AA(I),I=1,112) 212 3 FORMAT(10X,A4,1X,112A1) 213 110 CONTINUE 214 DO 120 I=3,NK 215 120 IF(I.NE.17.AND.I.NE.18.AND.I.NE.NK) AA(I)=H 216 130 IF (OUNIT .GT. 0) WRITE(OUNIT,4)(AA(I),I=1,112) 217 4 FORMAT('+',14X,112A1) 218 RETURN 219 END 220