要求：用EM演算法實現基於GMM的聚類演算法。

一、實驗資料

參考[1] 3.3.2章節。由兩個二維高斯分佈混合生成1000個數據，混合係數分別是0.4、0.6，均值和方差如下：

mu1=[-2,-2]

sigma1=[1.2, 0.5, 0.5, 1]

mean2=[2,2]

sigma2=[1.5, 0.7, 0.7, 1]

二、實驗過程、結果與分析

2.1 資料散點圖

2.2 用mclust包實現

R語言自帶mclust包可對混合高斯分佈實現EM聚類，cluster1有391個數據，cluster2有609個數據。

2.3 我的實現

演算法迭代47次後收斂，cluster1中有387個數據，cluster2中有613個數據。

2.4 分析

我的演算法得到的cluster1比mclust包得到的cluster1少了4個數據，故兩種演算法基本一致。將我的演算法得到的兩個高斯的均值、方差與預設值對比，也基本一致。

三、參考文獻

[1]Vipin Kumar. Data clustering algorithmsand applications

[2]Miin-Shen Yang. A robust EM clusteringalgorithm for Gaussian mixture models

四：原始碼（R程式 gmm2dim.R）

#------------------multivariate normaldistribution-------------------

library(MASS)

#gaussian 1

mean1<-c(-2,-2)

sigma1<-matrix(c(1.2,0.5,

                 0.5,1),nrow=2,ncol=2)

#gaussian 2

mean2<-c(2,2)

sigma2<-matrix(c(1.5,0.7,

                 0.7,1),nrow=2,ncol=2)

#The number of samples from the mixturedistribution

N = 1000           

#Sample N random uniforms U (均勻分佈 U值介於[0,1],N個U)

U = runif(N)

#matrix(dim:N,2) to store the samples fromthe mixture distribution                                         

samples = matrix(NA,nrow=N,ncol=2)

#Sampling from the mixture

for(i in 1:N){

   if(U[i]<.4){

       samples[i,] = mvrnorm(1,mean1,sigma1)

        }else{

       samples[i,] = mvrnorm(1,mean2,sigma2)

       }

}

plot(samples,xlim=c(-6,6),ylim=c(-6,6),xlab="",ylab="")

#----------------using the mclustpackage----------Gaussian finite mixture model fitted by EM algorithm

library(mclust)

mc<-Mclust(samples)

#----------------myalgorithm------------------------------

#-----initialization-----

 #themax iteration times

T = 100

 #means,attention:不能令mu1=mu2

mu1<-c(-1,-1)

mu2<-c(1,1)

 #covariances

sig1<-matrix(c(1,0,

               0,1),nrow=2,ncol=2)

sig2<-matrix(c(1,0,

               0,1),nrow=2,ncol=2)

 #mixing probabilities

pi1<- 0.5

pi2<- 1-pi1

 #probability matrix p(z(i)=j|x(i);pi,mu,sig)

prob<-matrix(rep(NA,N*M),nrow=N,ncol=M)

#-----EM------

 #dmvnorm(x,mu,sd) 生成多維正太分佈的密度值,包:mvtnorm

library(mvtnorm)

for(step in 1:T){

#----mu1賦值給oldmu1,在M step對mu1進行改變

 oldmu1<-mu1

 oldmu2<-mu2

 oldsig1<-sig1

 oldsig2<-sig2

 oldpi1<-pi1

 oldpi2<-pi2

#----E step

 for(i in 1:N){

   t1= dmvnorm(samples[i,],mu1,sig1) * pi1

   t2= dmvnorm(samples[i,],mu2,sig2) * pi2

  prob[i,1] = t1 / (t1+t2)

  prob[i,2] = t2 / (t1+t2)

 }

#----M step

#compute mu1,sig1,pi1,針對prob[,1]

 #--1--prob矩陣的第一列求和

  p1= sum(prob[,1])

 #--2--for i=1:N, prob[i,1]*samples[i,] [N*2]矩陣

  m1= prob[,1] * samples

 #c(m1的第一列求和,m1的第二列求和) [1*2]矩陣

  p2= c( sum(m1[,1]),sum(m1[,2]) )

 #--3--calculate sig1的分子

  p3= matrix(0,nrow=2,ncol=2)

 for(k in 1:N)

   p3 = p3 + prob[k,1] * (samples[k,] - mu1)%*%t(samples[k,] - mu1)

 #--4--

  mu1= p2 / p1

 sig1 = p3 / p1

  pi1= p1 / N

#compute mu2,sig2,pi2,針對prob[,2]

 #--1--prob矩陣的第二列求和

  p1= sum(prob[,2])

 #--2--for i=1:N, prob[i,2]*samples[i,] [N*2]矩陣

  m1= prob[,2] * samples

 #c(m1的第一列求和,m1的第二列求和) [1*2]矩陣

  p2= c( sum(m1[,1]),sum(m1[,2]) )

 #--3--calculate sig2的分子

  p3= matrix(0,nrow=2,ncol=2)

 for(k in 1:N)

   p3 = p3 + prob[k,2] * (samples[k,] - mu2)%*%t(samples[k,] - mu2)

 #--4--

  mu2= p2 / p1

 sig2 = p3 / p1

  pi2= p1 / N

#----判斷兩個cluster中點的數目

  c1= 0

 for( k in 1:N){

    if( prob[k,1] >= prob[k,2]){

       c1 = c1 + 1

    }

  }

#----whether get to convergence

 epsilo = 1e-10

  if(#均值 mu1,mu2每個量<epsilo

     abs(mu1[1]-oldmu1[1]) < epsilo &

     abs(mu1[2]-oldmu1[2]) < epsilo &

     abs(mu2[1]-oldmu2[1]) < epsilo &

      abs(mu2[2]-oldmu2[2]) < epsilo &

     #協方差矩陣 sig1,sig2每個量<epsilo

     abs(sig1[1,1]-oldsig1[1,1]) < epsilo &

     abs(sig1[1,2]-oldsig1[1,2]) < epsilo &

     abs(sig1[2,1]-oldsig1[2,1]) < epsilo &

     abs(sig1[2,2]-oldsig1[2,2]) < epsilo &

     abs(sig2[1,1]-oldsig2[1,1]) < epsilo &

     abs(sig2[1,2]-oldsig2[1,2]) < epsilo &

     abs(sig2[2,1]-oldsig2[2,1]) < epsilo &

     abs(sig2[2,2]-oldsig2[2,2]) < epsilo

     )break

 cat('Step',step,':mu1=',mu1,',mu2=',mu2,',sigma1=',sig1,',sigma2=',sig2,'c1=',c1,'c2=',N-c1,'\n')

}