# DATA
H <- c(5, 9, 8, 4, 7)   # jumlah Head
T <- c(5, 1, 2, 6, 3)   # jumlah Tail
n <- H + T

data_koin <- c("5H5T","9H1T","8H2T","4H6T","7H3T")
# INISIALISASI PARAMETER

theta_A <- 0.7   # peluang head koin A
theta_B <- 0.5   # peluang head koin B

# batas konvergensi
tolerance <- 0.0001

# iterasi maksimum
max_iter <- 100

# menyimpan hasil iterasi
hasil_iterasi <- data.frame()
# ITERASI EM

for(iter in 1:max_iter){

  cat("ITERASI :", iter, "\n")

  # E-STEP
  # likelihood masing-masing koin
  likelihood_A <- dbinom(H, size = n, prob = theta_A)
  likelihood_B <- dbinom(H, size = n, prob = theta_B)

  # probabilitas posterior
  gamma_A <- likelihood_A / (likelihood_A + likelihood_B)
  gamma_B <- likelihood_B / (likelihood_A + likelihood_B)

  # expected head dan tail
  expected_H_A <- gamma_A * H
  expected_T_A <- gamma_A * T

  expected_H_B <- gamma_B * H
  expected_T_B <- gamma_B * T

  # TAMPILKAN HASIL E-STEP
  tabel_estep <- data.frame(
    Data = data_koin,

    Rasio_Koin_A = round(gamma_A,4),
    Rasio_Koin_B = round(gamma_B,4),

    Est_H_A = round(expected_H_A,4),
    Est_T_A = round(expected_T_A,4),

    Est_H_B = round(expected_H_B,4),
    Est_T_B = round(expected_T_B,4)
  )

  print(tabel_estep)
  
  # M-STEP
  theta_A_new <- sum(expected_H_A) /
               sum(expected_H_A + expected_T_A)

  theta_B_new <- sum(expected_H_B) /
               sum(expected_H_B + expected_T_B)

  # HITUNG PERUBAHAN PARAMETER
  perubahan_A <- abs(theta_A_new - theta_A)
  perubahan_B <- abs(theta_B_new - theta_B)

  # SIMPAN HASIL ITERASI
  hasil_iterasi <- rbind(
    hasil_iterasi,
    data.frame(
    Iterasi = iter,

    Theta_A_Lama = theta_A,
    Theta_A_Baru = theta_A_new,
    Perubahan_A = perubahan_A,

    Theta_B_Lama = theta_B,
    Theta_B_Baru = theta_B_new,
    Perubahan_B = perubahan_B
    )
  )

  # CEK KONVERGENSI
  perubahan_A <- abs(theta_A_new - theta_A)
  perubahan_B <- abs(theta_B_new - theta_B)

  if(perubahan_A < tolerance &&
     perubahan_B < tolerance){

    theta_A <- theta_A_new
    theta_B <- theta_B_new

    break
  }

  # update parameter
  theta_A <- theta_A_new
  theta_B <- theta_B_new
}
## ITERASI : 1 
##   Data Rasio_Koin_A Rasio_Koin_B Est_H_A Est_T_A Est_H_B Est_T_B
## 1 5H5T       0.2949       0.7051  1.4744  1.4744  3.5256  3.5256
## 2 9H1T       0.9254       0.0746  8.3282  0.9254  0.6718  0.0746
## 3 8H2T       0.8416       0.1584  6.7327  1.6832  1.2673  0.3168
## 4 4H6T       0.1520       0.8480  0.6080  0.9119  3.3920  5.0881
## 5 7H3T       0.6948       0.3052  4.8639  2.0845  2.1361  0.9155
## ITERASI : 2 
##   Data Rasio_Koin_A Rasio_Koin_B Est_H_A Est_T_A Est_H_B Est_T_B
## 1 5H5T       0.1801       0.8199  0.9007  0.9007  4.0993  4.0993
## 2 9H1T       0.9316       0.0684  8.3840  0.9316  0.6160  0.0684
## 3 8H2T       0.8291       0.1709  6.6328  1.6582  1.3672  0.3418
## 4 4H6T       0.0726       0.9274  0.2905  0.4358  3.7095  5.5642
## 5 7H3T       0.6336       0.3664  4.4353  1.9008  2.5647  1.0992
## ITERASI : 3 
##   Data Rasio_Koin_A Rasio_Koin_B Est_H_A Est_T_A Est_H_B Est_T_B
## 1 5H5T       0.1339       0.8661  0.6697  0.6697  4.3303  4.3303
## 2 9H1T       0.9421       0.0579  8.4792  0.9421  0.5208  0.0579
## 3 8H2T       0.8356       0.1644  6.6849  1.6712  1.3151  0.3288
## 4 4H6T       0.0461       0.9539  0.1842  0.2763  3.8158  5.7237
## 5 7H3T       0.6134       0.3866  4.2940  1.8403  2.7060  1.1597
## ITERASI : 4 
##   Data Rasio_Koin_A Rasio_Koin_B Est_H_A Est_T_A Est_H_B Est_T_B
## 1 5H5T       0.1150       0.8850  0.5752  0.5752  4.4248  4.4248
## 2 9H1T       0.9479       0.0521  8.5310  0.9479  0.4690  0.0521
## 3 8H2T       0.8410       0.1590  6.7279  1.6820  1.2721  0.3180
## 4 4H6T       0.0364       0.9636  0.1457  0.2185  3.8543  5.7815
## 5 7H3T       0.6059       0.3941  4.2416  1.8178  2.7584  1.1822
## ITERASI : 5 
##   Data Rasio_Koin_A Rasio_Koin_B Est_H_A Est_T_A Est_H_B Est_T_B
## 1 5H5T       0.1076       0.8924  0.5381  0.5381  4.4619  4.4619
## 2 9H1T       0.9504       0.0496  8.5534  0.9504  0.4466  0.0496
## 3 8H2T       0.8436       0.1564  6.7490  1.6873  1.2510  0.3127
## 4 4H6T       0.0329       0.9671  0.1314  0.1971  3.8686  5.8029
## 5 7H3T       0.6031       0.3969  4.2219  1.8094  2.7781  1.1906
## ITERASI : 6 
##   Data Rasio_Koin_A Rasio_Koin_B Est_H_A Est_T_A Est_H_B Est_T_B
## 1 5H5T       0.1048       0.8952  0.5238  0.5238  4.4762  4.4762
## 2 9H1T       0.9514       0.0486  8.5624  0.9514  0.4376  0.0486
## 3 8H2T       0.8447       0.1553  6.7580  1.6895  1.2420  0.3105
## 4 4H6T       0.0315       0.9685  0.1261  0.1891  3.8739  5.8109
## 5 7H3T       0.6021       0.3979  4.2146  1.8063  2.7854  1.1937
## ITERASI : 7 
##   Data Rasio_Koin_A Rasio_Koin_B Est_H_A Est_T_A Est_H_B Est_T_B
## 1 5H5T       0.1037       0.8963  0.5183  0.5183  4.4817  4.4817
## 2 9H1T       0.9518       0.0482  8.5659  0.9518  0.4341  0.0482
## 3 8H2T       0.8452       0.1548  6.7616  1.6904  1.2384  0.3096
## 4 4H6T       0.0310       0.9690  0.1240  0.1861  3.8760  5.8139
## 5 7H3T       0.6017       0.3983  4.2120  1.8051  2.7880  1.1949
## ITERASI : 8 
##   Data Rasio_Koin_A Rasio_Koin_B Est_H_A Est_T_A Est_H_B Est_T_B
## 1 5H5T       0.1033       0.8967  0.5163  0.5163  4.4837  4.4837
## 2 9H1T       0.9519       0.0481  8.5673  0.9519  0.4327  0.0481
## 3 8H2T       0.8454       0.1546  6.7630  1.6908  1.2370  0.3092
## 4 4H6T       0.0308       0.9692  0.1233  0.1849  3.8767  5.8151
## 5 7H3T       0.6016       0.3984  4.2110  1.8047  2.7890  1.1953

Observasi 5 head dan 5 tail lebih cocok berasal dari koin B. Karena hasilnya relatif seimbang. Koin A cenderung menghasilkan lebih banyak head.

9 head dari 10 lemparan sangat mendukung bahwa koin yang digunakan adalah koin A. Karena koin A memiliki probabilitas head lebih tinggi. Hampir seluruh 9 head dianggap berasal dari koin A (8.5673 head). Koin B hanya mendapat kontribusi kecil (0.4327 head).

Banyaknya head masih lebih konsisten dengan karakteristik koin. Tetapi tidak sekuat observasi 9H1T. Mayoritas head diberikan ke koin A. Sebagian kecil tetap dialokasikan ke koin B.

Jumlah tail lebih banyak daripada head. Sangat tidak sesuai dengan karakteristik koin A. Sangat cocok dengan koin B. Hampir seluruh observasi dikaitkan ke koin B.

7 head cukup tinggi untuk mendukung koin A. Tetapi belum terlalu ekstrem. Karena itu, EM membagi kontribusi ke kedua koin.

# HASIL AKHIR
cat("HASIL AKHIR ALGORITMA EM\n")
## HASIL AKHIR ALGORITMA EM
cat("Estimasi akhir theta A :", round(theta_A,6), "\n")
## Estimasi akhir theta A : 0.796735
cat("Estimasi akhir theta B :", round(theta_B,6), "\n\n")
## Estimasi akhir theta B : 0.519613
# TABEL RIWAYAT ITERASI
print(round(hasil_iterasi,6))
##   Iterasi Theta_A_Lama Theta_A_Baru Perubahan_A Theta_B_Lama Theta_B_Baru
## 1       1     0.700000     0.756609    0.056609     0.500000     0.525635
## 2       2     0.756609     0.779862    0.023254     0.525635     0.525155
## 3       3     0.779862     0.789991    0.010129     0.525155     0.522390
## 4       4     0.789991     0.794152    0.004161     0.522390     0.520786
## 5       5     0.794152     0.795784    0.001631     0.520786     0.520070
## 6       6     0.795784     0.796408    0.000625     0.520070     0.519776
## 7       7     0.796408     0.796645    0.000237     0.519776     0.519659
## 8       8     0.796645     0.796735    0.000090     0.519659     0.519613
##   Perubahan_B
## 1    0.025635
## 2    0.000481
## 3    0.002765
## 4    0.001604
## 5    0.000716
## 6    0.000294
## 7    0.000117
## 8    0.000046