Lab6: Expected Maximization

1 Praktek 12: Ekspektasi Maksimalisasi

Dalam kelas ini kita akan menggunakan algoritma Ekspektasi Maksimalisasi untuk mengestimasi parameter Campuran Gaussian. Metode ini menyediakan pengklasifikasi yang tidak terjaga dan sangat berguna ketika distribusi Gaussian diasumsikan.

Misalkan kita ingin memodelkan parameter populasi yang diasumsikan menjadi salah satu dari dua populasi normal gaussian dengan probabilitas pencampuran \(\pi\). Artinya, \(x_i\sim N(\mu_1, \Sigma_1)\) dengan probabilitas \(\pi\), \(x_i\sim N(\mu_2, \Sigma_2)\) dengan probabilitas \((1-\pi)\). Dinotasikan \(\theta_1= (\mu_1, \Sigma_1)\) dan \(\theta_2= (\mu_2, \Sigma_2)\) parameter untuk tiap populasi dan \(\pi\) sebagai parameter pencampuran, \(\phi_1\) dan \(\phi_2\) untuk densitas tiap distribusi. Algoritma Ekspektasi Maksimalisasi dapat disimpulkan di langkah berikut.

1. Ambil tebakan awal untuk parameter \(\hat{\mu}_1, \hat{\mu}_2, \hat{\Sigma}_1, \hat{\Sigma}_2, \hat{\pi}\)

2. Langkah ekspektasi: Menghitung \[\hat{\gamma}_i=\dfrac{\hat{\pi}\phi_1(y_i)}{\hat{\pi}\phi_1(y_i)+(1-\hat{\pi})\phi_2(y_i)}\]

3. Langkah memaksimumkan: Menghitung Mean dan Varians \[\begin{align*} \hat{\mu_1} &= \dfrac{\sum_{i=1}^n\hat{\gamma}_ix_i}{\sum_{i=1}^n\hat{\gamma}_i} & \hat{\Sigma_1} = \dfrac{1}{n}\dfrac{\sum_{i=1}^n\hat{\gamma}_i(x_i-\mu_1)(x_i-\mu_1)^T}{\sum_{i=1}^n\hat{\gamma}_i}\\ \hat{\mu_2}& = \dfrac{\sum_{i=1}^n(1-\hat{\gamma}_i)x_i}{\sum_{i=1}^n(1-\hat{\gamma}_i)} & \hat{\Sigma_2} = \dfrac{1}{n}\dfrac{\sum_{i=1}^n\hat{\gamma}_i(x_i-\mu_2)(x_i-\mu_2)^T}{\sum_{i=1}^n(1-\hat{\gamma}_i)} \end{align*}\] dan probabilitas campuran \(\hat{\pi} = \dfrac{1}{n}\sum_{i=1}^n\hat{\gamma}_i\).

4. Iterasikan langkah 2 dan 3 hingga konvergen.

1.1 Latihan 1:

Simulasikan 300 sample dari Gaussian Mixture dengan probabilitas pencampuran setara dengan \(1/3\) sebagai berikut: \[\begin{equation*} Y_1 \sim N\bigg(\begin{bmatrix}1\\1\end{bmatrix},\begin{bmatrix}2&1\\1&1\end{bmatrix} \bigg) \qquad Y_2 \sim N\bigg(\begin{bmatrix}7\\7\end{bmatrix},\begin{bmatrix}2&2\\2&5\end{bmatrix} \bigg) \end{equation*}\]

1. Visualisasikan distribusinya dengan menggunakan scatter plot dan gunakan skala kontinu untuk memisahkan warna berdasarkan populasinya.

library(ggplot2)
library(mvtnorm)

Mu1 = c(1,1)
Mu2 = c(7,7)
Sigma1 = matrix(c(2, 1, 1, 1), 2,2)
Sigma2 = matrix(c(2, 2, 2, 5), 2,2)

pi = 1/3
n = 300

set.seed(2)
data = matrix(0, n, 2)
z = rep(0,n)
for (i in 1:n){
  z[i] = rbinom(1,1,pi)
  if (z[i] ==1){
    data[i,] = rmvnorm(1, Mu1,Sigma1)
  }else{
    data[i,] = rmvnorm(1, Mu2,Sigma2)
  }
}

to.plot = data.frame(x = data[,1], 
                  y = data[,2], 
                  class =  z)
ggplot(to.plot)+ aes(x, y, color = class)+
  geom_point()+geom_density_2d()

2. Mengimplementasikan algoritma EM untuk memperkirakan parameternya.

Bantuan: Buat tebakan untuk \(\mu_1\) dan \(\mu_2\) hanya dengan dua dari \(y_i\) dengan acak. Mulai \(\Sigma_1\) dan \(\Sigma_2\) sebagai sample matriks kovarians dan \(\pi\) dengan nilai \(0.5\). Kriteria berhenti dimana perbedaan antara dua nilai berturut-turut dari log ekspektasi lengkap kurang dari toleransi yang diinginkan \(|l_k-l_{k-1}|<\text{tol}\).

#initial values
pi = 0.5
mu1 = data[1,]
mu2 = data[150,]
sigma1 = cov(data)
sigma2 = cov(data)

library(SDMTools)

# EM algorithm
Alpha = NULL
tol = 10^-2
iter = 0
Q = 0
phi1 = dmvnorm(data, mu1, sigma1)
phi2 = dmvnorm(data, mu2, sigma2)
Q_ = sum(log(pi)+log(phi1)) +
  sum(log(pi)+log(phi2))

while (abs(Q-Q_)>=tol) {

  iter = iter+1
  Q = Q_

  # Expectation step
  alpha = pi*phi1/(pi*phi1+(1-pi)*phi2)
  Alpha = cbind(Alpha, alpha)

  # Maximization step
  mu1 = apply(data, 2, function(col) wt.mean(col, alpha))
  mu2 = apply(data, 2, function(col) wt.mean(col, 1-alpha))
  sigma1 = cov.wt(data, wt = alpha, method = "ML")$cov
  sigma2 = cov.wt(data, wt = 1-alpha, method = "ML")$cov
  pi = mean(alpha)

  phi1 = dmvnorm(data, mu1, sigma1)
  phi2 = dmvnorm(data, mu2, sigma2)
  Q_ = sum(log(pi)+log(phi1)) +
    sum(log(pi)+log(phi2))
}

3. Visualisasikan nilai yang diperoleh untuk probabilitas tiap poin iterasi terawal ke-8.

library(reshape)

## Warning: package 'reshape' was built under R version 4.0.2

to.plot = data.frame(x = data[,1],
                     y = data[,2],
                     iter = Alpha[,1:8])
to.plot = melt(to.plot, id.vars = c("x", "y"))

ggplot(to.plot)+aes(x,y, color = value)+
  geom_point()+facet_wrap(~variable, nrow = 2)

1.2 Latihan 2

Mari gunakan Algoritma EM untuk mengsegmentasikan sebuah gambar. Dalam kasus ini kita akan gunakan gambar melanoma.

library(OpenImageR)

## Warning: package 'OpenImageR' was built under R version 4.0.2

img = readImage("Melanoma.png")
imageShow(img)

Algoritme segmentasi gambar biasa mengubah larik tiga dimensi yaitu gambar menjadi matriks 3 kolom dengan jumlah baris sebanyak jumlah piksel.

img.vect = apply(img, 3, as.vector )
dim(img.vect)

## [1] 737280      3

Dalam kasus ini kita akan mencari proyeksi 1 dimensi dari matriks agar algoritma EM dapat digunakan untuk mengestimasi parameter distribusi yang menghasilkan proyeksi.

1. Implentasikan sebuah fungsi yang dapat memperkirakan parameter dari Gaussian dalam kasus 1 dimensi.

EM = function(data){
  
  set.seed(1)
  pi = 0.5
  mu1 = sample(data,1)
  mu2 = sample(data,1)
  sd1 = sd(data)
  sd2 = sd1
  
  tol = 10^-2
  iter = 0
  Q = 0
  phi1 = dnorm(data, mu1, sd1)
  phi2 = dnorm(data, mu2, sd2)
  Q_ = sum(log(pi)+log(phi1)) +
    sum(log(pi)+log(phi2))
  
  while (abs(Q-Q_)>=tol) {
    
    iter = iter+1
    Q = Q_
    
    # Expectation step 
    alpha = pi*phi1/(pi*phi1+(1-pi)*phi2)
    
    # Maximization step
    mu1 = wt.mean(data, alpha)
    mu2 = wt.mean(data, 1-alpha)
    sd1 = wt.sd(data, wt = alpha)
    sd2 = wt.sd(data, wt = 1-alpha)
    pi = mean(alpha)
    
    phi1 = dnorm(data, mu1, sd1)
    phi2 = dnorm(data, mu2, sd2)
    Q_ = sum(log(pi)+log(phi1)) +
      sum(log(pi)+log(phi2))
  }
    return(alpha)
 }

2. Diketahui dalam literatur bahwa pola melanoma biasanya diekspresikan dalam warna biru. Jalankan algoritma EM menggunakan informasi tersebut dan tampilkan hasil akhirnya.

# Blue array + EM
blue = img.vect[,3]
prob = EM(blue)
classification = (prob>0.5)+0
segmented.image = replicate(3,classification, simplify = T)
dim(segmented.image) = dim(img)  
imageShow(segmented.image)

3. Gunakan proyeksi yang mempertahankan variabilitas data sebanyak mungkin (komponen pertama di PCA) dan terapkan EM. Visualisasikan hasilnya.

# PCA + EM
s.img.vect = scale(img.vect)
sigma = cov(img.vect)
eig = eigen(sigma)
eigenvectors = eig$vectors[,1]
projection = s.img.vect%*% eigenvectors

prob = EM(projection)
classification = (prob>0.5)+0
segmented.image = replicate(3,classification, simplify = T)
dim(segmented.image) = dim(img)  
imageShow(segmented.image)

4. Gunakan proyeksi Pencahayaan \(L = (0.229, 0.588. 0.114)\) dan tampilkan hasil menerapkan algoritma EM.

# Luminance + EM
luminance = c(0.299, 0.587, 0.114)

projection = img.vect%*%luminance

prob = EM(projection)
classification = (prob>0.5)+0
segmented.image = replicate(3,classification, simplify = T)
dim(segmented.image) = dim(img)  
imageShow(segmented.image)

5. Pada akhirnya, kita akan gunakan varian dari algoritma Independent Histogram Pursuit. Algoritma ini menemukan proyeksi yang memaksimumkan amplitudo bimodal. sekali proyeksi ditemukan, kita lanjut ke contoh latihan sebelumnya.

# library(modes) # Could not find another package to replace it
# 
# # Independent Histogram Pursuit
# theta = expand.grid(seq(1,360,length.out = 50), 
#                     seq(1,360,length.out = 50))
# 
# w = matrix(0,nrow(theta),3)
# amplitude = rep(0, nrow(theta))
# for (i in 1:nrow(theta)){
#   w[i,] = c(sin(theta[i,1])*cos(theta[i,2]),
#             sin(theta[i,1])*sin(theta[i,2]),
#             cos(theta[i,1]))
#   amplitude[i] = bimodality_amplitude(img.vect%*%w[i,], fig = F)
# }
# 
# ind = which.max(amplitude)
# w.opt = w[ind,]
# projection = img.vect%*%w.opt
# 
# prob = EM(projection)
# classification = (prob>0.5)+0
# segmented.image = replicate(3,classification, simplify = T)
# dim(segmented.image) = dim(img)  
# imageShow(segmented.image)