@heading number right $$ INITIAL.CLUSPAC(SCLOVE) p.$$ -------------------------------------------------------------------------------- @end 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 C CALL SUBROUTINE TO COMPUTE DETERMINANT OF SAMPLE COVARIANCE C MATRIX: C IDET = 1 NRS1 = 0 C C COPY VARHAT INTO COVMX TO BE USED AS MATRIX INPUT TO MATEQ C BECAUSE MATRIX INPUT IS DESTROYED BY MATEQ AND VARHAT WILL C BE NEEDED LATER: DO 2402 JV = 1,IP DO 2402 JW = 1,IP COVMX(JV,JW) = VARHAT(JV,JW) 2402 CONTINUE CALL MATEQ(COVMX, IP,20,JFLG,DET,IDET,IV,NRS1,P,20) C C GENERAL FORM OF CALL IS: C CALL MATEQ(A,M,N,JFLG,DET,IDET,IV,NRS1,P,LL) C C THE INPUT MATRIX A IS DESTROYED IN THE COMPUTATION PROCESS. C ON RETURN, P CONTAINS THE INVERSE OF A. C SEE SUBROUTINE LISTING FOR FULLER EXPLANATION. C DET(VARHAT) = DET*10.0**IDET C IF (JFLG .GT. 0) WRITE (6,60000) JFLG C C C XIDET(1) = IDET C XLGDET(1) = DLOG(DET) + XIDET(1)*DLOG(10.0D+00) TRUDET(1) = DET*(10.0**IDET) C WRITE (6,62000) DET,IDET 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 WRITE(6,10900) XMN2LL 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 AIC = XMN2LL + 2.0*NOPARM WRITE(6,70000) AIC SCH = XMN2LL + ALOG(XN)*NOPARM WRITE (6,78000) 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 = 2,KUP 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 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 WRITE (6,20000) 20000 FORMAT(/ 1X,'INITIAL MEANS') C C SET INITIAL MEANS EQUAL TO K OF THE OBSERVATIONS:-- C DO 550 IG = 1,K I = INT( N/(2*K) ) + INT( (IG-1)*(N/K) ) DO 550 JV = 1,IP XBAR(IG,JV) = X(I,JV) 550 CONTINUE C EXAMPLE: IF N=150 AND K=3, THE INITIAL MEANS WILL BE C CASES 26, 76 AND 126. 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 DO 600 IG=1,K WRITE (6,21999) WRITE (6,22000) IG WRITE (6,DATFMT) ( XBAR(IG,JV), JV=1,IP ) 600 CONTINUE C C SET INITIAL COVARIANCE MATRICES:-- C (HERE SOME SLIGHT POSITIVE CORRELATION IN PROGRAMMED IN. C THIS COULD BE CHANGED.) C DO 630 IC = 1,K DO 610 JV=1,IP DO 610 JW=1,IP SIGMA(IC,JV,JW) = 0.40 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 THE WHOLE-SAMPLE VARIANCES: DO 611 JV=1,IP SDJV = DSQRT(VARHAT(JV,JV)) DO 611 JW=1,IP SDJW = DSQRT(VARHAT(JW,JW)) DO 611 IC = 1,K SIGMA(IC,JV,JW) = SDJV*SIGMA(IC,JV,JW)*SDJW 611 CONTINUE C COMPUTE INITIAL INVERSE COVARIANCE MATRIX: