Bagian 1: Pembentukan Data Pengamatan

Langkah paling awal adalah menyusun hasil eksperimen pelemparan koin ke dalam sebuah tabel data. Kita memasukkan jumlah kemunculan sisi angka (Head) dan sisi gambar (Tail) dari lima sesi percobaan yang berbeda. Perlu diingat bahwa kita sama sekali tidak mengetahui koin mana yang menghasilkan data pada masing-masing sesi tersebut.

coin_data <- data.frame(
  H = c(5, 9, 8, 4, 7),
  T = c(5, 1, 2, 6, 3)
)
print(coin_data)
##   H T
## 1 5 5
## 2 9 1
## 3 8 2
## 4 4 6
## 5 7 3

Bagian 2: Implementasi EM dengan Parameter Standar

Kita akan memulai proses pencarian nilai optimal dengan menebak probabilitas awal. Koin A diberi nilai tebakan probabilitas 0.6 dan koin B diberi nilai 0.5. Algoritma kemudian akan berputar melewati fase E-Step dan M-Step secara bergantian. Sistem akan terus memperbarui tebakannya hingga pergerakan nilai parameter tersebut menjadi sangat kecil atau sudah mencapai titik konvergen.

theta_A <- 0.60
theta_B <- 0.50
tol <- 0.0001
change <- 1
iter <- 1

history_1 <- data.frame(Iteration = numeric(), Theta_A = numeric(), Theta_B = numeric(), Change = numeric())

while(change > tol) {
  prob_A <- numeric(nrow(coin_data))
  prob_B <- numeric(nrow(coin_data))
  
  for(i in 1:nrow(coin_data)) {
    H <- coin_data$H[i]
    T <- coin_data$T[i]
    
    like_A <- (theta_A^H) * ((1 - theta_A)^T)
    like_B <- (theta_B^H) * ((1 - theta_B)^T)
    
    prob_A[i] <- like_A / (like_A + like_B)
    prob_B[i] <- like_B / (like_A + like_B)
  }
  
  expected_H_A <- sum(prob_A * coin_data$H)
  expected_T_A <- sum(prob_A * coin_data$T)
  expected_H_B <- sum(prob_B * coin_data$H)
  expected_T_B <- sum(prob_B * coin_data$T)
  
  new_theta_A <- expected_H_A / (expected_H_A + expected_T_A)
  new_theta_B <- expected_H_B / (expected_H_B + expected_T_B)
  
  change <- sqrt((new_theta_A - theta_A)^2 + (new_theta_B - theta_B)^2)
  
  history_1 <- rbind(history_1, data.frame(Iteration = iter, Theta_A = new_theta_A, Theta_B = new_theta_B, Change = change))
  
  theta_A <- new_theta_A
  theta_B <- new_theta_B
  iter <- iter + 1
}

print(round(history_1, 4))
##    Iteration Theta_A Theta_B Change
## 1          1  0.7130  0.5813 0.1392
## 2          2  0.7453  0.5693 0.0345
## 3          3  0.7681  0.5495 0.0301
## 4          4  0.7832  0.5346 0.0212
## 5          5  0.7911  0.5263 0.0115
## 6          6  0.7945  0.5224 0.0052
## 7          7  0.7959  0.5207 0.0022
## 8          8  0.7965  0.5200 0.0009
## 9          9  0.7967  0.5198 0.0003
## 10        10  0.7967  0.5197 0.0001
## 11        11  0.7968  0.5196 0.0001

Visualisasi di bawah ini menampilkan grafik pergerakan nilai parameter dari tebakan awal 0.6 dan 0.5 hingga akhirnya menemukan posisi yang stabil.

plot(history_1$Iteration, history_1$Theta_A, type = "b", ylim = c(0, 1),
     xlab = "Iterasi", ylab = "Nilai Parameter", main = "Konvergensi (Awal: 0.6 & 0.5)")
lines(history_1$Iteration, history_1$Theta_B, type = "b", col = "red")
legend("bottomright", legend = c("Theta A", "Theta B"), col = c("black", "red"), lty = 1)

Bagian 3: Eksperimen EM dengan Parameter Bebas

Sebagai bentuk analisis tambahan, kita akan menguji algoritma ini menggunakan angka tebakan awal yang cukup jauh berbeda. Kita menetapkan nilai 0.8 untuk probabilitas koin A dan 0.2 untuk koin B. Pengujian ini bertujuan untuk membuktikan apakah tebakan yang ekstrem akan mengubah hasil akhir atau sekadar memengaruhi durasi algoritma dalam mencapai kestabilan.

theta_A_var <- 0.80
theta_B_var <- 0.20
change_var <- 1
iter_var <- 1

history_var <- data.frame(Iteration = numeric(), Theta_A = numeric(), Theta_B = numeric(), Change = numeric())

while(change_var > tol) {
  prob_A_var <- numeric(nrow(coin_data))
  prob_B_var <- numeric(nrow(coin_data))
  
  for(i in 1:nrow(coin_data)) {
    H <- coin_data$H[i]
    T <- coin_data$T[i]
    
    like_A <- (theta_A_var^H) * ((1 - theta_A_var)^T)
    like_B <- (theta_B_var^H) * ((1 - theta_B_var)^T)
    
    prob_A_var[i] <- like_A / (like_A + like_B)
    prob_B_var[i] <- like_B / (like_A + like_B)
  }
  
  expected_H_A_var <- sum(prob_A_var * coin_data$H)
  expected_T_A_var <- sum(prob_A_var * coin_data$T)
  expected_H_B_var <- sum(prob_B_var * coin_data$H)
  expected_T_B_var <- sum(prob_B_var * coin_data$T)
  
  new_theta_A_var <- expected_H_A_var / (expected_H_A_var + expected_T_A_var)
  new_theta_B_var <- expected_H_B_var / (expected_H_B_var + expected_T_B_var)
  
  # Baris yang diperbaiki
  change_var <- sqrt((new_theta_A_var - theta_A_var)^2 + (new_theta_B_var - theta_B_var)^2)
  
  history_var <- rbind(history_var, data.frame(Iteration = iter_var, Theta_A = new_theta_A_var, Theta_B = new_theta_B_var, Change = change_var))
  
  theta_A_var <- new_theta_A_var
  theta_B_var <- new_theta_B_var
  iter_var <- iter_var + 1
}

print(round(history_var, 4))
##    Iteration Theta_A Theta_B Change
## 1          1  0.7513  0.4355 0.2405
## 2          2  0.7778  0.4862 0.0572
## 3          3  0.7893  0.5059 0.0228
## 4          4  0.7938  0.5140 0.0093
## 5          5  0.7956  0.5173 0.0038
## 6          6  0.7963  0.5186 0.0015
## 7          7  0.7966  0.5192 0.0006
## 8          8  0.7967  0.5194 0.0003
## 9          9  0.7968  0.5195 0.0001
## 10        10  0.7968  0.5196 0.0000

Grafik kedua ini menggambarkan laju perubahan parameter yang dimulai dari angka tebakan ekstrem. Anda bisa mengamati bentuk kurvanya dan membandingkannya secara langsung dengan grafik konvergensi pada bagian sebelumnya.

plot(history_var$Iteration, history_var$Theta_A, type = "b", ylim = c(0, 1),
     xlab = "Iterasi", ylab = "Nilai Parameter", main = "Konvergensi (Awal: 0.8 & 0.2)")
lines(history_var$Iteration, history_var$Theta_B, type = "b", col = "red")
legend("bottomright", legend = c("Theta A", "Theta B"), col = c("black", "red"), lty = 1)