GMM模型的R實現

預備知識：

友情提示：本程式碼配合GMM演算法原理中的步驟閱讀更佳哦！

本文分為一元高斯分佈的EM演算法以及多元高斯分佈的EM演算法，分別採用兩本書上的資料《統計學習方法》和《機器學習》。

一元高斯混合模型

步驟：

1、設定模型引數的初始值，以及給出測試的資料

data <- c(-67,-48,6,8,14,16,23,24,28,29,41,49,56,60,75)
a = c(0.5,0.5) #係數
u = c(0,1)     #初始均值
sigma2 =c(9,9) #初始方差

data資料來自《統計機器學習》

2：E步 迭代計算gamma值，也就是響應度

首先利用dnorm()計算高斯公式的密度值,sapply()來向量化帶入函式計算

p <- t(sapply(data,dnorm,u,sigma))
amatrix <- matrix(rep(a,N),nrow=N,byrow = T) #便於計算
r = (p * amatrix) / rowSums(p * amatrix)     
#r(gamma):分模型k對觀測資料data的響應度

M步 迭代更新均值、方差和模型係數

u <-  colSums(data*r)/colSums(r)
umatrix <- matrix(rep(u,N),nrow=N
 
,byrow = T)
sigma2 <- colSums(r*(mat-umatrix)^2)/colSums(r) 
a <- colSums(r)/length(data)

3：設定迭代停止條件

根據自己需要設定，這裡設定為均值和方差變化不大則停止（也可以設定迭代最大輪數）

4：根據每個觀測值對每個高斯分佈的響應度，響應度大則歸為該分佈（也就是簇）。

 cluster <- which(r == apply(r,1,max),arr.ind = T)

將上述過程寫成函式

GaussCluster <- function(data,k,a,u,sigma2
 
){
  #data:資料集,向量
  #k:聚類的個數，或高斯分佈的個數
  #a0:高斯分佈的先驗概率,選擇各個高斯分佈的概率,向量
  #u0:高斯分佈的初始均值
  #sigma2:一元高斯分佈即為方差,sigma為標準差

  mat = matrix(rep(data,k),ncol = k) #方便後續計算
  N = length(data)
  i=0
  while(TRUE){
    u0 <- u
    sigma  = sqrt(sigma2)
    sigma0 <- sigma
    umatrix <- matrix(rep(u,N),nrow=N,byrow = T)
    p <- t(sapply(data,dnorm,u,sigma)) 
    amatrix <- matrix(rep(a,N),nrow=N,byrow = T) 
    r = (p * amatrix) / rowSums(p * amatrix)     
    u <-  colSums(data*r)/colSums(r)
    sigma2 <- colSums(r*(mat-umatrix)^2)/colSums(r) 
    a <- colSums(r)/length(data)
    sumU <- sum((u-u0)^2)
    sumsigma <- sum((sigma-sigma0)^2)
    if((sumU <= 1e-4) & (sumsigma <= 1e-4)){
      break
    }
  }
  cluster <- which(r == apply(r,1,max),arr.ind = T)
  return(list(u = u,d = sigma2,a = a,cluster = cluster))
}

呼叫函式

GaussCluster(data,2,c(0.5,0.5),c(0,1),c(9,9))

最終得到的結果

多元高斯混合模型

多元高斯混合模型採用的是《機器學習》中的西瓜資料集。

多元高斯分佈中的密度值需要用mvtnorm包中的dmvnorm()來計算。

library(mvtnorm)
mulGaussCluster <- function(data,k,a,u,cov0){
  #data:資料集,向量或資料框
  #k:聚類的個數，或高斯分佈的個數
  #a:高斯分佈的先驗概率,選擇各個高斯分佈的概率,向量
  #u:高斯分佈的初始均值
  #cov0:多元高斯分佈的初始協方差陣
  N = nrow(data)
  covList = list()
  for(i in 1:k){
    covList[[i]] <- cov0
  }
  count=1
  while(TRUE){
    u0 <- u
    p = c()
    for(j in 1:k){
      prob = apply(data,1,dmvnorm,mean=u[j,],sigma = covList[[j]])
      p = cbind(p,prob)
    }
    amatrix <- matrix(rep(a,N),nrow=N,byrow = T)
    r = (p * amatrix) / rowSums(p * amatrix)  
    u <- t(r) %*% data / colSums(r) #求和可以轉化為向量相乘的形式，簡化計算
    for(j in 1:k){
      sigma = matrix(rep(0,4),ncol=2)
      for(i in 1:N){
        sigma = sigma + r[i,j] * (data[i,]-u0[j,]) %*% t(data[i,]-u0[j,])  #R中向量預設為列向量
      }
      covList[[j]]= sigma/sum(r[,j]) 
    }
    a <- colSums(r)/N
    count = count + 1
    if(count == 100){
      break
    }
  }
  cluster <- which(r == apply(r,1,max),arr.ind = T)
  cluster <- cluster[order(cluster[,1]),]
  return(list(u = u,covList = covList,a = a,cluster = cluster))
}

呼叫：

wamellondata <- read.csv('watermelon.csv')
data = as.matrix(wamellondata[,2:3])
a = c(1/3,1/3,1/3)
u = rbind(data[6,],data[22,],data[27,])
cov0 <- matrix(c(0.1,0,0,0.1),ncol = 2)
list=mulGaussCluster(data,k=3,a,u,cov0)
cluster = list$cluster

結果：

當然這時候還可以用ggplot2來繪製聚類圖

library(ggplot2)
ggplot(data=NULL,mapping=aes(x=data[,1],y=data[,2],group = cluster[,2])) + 
  geom_point(colour = cluster[,2]) + xlab('density')+ylab('sugarContent')+
  theme_minimal()

得到的圖：

參考網址：