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