Algoritma Expectation-Maximization (EM) adalah metode iteratif untuk menemukan estimasi maximum likelihood dari parameter dalam model statistik, di mana model tersebut bergantung pada variabel laten (tersembunyi) yang tidak teramati.
Dalam kasus ini, kita memiliki dua koin (Koin A dan Koin B) dengan probabilitas munculnya angka (Head) yang tidak diketahui (\(\theta_A\) dan \(\theta_B\)). Kita melakukan percobaan lemparan koin namun tidak mengetahui koin mana yang dipilih untuk setiap set lemparan.
Berdasarkan materi, kita melakukan 5 set percobaan, di mana setiap set terdiri dari 10 lemparan:
# Menyiapkan data observasi
data_H <- c(5, 9, 8, 4, 7)
n_lemparan <- 10
set_names <- paste("Set", 1:5)
Pada tahap ini, kita menggunakan tebakan awal sesuai materi: \(\theta_A = 0.6\) dan \(\theta_B = 0.5\).
Tahapan yang dilakukan:
# Inisialisasi
theta_A <- 0.6
theta_B <- 0.5
iterasi <- 10
log_ppt <- data.frame(Iterasi = 0:iterasi, Theta_A = NA, Theta_B = NA)
log_ppt[1, 2:3] <- c(theta_A, theta_B)
for(i in 1:iterasi) {
# E-Step: Menghitung responsibilitas (P(Z|E))
p_E_ZA <- dbinom(data_H, n_lemparan, theta_A)
p_E_ZB <- dbinom(data_H, n_lemparan, theta_B)
mu_A <- p_E_ZA / (p_E_ZA + p_E_ZB)
mu_B <- p_E_ZB / (p_E_ZA + p_E_ZB)
# M-Step: Update parameter
theta_A <- sum(data_H * mu_A) / sum(mu_A * n_lemparan)
theta_B <- sum(data_H * mu_B) / sum(mu_B * n_lemparan)
log_ppt[i+1, 2:3] <- c(round(theta_A, 3), round(theta_B, 3))
}
knitr::kable(log_ppt, caption = "Hasil Iterasi dengan Parameter Awal PPT")
| Iterasi | Theta_A | Theta_B |
|---|---|---|
| 0 | 0.600 | 0.500 |
| 1 | 0.713 | 0.581 |
| 2 | 0.745 | 0.569 |
| 3 | 0.768 | 0.550 |
| 4 | 0.783 | 0.535 |
| 5 | 0.791 | 0.526 |
| 6 | 0.795 | 0.522 |
| 7 | 0.796 | 0.521 |
| 8 | 0.796 | 0.520 |
| 9 | 0.797 | 0.520 |
| 10 | 0.797 | 0.520 |
Kesimpulan Simulasi 1: Algoritma mencapai konvergensi pada \(\theta_A \approx 0.796\) dan \(\theta_B \approx 0.520\). Ini menunjukkan Koin A cenderung berat sebelah ke angka, sementara Koin B mendekati koin adil.
Kita mencoba parameter awal yang lebih ekstrem: \(\theta_A = 0.9\) dan \(\theta_B = 0.1\).
# Inisialisasi Baru
theta_A_new <- 0.9
theta_B_new <- 0.1
log_baru <- data.frame(Iterasi = 0:iterasi, Theta_A = NA, Theta_B = NA)
log_baru[1, 2:3] <- c(theta_A_new, theta_B_new)
for(i in 1:iterasi) {
# E-Step
p_E_ZA <- dbinom(data_H, n_lemparan, theta_A_new)
p_E_ZB <- dbinom(data_H, n_lemparan, theta_B_new)
mu_A <- p_E_ZA / (p_E_ZA + p_E_ZB)
mu_B <- p_E_ZB / (p_E_ZA + p_E_ZB)
# M-Step
theta_A_new <- sum(data_H * mu_A) / sum(mu_A * n_lemparan)
theta_B_new <- sum(data_H * mu_B) / sum(mu_B * n_lemparan)
log_baru[i+1, 2:3] <- c(round(theta_A_new, 3), round(theta_B_new, 3))
}
knitr::kable(log_baru, caption = "Hasil Iterasi dengan Parameter Awal Baru (0.9 & 0.1)")
| Iterasi | Theta_A | Theta_B |
|---|---|---|
| 0 | 0.900 | 0.100 |
| 1 | 0.756 | 0.434 |
| 2 | 0.779 | 0.485 |
| 3 | 0.790 | 0.506 |
| 4 | 0.794 | 0.514 |
| 5 | 0.796 | 0.517 |
| 6 | 0.796 | 0.519 |
| 7 | 0.797 | 0.519 |
| 8 | 0.797 | 0.519 |
| 9 | 0.797 | 0.520 |
| 10 | 0.797 | 0.520 |
Berikut adalah beberapa poin penting dari perbandingan kedua simulasi: