Optimization Week 6: Expectation Maximization

1 Latihan 12 Exepektasi Maksimal

Di kelas ini kita akan menerapkan algoritma Expectation Maximization untuk memperkirakan parameter Campuran Gaussian. Metode ini menyediakan pengklasifikasi tanpa pengawasan yang sangat berguna ketika distribusi Gaussian diasumsikan.

Misalkan kita ingin memodelkan parameter populasi yang diasumsikan sebagai salah satu dari dua populasi normal gaussian dengan probabilitas pencampuran \(π\). Itu berarti \(x_i \sim N(μ_1,Σ_1)\) dengan probabilitas \(π,x_i \sim N(μ_2,Σ_2)\) dengan probabilitas \((1-π)\). Mari tunjukkan \(θ_1=(μ_1,Σ_1)\) dan \(θ_2=(μ_2,Σ_2)\) parameter setiap populasi dan \(π\) parameter pencampuran dan \(ϕ_1\) dan \(ϕ_2\) kepadatan setiap distribusi. Harapan Maksimalisasi Algoritma dapat diringkas dalam langkah-langkah berikut:

1.Ambil tebakan awal untuk parameter tersebut \(\hat{μ_1},\hat{μ_2},\hat{Σ_1},\hat{Σ_2},\hatπ\). 2.Langkah Harapan: Hitung \[\hat{γ_i}= \frac{\hatπ ϕ_1 (y_i)}{\hatπϕ_1 (y_i) + (1-\hatπ)ϕ_2(y_i)} \] 3.Langkah Maksimalisasi: Hitung mean dan varians tertimbang \[ \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 kemungkinan pencampuran \(\hat{π}=\frac{1}n ∑_i^n\hat{γ_i}\) 4.Ulangi langkah 2 dan 3 hingga konvergensi.

1.1 Latihan 1

Mensimulasikan 300 sampel Campuran Gaussian dengan probabilitas pencampuran sebesar \(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*} \]

a.Visualisasikan distribusinya dengan plot pencar dan gunakan skala berkelanjutan untuk memisahkan populasi berdasarkan warna.

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()

b.Ini mengimplementasikan algoritma EM untuk memperkirakan parameternya.

Petunjuk: Buat tebakan awal untuk \(μ_1\) dan \(μ_2\) cukup memilih dua \(y_1\) sembarang. Mulailah \(Σ_1\) dan \(Σ_2\) sebagai matriks kovarians sampel keseluruhan dan \(π\) pada nilai \(0.5\). Kriteria berhenti adalah bahwa perbedaan antara dua nilai berturut-turut dari kemungkinan log lengkap kurang dari toleransi yang diinginkan, \(|l_k-l_{k-1}|<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. Memvisualisasikan nilai yang diperoleh untuk probabilitas setiap titik dalam 8 iterasi pertama.

library(reshape)

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 melakukan segmentasi gambar. Dalam hal ini kita akan menggunakan gambar melanoma.

library(OpenImageR)
img = readImage("C://Users/IMELDA SIANTURI/Downloads/Melanomaa.jpg")
imageShow(img)

Algoritma 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 hal ini kita akan mencari proyeksi satu dimensi dari matriks tersebut sehingga algoritma EM dapat digunakan untuk mengestimasi parameter distribusi pada proyeksi yang dihasilkan.

1. Implement a function that estimates the parameters for a Gaussian blend in the one-dimensional case.

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. It is well known in the literature that the pattern of the melanoma is usually expressed in the color blue. Run the EM algorithm using that information and display the final result.

# 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 Luminance L = (0.229,0.588.0.114) dan tampilkan hasil penerapan 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. Terakhir, kami akan menggunakan varian dari algoritme Pengejaran Histogram Independen. Algoritma ini menemukan proyeksi yang memaksimalkan amplitudo bimodal. Setelah proyeksi itu ditemukan, kami melanjutkan ke latihan sebelumnya.

library(modes)

# 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)