C00000000111111111122222222223333333333444444444455555555556666666666777 C23456789012345678901234567890123456789012345678901234567890123456789012 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C CLUSPAC: Computer Programs for Mixture-Model Clustering C C C C COPYRIGHT (C) 199 Stanley Louis Sclove C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C CLUSTER ANALYSIS OF CASES BY MIXTURE MODEL C C C C PROGRAM MIXPDTA CLUSPAC: C C C C VARYING COVARIANCE MATRICES C C AUTOMATIC SETTING OF INITIAL PARAMETER ESTIMATES C C C C Programs are organized as C C programs within packages within the suite. C C PROGRAM: MIXPDTA PACKAGE: CLUSPAC SUITE: ISOPAC C C C C DOCUMENTATION: C C www.uic.edu/classes/idsc/ids594/ISOPAC/home.htm C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C MIXPDTA CLUSPAC C C VERSION 3.1 2003: Mar 24 C with IVERB for levels of verbosity in output C C VERSION 3.0 2003: Mar 05 C Subroutine INVS instead of MATEQ C C VERSION 2.99 2002: Jun 11 C VERSION 2.98 1994: Feb 15 C C C C PROGRAMMED BY: C C C C DR. STANLEY L. SCLOVE 312/996-2676 C C INFORMATION AND DECISION SCIENCES DEPT. M/C 294 C C UNIVERSITY OF ILLINOIS AT CHICAGO C C 601 S MORGAN ST C C CHICAGO, IL 60607-7124 C C C C 312-996-2676 C C slsclove@uic.edu C C www.uic.edu/~slsclove C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C ISOPAC consists of the CLUSPAC clustering programs, C C the TSPAC time-series segmentation programs and the C C IMPAC image-segmentation programs. C C C C CLUSPAC is a set of programs implementing clustering C C algorithms derived under the assumption of Gaussian C C class-conditional distributions. The ISDT* programs in C C CLUSPAC are based on the so-called "classification" C C likelihood. The MIX* programs are based on the mixture- C C model likelihood. C C C C Programs MIXP** in the CLUSPAC package are C C mixture-model programs for clustering multivariate data. C C (For univariate data the "MIX1**" programs may be used.) C C MIXPCM and MIXPCMA assume a common covariance matrix C C across distributions. MIXPDT and MIXPDTA allow different C C covariance matrices. C C MIXPCMA and MIXPDTA try a range of numbers of clusters, C C with automatic setting of initial values of the parameters. C C C C C C RESTRICTIONS (CAN BE MODIFIED): C C C N, SAMPLE SIZE, AT MOST 49999; C IP, NUMBER OF VARIABLES, AT MOST 20; C C K, NUMBER OF CLUSTERS, AT MOST 29; C C ITER, MAXIMUM NUMBER OF ITERATIONS, 50. C C C C SUBROUTINE(S) CALLED: C C INVS C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C CHARACTER*72 TITLE, DATFMT,INPFMT DIMENSION X(49999,20),SUM(29,20),ID(49999) C DIMENSION Y(N=IP+1,IP) DIMENSION Y(21,20) DIMENSION F(49999,29),PP(29,49999) DIMENSION ICLUS(49999),ICLSOL(49999) DIMENSION DENOM(49999),XMXPR(49999) DIMENSION SSD(29,20,20),SIGMA(29,20,20) DIMENSION PREC(29,20,20),CMPREC(20,20) DIMENSION DSQ(29) C DIMENSION TITLE(18) DIMENSION NG(29),XBAR(29,20) C DIMENSION DATFMT(18),INPFMT(18) C DIMENSION XIDET(29),XLGDET(29) DIMENSION TRUDET(29) C DIMENSION WGSS(20,20) DIMENSION VARHAT(20,20),COVMX(20,20) DIMENSION XMIN(20),XMAX(20) DIMENSION XJ(29,29) DIMENSION PR(29),XM2LL(29) DIMENSION AICVEC(29),SCHVEC(29),AICPPH(29),SCHPPH(29),IPARM(29) DIMENSION PPOLD(29,49999) DIMENSION PRIRLN(9) DIMENSION TOTAL(29) DIMENSION XMEAN(29) DIMENSION SUMSQS(20,20) DIMENSION SSDEVS(20,20) DIMENSION P(20,20),AMAT(20,20) C CHARACTER*6 ID C DOUBLE PRECISION SUM DOUBLE PRECISION WGSS,SSD,TERM,SSDEVS,SUMSQS,TOTAL DOUBLE PRECISION DEVV,DEVW,PREC,CMPREC DOUBLE PRECISION VARHAT,COVMX,SDJV,SDJW,AMAT DOUBLE PRECISION P, SIGMA DOUBLE PRECISION DET C DOUBLE PRECISION XLGDET DOUBLE PRECISION DSQ,ZSQ DOUBLE PRECISION XBAR DOUBLE PRECISION F,PR DOUBLE PRECISION SUMPXF, SUMLNF DOUBLE PRECISION SCHPPH,AICPPH,SPPHA,SPPHS,PP DOUBLE PRECISION DENOM DOUBLE PRECISION PPOLD DOUBLE PRECISION PRIRLN DOUBLE PRECISION TRUDET,COMDET DOUBLE PRECISION TEST C DOUBLE PRECISION XLIKLY C DOUBLE PRECISION PI C C for Subroutine INVS: DIMENSION AY(20),W(20) DOUBLE PRECISION AY, W C C C00000000111111111122222222223333333333444444444455555555556666666666777 C23456789012345678901234567890123456789012345678901234567890123456789012 C C C .................................................................. C C CONTROL STATEMENTS: C ------------------ C VARIABLE FORMAT DESCRIPTION C -------- ------ ----------- C TITLE (A80) dataset / job title C N (2X,I5) sample size (number of cases) C IP (3X,I2) number of variables C KLO (4X,I2) min number of clusters to be fit C KHI (4X,I2) max number of clusters to be fit C IVERB (6X,I1) level of verbosity of output C IVERB = 1: terse output C upon convergence, values of BIC only C IVERB = 2: medium output C above, plus means C IVERB = 3: verbose output C above, plus covariance matrices C IVERB = 4: complete output C above, for every iteration C C If KLO=KHI, LABELS will be written. C C0000000011111111112222222222333333333344444444445555555555666666666677777777778 C2345678901234567890123456789012345678901234567890123456789012345678901234567890 C C INPFMT (A80) input format, e.g., (4F4.1,1X,A6) C Allows for a last variable which is C an alphanumeric case label. C DATFMT (A80) data format, e.g., (4F4.1) C "DATFMT" WILL ALSO BE USED FOR OUTPUT, SO ALLOW C AT LEAST ONE BLANK AT THE BEGINNING FOR CARRIAGE C CONTROL. C X data, a case at a time, in format specified by DATFMT C ...................................................................... C C READ (5,15000) TITLE C C WRITE PROGRAM INFORMATION: WRITE (6,10000) WRITE (6,10050) WRITE (6,10999) WRITE (6,11001) TITLE C C THE FOLLOWING ARE IMSL DATE AND TIME FUNCTIONS: C CALL TDATE(IDAY,MONTH,IYEAR) C CALL DTIME(IHOUR,MINUTE,ISEC) C WRITE (6,11100) MONTH,IDAY,IYEAR,IHOUR,MINUTE,ISEC C C READ SAMPLE SIZE, N. READ (5,12000) N C N is entered in format (2X,I5) C READ NUMBER OF VARIABLES, IP. READ (5,54000) IP WRITE (6,11400) IP WRITE (6,13000) N C READ (5,24000) KLO READ (5,24000) KHI C READ (5,23000) IVERB WRITE (6,*) 'Level of verbosity is ',IVERB, ' .' WRITE (6,*) 'IVERB = 1: terse output' WRITE (6,*) 'upon convergence, values of BIC only' WRITE (6,*) 'IVERB = 2: medium output' WRITE (6,*) 'above, plus means' WRITE (6,*) 'IVERB = 3: verbose output' WRITE (6,*) 'above, plus covariance matrices' WRITE (6,*) 'IVERB = 4: complete output' WRITE (6,*) 'above, for every iteration' WRITE (6,*) 'If KLO=KHI, LABELS will be written.' C C C READ INPUT FORMAT (INPUT INCLUDES DATA AND ID.) READ (5,15000) INPFMT C READ DATA FORMAT. READ (5,15000) DATFMT C C READ DATA. DO 500 I = 1,N READ (5,INPFMT) (X(I,JV), JV = 1,IP),ID(I) IF (I .EQ. 1) GO TO 100 GO TO 300 100 CONTINUE DO 200 JV = 1,IP XMAX(JV) = X(1,JV) XMIN(JV) = X(1,JV) 200 CONTINUE 300 CONTINUE DO 400 JV = 1,IP IF (X(I,JV) .LT. XMIN(JV)) XMIN(JV)=X(I,JV) IF (X(I,JV) .GT. XMAX(JV)) XMAX(JV)=X(I,JV) 400 CONTINUE 500 CONTINUE WRITE (6,19000) WRITE (6,DATFMT) (XMIN(JV),JV=1,IP) WRITE (6,21000) WRITE (6,DATFMT) (XMAX(JV),JV=1,IP) C WRITE (6,11111) WRITE (6,11112) DO 111 I=1,4 WRITE (6,INPFMT) (X(I,JV),JV=1,IP),ID(I) C KUP = MIN(KHI,N) KREACH = KUP 111 CONTINUE WRITE (6,11096) KUP C C SET CONSTANTS: C XN = N PI = .314159265358979324D+01 C C C SET INITIAL VALUES OF MIXING PROBABILITIES:-- C C XJ(IC,K), K=2 TO 9, IC = 1 TO K C OPTIMAL CLASS PROBABILITIES - NORMAL DISTRIBUTION C JOHARI & SCLOVE (1975). "PARTITIONING A DISTRIBUTION." C COMMUNICATIONS IN STATISTICS. C XJ(1,2)=.5 XJ(2,2)=.5 XJ(1,3)=.2703 XJ(2,3)=.4594 XJ(3,3)=.2703 XJ(1,4)=.1631 XJ(2,4)=.3369 XJ(3,4)=.3369 XJ(4,4)=.1631 XJ(1,5)=.1068 XJ(2,5)=.2444 XJ(3,5)=.2976 XJ(4,5)=.2444 XJ(5,5)=.1068 XJ(1,6)=.0739 XJ(2,6)=.1810 XJ(3,6)=.2451 XJ(4,6)=.2451 XJ(5,6)=.1810 XJ(6,6)=.0739 XJ(1,7)=.0536 XJ(2,7)=.1375 XJ(3,7)=.1986 XJ(4,7)=.2106 XJ(5,7)=.1986 XJ(6,7)=.1375 XJ(7,7)=.0536 XJ(1,8)=.0402 XJ(2,8)=.1067 XJ(3,8)=.1613 XJ(4,8)=.1918 XJ(5,8)=.1918 XJ(6,8)=.1613 XJ(7,8)=.1067 XJ(8,8)=.0402 XJ(1,9)=.0310 XJ(2,9)=.0845 XJ(3,9)=.1324 XJ(4,9)=.1643 XJ(5,9)=.1756 XJ(6,9)=.1643 XJ(7,9)=.1324 XJ(8,9)=.0845 XJ(9,9)=.0310 C C EVALUATE MODEL-SELECTION CRITERIA FOR K=1 (NO CLUSTERING): C C COMPUTE SUM VECTOR AND SUM-OF-PRODUCTS MATRIX: DO 105 JV = 1,IP TOTAL(JV) = 0.0 DO 105 JW = 1,IP SUMSQS(JV,JW) = 0.0 105 CONTINUE DO 110 JV = 1,IP DO 110 I = 1,N TOTAL(JV) = TOTAL(JV) + X(I,JV) 110 CONTINUE DO 120 JV = 1,IP DO 120 JW = 1,IP DO 120 I = 1,N SUMSQS(JV,JW) = SUMSQS(JV,JW) + X(I,JV)*X(I,JW) 120 CONTINUE DO 130 JV = 1,IP XMEAN(JV) = TOTAL(JV)/XN 130 CONTINUE DO 140 JV = 1,IP DO 140 JW = 1,IP SSDEVS(JV,JW) = SUMSQS(JV,JW)- TOTAL(JV)*TOTAL(JW)/XN VARHAT(JV,JW) = SSDEVS(JV,JW)/XN 140 CONTINUE C WRITE SUMMARY STATISTICS FOR WHOLE SAMPLE: WRITE(6,10500) WRITE(6,10100) ( XMEAN(JV), JV = 1,IP ) WRITE(6,10200) DO 150 JV = 1,IP WRITE(6,48000) ( VARHAT(JV,JW), JW=1,IP ) 150 CONTINUE C C Call subroutine INVS to compute determinant of sample cov.mx.: C C ********************** INVS001 C * SUBROUTINE INVS * INVS002 C ********************** INVS003 C INVS004 C INVERT A SYMMETRIC MATRIX (MS=1) IN PLACE AND CALCULATE THE INVS005 C DETERMINANT INVS006 C INVS007 C CALL INVS (A,N,DET,W) INVS008 C INVS009 C A .......... INPUT-OUTPUT MATRIX, N BY N,SYMMETRIC (MS=1) INVS010 C N .......... NUMBER OF ROWS IN A, EQUAL TO NUMBER OF COLUMNS INVS011 C DET ........ OUTPUT SCALAR, DETERMINANT OF A INVS012 C W .......... WORKING VECTOR OF LENGTH N INVS013 C C C Pack matrix into vector used by subroutine INVS: C C I use "AY" rather than "A" for the matrix of INVS C because "A" was used elsewhere. C Pack VARHAT into AY: DO 2610 IVAR = 1,IP DO 2610 JVAR = 1,IVAR KAY = IVAR*(IVAR-1)/2+JVAR AY(KAY) = VARHAT(IVAR,JVAR) 2610 CONTINUE C C C Form of CALL: C CALL INVS (A,N,DET,W) INVS008 C INVS009 C A .......... INPUT-OUTPUT MATRIX, N BY N,SYMMETRIC (MS=1) INVS010 C N .......... NUMBER OF ROWS IN A, EQUAL TO NUMBER OF COLUMNS INVS011 C DET ........ OUTPUT SCALAR, DETERMINANT OF A INVS012 C W .......... WORKING VECTOR OF LENGTH N INVS013 C CALL INVS (AY,IP,DET,W) C TRUDET(1) = DET C C MAX LIKELIHOOD = C (2*PI)**(-N*P/2)*(DET OF COV MX)**(-N/2)*EXP(-N*P/2) C C COMPUTE LOG OF MAX LIKELIHOOD: C XLL=-(XN*IP/2.0)*DLOG(2.0*PI)-(XN/2.0)*DLOG(TRUDET(1))-XN*IP/2.0 XMN2LL = (-2.0)*XLL XM2LL(1) = XMN2LL C C NO. OF PARAMETERS FOR UNCLUSTERED SAMPLE IS: C 1 MEAN VECTOR + 1 COV MX. C NOPARM = IP + IP*(IP+1)/2 IPARM(1) = NOPARM C WRITE (6,72000) NOPARM WRITE(6,70900) XMN2LL AIC = XMN2LL + 2.0*NOPARM WRITE(6,70000) AIC C Note SCH is BIC. SCH = XMN2LL + ALOG(XN)*NOPARM WRITE (6,71000) SCH AICVEC(1) = AIC SCHVEC(1) = SCH C WRITE (6,10050) C C C PERFORM CLUSTERING FOR VARIOUS NUMBERS OF CLUSTERS, K: C DO 850 K = KLO,KHI WRITE (6,26000) K C C C SET CONSTANTS FOR THIS VALUE OF K: C C COUNT NUMBER OF INDEPENDENT PARAMETERS:-- C K MEAN VECTORS OF DIMENSION P AND K P-BY-P COVARIANCE MATRICES, C WHERE P IS THE NUMBER OF VARIABLES, AND K-1 INDEPENDENT C MIXING PROBABILITIES. C NOPARM = K*IP + K*IP*(IP+1)/2 + (K-1) IPARM(K) = NOPARM C IF(IVERB .GE. 2) WRITE (6,72000) NOPARM C C SET INITIAL VALUES OF PARAMETER ESTIMATES FOR THIS VALUE OF K: C C SET INITIAL VALUES OF PRIOR PROBABILITIES EQUAL TO THE C OPTIMAL VALUES FOR A NORMAL DISTRIBUTION:-- C DO 360 IC = 1,K PR(IC) = XJ(IC,K) 360 CONTINUE C C TO SET INITIAL VALUES OF PRIOR PROBABILITIES EQUAL TO C BINOMIAL(X=IC-1;N=K-1,P=1/2). C FOR IC = 1,2,...,K, C REMOVE THE "C" FROM CC1 OF THE NEXT LINE: C PR(IC) = GAMMA(XK)/(GAMMA(XC)*GAMMA(XK-XC+1)*(2**K)) C C FOR EQUAL INITIAL VALUES OF PRIOR PROBABILITIES, C REMOVE THE "C" FROM CC1 OF THE NEXT LINE: C PR(IC) = 1.0/K . C C INITIALIZE POSTERIOR PROBABILITIES: DO 549 I = 1,N DO 549 IC = 1,K PPOLD(IC,I) = 1.0/K 549 CONTINUE C IF (IVERB .GE. 2) WRITE (6,20000) C C Set initial means equal to K of the observations:-- C DO 550 IG = 1,K C 1st observation for IG=1, N-th observation for IG=K: I = INT( 1 + (IG-1)*(N-1)/(K-1) ) C For more centered means, use instead: C I = INT( N/(2*K) ) + INT( (IG-1)*(N/K) ) DO 550 JV = 1,IP XBAR(IG,JV) = X(I,JV) 550 CONTINUE C C EXAMPLE: If N=150 and K=3, the initial means will be C cases 26, 76 AND 126. C C This initialization scheme may be good when the data are C ordered with respect to an important variable. C C C ALTERNATIVELY, HERE IS ANOTHER SCHEME FOR COMPUTING INITIAL MEANS: C (IN 2 DIMENSIONS THIS SCHEME RESULTS IN MEANS ORIENTED FROM C SOUTHWEST TO NORTHEAST): C DO 601 JV = 1,IP C DO 601 IG = 1,K C W = IG/(K+1.0) C XBAR(IG,JV) = (1-W)*XMIN(JV) + W*XMAX(JV) C 601 CONTINUE C C IF (IVERB .GE. 2) WRITE (6,21999) DO 600 IG=1,K IF (IVERB .GE. 2) WRITE(6,22000) IG IF (IVERB. GE. 2) WRITE(6,DATFMT) (XBAR(IG,JV), JV=1,IP) 600 CONTINUE C C SET INITIAL COVARIANCE MATRICES:-- C C First recall how the initial were C set equal to K of the observations:-- C C DO 550 IG = 1,K C I = INT( N/(2*K) ) + INT( (IG-1)*(N/K) ) C DO 550 JV = 1,IP C XBAR(IG,JV) = X(I,JV) C 550 CONTINUE C C The K initial covariance matrices will be set C using K blocks of P+1 observations. C Each block starts with one of the seed points. C C COMPUTE SUM VECTOR AND SUM-OF-PRODUCTS MATRIX: C C BY-PASS THIS METHOD OF INITIALIZING COV MCS: C GO TO 113 C DO 142 IC = 1,K C C Get X's for this value of IC: C Store them as vectors Y(1),Y(2),...,Y(P),Y(P+1): C C Compute first subscript for this block of P+1: C ICN = INT( N/(2*K) ) + INT( (IC-1)*(N/K) ) C DO 107 I=1,IP+1 C DO 107 JV = 1,IP C Y(I,JV) = X(ICN+I-1,JV) C 107 CONTINUE C DO 106 JV = 1,IP C TOTAL(JV) = 0.0 C DO 106 JW = 1,IP C SUMSQS(JV,JW) = 0.0 C 106 CONTINUE C DO 112 JV = 1,IP C DO 112 I = 1,N C TOTAL(JV) = TOTAL(JV) + Y(I,JV) C 112 CONTINUE C DO 121 JV = 1,IP C DO 121 JW = 1,IP C DO 121 I = 1,N C SUMSQS(JV,JW) = SUMSQS(JV,JW) + Y(I,JV)*Y(I,JW) C 121 CONTINUE C DO 131 JV = 1,IP C XMEAN(JV) = TOTAL(JV)/XN C 131 CONTINUE C DO 141 JV = 1,IP C DO 141 JW = 1,IP C SSDEVS(JV,JW) = SUMSQS(JV,JW)- TOTAL(JV)*TOTAL(JW)/XN C VARHAT(JV,JW) = SSDEVS(JV,JW)/XN C SIGMA(IC,JV,JW) = VARHAT(JV,JW) C 141 CONTINUE C C 142 CONTINUE C C C C Here is an alternative way of setting initial cov. mcs.: C Set initial covariance matrices, and scale them up to C covariance matrices. C Some correlation could be programmed in in stmt no. 114: C 113 CONTINUE C DO 630 IC = 1,K DO 610 JV=1,IP DO 610 JW=1,IP 114 SIGMA(IC,JV,JW) = 0.00 610 CONTINUE DO 620 JV = 1,IP SIGMA(IC,JV,JV) = 1.0 620 CONTINUE 630 CONTINUE C SCALE UP THIS CORRELATION MATRIX TO A COVARIANCE MATRIX C USING half of THE WHOLE-SAMPLE VARIANCES C (as if half the variance is between-group and half is C within-group): DO 611 JV=1,IP SDJV = DSQRT(VARHAT(JV,JV)/2.0) DO 611 JW=1,IP SDJW = DSQRT(VARHAT(JW,JW)/2.0) DO 611 IC = 1,K SIGMA(IC,JV,JW) = SDJV*SIGMA(IC,JV,JW)*SDJW 611 CONTINUE C C COMPUTE DETERMINANT AND INVERSE OF C INITIAL COVARIANCE MATRICES: C C DO 640 IC=1,K C C C COPY SIGMA(IC,.,.) INTO COVMX TO BE USED AS MATRIX INPUT TO SUBROUTINE C BECAUSE MATRIX INPUT IS DESTROYED BY SUBROUTINE AND SIGMA(IC,.,.) C WILL BE NEEDED LATER: DO 2404 JV = 1,IP DO 2404 JW = 1,IP COVMX(JV,JW) = SIGMA(IC,JV,JW) 2404 CONTINUE C C C I use "AY" rather than "A" for the matrix of INVS C because "A" was used elsewhere. C C Pack COVMX into AY: DO 2612 IVAR = 1,IP DO 2612 JVAR = 1,IVAR KAY = IVAR*(IVAR-1)/2+JVAR AY(KAY) = COVMX(IVAR,JVAR) 2612 CONTINUE C C CALL INVS (AY,IP,DET,W) C C C Unpack AY into P: DO 2622 IVAR = 1,IP DO 2622 JVAR = 1,IVAR KAY = IVAR*(IVAR-1)/2+JVAR P(IVAR,JVAR) = AY(KAY) 2622 CONTINUE C Fill other triangular portion of matrix P: DO 2632 IVAR=1,IP IP1 = IVAR+1 DO 2632 JVAR = IP1,IP P(IVAR,JVAR) = P(JVAR,IVAR) 2632 CONTINUE C C C C SET INITIAL PRECISION MATRICES FOR THIS K: DO 641 JV = 1,IP DO 641 JW = 1,IP PREC(IC,JV,JW) = P(JV,JW) 641 CONTINUE TRUDET(IC) = DET C C END PARAMETER INITIALIZATION FOR THIS K: C 640 CONTINUE C C C BEGIN ITERATIONS: C ---------------- C C C Set max number of iterations: MAXITER = 50 ITER = 1 700 CONTINUE IF (IVERB .GE. 4) WRITE(6,32000) ITER IF (ITER .EQ. 1) GO TO 901 DO 800 I = 1,N C Store old labels: ICLSOL(I) = ICLUS(I) 800 CONTINUE C 901 CONTINUE C C COMMENCE COMPUTATION OF F(I,IC), I=1,...,N, IC=1,...,K: C C COMPUTE MAHALANOBIS D-SQUARE BETWEEN THE I-TH OBSERVATION AND C THE L-TH MEAN, L=1,2,...,K, I=1,2,...,N:-- C DO 1200 I = 1,N DO 1000 L=1,K DSQ(L) = 0.0D+00 DO 1310 JV=1,IP DEVV = XBAR(L,JV) - X(I,JV) DO 1310 JW=1,IP DEVW = XBAR(L,JW) - X(I,JW) DSQ(L) = DSQ(L) + DEVV*PREC(L,JV,JW)*DEVW 1310 CONTINUE ZSQ = DSQ(L) C IF D-SQ IS INORDINATELY LARGE, SET VALUE OF PDF TO ZERO C (IT IS EXTREMELY SMALL ANYWAY, AND THIS AVOIDS UNDERFLOW): IF ( ZSQ/2.0 .LE. 174.673D+00 ) GO TO 1090 F(I, L) = 0.0D+00 GO TO 1100 C 1090 CONTINUE C Compute C F(I,L) = DEXP(-ZSQ/2.0)/((2*PI)**(IP/2)*DSQRT(TRUDET(L))): C XIP = IP F(I,L) = (-ZSQ/2.0)-(XIP/2.0)*DLOG(2*PI)-0.5*DLOG(TRUDET(L)) IF(F(I,L) .GT. 174.673) WRITE(6,87000) ZSQ,TRUDET(L) F(I,L) = DEXP( F(I,L) ) 1100 CONTINUE 1000 CONTINUE 1200 CONTINUE C C COMPUTE LOG LIKELIHOOD AND VALUES OF MODEL-SELECTION CRITERIA: C C XLIKLY = 1.0 SUMLNF = 0.0 DO 2200 I=1,N SUMPXF = 0.0 DO 2100 IC=1,K SUMPXF = SUMPXF + PR(IC)*F(I,IC) 2100 CONTINUE C C XLIKLY = XLIKLY*SUMPXF IF ( SUMPXF .LE. 0.0 ) GO TO 2190 SUMLNF = SUMLNF + DLOG(SUMPXF) 2190 CONTINUE 2200 CONTINUE C C SUMLNF = DLOG(XLIKLY) XMN2LL = -2.0*SUMLNF XM2LL(K) = XMN2LL C IF (IVERB .GE. 4) WRITE (6,10900) XMN2LL C C C COMPUTE MODEL SELECTION CRITERIA AIC = XMN2LL + 2.0*NOPARM SCH = XMN2LL + ALOG(XN)*NOPARM IF (IVERB .GE. 4) WRITE (6,70000) AIC IF (IVERB .GE. 4) WRITE (6,71000) SCH AICVEC(K) = AIC SCHVEC(K) = SCH C C C COMPUTE POSTERIOR PROBABILITIES OF GROUP MEMBERSHIP, C PP(IC,I), THE CONDITIONAL PROBABILITY OF POP'N IC, GIVEN X(I), C AS PR(IC)*F(I,IC)/DENOM(I): DO 1350 I = 1,N DENOM(I) = 0.0 DO 1350 IC=1,K DENOM(I) = DENOM(I) + PR(IC)*F(I,IC) 1350 CONTINUE C C IF DENOM(I)=0, PP(IC,I) IS SET EQUAL TO OLD PP(IC,I) SO THAT C THE I-TH CASE WILL STILL HAVE WEIGHT IN SUCCEEDING C COMPUTATIONS: C C DO 1400 I = 1,N IF ( DENOM(I) .GT. 0.0D+00 ) GO TO 1410 IF (IVERB .GE. 4) WRITE (6,14560) ID(I) DO 1456 IC = 1,K PP(IC,I) = PPOLD(IC,I) 1456 CONTINUE GO TO 1401 C 1410 CONTINUE DO 1402 IC=1,K PP(IC,I) = PR(IC)*F(I,IC)/DENOM(I) 1402 CONTINUE C 1401 CONTINUE 1400 CONTINUE C IF (IVERB .GE. 4) WRITE (6,55555) IF ( N-15 ) 1421,1421,1422 1422 CONTINUE IF (IVERB .GE. 4) WRITE(6,76000) DO 1420 I = 1,4 IF (IVERB .GE. 4) WRITE(6,82000) ( PP(IC,I), IC = 1,K ) 1420 CONTINUE GO TO 1423 1421 CONTINUE DO 1424 I = 1,N IF (IVERB .GE. 4) WRITE (6,25000) ID(I), (PP(IC,I),IC=1,K) 1424 CONTINUE 1423 CONTINUE C C C UPDATE PARAMETER ESTIMATES: C C UPDATE CLUSTER PRIOR PROBABILITIES PR(IC):-- DO 3800 IC = 1,K PR(IC) = 0.0 DO 3800 I = 1,N PR(IC) = PR(IC) + PP(IC,I) 3800 CONTINUE DO 3900 IC = 1,K PR(IC) = PR(IC)/N 3900 CONTINUE C IF MIXING PROBABILITY IS ZERO FOR ANY CLUSTER, STOP. DO 3901 IC = 1,K IF ( PR(IC) .GT. 0.0D+00 ) GO TO 3901 WRITE (6,40040) K,IC 40040 FORMAT(/' IN FITTING ',I1,' CLUSTERS, PR(',I1,')=0. STOP.'/) GO TO 4004 3901 CONTINUE DO 1750 IG = 1,K DO 1750 JV = 1,IP SUM(IG,JV) = 0.0 DO 1750 JW = 1,IP SSD(IG,JV,JW) = 0.0 1750 CONTINUE C C UPDATE MEANS:-- DO 1875 IC = 1,K DO 1875 JV = 1,IP DO 1875 I = 1,N SUM(IC,JV) = SUM(IC,JV) + PP(IC,I)*X(I,JV) 1875 CONTINUE DO 1900 IG = 1,K TEST = XN*PR(IG) IF ( TEST .EQ. 0.0D+00 ) GO TO 2050 C IF (N*PR(IG) .LT. 0.5) GO TO 2050 GO TO 2150 2050 WRITE (6,74000) IG NGFLAG = 1 C IF EXP.NO. OF OBSERVATIONS FOR A CLUSTER IS LESS THAN 0.5, C END COMPUTATION FOR THIS NUMBER OF CLUSTERS C AND GO ON TO NEXT VALUE OF K: GO TO 850 2150 CONTINUE DO 1900 JV = 1,IP XBAR(IG,JV) = SUM(IG,JV)/(XN*PR(IG)) 1900 CONTINUE C C C CALL XUFLOW(0) C C Compute separate covariance matrices SIGMA C with elements SIGMA(IC,JV,JW): DO 3600 IC=1,K DO 3600 JV=1,IP DO 3600 JW=JV,IP DO 3600 I=1,N DEVV=X(I,JV)-XBAR(IC,JV) DEVW=X(I,JW)-XBAR(IC,JW) IF ( PP(IC,I) .EQ. 0.0 ) GO TO 3600 IF ( DEVV .EQ. 0.0 ) GO TO 3600 IF ( DEVW .EQ. 0.0 ) GO TO 3600 SSD(IC,JV,JW)=SSD(IC,JV,JW) + PP(IC,I)*DEVV*DEVW 3600 CONTINUE C DO 3700 IC=1,K DO 3700 JV=2,IP DO 3700 JW=1,JV-1 SSD(IC,JV,JW)=SSD(IC,JW,JV) 3700 CONTINUE C DO 3710 IC = 1,K DO 3710 JV = 1,IP DO 3710 JW = 1,IP SIGMA(IC,JV,JW) = SSD(IC,JV,JW)/(N*PR(IC)) 3710 CONTINUE C C C Pool to compute VARHAT, estimate of common cov.mx.: C DO 1950 JV = 1,IP DO 1950 JW = 1,IP WGSS(JV,JW) = 0.0 1950 CONTINUE C DO 2300 JV = 1,IP DO 2300 JW = 1,IP DO 2300 IC = 1,K WGSS(JV,JW) = WGSS(JV,JW) + PR(IC)*SSD(IC,JV,JW) 2300 CONTINUE C COMPUTE VARHAT, MLE OF COMMON COVARIANCE MATRIX: DO 2400 JV = 1,IP DO 2400 JW = 1,IP VARHAT(JV,JW) = WGSS(JV,JW)/N 2400 CONTINUE C C UPDATE PRECISION MATRICES: C C NRS1 = 0 C IDET = 1 C C Note VARHAT is pooled cov.mx. C Copy VARHAT into COVMX to be used as matrix input to C subroutine because C matrix input is destroyed by the subroutine C and the common covariance matrix VARHAT will C be needed later: DO 2406 JV = 1,IP DO 2406 JW = 1,IP KAY = JV*(JV-1)/2+JW AY(KAY) = VARHAT(JV,JW) 2406 CONTINUE C COMPUTE INVERSE OF POOLED COVARIANCE MATRIX:-- C C CALL INVS (AY,IP,DET,W) C C Unpack AY into P: DO 3811 IVAR = 1,IP DO 3811 JVAR = 1,IVAR KAY = IVAR*(IVAR-1)/2+JVAR P(IVAR,JVAR) = AY(KAY) 3811 CONTINUE C Fill other triangular portion of matrix P: DO 2630 IVAR=1,IP IP1 = IVAR+1 DO 2630 JVAR = IP1,IP P(IVAR,JVAR) = P(JVAR,IVAR) 2630 CONTINUE C C DO 3813 JV = 1,IP DO 3813 JW = 1,IP CMPREC(JV, JW)=P(JV,JW) 3813 CONTINUE C COMDET = DET*(10.0**IDET) COMDET = DET C C C UPDATE PRECISION MATRICES: C Copy SIGMA into AMAT, a matrix temp: DO 3711 IC = 1,K DO 3712 JV = 1,IP DO 3712 JW = 1,IP AMAT(JV,JW) = SIGMA(IC,JV,JW) 3712 CONTINUE C Pack VARHAT into AY: DO 4100 IVAR = 1,IP DO 4100 JVAR = 1,IVAR KAY = IVAR*(IVAR-1)/2+JVAR AY(KAY) = AMAT(IVAR,JVAR) 4100 CONTINUE C C C CALL INVS (AY,IP,DET,W) C C AY now contains the elements of the inverse. C Unpack AY into P: DO IVAR = 1,IP DO JVAR = 1,IVAR KAY = IVAR*(IVAR-1)/2+JVAR P(IVAR,JVAR) = AY(KAY) ENDDO ENDDO C Fill other triangular portion of matrix P: DO IVAR=1,IP IP1 = IVAR+1 DO JVAR = IP1,IP P(IVAR,JVAR) = P(JVAR,IVAR) ENDDO ENDDO C C Load P as precision matrix of the IC-th group: DO 3713 JV = 1,IP DO 3713 JW = 1,IP PREC(IC,JV,JW)=P(JV,JW) 3713 CONTINUE C TRUDET(IC) = DET*(10.0**IDET) C IF DET=0, REPLACE COV MX WITH COMMON COV MX: IF ( TRUDET(IC) .GT. 0.0D+00 ) GO TO 3711 TRUDET(IC) = COMDET IF (IVERB .GE. 4) WRITE (6, 17000) IC DO 4000 JV = 1,IP DO 4000 JW = 1,IP PREC(IC,JV,JW) = CMPREC(JV,JW) 4000 CONTINUE 3711 CONTINUE C C C (END PARAMETER-ESTIMATE UPDATE SEQUENCE) C C IF (ITER .EQ. 1) GO TO 900 C STORE OLD LABELS: DO 825 I = 1,N ICLSOL(I) = ICLUS(I) 825 CONTINUE C 900 CONTINUE C C COMPUTE NEW LABELS BY MAX POSTERIOR PROBABILITY: C DO 1600 I = 1,N XMXPR(I) = PP(1,I) ICLUS(I) = 1 DO 1600 IC = 2,K IF ( PP(IC,I) .GT. XMXPR(I) ) GO TO 1500 GO TO 1600 1500 XMXPR(I) = PP(IC,I) ICLUS(I) = IC 1600 CONTINUE C DO 1601 I = 1,N DO 1601 IC = 1,K PPOLD(IC,I) = 1.0/K PPOLD(IC,I) = PP(IC,I) 1601 CONTINUE C IF (N .GE. 31) GO TO 250 IF (IVERB .GE. 4) WRITE (6,14000) C IF (IVERB .GE. 4) WRITE (6,16000) IF (IVERB .GE. 4) WRITE (6,18000) (ID(I),ICLUS(I), I=1,N) 250 CONTINUE C C C IF (ITER .EQ. 1) GO TO 3000 C IF NEW CLUSTERING IS SAME AS OLD, NO FURTHER ITERATION IS C NECESSARY AND RESULTS FOR THIS VALUE OF K WILL BE PRINTED. C OTHERWISE, FURTHER ITERATIONS WILL BE PERFORMED. C DO 2900 I = 1,N IF (ICLUS(I) .EQ. ICLSOL(I)) GO TO 2900 GO TO 3000 2900 CONTINUE GO TO 3300 3000 CONTINUE C C COUNT NUMBERS IN CLUSTERS: DO 3311 IC = 1,K NG(IC) = 0 3311 CONTINUE DO 3321 I = 1,N IC = ICLUS(I) NG(IC) = NG(IC) + 1 3321 CONTINUE IF (IVERB .GE. 2) WRITE (6,36000) (NG(IC),IC=1,K) C ITER = ITER + 1 C IF ( ITER .GE. MAXITER+1 ) GO TO 3100 C C Unless the maximum number of iterations has already been C performed, go to stmt no. 700 and begin another: GO TO 700 3100 WRITE (6,88000) C C 3300 CONTINUE C C IF (IVERB .GE. 2) WRITE (6,40000) (PR(IC),IC=1,K) DO 2800 IC = 1,K IF (IVERB .GE. 2) XWRITE (6,30000) IC, (XBAR(IC,JV),JV=1,IP) 2800 CONTINUE C IF (IVERB .GE. 3) WRITE (6,42000) DO 2500 JV=1,IP IF (IVERB .GE. 3) WRITE (6,48000) (VARHAT(JV,JW),JW=1,IP) 2500 CONTINUE C IF (N .GE. 31) GO TO 3510 C LIST CLUSTER MEMBERSHIP OF EACH CASE: IF (KLO .EQ. KHI) WRITE (6,52000) DO 3500 I = 1,N IF (KLO .EQ. KHI) 1WRITE (6,50000) ID(I), ICLUS(I), ( X(I,JV),JV=1,IP ) 3500 CONTINUE 3510 CONTINUE C C C IF (IVERB .GE. 4) WRITE (6,58000) ITER C C C C COMPUTE VALUES OF MODEL-SELECTION CRITERIA: C C AIC = XMN2LL + 2.0*NOPARM SCH = XMN2LL + ALOG(XN)*NOPARM C IF (IVERB .LE. 3) GO TO 852 IF ( N .GE. 301 ) GO TO 852 DO 851 IC=1,K IF (IVERB .GE. 4) WRITE (6,44444) IC DO 851 I = 1,N IF ( ICLUS(I) .EQ. IC ) WRITE (6,44445) ID(I) 851 CONTINUE 852 CONTINUE C C C GO ON TO NEXT NUMBER OF CLUSTERS (NEXT VALUE OF K), C if KHI has not already been reached. C WRITE (6,10050) C C End DO-LOOP on K, the number of clusters: 850 CONTINUE C C C 4004 CONTINUE KREACH = K-1 C C Posterior probabilities of alternative models in terms of MSCs: C SCHPPH(K) = PR(H(K))* C CONST*(2*PI)**(-IPARM(K)/2.0)*EXP(-SCHVEC(K)/2.0)*PRIOR PDF(K); C Here, for SCH, i.e., BIC, PR(H(K)) IS TAKEN TO BE UNIFORM. C AICPPH(K) = PR(H(K))* C CONST*(2*PI)**(-IPARM(K)/2.0)*EXP(-SCHVEC(K)/2.0)*PRIOR PDF(K); C Here, for AIC, PR(H(K)) IS TAKEN TO BE Const*EXP(-NOPARM(K)). SCHPPH(1) = (IPARM(1)/2.0)*DLOG(2.0*PI) - SCHVEC(1)/2.0 AICPPH(1) = (IPARM(1)/2.0)*DLOG(2.0*PI) - SCHVEC(1)/2.0 C Add in log of prior Pr[H(k)] for AIC: AICPPH(1) = AICPPH(1) - IPARM(1) C Add in log of prior Pr[H(k)] for Schwarz criterion: SCHPPH(1) = SCHPPH(1) + ALOG(1.0/KUP) C Approximate log of PRIOR PDF: C PRIRLN(1) = -(1*IP/2.0)*DLOG(2*PI) -IP/2.0 C Note that this assumes unit variance: It will have to be C corrected to use a variance more reflective of that in the C observations. C Meantime, set it equal to zero. PRIRLN(1) = 0.0 C Add in log of approximate prior pdf: SCHPPH(1) = SCHPPH(1) + PRIRLN(1) AICPPH(1) = AICPPH(1) + PRIRLN(1) C C ADD CONSTANT SCHVEC(1)/2 OR AICVEC(1)/2 TO PREVENT UNDERFLOW: C (This is okay because the added term cancels out below.) C SCHPPH(1) = SCHPPH(1) + SCHVEC(1)/2 C AICPPH(1) = AICPPH(1) + AICVEC(1)/2 C DO 6110 K = 2,KREACH SCHPPH(K) = (IPARM(K)/2.0)*DLOG(2*PI) - SCHVEC(K)/2.0 AICPPH(K) = (IPARM(K)/2.0)*DLOG(2*PI) - SCHVEC(K)/2.0 C Add in log Pr[H(k)]: AICPPH(K) = AICPPH(K) - IPARM(K) SCHPPH(K) = SCHPPH(K) + ALOG(1.0/KUP) C Approximate log of PRIOR PDF: C PRIRLN(K) = -(K*IP/2.0)*DLOG(2*PI) -K*IP/2.0 C Note that this assumes unit variance: It will have to be C corrected to use a variance more reflective of that in the C observations. C Meantime, set it equal to zero. PRIRLN(K) = 0.0 AICPPH(K) = AICPPH(K) + PRIRLN(K) SCHPPH(K) = SCHPPH(K) + PRIRLN(K) C C ADD CONSTANT SCHVEC(1)/2 OR AICVEC(1)/2 TO PREVENT UNDERFLOW: C SCHPPH(K) = SCHPPH(K) + SCHVEC(1)/2 C AICPPH(K) = AICPPH(K) + AICVEC(1)/2 C 6110 CONTINUE C AT THIS POINT, WE HAVE THE LOGS OF THE PPHS; WE NOW CONVERT C THEM TO THE PPHS: SCHPPH(1) = DEXP(SCHPPH(1)) AICPPH(1) = DEXP(AICPPH(1)) DO 6119 K = 2,KREACH C If an argument is less or = -174.673, there will be underflow: IF ( SCHPPH(K) .LE. -174. ) GO TO 6115 SCHPPH(K) = DEXP(SCHPPH(K)) GO TO 6119 6115 SCHPPH(K) = 0.0 6119 CONTINUE DO 6121 K = 2,KREACH IF ( AICPPH(K) .LE. -174.673 ) GO TO 6126 AICPPH(K) = DEXP(AICPPH(K)) GO TO 6121 6126 AICPPH(K) = 0.0 6121 CONTINUE SPPHS = 0.0 SPPHA = 0.0 DO 6120 K = 1,KREACH IF ( SCHPPH(K) .GT. 0.0D+00 ) SPPHS = SPPHS + SCHPPH(K) IF ( AICPPH(K) .GT. 0.0D+00 ) SPPHA = SPPHA + AICPPH(K) 6120 CONTINUE DO 6130 K = 1,KREACH SCHPPH(K) = SCHPPH(K)/SPPHS AICPPH(K) = AICPPH(K)/SPPHA 6130 CONTINUE C C C WRITE VALUES OF MSC'S FOR VARIOUS K: WRITE (6,11001) TITLE WRITE(6,34000) DO 605 K = 1,KREACH WRITE(6,46000) K,XM2LL(K),IPARM(K),AICVEC(K),AICPPH(K), 1SCHVEC(K),SCHPPH(K) 605 CONTINUE C C WRITE (6,10050) WRITE (6,11001) TITLE WRITE (6,96666) 96666 FORMAT(' CONDITION CODE FOR THIS RUN: ') IF ( NGFLAG .EQ. 0 ) GO TO 4002 WRITE (6,96000) GO TO 4001 4002 CONTINUE WRITE (6,98000) 4001 CONTINUE C WRITE (6,10050) WRITE (6,11097) C WRITE (6,86000) C WRITE (6,11100) MONTH,IDAY,IYEAR,IHOUR,MINUTE,ISEC C CALL DTIME(IHOUR,MINUTE,ISEC) C CALL TDATE(IDAY,MONTH,IYEAR) C WRITE (6,89000) C WRITE (6,11100) MONTH,IDAY,IYEAR,IHOUR,MINUTE,ISEC STOP C C 10000 FORMAT(/' ...................................................'/ X1X,' ',/ X1X,'MIXPDTA CLUSPAC - MIXTURE-MODEL CLUSTERING OF CASES',// X' VARYING COVARIANCE MATRICES '/ X' AUTOMATIC SETTING OF INITIAL PARAMETER ESTIMATES '// X' Copyright (C) 1992 Stanley Louis Sclove '// X' Developed and programmed by: '// X19X,'Prof. Stanley L. Sclove, Ph.D. '/ X19X,'Information & Decision Sciences Dept. (MC 294) '/ X19X,'University of Illinois at Chicago '/ X19X,'601 S Morgan St. '/ X19X,'Chicago, IL 60607-7124 '/ X19X,' '/ X19X,'312-996-2676 '/ X19X,'slsclove@uic.edu '/ X19X,'www.uic.edu/~slsclove '/ X' ...................................................'/ X' '/ X' MIXPDTA VERSION 3.02 2003: Mar 23 '/) C C C00000000111111111122222222223333333333444444444455555555556666666666777 C23456789012345678901234567890123456789012345678901234567890123456789012 10050 FORMAT(/' .....................................................'/) 10100 FORMAT(' MEAN VECTOR: ',(9F8.2/)/) 10999 FORMAT( ' PROBLEM TITLE IS' ) 10200 FORMAT(/' COVARIANCE MATRIX: ') 10500 FORMAT(/,' STATISTICS FOR WHOLE (UNCLUSTERED) SAMPLE: ') 11001 FORMAT(1X,A72/) C C11002 FORMAT(/1X,18A4/) 11112 FORMAT(' ----- ---- ---- ------'/) 11096 FORMAT(/' NUMBER OF CLUSTERS TO REPORT. . . . . . . . 1 TO', XI2/) 11111 FORMAT(/' FIRST FOUR DATA POINTS') 11097 FORMAT(' MIXPDTA CLUSPAC - MIXTURE-MODEL CLUSTERING OF CASES',/) 11100 FORMAT(1X,I2,'/',I2,'/',I5,3X,'AT ',I2,':',I2,':',I2/) 19000 FORMAT( 1X,'MINIMUM OF EACH VARIABLE: ') 20000 FORMAT(1X,'INITIAL MEANS') 23000 FORMAT(6X,I1) 24000 FORMAT(4X,I2) 25000 FORMAT(1X, A6, 9F8.3) 55555 FORMAT(/' CASE POSTERIOR PROBS OF CLUSTER MEMBERSHIP: ') 44444 FORMAT(' CLUSTER ',I2) 44445 FORMAT(1X,A6) 10900 FORMAT(/' MINUS 2 LOG LIKELIHOOD = ',F13.1) C11000 FORMAT(/,1X,18A4 ) 11500 FORMAT(' DET = ',E15.7, ' IDET = ',I2/) 12000 FORMAT(2X,I5) 13000 FORMAT(1X,'NUMBER OF OBSERVATIONS (SAMPLE SIZE), N . . ', I5) 11400 FORMAT(1X,'NUMBER OF VARIABLES . . . . . . . . . . . . ', I5) 14000 FORMAT(1X, 'CLUSTERING:') C15000 FORMAT(18A4) 15000 FORMAT(A72) C15000 FORMAT(80A1) 14560 FORMAT(1X,'THE DENOMINATOR OF THE POSTERIOR PROBS FOR ', X1X,A6, ' IS 0, SO THE POSTERIOR PROBS PP(IC,I), IC=1,...,K,'/ X' HAVE BEEN SET EQUAL TO THOSE FROM THE ', X' PRECEDING ITERATION SO THAT THE I-TH CASE WILL STILL '/ X' HAVE WEIGHT IN SUCCEEDING COMPUTATIONS.') 16000 FORMAT(/,1X,'CASES AND LABELS:--'/) 17000 FORMAT(' PREC.MX. OF CLUSTER ',I2,' HAS BEEN SET TO POOLED ', X'PREC.MX.') 18000 FORMAT(1X, (7(1X,A6,': ',I2)/) ) 21999 FORMAT(' '//) 22000 FORMAT(1X,'MEAN VECTOR FOR CLUSTER ',I2,': ') 26000 FORMAT(1X,'K = ',I2,' CLUSTERS') C 30000 FORMAT(1X,'MEAN VECTOR FOR CLUSTER ',I2,': ',/,(9F13.3/)) 30000 FORMAT(1X,'MEAN VECTOR FOR CLUSTER ',I2,': ',(9F13.3/)) 32000 FORMAT(1X,'ITERATION ', I2) 21000 FORMAT(1X, 'MAXIMUM OF EACH VARIABLE: ') 36000 FORMAT(/,1X,'NUMBERS IN CLUSTERS:',3X,9(I5,3X)/) 40000 FORMAT(1X,'MIXING PROBABILITIES:',3X,9(F4.2,3X)/) 42000 FORMAT(1X,'COMMON COVARIANCE MATRIX (MLE):',/ ) 34000 FORMAT(1X,'MODEL SELECTION CRITERIA'// X' MODEL-SELECTION CRITERIA '/ X' NUMBER MINUS 2 NUMBER (PPH=POST.PROB.OF MODEL):'/ X' OF TIMES OF AKAIKE SCHWARZ '/ X' CLUSTERS, LOG OF PARAM- ------------ ------------ '/ X' K LIKELIHOOD ETERS VALUE PPH VALUE PPH '/ X' --------- ---------- ----- ------- --- ------- ---'//) C 9 XXXX.X 143 XXXX.X .XX XXXX.X .XX C WRITE(6,46000) K,XM2LL(K),IPARM(K),AICVEC(K),AICPPH(K), C XSCHVEC(K),SCHPPH(K) 46000 FORMAT(5X,I1,6X,F8.1, 4X,I5,1X,F9.1,1X,F4.2,1X,F9.1,1X,F4.2) 48000 FORMAT(1X,5F13.3) 50000 FORMAT(1X,A6,': ',I2,(8F13.3/)) 52000 FORMAT(//,1X,'CLUSTER MEMBERSHIP OF EACH CASE:'/ X1X,'CASE LABEL DATA'/ ) 54000 FORMAT(3X,I2) 58000 FORMAT(1X, 'CONVERGENCE: NO CASE CHANGED CLUSTERS AFTER ', X'ITERATION ',I2,'.'/ ) 60000 FORMAT(/,1X,'JFLG = ',I2,'. IF JFLG=0, COMPUTATION OF DET', X' WENT WELL; OTHERWISE, THERE WAS TROUBLE OR MATRIX WAS ', X'ILL-CONDITIONED.'//) 62000 FORMAT(/1X,'DET = ',E13.5,' IDET = ',I3,5X, X'ACTUAL DET. = DET*10**IDET',//) 70900 FORMAT(1X,' MINUS 2 LOG LIKELIHOOD = ',F13.1) 70000 FORMAT(1X,' AIC = ',F13.1) 71000 FORMAT(1X,' BIC = ',F13.1) 72000 FORMAT(1X,' NUMBER OF PARAMETERS = ',I5) 74000 FORMAT(/1X,'EXPECTED NO. OF OBSERVATIONS IN GROUP ',I3,'IS ', X' LESS THAN 0.5. ', X' END COMPUTATION FOR THIS NUMBER OF CLUSTERS. '/) 76000 FORMAT(/' POSTERIOR PROB OF GROUP MEMBERSHIP FOR 1ST 4 CASES:') 82000 FORMAT(1X,29F5.2) 87000 FORMAT(1X,2E15.7) 86000 FORMAT(/' TIME BEGAN: ') 89000 FORMAT(/' TIME FINISHED: ') 88000 FORMAT(' ALGORITHM HAS NOT CONVERGED IN MAX NUMBER OF ', X'ITERATIONS. STOP') 96000 FORMAT(/' NGFLAG=1: FOR SOME K, AT LEAST ONE CLUSTER HAS '/ X' AN EXP.NO. OF OBSERVATIONS LESS THAN 0.5 .') 98000 FORMAT(//1X,'PROGRAM ENDED NORMALLY.') END C C INVS000 C ********************** INVS001 C * SUBROUTINE INVS * INVS002 C ********************** INVS003 C INVS004 C INVERT A SYMMETRIC MATRIX (MS=1) IN PLACE AND CALCULATE THE INVS005 C DETERMINANT INVS006 C INVS007 C CALL INVS (A,N,DET,W) INVS008 C INVS009 C A .......... INPUT-OUTPUT MATRIX, N BY N,SYMMETRIC (MS=1) INVS010 C N .......... NUMBER OF ROWS IN A, EQUAL TO NUMBER OF COLUMNS INVS011 C DET ........ OUTPUT SCALAR, DETERMINANT OF A INVS012 C W .......... WORKING VECTOR OF LENGTH N INVS013 C INVS014 SUBROUTINE INVS(A,N,C,W) INVS015 DOUBLE PRECISION U,X,Y,Z,D,A,C,W INVS016 DIMENSION A(1), W(1) INVS017 INTEGER DIAGMK INVS018 INTEGER DIAG,DIAG2,ROWNO,ROWCOL INVS019 INTEGER COLNO INVS020 D=A(1) INVS021 IF(D) 10,20,10 INVS022 10 A(1)=1.0D0/D INVS023 IF(N-1) 20,20,200 INVS024 200 DIAG=1 INVS025 DO 140 K=2,N INVS026 KM1=K-1 INVS027 DIAGMK=DIAG INVS028 DIAG=DIAG+K INVS029 U=A(DIAG) INVS030 COLNO=DIAGMK INVS031 DIAG2=0 INVS032 DO 90 I=1,KM1 INVS033 X=0.0D0 INVS034 COLNO=COLNO+1 INVS035 ROWNO=DIAGMK INVS036 J=1 INVS037 ROWCOL=DIAG2 INVS038 300 IF(J-I) 310,320,320 INVS039 310 ROWCOL=ROWCOL+1 INVS040 ROWNO=ROWNO+1 INVS041 Y=A(ROWCOL) INVS042 Z=A(ROWNO) INVS043 X=X+Y*Z INVS044 J=J+1 INVS045 GO TO 300 INVS046 320 ROWCOL=ROWCOL+1 INVS047 60 ROWNO=ROWNO+1 INVS048 Y=A(ROWCOL) INVS049 Z=A(ROWNO) INVS050 X=X+Y*Z INVS051 ROWCOL=ROWCOL+J INVS052 J=J+1 INVS053 IF(J-K) 60,70,70 INVS054 70 W(I)=-X INVS055 Y=A(COLNO) INVS056 U = U-X*Y INVS057 DIAG2=DIAG2+I INVS058 90 CONTINUE INVS059 D=D*U INVS060 IF(U) 100,20,100 INVS061 100 A(DIAG)=1.0D0/U INVS062 ROWNO=DIAGMK INVS063 DIAG2=0 INVS064 DO 140 I=1,KM1 INVS065 ROWNO=ROWNO+1 INVS066 DIAG2=DIAG2+I INVS067 X=W(I) INVS068 X=X/U INVS069 A(ROWNO)=X INVS070 ROWCOL=DIAG2 INVS071 DO 140 J=I,KM1 INVS072 Y=W(J) INVS073 Z=A(ROWCOL) INVS074 A(ROWCOL)=Z+X*Y INVS075 ROWCOL=ROWCOL+J INVS076 140 CONTINUE INVS077 20 C=D INVS078 C WRITE (*,*) ' C=', C C IF(C)6665,99,6665 INVS079 IF(C) 6665,99,6665 INVS079 99 WRITE(6,97) INVS080 97 FORMAT(' ERROR MESSAGE FROM SUBROUTINE INVS ......... MATRIX IS SI INVS081 1NGULAR'/'1') INVS082 C CALL EXIT STOP INVS083 6665 CONTINUE C6665 CALL SLITET(4,LTEST) INVS084 C IF(LTEST-1)6666,6667,6666 INVS085 C6667 CALL SLITE(4) INVS086 C CALL DPRNT(A,N,N,1,6HINVS ) INVS087 6666 RETURN INVS088 END INVS089