	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,' AND PARTITION',I2,
     :		  ' 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,' 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,' 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(1X,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,' 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,' 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(/,1X,'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(/,1X,'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(1X,I2,1X,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
