###1. Menentukan Data Percobaan

# Jumlah Head dan Tail pada setiap set
H <- c(5, 9, 8, 4, 7)
T <- c(5, 1, 2, 6, 3)

set <- 1:5

data.frame(Set=set, Head=H, Tail=T)
##   Set Head Tail
## 1   1    5    5
## 2   2    9    1
## 3   3    8    2
## 4   4    4    6
## 5   5    7    3

###2. Parameter

theta_A <- 0.6
theta_B <- 0.4

theta_A
## [1] 0.6
theta_B
## [1] 0.4

###3. E-step

# Probabilitas hasil muncul jika menggunakan koin A
prob_A <- theta_A^H * (1-theta_A)^T

# Probabilitas hasil muncul jika menggunakan koin B
prob_B <- theta_B^H * (1-theta_B)^T

data.frame(
  Set=set,
  H=H,
  T=T,
  Prob_A=prob_A,
  Prob_B=prob_B
)
##   Set H T       Prob_A       Prob_B
## 1   1 5 5 0.0007962624 0.0007962624
## 2   2 9 1 0.0040310784 0.0001572864
## 3   3 8 2 0.0026873856 0.0002359296
## 4   4 4 6 0.0005308416 0.0011943936
## 5   5 7 3 0.0017915904 0.0003538944

###4. Rasio koin A dan B memunculkan E

gamma_A <- prob_A/(prob_A + prob_B)
gamma_B <- prob_B/(prob_A + prob_B)

data.frame(
  Set=set,
  P_A_given_E=gamma_A,
  P_B_given_E=gamma_B
)
##   Set P_A_given_E P_B_given_E
## 1   1   0.5000000  0.50000000
## 2   2   0.9624468  0.03755318
## 3   3   0.9192938  0.08070618
## 4   4   0.3076923  0.69230769
## 5   5   0.8350515  0.16494845

###5. Estimasi H dan T setiap koin

HA <- gamma_A * H
TA <- gamma_A * T

HB <- gamma_B * H
TB <- gamma_B * T

data.frame(
  Set=set,
  HA=HA,
  TA=TA,
  HB=HB,
  TB=TB
)
##   Set       HA        TA        HB         TB
## 1   1 2.500000 2.5000000 2.5000000 2.50000000
## 2   2 8.662021 0.9624468 0.3379786 0.03755318
## 3   3 7.354351 1.8385876 0.6456494 0.16141236
## 4   4 1.230769 1.8461538 2.7692308 4.15384615
## 5   5 5.845361 2.5051546 1.1546392 0.49484536

###6. Menjumlahkan Estimasi H dan T

sum_HA <- sum(HA)
sum_TA <- sum(TA)

sum_HB <- sum(HB)
sum_TB <- sum(TB)

sum_HA
## [1] 25.5925
sum_TA
## [1] 9.652343
sum_HB
## [1] 7.407498
sum_TB
## [1] 7.347657

###7. M-STEP

theta_A_new <- sum_HA / (sum_HA + sum_TA)
theta_B_new <- sum_HB / (sum_HB + sum_TB)

theta_A_new
## [1] 0.7261346
theta_B_new
## [1] 0.5020278

###8. Iterasi EM

theta_A <- 0.6
theta_B <- 0.4

for(i in 1:10){

  # E-Step
  prob_A <- theta_A^H * (1-theta_A)^T
  prob_B <- theta_B^H * (1-theta_B)^T
  
  gamma_A <- prob_A/(prob_A + prob_B)
  gamma_B <- prob_B/(prob_A + prob_B)

  # Estimasi
  HA <- sum(gamma_A * H)
  TA <- sum(gamma_A * T)
  
  HB <- sum(gamma_B * H)
  TB <- sum(gamma_B * T)

  # M-Step
  theta_A_new <- HA/(HA + TA)
  theta_B_new <- HB/(HB + TB)

  cat("Iterasi:", i,"\n")
  cat("Theta A =", theta_A_new,"\n")
  cat("Theta B =", theta_B_new,"\n\n")

  theta_A <- theta_A_new
  theta_B <- theta_B_new
}
## Iterasi: 1 
## Theta A = 0.7261346 
## Theta B = 0.5020278 
## 
## Iterasi: 2 
## Theta A = 0.767965 
## Theta B = 0.5193574 
## 
## Iterasi: 3 
## Theta A = 0.7851831 
## Theta B = 0.5207994 
## 
## Iterasi: 4 
## Theta A = 0.7922579 
## Theta B = 0.5202544 
## 
## Iterasi: 5 
## Theta A = 0.7950516 
## Theta B = 0.5198623 
## 
## Iterasi: 6 
## Theta A = 0.7961281 
## Theta B = 0.5196892 
## 
## Iterasi: 7 
## Theta A = 0.7965384 
## Theta B = 0.5196219 
## 
## Iterasi: 8 
## Theta A = 0.7966941 
## Theta B = 0.5195971 
## 
## Iterasi: 9 
## Theta A = 0.7967531 
## Theta B = 0.5195881 
## 
## Iterasi: 10 
## Theta A = 0.7967754 
## Theta B = 0.5195848

###9. Hasil

data.frame(
Set=set,
H=H,
T=T,
P_A=gamma_A,
P_B=gamma_B
)
##   Set H T       P_A        P_B
## 1   1 5 5 0.1030700 0.89693000
## 2   2 9 1 0.9519996 0.04800044
## 3   3 8 2 0.8454855 0.15451445
## 4   4 4 6 0.0307301 0.96926990
## 5   5 7 3 0.6015416 0.39845843