EM Algorithm

h <- c(5, 9, 8, 4, 7) # jumlah gambar/head
n <- c(10, 10, 10, 10, 10) # total lemparan per set

run_EM <- function(theta_A_init, theta_B_init, max_iter = 50, tol = 1e-5) {
  theta_A <- theta_A_init
  theta_B <- theta_B_init

# menyimpan riwayat iterasi
  history <- data.frame(
    Iterasi = 0,
    theta_A = theta_A,
    theta_B = theta_B,
    Perubahan = NA
  )
  
  for (iter in 1:max_iter) {
    # E-STEP
    likelihood_A <- (theta_A^h) * ((1 - theta_A)^(n - h))
    likelihood_B <- (theta_B^h) * ((1 - theta_B)^(n - h))
    
    p_Z_A <- likelihood_A / (likelihood_A + likelihood_B)
    p_Z_B <- likelihood_B / (likelihood_A + likelihood_B)
    
    # M-STEP
    theta_A_new <- sum(h * p_Z_A) / sum(n * p_Z_A)
    theta_B_new <- sum(h * p_Z_B) / sum(n * p_Z_B)
    
    # menghitung perubahan menggunakan Euclidean Distance
    change <- sqrt((theta_A - theta_A_new)^2 + (theta_B - theta_B_new)^2)
    
    # hasil iterasi
    history <- rbind(history, data.frame(
      Iterasi = iter,
      theta_A = theta_A_new,
      theta_B = theta_B_new,
      Perubahan = change
    ))
    
    # memperbarui nilai untuk iterasi berikutnya
    theta_A <- theta_A_new
    theta_B <- theta_B_new
    
    if (change < tol) {
      break
    }
  }
  
  return(history)
}

# Set 1: theta_A = 0.7, theta_B = 0.6
cat("Inisialisasi (0.7, 0.6) \n")

## Inisialisasi (0.7, 0.6)

hasil_set1 <- run_EM(theta_A_init = 0.7, theta_B_init = 0.6)
print(hasil_set1, digits = 4)

##    Iterasi theta_A theta_B Perubahan
## 1        0  0.7000  0.6000        NA
## 2        1  0.7289  0.5883 3.120e-02
## 3        2  0.7529  0.5660 3.275e-02
## 4        3  0.7730  0.5457 2.855e-02
## 5        4  0.7858  0.5322 1.867e-02
## 6        5  0.7923  0.5251 9.579e-03
## 7        6  0.7950  0.5219 4.231e-03
## 8        7  0.7961  0.5205 1.737e-03
## 9        8  0.7965  0.5200 6.927e-04
## 10       9  0.7967  0.5197 2.735e-04
## 11      10  0.7968  0.5196 1.077e-04
## 12      11  0.7968  0.5196 4.248e-05
## 13      12  0.7968  0.5196 1.679e-05
## 14      13  0.7968  0.5196 6.655e-06

# Set 2: theta_A = 0.7, theta_B = 0.4
cat("Inisialisasi (0.7, 0.4) \n")

## Inisialisasi (0.7, 0.4)

hasil_set2 <- run_EM(theta_A_init = 0.7, theta_B_init = 0.4)
print(hasil_set2, digits = 4)

##    Iterasi theta_A theta_B Perubahan
## 1        0  0.7000  0.4000        NA
## 2        1  0.7571  0.4759 9.500e-02
## 3        2  0.7813  0.5025 3.596e-02
## 4        3  0.7908  0.5126 1.389e-02
## 5        4  0.7945  0.5167 5.486e-03
## 6        5  0.7959  0.5184 2.193e-03
## 7        6  0.7964  0.5191 8.822e-04
## 8        7  0.7967  0.5194 3.561e-04
## 9        8  0.7967  0.5195 1.441e-04
## 10       9  0.7968  0.5195 5.844e-05
## 11      10  0.7968  0.5196 2.373e-05
## 12      11  0.7968  0.5196 9.652e-06

# Set 3: theta_A = 0.6, theta_B = 0.4
cat("Inisialisasi (0.6, 0.4) \n")

## Inisialisasi (0.6, 0.4)

hasil_set3 <- run_EM(theta_A_init = 0.6, theta_B_init = 0.4)
print(hasil_set3, digits = 4)

##    Iterasi theta_A theta_B Perubahan
## 1        0  0.6000  0.4000        NA
## 2        1  0.7261  0.5020 1.622e-01
## 3        2  0.7680  0.5194 4.528e-02
## 4        3  0.7852  0.5208 1.728e-02
## 5        4  0.7923  0.5203 7.096e-03
## 6        5  0.7951  0.5199 2.821e-03
## 7        6  0.7961  0.5197 1.090e-03
## 8        7  0.7965  0.5196 4.158e-04
## 9        8  0.7967  0.5196 1.577e-04
## 10       9  0.7968  0.5196 5.969e-05
## 11      10  0.7968  0.5196 2.258e-05
## 12      11  0.7968  0.5196 8.538e-06

# Insight: Mengubah parameter awal pada Algoritma EM tidak mengubah hasil kesimpulan akhir (nilai estimasi parameter tetap sama), namun sangat memengaruhi efisiensi komputasi (jumlah iterasi). Semakin dekat dan semakin tegas perbedaan tebakan awal kita terhadap karakteristik asli data, maka algoritma akan konvergen semakin cepat.