Mixture 응용 - 클러스터링

Clustering

먼저 원 데이터의 산점도를 확인해보자.

## Using mixture models for clustering in the iris dataset
## Compare k-means clustering and a location and scale mixture model with K normals

### Loading data and setting up global variables
rm(list=ls())
library(mclust)

## Package 'mclust' version 5.4.10
## Type 'citation("mclust")' for citing this R package in publications.

library(mvtnorm)

## 
## 다음의 패키지를 부착합니다: 'mvtnorm'

## The following object is masked from 'package:mclust':
## 
##     dmvnorm

### Defining a custom function to create pair plots
### This is an alternative to the R function pairs() that allows for 
### more flexibility. In particular, it allows us to use text to label 
### the points
pairs2 = function(x, col="black", pch=16, labels=NULL, names = colnames(x)){
  n = dim(x)[1]
  p = dim(x)[2]
  par(mfrow=c(p,p))
  for(k in 1:p){
    for(l in 1:p){
      if(k!=l){
        par(mar=c(3,3,1,1)+0.1)
        plot(x[,k], x[,l], type="n", xlab="", ylab="")
        if(is.null(labels)){
          points(x[,k], x[,l], pch=pch, col=col)
        }else{
          text(x[,k], x[,l], labels=labels, col=col)
        }
      }else{
        plot(seq(0,5), seq(0,5), type="n", xlab="", ylab="", axes=FALSE)
        text(2.5,2.5,names[k], cex=1.2)
      }
    }
  }
}

## Setup data
# 산점도를 그려서 데이터를 확인한다.
data("diabetes")
x       = as.matrix(diabetes[,-1])
n       = dim(x)[1]
p       = dim(x)[2]       # Number of features
KK      = 3
epsilon = 0.0000001
par(mfrow=c(1,1))
par(mar=c(4,4,1,1))
colscale = c("black","blue","red")
# Chemical: black, c
# Normal: blue, n
# Overt: red, o

shortnam  = c("c","n","o")
pairs2(x, col=colscale[diabetes[,1]], labels=shortnam[as.numeric(diabetes[,1])])

다음으로 EM 알고리즘을 사용해서 Mixture 모델의 파라미터를 추정하고, 최종 Mixture 모델의 cc 값을 이용해서 군집 결과를 플롯한다.

# Initialize the parameters of the algorithm
set.seed(63252)
numruns = 15
v.sum   = array(0, dim=c(numruns, n, KK))
QQ.sum  = rep(0, numruns)

for(ss in 1:numruns){
  w   = rep(1,KK)/KK  #Assign equal weight to each component to start with
  mu  = rmvnorm(KK, apply(x,2,mean), 3*var(x))   #Cluster centers randomly spread over the support of the data
  Sigma      = array(0, dim=c(KK,p,p))  #Initial variances are assumed to be the same
  Sigma[1,,] = var(x)
  Sigma[2,,] = var(x)
  Sigma[3,,] = var(x)
  
  sw     = FALSE
  QQ     = -Inf
  QQ.out = NULL
  s      = 0
  
  while(!sw){
    ## E step
    v = array(0, dim=c(n,KK))
    for(k in 1:KK){  #Compute the log of the weights
      v[,k] = log(w[k]) + mvtnorm::dmvnorm(x, mu[k,], Sigma[k,,], log=TRUE) 
    }
    for(i in 1:n){
      v[i,] = exp(v[i,] - max(v[i,]))/sum(exp(v[i,] - max(v[i,])))  #Go from logs to actual weights in a numerically stable manner
    }
    
    ## M step
    w = apply(v,2,mean)
    mu = matrix(0, nrow=KK, ncol=p)
    for(k in 1:KK){
      for(i in 1:n){
        mu[k,]    = mu[k,] + v[i,k]*x[i,]
      }
      mu[k,] = mu[k,]/sum(v[,k])
    }
    Sigma = array(0,dim=c(KK, p, p))
    for(k in 1:KK){
      for(i in 1:n){
        Sigma[k,,] = Sigma[k,,] + v[i,k]*(x[i,] - mu[k,])%*%t(x[i,] - mu[k,])
      }
      Sigma[k,,] = Sigma[k,,]/sum(v[,k])
    }
    
    ##Check convergence
    QQn = 0
    for(i in 1:n){
      for(k in 1:KK){
        QQn = QQn + v[i,k]*(log(w[k]) + mvtnorm::dmvnorm(x[i,],mu[k,],Sigma[k,,],log=TRUE))
      }
    }
    if(abs(QQn-QQ)/abs(QQn)<epsilon){
      sw=TRUE
    }
    QQ = QQn
    QQ.out = c(QQ.out, QQ)
    s = s + 1
  }
  
  v.sum[ss,,] = v
  QQ.sum[ss]  = QQ.out[s]
  print(paste("ss =", ss))
}

## [1] "ss = 1"
## [1] "ss = 2"
## [1] "ss = 3"
## [1] "ss = 4"
## [1] "ss = 5"
## [1] "ss = 6"
## [1] "ss = 7"
## [1] "ss = 8"
## [1] "ss = 9"
## [1] "ss = 10"
## [1] "ss = 11"
## [1] "ss = 12"
## [1] "ss = 13"
## [1] "ss = 14"
## [1] "ss = 15"

# --------------------------------------------------------
## Cluster reconstruction under the mixture model
# 지시자를 결정하기 위해서 가장 큰 확률값을 고른다.
# v 는 관측치가 각 성분에 포함될 확률이다.
# v.sum 은 15번 반복에 대해서 각 관측치의 성분에 포함될 확률이다.
# QQ.sum 은 15개의 각 반복단계에서의 Q 함수값이다. 
# cc 는 어떤 성분이 가장 큰 확률값을 가지는지에 대한 지시자이다.
cc = apply(v.sum[which.max(QQ.sum),,], 1 ,which.max) 

# 관측치 플롯
# colscale = c("black","blue","red")
pairs2(x, col=colscale[cc], labels=cc)

adjustedRandIndex 는 특별한 함수이다. 군집 결과를 비교하고 원 데이터의 라벨과 모두 일치하면 1의 값을 반환한다.

# 서로 다른 값들로 구성된 벡터가 얼마나 유사한지를 평가하는 함수가 adjustedRandIndex 함수이다.
# 우리가 생성한 데이터 파티션 cc 가 진짜 데이터 파티션과 비교해서
# 얼마나 가까운지를 비교한다.
# adjustedRandIndex 는 1까지 값을 취할 수 있다.
ARImle = adjustedRandIndex(cc, as.numeric(diabetes[,1]))  # Higher values indicate larger agreement
ARImle

## [1] 0.6640181

다음은, K-Means clustering 알고리즘을 수행하고 플롯과 adjustedRandIndex 를 사용해서 군집결과를 비교한다.

# --------------------------------------------------------
## Cluster reconstruction under the K-means algorithm
# 위의 결과를 K-means clustering 결과와 비교한다.
# 아래 kmeans 함수 호출 시 클러스터 개수인 3을 전달한다.
# nstart: 알고리즘 반복 횟수
irisCluster <- kmeans(x, 3, nstart = numruns)
# colscale = c("black","blue","red")
pairs2(x, col=colscale[irisCluster$cluster], labels=irisCluster$cluster)

ARIkmeans = adjustedRandIndex(irisCluster$cluster, as.numeric(diabetes[,1]))
ARIkmeans

## [1] 0.3764845

결과를 보면, EM for mixture 의 값은 0.66 이고, K-Means 의 값은 0.38 가 된다.

Mixture 응용 - 클러스터링 - diabetes

Kong, Seokkyu

2022-09-09

Summary

Clustering