zmes<-function(x,ppi,mui,sigi,N)
{
# Implementacia EM algoritmu pre odhad parametrov
# zmesi k jednorozmernych normalnych rozdeleni
# ppi ... inicialne vahy zmesi
# mui,sigi ... inicialne parametre komponent zmesi
# N je pocet iteracii algoritmu

f<-function(x,mu,sig){exp(-(x-mu)^2/(2*sig^2))/(sqrt(2*pi)*sig)}

pp<-ppi; mu<-mui; sig<-sigi
n<-length(x); k<-length(pp)
pie<-matrix(0,ncol=k,nrow=n)

for(t in 1:N){
  for(i in 1:n){
    s<-0
    for(l in 1:k) s<-s+pp[l]*f(x[i],mu[l],sig[l])
    for(l in 1:k) pie[i,l]<-pp[l]*f(x[i],mu[l],sig[l])/s
    }
  for(l in 1:k){
    pp[l]<-mean(pie[,l])
    mu[l]<-weighted.mean(x,pie[,l])
    sig[l]<-sqrt(weighted.mean((x-mu[l])^2,pie[,l]))
    }
 
  #Nasledujuci odsek je len na ilustraciu priebehu algoritmu
  print("New estimates"); print(pp); print(mu); print(sig)
  xx<-seq(min(x),max(x),by=0.1);yy<-0*xx
  for(ii in 1:length(xx)){
     for(l in 1:k) yy[ii]<-yy[ii]+pp[l]*f(xx[ii],mu[l],sig[l])
     }
  hist(x,freq=FALSE,ylim=c(0,0.4))
  lines(xx,yy,xlab="x",ylab="y",col="red",lwd=2)
  }

list(p=pp,mu=mu,sig=sig)
}

# Demo
mu<-c(-1.5,1.2);dev<-c(0.7,1.1);n<-10000
choice<-sample(1:2,n,replace=TRUE,prob=c(0.3,0.7))
x<-rnorm(n,mean=mu[choice],sd=dev[choice])
zmes(x,c(0.7,0.3),c(-3,3),c(3,3),80)

mu<-c(-2,1,5);dev<-c(0.6,1.4,0.7);n<-10000
choice<-sample(1:3,n,replace=TRUE,prob=c(0.2,0.3,0.5))
x<-rnorm(n,mean=mu[choice],sd=dev[choice])
zmes(x,rep(1/3,3),c(-1,1,3),c(1,2,3),30)