EM Algorithm

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Data paparan (Zij) - 8 individu, 10 tahun paparan, kolom ke-11 adalah status infeksi akhir

data <- matrix(c(
  0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1,
  0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1,
  0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1,
  0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0,
  0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,
  1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
  0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0
), nrow = 8, byrow = TRUE)

Pemisahan data

Zij <- data[, 1:10]  # Data paparan (8 individu, 10 tahun)
H <- data[, 11]      # Status infeksi akhir

Inisialisasi probabilitas infeksi awal (P_j)

J <- ncol(Zij)  # Jumlah tahun paparan
n <- nrow(Zij)  # Jumlah individu
P_j <- rep(0.1, J)  # Probabilitas infeksi awal 
max_iter <- 100  # Jumlah iterasi maksimum
tol <- 1e-6      # Toleransi konvergensi

Algoritma EM

for (iter in 1:max_iter) {
  # E-step: Menghitung ekspektasi E(Xij | Hi, P_j)
  E_Xij <- matrix(0, n, J)  # Inisialisasi E(Xij) dengan nol
  
  for (i in 1:n) {
    if (H[i] == 1) {  # Hanya untuk individu yang akhirnya terinfeksi
      denom <- 1 - prod(1 - Zij[i, ] * P_j, na.rm = TRUE)  # Penyebut
      
      if (!is.na(denom) && denom > 0) {  # Pastikan denom bukan NaN atau nol
        E_Xij[i, ] <- Zij[i, ] * P_j / denom
      } else {
        E_Xij[i, ] <- 0  # Jika denom nol, tetapkan ke nol untuk menghindari NaN
      }
    }
  }
  
  # M-step: Update probabilitas infeksi P_j
  sum_Zij <- colSums(Zij, na.rm = TRUE)  # Menghitung jumlah paparan per tahun
  sum_E_Xij <- colSums(E_Xij, na.rm = TRUE)  # Menghitung jumlah ekspektasi
  
  new_P_j <- ifelse(sum_Zij > 0, sum_E_Xij / sum_Zij, 0)  # Update dengan kondisi valid
  
  # Pastikan nilai P_j tetap dalam rentang [0, 1]
  new_P_j <- pmax(0, pmin(new_P_j, 1))
  
   # Cetak hasil iterasi
  cat("Iterasi:", iter, "P_j:", round(new_P_j, 6), "\n")
  
  # Cek konvergensi
  if (sum(abs(new_P_j - P_j), na.rm = TRUE) < tol) {
    cat("Konvergen pada iterasi:", iter, "\n")
    break
  }
  
  P_j <- new_P_j  # Update nilai P_j untuk iterasi selanjutnya
}

## Iterasi: 1 P_j: 0 0.123001 0 0.184502 0.246002 0.5 0.369004 0 0.123001 0 
## Iterasi: 2 P_j: 0 0.097601 0 0.150734 0.329189 0.5 0.602937 0 0.097601 0 
## Iterasi: 3 P_j: 0 0.071701 0 0.097399 0.383639 0.5 0.779195 0 0.071701 0 
## Iterasi: 4 P_j: 0 0.050976 0 0.05552 0.418536 0.5 0.888316 0 0.050976 0 
## Iterasi: 5 P_j: 0 0.035674 0 0.029574 0.441532 0.5 0.946361 0 0.035674 0 
## Iterasi: 6 P_j: 0 0.024739 0 0.01523 0.457778 0.5 0.974695 0 0.024739 0 
## Iterasi: 7 P_j: 0 0.017029 0 0.007719 0.469778 0.5 0.988045 0 0.017029 0 
## Iterasi: 8 P_j: 0 0.011639 0 0.003884 0.47868 0.5 0.994299 0 0.011639 0 
## Iterasi: 9 P_j: 0 0.007906 0 0.001948 0.485172 0.5 0.997251 0 0.007906 0 
## Iterasi: 10 P_j: 0 0.005342 0 0.000975 0.489807 0.5 0.998662 0 0.005342 0 
## Iterasi: 11 P_j: 0 0.003596 0 0.000488 0.493055 0.5 0.999343 0 0.003596 0 
## Iterasi: 12 P_j: 0 0.002413 0 0.000244 0.495297 0.5 0.999676 0 0.002413 0 
## Iterasi: 13 P_j: 0 0.001616 0 0.000122 0.49683 0.5 0.999839 0 0.001616 0 
## Iterasi: 14 P_j: 0 0.001081 0 6.1e-05 0.49787 0.5 0.99992 0 0.001081 0 
## Iterasi: 15 P_j: 0 0.000722 0 3.1e-05 0.498572 0.5 0.99996 0 0.000722 0 
## Iterasi: 16 P_j: 0 0.000482 0 1.5e-05 0.499044 0.5 0.99998 0 0.000482 0 
## Iterasi: 17 P_j: 0 0.000322 0 8e-06 0.499361 0.5 0.99999 0 0.000322 0 
## Iterasi: 18 P_j: 0 0.000215 0 4e-06 0.499573 0.5 0.999995 0 0.000215 0 
## Iterasi: 19 P_j: 0 0.000143 0 2e-06 0.499715 0.5 0.999998 0 0.000143 0 
## Iterasi: 20 P_j: 0 9.5e-05 0 1e-06 0.49981 0.5 0.999999 0 9.5e-05 0 
## Iterasi: 21 P_j: 0 6.4e-05 0 0 0.499873 0.5 0.999999 0 6.4e-05 0 
## Iterasi: 22 P_j: 0 4.2e-05 0 0 0.499915 0.5 1 0 4.2e-05 0 
## Iterasi: 23 P_j: 0 2.8e-05 0 0 0.499943 0.5 1 0 2.8e-05 0 
## Iterasi: 24 P_j: 0 1.9e-05 0 0 0.499962 0.5 1 0 1.9e-05 0 
## Iterasi: 25 P_j: 0 1.3e-05 0 0 0.499975 0.5 1 0 1.3e-05 0 
## Iterasi: 26 P_j: 0 8e-06 0 0 0.499983 0.5 1 0 8e-06 0 
## Iterasi: 27 P_j: 0 6e-06 0 0 0.499989 0.5 1 0 6e-06 0 
## Iterasi: 28 P_j: 0 4e-06 0 0 0.499993 0.5 1 0 4e-06 0 
## Iterasi: 29 P_j: 0 2e-06 0 0 0.499995 0.5 1 0 2e-06 0 
## Iterasi: 30 P_j: 0 2e-06 0 0 0.499997 0.5 1 0 2e-06 0 
## Iterasi: 31 P_j: 0 1e-06 0 0 0.499998 0.5 1 0 1e-06 0 
## Iterasi: 32 P_j: 0 1e-06 0 0 0.499999 0.5 1 0 1e-06 0 
## Iterasi: 33 P_j: 0 0 0 0 0.499999 0.5 1 0 0 0 
## Konvergen pada iterasi: 33

Hasil estimasi probabilitas infeksi per tahun

cat("Probabilitas infeksi final per tahun:\n")

## Probabilitas infeksi final per tahun:

print(round(P_j, 6))

##  [1] 0.000000 0.000001 0.000000 0.000000 0.499999 0.500000 1.000000 0.000000
##  [9] 0.000001 0.000000

EM Algorithm

Siskarossa Ika Oktora

2025-02-18