PROGRAM MULTIMIX * This program fits a mixture of multivariate distributions using * the EM algorithm. The data file contains both categorical and * continuous variables. * If the program does not converge after iter=200 iterations, the * estimates of the parameters will be entered into * EMPARAMEST.OUT. This file can then be used as the parameter * input file for PROGRAM MULTIMIX if desired. * The assignment of the observations to groups (IGP(i)) and the * posterior probabilities (Zij's) are entered into GROUPS.OUT. * This file can be used in MINITAB etc. for further analysis. * NOTE: * (1) THIS PROGRAM REQUIRES VARIABLES IN A PARTITION TO BE STORED * CONTIGUOUSLY. HENCE THE DATA IS READ IN WITH THE VARIABLE ORDER * BEING SPECIFIED BY JP(J). INTYPE(J) AND NCAT(J) BOTH REFER TO * THE REARRANGED DATA * (2) IF SDENS=0, THEN Z(II,K) IS SET TO 0.01 * (3) THE PROGRAM CURRENTLY HAS A MAXIMUM OF * 1500 OBSERVATIONS (iob=1500) * 6 GROUPS (ik6=6) * 15 ATTRIBUTES & 15 PARTITIONS (ip15=15) * 10 LEVELS OF CATEGORIES (im10=10) * 200 ITERATIONS FOR CONVERGENCE (iter=200) * ******ALTER IF REQUIRED****** * (REMEMBER TO ALTER PARAMETERS IN DETINV ALSO) * The parameter file contains:- * NG - the number of groups * NOBS - the number of observations * NVAR - the number of variables * NPAR - the number of partitions * ISPEC - an indicator variable for a specified grouping of the * observations (1 = observations are not specified into * groups, 2 = observations are specified into groups) * JP(j) - column in the data array in which the jTH variable of * the file will be stored * IP(l) - number of variables in the lTH partition, l=1,NPAR * IPC(l) - number of continuous variables in partition l, l=1,NPAR * ISV(l) - indicator starting value for the partition, l=1,NPAR * IEV(l) - indicator end value for the partition, l=1,NPAR * ITYPE(l) - indicator giving the type of model each partition is * (1 = categorical, 2 = MVN, 3 = location models), l=1,NPAR * INTYPE(j) - an indicator variable giving type of variable, j=1,NVAR * (1 = categorical variable, 2 = continuous variable, * 3 = categorical variable involved in location model, * 4 = continuous variable involved in the location model) * NCAT(j) - number of categories of jTH variable. * (For continuous variables, set NCAT(J)=0) * PI(k) - estimated mixing proportions for each group, k=1,NG * THETA(K,J,M) - estimated probability that the jTH categorical * variable is at level M, given that in group k * EMUL(k,l,j,m)- estimated mean vector for each group, each * partition, and each location of the location model variables * EMU(K,L,J) -estimated mean vector for each group and each * partition of continuous variables * VARIX(K,L,I,J) -estimated covariance matrices for each group IMPLICIT DOUBLE PRECISION(A-H,O-Z) CHARACTER*16 infile, datafile PARAMETER (PIE=3.141592653589792, iob=2000, ip15=15, ik6=6, : im10=10, iter=200) DIMENSION EMU(ik6,ip15,ip15),VAR(ik6,ip15,ip15,ip15), : Z(iob,ik6), PI(ik6), DENS(ik6,ip15), VARIN(ik6,ip15,ip15,ip15), : ZSUM(ik6), XSUM(ik6,ip15), ADET(ik6,ip15), CLOGLI(iter), : VARIX(ik6,ip15,ip15,ip15), APRODENS(ik6), IP(ip15), IPC(ip15), : ISV(ip15), IEV(ip15), ITYPE(ip15), EMUL(ik6,ip15,ip15,im10), : THETA(ik6,ip15,0:im10), PROB(ik6,iob),INTYPE(ip15),NCAT(ip15), : IM(ip15), IGP(iob), NUM(ik6), XSUM2(ik6,ip15,ip15),ZM(ik6,ip15), : IGRP(iob), JP(ip15), IX(iob,ip15), X(iob,ip15) PRINT *, ' MIXTURE ESTIMATION BY EM' PRINT *, '-------------------------------' PRINT *, 'Data file: ' READ (*,10) datafile 10 FORMAT (A16) PRINT *, 'Parameter file: ' READ (*,10) infile OPEN (7, FILE='GENERAL.OUT', STATUS='NEW') OPEN (8, FILE=datafile, STATUS='OLD') OPEN (9, FILE=infile, STATUS='OLD') OPEN(12, FILE='GROUPS.OUT',STATUS='NEW') READ (9,*) NG, NOBS, NVAR, NPAR,ISPEC READ (9,*) (JP(J),J=1,NVAR) READ (9,*) (IP(L),L=1,NPAR) READ (9,*) (IPC(L), L=1,NPAR) READ (9,*) (ISV(L),L=1,NPAR) READ (9,*) (IEV(L),L=1,NPAR) READ (9,*) (ITYPE(L),L=1,NPAR) READ (9,*) (INTYPE(J),J=1,NVAR) READ (9,*) (NCAT(J),J=1,NVAR) * read in estimates of the parameters if grouping not specified IF(ISPEC.EQ.1) THEN READ(9,*) (PI(K),K=1,NG) DO 12 K=1,NG DO 12 L=1,NPAR IF (ITYPE(L).EQ.1) THEN DO 13 J=ISV(L),IEV(L) READ(9,*) (THETA(K,J,M),M=1,NCAT(J)) 13 CONTINUE ELSE IF(ITYPE(L).EQ.2) THEN READ(9,*)(EMU(K,L,J),J=1,IPC(L)) ELSE IF(ITYPE(L).EQ.3) THEN DO 14 J=ISV(L),IEV(L) IF(INTYPE(J).EQ.3) THEN READ(9,*)(THETA(K,J,M),M=1,NCAT(J)) IM(L)=NCAT(J) END IF 14 CONTINUE J1=0 DO 15 J=ISV(L),IEV(L) IF(INTYPE(J).EQ.4) THEN J1=J1+1 READ(9,*) : (EMUL(K,L,J1,M),M=1,IM(L)) END IF 15 CONTINUE END IF 12 CONTINUE DO 16 K=1,NG DO 16 L=1,NPAR IF(ITYPE(L).NE.1) : READ(9,*)((VARIX(K,L,I,J),J=1,IPC(L)),I=1,IPC(L)) 16 CONTINUE ELSE READ(9,*) (IGRP(I),I=1,NOBS) DO 18 I=1,NOBS DO 18 K=1,NG Z(I,K)=0.0 18 CONTINUE DO 31 I=1,NOBS IK=IGRP(I) Z(I,IK)=1.0 31 CONTINUE DO 33 L=1,NPAR IF(ITYPE(L).EQ.3) THEN DO 32 J=ISV(L),IEV(L) IF (INTYPE(J).EQ.3) IM(L)=NCAT(J) 32 CONTINUE END IF 33 CONTINUE END IF DO 17 I=1,NOBS READ(8,* ) (X(I,JP(J)),J=1,NVAR) 17 CONTINUE * Check the categorical variable to see if 0 is a category. * Add 1 to the variables if it is so that Xij=0 is taken as * level 1, Xij=1 is level 2 etc. DO 80 J=1,NVAR IF ((INTYPE(J).EQ.1).OR.(INTYPE(J).EQ.3)) THEN IMIN=NINT(X(1,J)) DO 81 I=2,NOBS INEXT=NINT(X(I,J)) IF(INEXT.LT.IMIN) IMIN=INEXT 81 CONTINUE IF (IMIN.EQ.0) THEN DO 82 I=1,NOBS IX(I,J)=NINT(X(I,J)) + 1 82 CONTINUE ELSE DO 83 I=1,NOBS IX(I,J)=NINT(X(I,J)) 83 CONTINUE END IF END IF 80 CONTINUE * print the input values WRITE(7,50)NG,NOBS,NVAR,NPAR 50 FORMAT(/,' NO OF GROUPS IS ',I3,/,' NO OF OBSERVATIONS IS ', : I5,/,' NO OF VARIABLES ',I3,/,' NO OF PARTITIONS', I3) WRITE(7,51)(IP(L),L=1,NPAR) 51 FORMAT(/,' THE NO OF VARIABLES IN EACH PARTITION IS',/,10I3) * Send to the M step if ISPEC = 2 (groups specified) IF(ISPEC.EQ.2) GO TO 100 * otherwise, print the parameter estimates WRITE(7,52) 52 FORMAT(/,' MIXING PROPORTIONS') WRITE(7,53)(PI(K),K=1,NG) 53 FORMAT(/,10F6.3) DO 56 K=1,NG DO 56 L=1,NPAR IF(ITYPE(L).EQ.1) THEN WRITE(7,54) K 54 FORMAT(/,2X,' THETA(K,J,M) FOR GROUP ',I2) DO 64 J=ISV(L),IEV(L) WRITE(7,55)(THETA(K,J,M),M=1,NCAT(J)) 55 FORMAT(10F8.4) 64 CONTINUE ELSE IF(ITYPE(L).EQ.2) THEN WRITE(7,57) K,L 57 FORMAT(/,' FOR GROUP ',I2, X, ' AND PARTITION',I2,X, : ' THE MEAN IS') WRITE(7,58)(EMU(K,L,J),J=1,IPC(L)) 58 FORMAT(10F12.6) ELSE WRITE(7,57) K,L WRITE(7,58)((EMUL(K,L,J,M),M=1,IM(L)),J=1,IPC(L)) DO 62 J=ISV(L),IEV(L) IF (INTYPE(J).EQ.3) THEN WRITE(7,54) WRITE(7,55)(THETA(K,J,M),M=1,NCAT(J)) END IF 62 CONTINUE END IF 56 CONTINUE DO 59 K=1,NG DO 59 L=1,NPAR IF (ITYPE(L).NE.1) THEN WRITE(7,60) L,K 60 FORMAT(/,2X,'VARIANCE FOR PARTITION',I2,X,' AND GROUP',I2) DO 63 I=1,IPC(L) WRITE(7,61)(VARIX(K,L,I,J),J=1,IPC(L)) 61 FORMAT(10F13.6) 63 CONTINUE END IF 59 CONTINUE DO 74 L=1,NPAR IF (ITYPE(L).EQ.1) THEN WRITE(7,75) L,IP(L),ITYPE(L) 75 FORMAT(' PARTITION',I3,' HAS',I2,' VARIABLES',/,' ITYPE', : ' IS',I2,' HENCE A CATEGORICAL MODEL FOR THIS PARTITION') ELSE IF (ITYPE(L).EQ.2) THEN WRITE(7,76) L,IP(L),ITYPE(L) 76 FORMAT(' PARTITION',I3,' HAS',I2,' VARIABLES',/,' ITYPE', : ' IS',I2,' HENCE A MVN MODEL FOR THIS PARTITION') ELSE WRITE(7,77) L,IP(L),ITYPE(L) 77 FORMAT(' PARTITION',I3,' HAS',I2,'VARIABLES',/,' ITYPE', : ' IS',I2,' HENCE A LOCATION MODEL FOR THIS PARTITION') END IF 74 CONTINUE * E step of the EM algorithm. * Estimate the complete data sufficient statistics given the * data, & current values of means,variances and mixing proportions. * call DETINV to calculate the determinants and inverses for each * covariance matrix. This subroutine uses NAG to calculate * the determinants and inverses. ICNT=1 99999 DO 19 L=1,NPAR IF (ITYPE(L).NE.1) THEN ITYPE_CURR=ITYPE(L) IPC_CURR=IPC(L) L_CURR=L CALL DETINV(NG,L_CURR,IPC_CURR,ADET,VARIX,VARIN) END IF 19 CONTINUE ALL=0.0 DO 20 II = 1,NOBS SDENS = 0.0 DO 21 K=1,NG PROB(K,II)=1.0 APRODENS(K)=1.0 * evaluate the discrete variables contribution to the densities DO 22 J=1,NVAR IF((INTYPE(J).EQ.1).OR.(INTYPE(J).EQ.3)) : PROB(K,II) = THETA(K,J,IX(II,J))*PROB(K,II) 22 CONTINUE * evaluate the continuous variables contribution to the densities DO 23 L=1,NPAR IF (ITYPE(L).NE.1) THEN DENS(K,L)=0.0 I1=0 * (i) evaluate the MVN contribution IF(ITYPE(L).EQ.2) THEN DO 24 I=ISV(L),IEV(L) J1=0 I1=I1+1 DO 24 J=ISV(L),IEV(L) J1=J1+1 DENS(K,L) = DENS(K,L) : + (X(II,I)-EMU(K,L,I1))*VARIN(K,L,I1,J1)*(X(II,J)-EMU(K,L,J1)) 24 CONTINUE * (ii) evaluate the continuous location variables contribution ELSE IF(ITYPE(L).EQ.3) THEN DO 25 I=ISV(L),IEV(L) IF(INTYPE(I).EQ.3) M=IX(II,I) 25 CONTINUE DO 26 I=ISV(L),IEV(L) IF (INTYPE(I).EQ.4) THEN J1=0 I1=I1+1 DO 30 J=ISV(L),IEV(L) IF (INTYPE(J).EQ.4) THEN J1=J1+1 DENS(K,L)=DENS(K,L) + (X(II,I)-EMUL(K,L,I1,M)) : *VARIN(K,L,I1,J1)*(X(II,J)-EMUL(K,L,J1,M)) END IF 30 CONTINUE END IF 26 CONTINUE END IF DENS(K,L) = DEXP(-0.5*DENS(K,L)) A=0.5*FLOAT(IPC(L)) DENS(K,L)=DENS(K,L)/((2.0*PIE)**(A)*DSQRT(ADET(K,L))) APRODENS(K)=APRODENS(K)*DENS(K,L) END IF 23 CONTINUE APRODENS(K)=APRODENS(K)*PROB(K,II)*PI(K) SDENS=SDENS+APRODENS(K) 21 CONTINUE IF (SDENS.NE.0.0) THEN DO 27 K=1,NG Z(II,K)=APRODENS(K)/SDENS 27 CONTINUE ALL=ALL+DLOG(SDENS) ELSE DO 28 K=1,NG Z(II,K)=0.01 28 CONTINUE WRITE(7,29)II 29 FORMAT(//'SUM OF DENSITY FUNCTIONS IS ZERO', : 'FOR OBSERVATION',I4) END IF 20 CONTINUE * Check on convergence - look at the likelihood function * If the absolute value of the difference in 2 likelihoods is * less than a tolerance value the estimates are written out. * A check is made on the number of iterations (200 max). * statement 100 - the M step * statement 500 - print out the current estimates (algorithm * has converged) * statement 501 - algorithm hasn't converged, current estimates * are printed out. CLOGLI(ICNT) = ALL WRITE(7,888) ICNT,CLOGLI(ICNT) 888 FORMAT(1X,'FOR LOOP',I5,'LOGLIKELIHOOD IS',F16.8) IF (ICNT.LE.10) GO TO 100 C=1.D-06 TOL=ABS(CLOGLI(ICNT) - CLOGLI(ICNT-10)) IF(TOL.LE.C) GO TO 500 IF(ICNT.GE.iter) GO TO 501 * M step of the EM algorithm * Calculate an updated estimate of the parameters * means and covariances. * (i) the mixing proportions, 100 DO 101 K=1,NG ZSUM(K)=0.0 DO 102 II=1,NOBS ZSUM(K)=ZSUM(K) + Z(II,K) 102 CONTINUE PI(K) = ZSUM(K)/NOBS 101 CONTINUE * (ii) the conditional probabilities DO 103 K=1,NG DO 103 L=1,NPAR IF(ITYPE(L).EQ.1) THEN DO 104 J=ISV(L),IEV(L) DO 104 M=1,NCAT(J) THETA(K,J,M) = 0.0 DO 105 II=1,NOBS IF (IX(II,J).EQ.M) THETA(K,J,M) : = THETA(K,J,M) + Z(II,K) 105 CONTINUE THETA(K,J,M)=THETA(K,J,M)/ZSUM(K) 104 CONTINUE ELSE IF(ITYPE(L).EQ.2) THEN * (iii) the means (EMU(K,L,J)) J1=0 DO 106 J=ISV(L),IEV(L) J1=J1+1 XSUM(K,J1)=0.0 DO 107 II=1,NOBS XSUM(K,J1)=XSUM(K,J1) + X(II,J)*Z(II,K) 107 CONTINUE EMU(K,L,J1)=XSUM(K,J1)/ZSUM(K) 106 CONTINUE * (iv) the location model parameters * (i)the discrete variable ELSE IF(ITYPE(L).EQ.3) THEN DO 108 J=ISV(L),IEV(L) IF(INTYPE(J).EQ.3) THEN DO 109 M=1,NCAT(J) THETA(K,J,M)=0.0 DO 110 II=1,NOBS IF(IX(II,J).EQ.M) THETA(K,J,M) : = THETA(K,J,M) + Z(II,K) 110 CONTINUE THETA(K,J,M)=THETA(K,J,M)/ZSUM(K) 109 CONTINUE END IF 108 CONTINUE * (ii) the continuous variables the means EMUL(K,L,J,M) DO 111 J=ISV(L),IEV(L) IF (INTYPE(J).EQ.3) THEN DO 112 M=1,IM(L) J1=0 DO 112 JJ=ISV(L),IEV(L) IF (INTYPE(JJ).EQ.4) THEN J1=J1+1 XSUM2(K,J1,M)=0.0 ZM(K,M)=0.0 DO 113 II=1,NOBS IF(IX(II,J).EQ.M) THEN XSUM2(K,J1,M)=XSUM2(K,J1,M) + : X(II,JJ)*Z(II,K) ZM(K,M)=ZM(K,M)+Z(II,K) END IF 113 CONTINUE EMUL(K,L,J1,M)=XSUM2(K,J1,M)/ZM(K,M) END IF 112 CONTINUE END IF 111 CONTINUE WRITE(7,601)K,((EMUL(K,L,J1,M),M=1,IM(L)),J1=1,IPC(L)) 601 FORMAT(/2X,'FOR GROUP',I2,X,'EMUL',10F10.4) END IF 103 CONTINUE * Calculate updated estimates of the variances * (i) the multivariate normal data DO 115 K=1,NG DO 116 L=1,NPAR IF (ITYPE(L).NE.1) THEN DO 117 J=1,IPC(L) DO 117 I=1,J VAR(K,L,I,J)=0.0 117 CONTINUE END IF 116 CONTINUE DO 118 II=1,NOBS DO 118 L=1,NPAR IF (ITYPE(L).EQ.2) THEN J1=0 DO 119 J=ISV(L),IEV(L) I1=0 J1=J1+1 DO 119 I=ISV(L),ISV(L)+J1 I1=I1+1 VAR(K,L,I1,J1) = VAR(K,L,I1,J1) + (X(II,J) : -EMU(K,L,J1))*(X(II,I)-EMU(K,L,I1))*Z(II,K) 119 CONTINUE END IF 118 CONTINUE * (ii) the continuous location data DO 120 L=1,NPAR IF (ITYPE(L).EQ.3) THEN DO 121 J=ISV(L),IEV(L) IF (INTYPE(J).EQ.3) THEN DO 122 II=1,NOBS DO 122 M=1,IM(L) IF(IX(II,J).EQ.M) THEN J1=0 DO 123 JJ=ISV(L),IEV(L) IF (INTYPE(JJ).EQ.4) THEN I1=0 J1=J1+1 DO 124 I=ISV(L),ISV(L)+J1 IF (INTYPE(I).EQ.4) THEN I1=I1+1 VAR(K,L,I1,J1) = VAR(K,L,I1,J1)+(X(II,JJ) : -EMUL(K,L,J1,M))*(X(II,I)-EMUL(K,L,I1,M)) : *Z(II,K) END IF 124 CONTINUE END IF 123 CONTINUE END IF 122 CONTINUE END IF 121 CONTINUE END IF 120 CONTINUE DO 126 L=1,NPAR IF (ITYPE(L).NE.1) THEN DO 125 J=1,IPC(L) DO 125 I=1,J VAR(K,L,I,J)=VAR(K,L,I,J)/ZSUM(K) VAR(K,L,J,I)=VAR(K,L,I,J) 125 CONTINUE END IF 126 CONTINUE 115 CONTINUE * Make a copy of the covariance matrix before we use NAG DO 130 K=1,NG DO 130 L=1,NPAR IF(ITYPE(L).NE.1) THEN DO 131 J=1,IPC(L) DO 131 I=1,IPC(L) VARIX(K,L,I,J)=VAR(K,L,I,J) 131 CONTINUE END IF 130 CONTINUE ICNT=ICNT+1 * send back to the E step GO TO 99999 * Write out the current estimates of the parameters. * (1) If the algorithm has not converged:- 501 WRITE(7,502) 502 FORMAT(//'----------------------------------------------------') WRITE(7,503) iter 503 FORMAT(//' THE EM ALGORITHM HAS NOT CONVERGED AFTER ',I3, : /,'ITERATIONS BUT THE CURRENT ESTIMATES OF THE PARAMETERS ', : /,'WILL BE PRINTED OUT.') WRITE(7,502) * The parameters are to be written to 'EMPARAMEST.DAT' to be * used as input for the PROGRAM MULTIMIX. ISPEC is set to 1. OPEN(11, FILE='EMPARAMEST.OUT',STATUS='NEW') ISPEC=1 WRITE(11,504) NG, NOBS, NVAR, NPAR, ISPEC 504 FORMAT(X,5I6) WRITE(11,505) (IP(L),L=1,NPAR) WRITE(11,505) (IPC(L),L=1,NPAR) WRITE(11,505) (ISV(L),L=1,NPAR) WRITE(11,505) (IEV(L),L=1,NPAR) WRITE(11,505) (ITYPE(L),L=1,NPAR) WRITE(11,505) (INTYPE(J),J=1,NVAR) WRITE(11,505) (NCAT(J),J=1,NVAR) 505 FORMAT(10I4) WRITE(11,506) (PI(K),K=1,NG) 506 FORMAT(10F10.6) DO 507 K=1,NG DO 507 L=1,NPAR IF (ITYPE(L).EQ.1) THEN DO 508 J=ISV(L),IEV(L) WRITE(11,509)(THETA(K,J,M),M=1,NCAT(J)) 509 FORMAT(10F10.6) 508 CONTINUE ELSE IF(ITYPE(L).EQ.2) THEN WRITE(11,510) (EMU(K,L,J),J=1,1,IPC(L)) 510 FORMAT(10F13.6) ELSE DO 511 J=ISV(L),IEV(L) IF (INTYPE(J).EQ.3) THEN WRITE(11,509)(THETA(K,J,M),M=1,NCAT(J)) END IF 511 CONTINUE DO 512 J=1,IPC(L) WRITE(11,510) (EMUL(K,L,J,M),M=1,IM(L)) 512 CONTINUE END IF 507 CONTINUE DO 513 K=1,NG DO 513 L=1,NPAR IF (ITYPE(L).NE.1) THEN DO 514 I=1,IPC(L) WRITE(11,515)(VARIX(K,L,I,J),J=1,IPC(L)) 515 FORMAT(10F13.6) 514 CONTINUE END IF 513 CONTINUE * (2) the current estimates of the parameters are printed out * Estimates of the proportions in each group and loglikelihood 500 WRITE(7,888) ICNT,CLOGLI(ICNT) DO 540 K=1,NG WRITE(7,541) K,PI(K) 541 FORMAT(//'THE ESTIMATE OF THE MIXING PROPORTION IN GROUP ',I3, : 'IS ', F10.8) 540 CONTINUE * estimates of the probabilities for each group DO 525 K=1,NG DO 525 L=1,NPAR WRITE(7,524)K,L 524 FORMAT(//,'THE CURRENT ESTIMATES FOR GROUP',I3,' AND' : ' PARTITION ',I3, ' ARE ') IF (ITYPE(L).EQ.1) THEN DO 528 J=ISV(L),IEV(L) WRITE(7,529)J 529 FORMAT(/,' FOR VARIABLE ',I3,X,'THETA(K,J,M) IS') WRITE(7,530) (THETA(K,J,M),M=1,NCAT(J)) 530 FORMAT(10F10.6) 528 CONTINUE ELSE IF(ITYPE(L).EQ.2) THEN * Estimate of the means for each group. WRITE(7,531) (EMU(K,L,J),J=1,IPC(L)) 531 FORMAT(/,' THE MEAN FOR THIS PARTITION IS ',/,10F13.6) * estimates of the location model parameters ELSE IF(ITYPE(L).EQ.3) THEN DO 533 J=ISV(L),IEV(L) IF (INTYPE(J).EQ.3) THEN WRITE(7,534)J,(THETA(K,J,M),M=1,NCAT(J)) 534 FORMAT(/,' FOR VARIABLE',I3,X,'THETA(K,J,M)' : ' IS',10F10.6) END IF 533 CONTINUE J1=0 DO 535 J=ISV(L),IEV(L) IF(INTYPE(J).EQ.4) THEN J1=J1+1 WRITE(7,536) (EMUL(K,L,J1,M),M=1,IM(L)) 536 FORMAT(/,X,'THE MEAN FOR THE CONTINUOUS LOCATION' : ' VARIABLES IS',/,10F13.6) END IF 535 CONTINUE END IF 525 CONTINUE * Estimates of the variances for each group. DO 523 K=1,NG WRITE(7,526)K 526 FORMAT(/,' THE CURRENT ESTIMATE OF THE COVARIANCE MATRIX FOR', : ' GROUP',I3) DO 523 L=1,NPAR IF(ITYPE(L).NE.1) THEN WRITE(7,522)L 522 FORMAT(' AND PARTITION ',I3) DO 543 I=1,IPC(L) WRITE(7,527)(VAR(K,L,I,J),J=1,IPC(L)) 527 FORMAT(10F15.6) 543 CONTINUE END IF 523 CONTINUE * Determine the assignment of the observations to groups DO 516 I=1,NOBS IMAX = 1 DO 517 K=2,NG IF(Z(I,K).GT.Z(I,IMAX)) THEN IMAX=K END IF 517 CONTINUE IGP(I)=IMAX 516 CONTINUE WRITE(7,542) 542 FORMAT(/,X,'THE ASSIGNMENT OF OBSERVATIONS TO GROUPS',/) WRITE(7,518) (IGP(I),I=1,NOBS) 518 FORMAT(10I3) DO 537 K=1,NG NUM(K)=0 537 CONTINUE DO 538 I=1,NOBS DO 538 K=1,NG IF(IGP(I).EQ.K) NUM(K)=NUM(K)+1 538 CONTINUE WRITE(7,539)(NUM(K),K=1,NG) 539 FORMAT(1X,/,' TOTAL NUMBERS IN EACH GROUP',/,10I5) * have a look at the Zij,s, and write out the assigned groups and * the Zij's to GROUPS.OUT WRITE(7,521) 521 FORMAT(//'THE ESTIMATES OF THE POSTERIOR PROBALITIES') DO 519 I=1,NOBS WRITE(7,520)I,(Z(I,K),K=1,NG) 520 FORMAT('OBSERVATION',I4,2X,10F10.6) WRITE(12,532) IGP(I), (Z(I,K),K=1,NG) 532 FORMAT(X,I2,X,10F9.6) 519 CONTINUE END SUBROUTINE DETINV(NG,L_CURR,IPC_CURR,ADET,VARIX,VARIN) IMPLICIT DOUBLE PRECISION(A-H,O-Z) PARAMETER (ik6=6,ip15=15,IIP=ip15+1) DIMENSION VARIX(ik6,ip15,ip15,ip15),VARIN(ik6,ip15,ip15,ip15), : ADET(ik6,ip15) INTEGER IA,IT,NNULL REAL*8 TEMP(100),WKSPCE(100),B(100),TOL,CMY(100),PROD,RMAX IA=IIP IB=IIP NNULL = 0 DO 201 K=1,NG IT=0 TOL=0.0 N=IPC_CURR DO 202 J=1,IPC_CURR TOL=TOL+SQRT(VARIX(K,L_CURR,J,J)) DO 202 I=1,J IT=IT+1 TEMP(IT) = VARIX(K,L_CURR,I,J) 202 CONTINUE TOL=(TOL/FLOAT(IPC_CURR))*0.000001 CALL SYMINV(TEMP,N,B,WKSPCE,NULL,IER,RMAX,CMY,TOL) IF (IER.NE.0) THEN WRITE(7,555) K,IER 555 FORMAT (/2X,'Terminal error in SYMINV for matrix',I3, : ' as IFAULT is ',I3) RETURN ELSE IF (NULL.NE.0) THEN WRITE (7,556) K,NULL 556 FORMAT (/2X,'Rank deficiency of covariance matrix' : ,I3,' is',I3) WRITE (*,556) K,NULL NNULL=NNULL+1 ENDIF IT=0 PROD=1.0 DO 220 J=1,IPC_CURR JJ=J*(J+1)/2 PROD=PROD*CMY(JJ)*CMY(JJ) DO 220 I=1,J IT=IT+1 VARIN(K,L_CURR,I,J)=B(IT) VARIN(K,L_CURR,J,I)=VARIN(K,L_CURR,I,J) 220 CONTINUE ADET(K,L_CURR)=PROD ENDIF 201 CONTINUE NULL=NNULL RETURN END subroutine syminv(a,n,c,w,nullty,ifault,rmax,CMY,TOL) implicit double precision (a-h,o-z) c c algorithm as7, applied statistics, vol.17, 1968. c c arguments:- c a() = input, the symmetric matrix to be inverted, stored in c lower triangular form c n = input, order of the matrix c c() = output, the inverse of a (a generalized inverse if c is c singular), also stored in lower triangular. c c and a may occupy the same locations. c w() = workspace, dimension at least n. c nullty = output, the rank deficiency of a. c ifault = output, error indicator c = 1 if n < 1 c = 2 if a is not +ve semi-definite c = 0 otherwise c rmax = output, approximate bound on the accuracy of the diagonal c elements of c. e.g. if rmax = 1.e-04 then the diagonal c elements of c will be accurate to about 4 dec. digits. c c latest revision - 18 april 1981 c c*************************************************************************** c dimension a(1),c(1),w(n),CMY(100) nrow=n ifault=1 if(nrow.le.0) go to 100 ifault=0 c c cholesky factorization of a, result in c c call chola(a,nrow,c,nullty,ifault,rmax,w,TOL) if(ifault.ne.0) go to 100 c c invert c & form the product (cinv)'*cinv, where cinv is the inverse c of c, row by row starting with the last row. c irow = the row number, ndiag = location of last element in the row. c nn=nrow*(nrow+1)/2 do 200 imy=1,nn 200 cmy(imy)=c(imy) irow=nrow ndiag=nn 10 if(c(ndiag).eq.0.0) go to 60 l=ndiag do 20 i=irow,nrow w(i)=c(l) l=l+i 20 continue icol=nrow jcol=nn mdiag=nn 30 l=jcol x=0.0 if(icol.eq.irow) x=1.0/w(irow) k=nrow 40 if(k.eq.irow) go to 50 x=x-w(k)*c(l) k=k-1 l=l-1 if(l.gt.mdiag) l=l-k+1 go to 40 50 c(l)=x/w(irow) if(icol.eq.irow) go to 80 mdiag=mdiag-icol icol=icol-1 jcol=jcol-1 go to 30 60 l=ndiag do 70 j=irow,nrow c(l)=0.0 l=l+j 70 continue 80 ndiag=ndiag-irow irow=irow-1 if(irow.ne.0) go to 10 100 return end c c subroutine chola(a,n,u,nullty,ifault,rmax,r,TOL) implicit double precision (a-h,o-z) c c algorithm as6, applied statistics, vol.17, 1968, with c modifications by a.j.miller c c arguments:- c a() = input, a +ve definite matrix stored in lower-triangular c form. c n = input, the order of a c u() = output, a lower triangular matrix such that u*u' = a. c a & u may occupy the same locations. c nullty = output, the rank deficiency of a. c ifault = output, error indicator c = 1 if n < 1 c = 2 if a is not +ve semi-definite c = 0 otherwise c rmax = output, an estimate of the relative accuracy of the c diagonal elements of u. c r() = output, array containing bounds on the relative accuracy c of each diagonal element of u. c c latest revision - 18 april 1981 c c*************************************************************************** c dimension a(1),u(1),r(n) c c eta should be set equal to the smallest +ve value such that c 1.0 + eta is calculated as being greater than 1.0 in the accuracy c being used. c C data eta/1.e-07/ ETA=TOL ifault=1 if(n.le.0) go to 100 ifault=2 nullty=0 rmax=eta r(1)=eta j=1 k=0 c c factorize column by column, icol = column no. c do 80 icol=1,n l=0 c c irow = row number within column icol c do 40 irow=1,icol k=k+1 w=a(k) if(irow.eq.icol) rsq=(w*eta)**2 m=j do 10 i=1,irow l=l+1 if(i.eq.irow) go to 20 w=w-u(l)*u(m) if(irow.eq.icol) rsq=rsq+(u(l)**2*r(i))**2 m=m+1 10 continue 20 if(irow.eq.icol) go to 50 if(u(l).eq.0.0) go to 30 u(k)=w/u(l) go to 40 30 u(k)=0.0 if(abs(w).gt.abs(rmax*a(k))) go to 100 40 continue c c end of row, estimate relative accuracy of diagonal element. c 50 rsq=sqrt(rsq) if(abs(w).le.5.*rsq) go to 60 if(w.lt.0.0) go to 100 u(k)=sqrt(w) r(i)=rsq/w if(r(i).gt.rmax) rmax=r(i) go to 70 60 u(k)=0.0 nullty=nullty+1 70 j=j+icol 80 continue ifault=0.0 100 return end