[R][源码]EM算法实现基于高斯混合模型（GMM）的聚类

要求：用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')

}