# 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