MACRO	# minitab local macro for versions 11-13

MIX13 X GUESS MOD; 
PARAM M S P;
POSTPROB PP.1-PP.10.

#______________________________________________________
# macro uses em algorithm to find maximum likelihood 
# solution to fitting a mixture of n normals 
# to a given univariate sample.

# macro provides lists of the estimated means, standard deviations,
# and proportions in the session window. Also provides a plot of the mixture
# model and its component normal densities. The macro optionally
# provides output of the estimated parameters and posterior
# probabilities if the subcommands are used to specify output columns.
# for param specify 3 columns. For postprob specify 10 columns.

# last modified: 07/02/06

# author: michael thompson

# further information: amc technical brief no 22

#______________________________________________________

#input declarations

MCOLUMNS X 	# input data in a column, not changed by procedure.
 
MCOLUMNS GUESS 	# starting estimates of component means in a column, 
		# unchanged by procedure.

MCONSTANTS MOD	# if mod is set to 1, the algorithm will calculate a 
		# pooled standard deviation, otherwise each componment 
		# will have a separate standard deviation.
#_____________________________________________________________

#OUTPUT DECLARATIONS

MCOLUMNS M S P 	# estimated component parameters
MCOLUMNS PP.1-PP.10	# posterior probabilities of membership
#_______________________________________________________________

# workspace 

MCOLUMNS Y	
MCOLUMNS SUMPF	
MCOLUMNS PF.1-PF.10	
MCOLUMNS F.1-F.10	
MCOLUMNS NIT	
MCOLUMNS CLOGL	
MCOLUMNS SUMP
MCOLUMNS WSS
MCOLUMNS XPOINT

MCONSTANTS NDATA 
MCONSTANTS I INDEX
MCONSTANTS LOW HIGH INC TOP
MCONSTANTS PREVIOUS
MCONSTANTS NCOMP
MCONSTANTS LOGL
MCONSTANTS POOL
MCONSTANTS TEMP TEMPA

#____________________________________________________
NOTITLE
COPY GUESS M
LET NCOMP=N(M)	
LET NDATA=N(X)

ERASE S P

# define prior probabilities for group membership

SORT X Y 
LET LOW=Y(1)

DO I=1:NCOMP
	IF I<NCOMP
		LET HIGH=(M(I)+M(I+1))/2
	ELSE
		LET HIGH=Y(NDATA)
	ENDIF

	CODE (LOW:HIGH)99999 Y PP.I
	CODE (-99999:99998)0 (99999)1 PP.I PP.I

	LET LOW=HIGH*1.00001
ENDDO

# start em algorithm

LET INDEX=0
LET PREVIOUS=0

WHILE INDEX<200
LET INDEX=INDEX+1

# the e-step (estimates parameters)



DO I=1:NCOMP
	LET SUMP(I)=SUM(PP.I)
	LET P(I)=SUMP(I)/NDATA
	LET M(I)=SUM(Y*PP.I)/SUMP(I)	
	LET WSS(I)=SUM((Y-M(I))**2*PP.I)
	LET S(I)=SQRT(WSS(I)/SUMP(I))
ENDDO

IF MOD=1
	LET POOL=SQRT(SUM(WSS)/SUM(SUMP))
	DO I=1:NCOMP
		LET S(I)=POOL
	ENDDO
ENDIF
	
	

# the m-step (estimates posterior probabilities)

DO I=1:NCOMP
	LET F.I=0.398942*EXPO(-(Y-M(I))**2/(2*S(I)**2))/S(I)
	LET PF.I=P(I)*F.I
ENDDO

RSUM PF.1-PF.NCOMP SUMPF

DO I=1:NCOMP
	LET PP.I=PF.I/SUMPF
ENDDO

LET LOGL=SUM(LOGE(SUMPF))

# stopping criterion

LET NIT(1)=INDEX
IF ABSOLUTE((LOGL-PREVIOUS)/LOGL)<1.0E-9
	LET INDEX=201
ENDIF

LET PREVIOUS=LOGL

ENDWHILE

# set up of session window output

NOTE
NOTE
NOTE NORMAL MIXTURE MODEL
NOTE

NOTE
NOTE ESTIMATED MEANS (M), STANDARD DEVIATIONS (S), AND PROPORTIONS (P)
print M S P 
NOTE
LET TEMP="LOG LIKELIHOOD = "
LET CLOGL(1)=LOGL
NTOA CLOGL CLOGL
LET TEMPA=CLOGL(1)
KKCAT TEMP TEMPA TEMP
NAME TEMP '-'
PRINT TEMP
LET TEMP="NUMBER OF ITERATIONS = "
NTOA NIT NIT
LET TEMPA=NIT(1)
KKCAT TEMP TEMPA TEMP
PRINT TEMP
# set up plot of model and components

# define plot points

LET LOW=Y(1)*0.9
LET HIGH=Y(NDATA)*1.1
LET INC=(HIGH-LOW)/99

SET XPOINT
(LOW:HIGH/INC)
END

# calculate component densities at plot points

DO I=1:NCOMP
	LET PF.I=P(I)*0.398942*EXPO(-(XPOINT-M(I))**2/(2*S(I)**2))/S(I)
ENDDO

RSUM PF.1-PF.NCOMP SUMPF
LET TOP=1.1*MAXIMUM(SUMPF)

BRIEF 0
LAYOUT

DO I=1:NCOMP

PLOT PF.I*XPOINT;
	CONNECT;
		COLOR 4;
	MAXIMUM 1 HIGH;
	MINIMUM 1 LOW;
	MAXIMUM 2 TOP;
	MINIMUM 2 0;
	AXIS 1;
	LABEL " ";
	AXIS 2;
	LABEL " ".
ENDDO

PLOT SUMPF*XPOINT;
	TITLE "Mixture model and normal components";
	CONNECT;
		COLOR 2;
	MAXIMUM 1 HIGH;
	MINIMUM 1 LOW;
	MAXIMUM 2 TOP;
	MINIMUM 2 0;
	AXIS 1;
	LABEL "Result";
	AXIS 2;
	LABEL "Density".

ENDLAYOUT
BRIEF 2
TITLE
ENDMACRO