Algoritma Expectation-Maximization (EM) merupakan metode iteratif yang digunakan untuk melakukan estimasi parameter pada data yang memiliki nilai hilang (missing values) atau variabel laten (latent variables). Algoritma ini banyak digunakan dalam bidang statistika, machine learning, pengolahan citra, analisis klaster, hingga pemodelan campuran distribusi (mixture models).
Secara umum, algoritma EM bekerja dengan cara memperkirakan informasi yang tidak diketahui, kemudian menggunakan hasil perkiraan tersebut untuk memperbarui estimasi parameter model secara berulang hingga mencapai konvergensi. Algoritma EM terdiri dari dua tahap utama yang dilakukan secara iteratif, yaitu:
Pada tahap ini dilakukan perhitungan nilai harapan dari variabel yang tidak diketahui berdasarkan parameter sementara yang dimiliki saat ini. Dengan kata lain, algoritma mencoba “mengisi” atau memperkirakan informasi yang hilang menggunakan estimasi parameter sebelumnya.
Setelah nilai harapan diperoleh, tahap berikutnya adalah memperbarui estimasi parameter dengan memaksimalkan fungsi likelihood berdasarkan hasil pada tahap E-Step.
Proses E-Step dan M-Step dilakukan berulang kali hingga perubahan nilai parameter sangat kecil atau telah mencapai kondisi konvergen.
Misalkan terdapat dua buah koin, yaitu Koin A dan Koin B. Kedua koin memiliki probabilitas muncul sisi Head (H) yang berbeda, tetapi probabilitas tersebut belum diketahui. Dilakukan percobaan pelemparan koin sebanyak 5 set. Pada setiap set, salah satu koin dipilih secara acak, kemudian dilempar sebanyak 10 kali. Akan tetapi, identitas koin yang digunakan pada setiap set tidak dicatat sehingga menjadi variabel laten.
Hasil percobaan yang diperoleh adalah sebagai berikut:
Percobaan <- data.frame(Set = 1:5,
Banyak_Head = c(5,9,8,4,7),
Banyak_Tail = c(5,1,2,6,3))
Percobaan## Set Banyak_Head Banyak_Tail
## 1 1 5 5
## 2 2 9 1
## 3 3 8 2
## 4 4 4 6
## 5 5 7 3Jika diketahui parameter awal (probabilitas munculnya head pada setiap koin): \[\theta{}_A = 0.6\] \[\theta{}_B = 0.5\] Tentukan estimasi probabilitas munculnya sisi Head pada masing-masing koin menggunakan algoritma Expectation-Maximization (EM) hingga mencapai konvergensi.
Dari soal, kita telah mendapatkan informasi parameter awal, yaitu \(\theta{}_A = 0.6\) dan \(\theta{}_B = 0.5\).
# Parameter Awal
Theta_A = 0.6
Theta_B = 0.5Pada step ini digunakan teori probabilitas binomial karena kasus pelemparan koin merupakan kasus binomial sehingga: \[P(E|Z_x)=\frac{n!}{h!\cdot{}(n-h)!} \cdot \theta_x^h \cdot (1-\theta_x)^{n-h}\] Ket:
\(P(E|Z_x):\) Probabilitas munculnya E dalam setipa koin x
\(n:\) Total lemparan dalam satu set E
\(h:\) Total munculnya h dalam satu set E
\(\theta_x:\) Probabilitas munculnya h menggunakan koin x
\[P(Z_x|E)=\frac{P(E|Z_x)\cdot{}P(Z_x)}{P(E|Z_x)\cdot{}P(Z_x)+P(E|Z_y)\cdot{}P(Z_y)}\] Ket:
\(P(Z_x|E):\) Probabilitas pemilihan koin x dengan syarat E
\(P(Z_X):\) Probabilitas pemilihan koin x
Maka:
\[P(Z_A|E)=\frac{P(E|Z_A)\cdot{}P(Z_A)}{P(E|Z_A)\cdot{}P(Z_A)+P(E|Z_B)\cdot{}P(Z_B)}\] Karena pada soal hanya terdapat dua koin, maka probabilitas adil pada \(P(Z_A)\) atau \(P(Z_B)\) adalah 0.5 sehingga: \[P(Z_A|E)=\frac{P(E|Z_A)}{P(E|Z_A)+P(E|Z_B)}\]
Percobaan## Set Banyak_Head Banyak_Tail
## 1 1 5 5
## 2 2 9 1
## 3 3 8 2
## 4 4 4 6
## 5 5 7 3# Peluang Za dengan syarat (E = 2)
Pza_E1 = dbinom(x = 9, size = 10, prob = 0.6)/(dbinom(x = 9, size = 10, prob = 0.6)+dbinom(x = 9, size = 10, prob = 0.5))
Pza_E1## [1] 0.8049855# Tampilkan hasil
for (i in 1:length(Percobaan$Set)) {
a <- Percobaan$Banyak_Head[i]
# Peluang Za dengan syarat E
Percobaan$Peluang_Koin_A[i] = dbinom(x = a, size = 10, prob = 0.6)/
(dbinom(x = a, size = 10, prob = 0.6)+dbinom(x = a, size = 10, prob = 0.5))
# Peluang Zb dengan syarat E
Percobaan$Peluang_Koin_B[i] = dbinom(x = a, size = 10, prob = 0.5)/
(dbinom(x = a, size = 10, prob = 0.6)+dbinom(x = a, size = 10, prob = 0.5))
}
Percobaan## Set Banyak_Head Banyak_Tail Peluang_Koin_A Peluang_Koin_B
## 1 1 5 5 0.4491489 0.5508511
## 2 2 9 1 0.8049855 0.1950145
## 3 3 8 2 0.7334672 0.2665328
## 4 4 4 6 0.3521561 0.6478439
## 5 5 7 3 0.6472151 0.3527849# Estimasi H setiap koin
Percobaan$Est_H_K_A <- Percobaan$Banyak_Head * Percobaan$Peluang_Koin_A
Percobaan$Est_H_K_B <- Percobaan$Banyak_Head * Percobaan$Peluang_Koin_B
Percobaan## Set Banyak_Head Banyak_Tail Peluang_Koin_A Peluang_Koin_B Est_H_K_A Est_H_K_B
## 1 1 5 5 0.4491489 0.5508511 2.245745 2.754255
## 2 2 9 1 0.8049855 0.1950145 7.244870 1.755130
## 3 3 8 2 0.7334672 0.2665328 5.867737 2.132263
## 4 4 4 6 0.3521561 0.6478439 1.408625 2.591375
## 5 5 7 3 0.6472151 0.3527849 4.530506 2.469494\[\theta_A^{'}=\frac{\text{Total Estimasi H Koin A}}{\text{Banyak pelemparan setiap set} \cdot{} \text{Total peluang koin A}}\]
\[\theta_B^{'}=\frac{\text{Total Estimasi H Koin B}}{\text{Banyak pelemparan setiap set} \cdot{} \text{Total peluang koin B}}\]
# Theta A Aksen
Theta_A_Aksen <- sum(Percobaan$Est_H_K_A)/(10*sum(Percobaan$Peluang_Koin_A))
# Theta B Aksen
Theta_B_Aksen <- sum(Percobaan$Est_H_K_B)/(10*sum(Percobaan$Peluang_Koin_B))
cat("Theta A Aksen =", Theta_A_Aksen, "\n")## Theta A Aksen = 0.7130122cat("Theta B Aksen =", Theta_B_Aksen, "\n")## Theta B Aksen = 0.5813393Kita kembali lagi ke expectation step namun kita gunakan parameter yang baru didapatkan saat maximization step. Kita lakukan itu secara berulang (\(E \to M \to E \to M \to....\)) sampai batas iterasi atau kondisi konvergen (Hasil M tidak banyak berubah)
# Pembuatan Data
Iterasi = data.frame(Iterasi = 0,
Theta_A = 0.6,
Theta_B = 0.5)Theta_A <- numeric()
Theta_B <- numeric()
Theta_A <- 0.6
Theta_B <- 0.5
for(I in 1:20){
# Buat Data Frame per Iterasi
percobaan <- data.frame(Set = 1:5,
Banyak_Head = c(5,9,8,4,7),
Banyak_Tail = c(5,1,2,6,3))
# Membuat kolom Peluang setiap koin
for (i in 1:length(percobaan$Set)) {
a <- percobaan$Banyak_Head[i]
# Peluang Za dengan syarat E
percobaan$Peluang_Koin_A[i] = dbinom(x = a, size = 10, prob = Theta_A)/
(dbinom(x = a, size = 10, prob = Theta_A)+dbinom(x = a, size = 10, prob = Theta_B))
# Peluang Zb dengan syarat E
percobaan$Peluang_Koin_B[i] = dbinom(x = a, size = 10, prob = Theta_B)/
(dbinom(x = a, size = 10, prob = Theta_A)+dbinom(x = a, size = 10, prob = Theta_B))
}
# Estimasi H setiap koin
percobaan$Est_H_K_A <- percobaan$Banyak_Head * percobaan$Peluang_Koin_A
percobaan$Est_H_K_B <- percobaan$Banyak_Head * percobaan$Peluang_Koin_B
# Theta A Aksen
Theta_A_Aksen <- sum(percobaan$Est_H_K_A)/(10*sum(percobaan$Peluang_Koin_A))
# Theta B Aksen
Theta_B_Aksen <- sum(percobaan$Est_H_K_B)/(10*sum(percobaan$Peluang_Koin_B))
# Masukan ke dalam tabel
Iterasi = rbind(Iterasi, data.frame(Iterasi = I,
Theta_A = Theta_A_Aksen,
Theta_B = Theta_B_Aksen))
# Mengganti Theta dengan Theta baru
Theta_A <- Theta_A_Aksen
Theta_B <- Theta_B_Aksen
}Iterasi## Iterasi Theta_A Theta_B
## 1 0 0.6000000 0.5000000
## 2 1 0.7130122 0.5813393
## 3 2 0.7452920 0.5692558
## 4 3 0.7680988 0.5495359
## 5 4 0.7831646 0.5346175
## 6 5 0.7910552 0.5262812
## 7 6 0.7945325 0.5223904
## 8 7 0.7959287 0.5207299
## 9 8 0.7964656 0.5200472
## 10 9 0.7966683 0.5197704
## 11 10 0.7967441 0.5196587
## 12 11 0.7967724 0.5196136
## 13 12 0.7967829 0.5195954
## 14 13 0.7967868 0.5195881
## 15 14 0.7967882 0.5195851
## 16 15 0.7967888 0.5195839
## 17 16 0.7967890 0.5195835
## 18 17 0.7967890 0.5195833
## 19 18 0.7967891 0.5195832
## 20 19 0.7967891 0.5195831
## 21 20 0.7967891 0.5195831Perubahan setiap iterasinya dapat dihitung menggunakan persamaan umum seperti Eucludian. \[d(q,p)=\sqrt{(q_1-p_1)^2+(q_2-p_2)^2+...+(q_n-p_n)^2}\]
# Buat kolom perubahan
Iterasi$Perubahan <- 0
for(a in 2:length(Iterasi$Iterasi)){
Iterasi$Perubahan[a] <- sqrt((Iterasi$Theta_A[a-1]-Iterasi$Theta_A[a])^2 + (Iterasi$Theta_B[a-1]-Iterasi$Theta_B[a])^2)
}
Iterasi## Iterasi Theta_A Theta_B Perubahan
## 1 0 0.6000000 0.5000000 0.000000e+00
## 2 1 0.7130122 0.5813393 1.392403e-01
## 3 2 0.7452920 0.5692558 3.446735e-02
## 4 3 0.7680988 0.5495359 3.014999e-02
## 5 4 0.7831646 0.5346175 2.120229e-02
## 6 5 0.7910552 0.5262812 1.147851e-02
## 7 6 0.7945325 0.5223904 5.218174e-03
## 8 7 0.7959287 0.5207299 2.169478e-03
## 9 8 0.7964656 0.5200472 8.685631e-04
## 10 9 0.7966683 0.5197704 3.430641e-04
## 11 10 0.7967441 0.5196587 1.350369e-04
## 12 11 0.7967724 0.5196136 5.318122e-05
## 13 12 0.7967829 0.5195954 2.098702e-05
## 14 13 0.7967868 0.5195881 8.303624e-06
## 15 14 0.7967882 0.5195851 3.294411e-06
## 16 15 0.7967888 0.5195839 1.310625e-06
## 17 16 0.7967890 0.5195835 5.227978e-07
## 18 17 0.7967890 0.5195833 2.090695e-07
## 19 18 0.7967891 0.5195832 8.380802e-08
## 20 19 0.7967891 0.5195831 3.367047e-08
## 21 20 0.7967891 0.5195831 1.355530e-08Kita bisa lihat bahwa semakin dilakukan interaksi, semakin kecil pula perubahannya. Kita ambil batas iterasinya aja. Kita mendapatkan \(\theta_A = 0.796781\) dan \(\theta_B = 0.5195831\). Maka kemungkinan koin setiap set:
Percobaan$Peluang_H <- Percobaan$Banyak_Head/(Percobaan$Banyak_Head + Percobaan$Banyak_Tail)
Percobaan$Kemungkinan_Koin <- c("B","A","A","B","A")
library(dplyr)## Warning: package 'dplyr' was built under R version 4.4.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionPercobaan %>% select(Set, Banyak_Head, Banyak_Tail, Peluang_H, Kemungkinan_Koin)## Set Banyak_Head Banyak_Tail Peluang_H Kemungkinan_Koin
## 1 1 5 5 0.5 B
## 2 2 9 1 0.9 A
## 3 3 8 2 0.8 A
## 4 4 4 6 0.4 B
## 5 5 7 3 0.7 AKita coba bagaimana jadinya jika \(\theta_A = 0.3\) dan \(\theta_B = 0.6\).
# Data seperti contoh
Percobaan2 <- data.frame(Set = 1:5,
Banyak_Head = c(5,9,8,4,7),
Banyak_Tail = c(5,1,2,6,3))
Percobaan2## Set Banyak_Head Banyak_Tail
## 1 1 5 5
## 2 2 9 1
## 3 3 8 2
## 4 4 4 6
## 5 5 7 3# Pembuatan Data
Iterasi2 = data.frame(Iterasi = 0,
Theta_A = 0.3,
Theta_B = 0.6)
Theta_A <- numeric()
Theta_B <- numeric()
Theta_A <- 0.3
Theta_B <- 0.6
for(I in 1:20){
# Buat Data Frame per Iterasi
percobaan <- data.frame(Set = 1:5,
Banyak_Head = c(5,9,8,4,7),
Banyak_Tail = c(5,1,2,6,3))
# Membuat kolom Peluang setiap koin
for (i in 1:length(percobaan$Set)) {
a <- percobaan$Banyak_Head[i]
# Peluang Za dengan syarat E
percobaan$Peluang_Koin_A[i] = dbinom(x = a, size = 10, prob = Theta_A)/
(dbinom(x = a, size = 10, prob = Theta_A)+dbinom(x = a, size = 10, prob = Theta_B))
# Peluang Zb dengan syarat E
percobaan$Peluang_Koin_B[i] = dbinom(x = a, size = 10, prob = Theta_B)/
(dbinom(x = a, size = 10, prob = Theta_A)+dbinom(x = a, size = 10, prob = Theta_B))
}
# Estimasi H setiap koin
percobaan$Est_H_K_A <- percobaan$Banyak_Head * percobaan$Peluang_Koin_A
percobaan$Est_H_K_B <- percobaan$Banyak_Head * percobaan$Peluang_Koin_B
# Theta A Aksen
Theta_A_Aksen <- sum(percobaan$Est_H_K_A)/(10*sum(percobaan$Peluang_Koin_A))
# Theta B Aksen
Theta_B_Aksen <- sum(percobaan$Est_H_K_B)/(10*sum(percobaan$Peluang_Koin_B))
# Masukan ke dalam tabel
Iterasi2 = rbind(Iterasi2, data.frame(Iterasi = I,
Theta_A = Theta_A_Aksen,
Theta_B = Theta_B_Aksen))
# Mengganti Theta dengan Theta baru
Theta_A <- Theta_A_Aksen
Theta_B <- Theta_B_Aksen
}
# Buat kolom perubahan
Iterasi2$Perubahan <- 0
for(a in 2:length(Iterasi2$Iterasi)){
Iterasi2$Perubahan[a] <- sqrt((Iterasi2$Theta_A[a-1]-Iterasi2$Theta_A[a])^2 +
(Iterasi2$Theta_B[a-1]-Iterasi2$Theta_B[a])^2)
}
# Cek Hasil
Iterasi2## Iterasi Theta_A Theta_B Perubahan
## 1 0 0.3000000 0.6000000 0.000000e+00
## 2 1 0.4505366 0.7147889 1.893086e-01
## 3 2 0.4960139 0.7647872 6.758708e-02
## 4 3 0.5106875 0.7843936 2.448927e-02
## 5 4 0.5160066 0.7920190 9.297270e-03
## 6 5 0.5181067 0.7949582 3.612399e-03
## 7 6 0.5189686 0.7960865 1.419868e-03
## 8 7 0.5193269 0.7965193 5.619025e-04
## 9 8 0.5194764 0.7966854 2.234263e-04
## 10 9 0.5195387 0.7967492 8.917514e-05
## 11 10 0.5195647 0.7967737 3.570773e-05
## 12 11 0.5195755 0.7967831 1.433970e-05
## 13 12 0.5195799 0.7967868 5.773719e-06
## 14 13 0.5195818 0.7967882 2.330222e-06
## 15 14 0.5195826 0.7967887 9.424513e-07
## 16 15 0.5195829 0.7967889 3.818903e-07
## 17 16 0.5195830 0.7967890 1.550030e-07
## 18 17 0.5195831 0.7967890 6.300482e-08
## 19 18 0.5195831 0.7967891 2.564232e-08
## 20 19 0.5195831 0.7967891 1.044759e-08
## 21 20 0.5195831 0.7967891 4.260723e-09Dari hasil di atas, kita mendapatkan \(\theta_A = 0.5195831\) dan \(\theta_B = 0.796781\).
Perbandingan = data.frame(Parameter_Awal = c("Sama seperti contoh", "Berbeda"),
Theta_A = c(Iterasi$Theta_A[21], Iterasi2$Theta_A[21]),
Theta_B = c(Iterasi$Theta_B[21], Iterasi2$Theta_B[21]))
Perbandingan## Parameter_Awal Theta_A Theta_B
## 1 Sama seperti contoh 0.7967891 0.5195831
## 2 Berbeda 0.5195831 0.7967891Dari tabel di atas terlihat bahwa ketika kita mengganti parameter awal dengan probabilitas yang berbeda, khususnya berbalik dari contoh (\(\theta_A < \theta_B\)) maka akan menyebabkan estimasi parameter yang terbalik. Jadi, algoritma ini sensitif terhadap nilai parameter awal. Kemungkinan koin dari ini adalah sebagai berikut.
Percobaan2$Peluang_H <- Percobaan2$Banyak_Head/(Percobaan2$Banyak_Head + Percobaan2$Banyak_Tail)
Percobaan2$Kemungkinan_Koin <- c("A","B","B","A","B")
data.frame(Kemungkinan_Koin_Awal = Percobaan$Kemungkinan_Koin,
Kemungkinan_Koin_Dengan_Parameter_Berbeda = Percobaan2$Kemungkinan_Koin)## Kemungkinan_Koin_Awal Kemungkinan_Koin_Dengan_Parameter_Berbeda
## 1 B A
## 2 A B
## 3 A B
## 4 B A
## 5 A B