normem=function(x, mu, s, p)
{
#In addition to the data set x, this program requires initial estimates
# of the mean vector, standard deviation vector, and mixing proportions
#vector to be provided.  A typical call would be:
#          normem(x,c(0,1),c(1,1),c(.5,.5))
	nt <- length(mu)
	n <- length(x)
	diff <- 1
	iter <- 0
	onet <- rep(1, nt)
	onen <- rep(1, n)
	onetn <- onet %o% onen	
#g is a ntxn matrix containing nt normal pdfs evaluated at n data
#points. F is a 1xn matrix containing the mixture pdf evaluated at
#n data points.  A is a txn matrix containing the conditional
#probabilities of group membership for each data point.  Much of the
#processing is done using ntxn arrays and operations
	gexp <-  - (onet %o% x - mu %o% onen)^2/(2 * s^2 %o% onen)
	g <- exp(gexp)/((sqrt(2 * pi) * onetn) * (s %o% onen))
	f <- p %*% g
	l <- sum(log(f))
	while(abs(diff) > 1e-06 & iter < 1000) {
     
	a <- ((p %o% onen) * g)/(onet %*% f)
		p <- apply(a, 1, sum)/n
		mu <- as.numeric(a %*% x/(n * p))
		s <- apply(a * ((onet %o% x - mu %o% onen)^2), 1, sum)
		s <- sqrt(s/(n * p))
		gexp <-  - (onet %o% x - mu %o% onen)^2/(2 * s^2 %o% onen)
		g <- exp(gexp)/((sqrt(2 * pi) * onetn) * (s %o% onen))
		f <- p %*% g
		lnew <- sum(log(f))
		diff <- abs(lnew - l)/max(abs(lnew), 1)
		l <- lnew
		iter <- iter + 1
	}
	out <- list(mu, s, p, iter)
        names(out)=c("Mu Vector","Sigma Vector","P Vector","No. Iterations")
        gradient.em(f,x)
	out
}
gradient.em=function(f,x){
#gradient.em implements a diagnostic for including additional support
#points in the finite mixing distribution.  It is usually presented as 
#a 1-dim procedure, but has been implemented here as a 2-d procedure for
#both the mean and variance of a mixed normal density
n=length(x)
m.g=seq(min(x),max(x),length=15) 
s.g=seq(diff(range(x))/8,diff(range(x))/4,length=15)
one15=rep(1,15)
onen=rep(1,n)
fmix=one15%o%one15%o%as.numeric(f)
fexp=-((one15%o%one15%o%x-m.g%o%one15%o%onen)^2)/((one15%o%s.g%o%onen)^2)
fdelt=(1/sqrt(2*pi))*(1/(one15%o%s.g%o%onen))*exp(fexp/2)
grad=apply(fdelt/fmix,c(1,2),sum)-n
win.graph()
#Simple version of filled.contour that could be a default:
#filled.contour(m.g,s.g,grad,main="contour plot of gradient function",xlab="Mean",ylab="Sigma")
pgpos=pretty(grad[grad>0],10)
np=length(pgpos)
pgneg=pretty(grad[grad<=0],10)
ng=length(pgneg)
if(np==0){
y=filled.contour(m.g,s.g,grad,levels=pgneg,col=c(rainbow(ng,start=.7,end=.8)[ng:1]),main="Contour plot of gradient function",xlab="Mean",ylab="Sigma")
}
else {if(pgpos[1]==0){pgpos=pgpos+1}
y=filled.contour(m.g,s.g,grad,levels=c(pgneg,pgpos),col=c(rainbow(ng,start=.7,end=.8)[ng:1],terrain.colors(np)),main="Contour plot of gradient function",xlab="Mean",ylab="Sigma")
}
y
}