Praktikum Pemrograman Statistika

Expectation–Maximization: Dari Data Hilang hingga Estimasi Optimal

Pendahuluan

Dalam banyak permasalahan statistika modern, keberadaan data yang tidak lengkap—baik karena hilang, tidak teramati, maupun secara intrinsik bersifat laten—menjadi tantangan penting dalam melakukan pendugaan parameter. Metode klasik seperti Maximum Likelihood pada umumnya mensyaratkan ketersediaan data lengkap agar fungsi likelihood dapat dievaluasi dan dioptimasi secara langsung. Akan tetapi, ketika sebagian komponen data tidak tersedia, proses optimasi dapat menjadi sangat kompleks atau bahkan tidak dapat dilakukan secara analitik.

Algoritma Expectation–Maximization (EM) hadir sebagai pendekatan sistematis untuk menyelesaikan permasalahan tersebut. EM memungkinkan proses pendugaan parameter dilakukan meskipun sebagian informasi tidak terobservasi, dengan memanfaatkan struktur probabilistik antara data teramati (observed data) dan data tidak teramati (missing atau latent data). Intuisi utama metode ini adalah bahwa pendugaan parameter dapat disederhanakan apabila kita dapat “melengkapi” data yang hilang—bukan secara langsung, tetapi melalui ekspektasi berdasarkan parameter sementara. Pendekatan inilah yang membuat EM sangat fleksibel dan relevan dalam berbagai konteks.

Algoritma EM bekerja melalui dua langkah berulang. Langkah Ekspektasi (E-step) membentuk nilai harapan dari log-likelihood data lengkap dengan kondisi data yang teramati dan nilai parameter saat ini. Selanjutnya, Langkah Maksimisasi (M-step) memperbarui parameter dengan memaksimalkan nilai harapan tersebut. Proses ini menghasilkan urutan estimasi parameter yang secara teoretis dijamin tidak menurunkan nilai log-likelihood (monotonicity property), dan pada akhirnya mengarah pada nilai konvergen yang stabil.

Sejak diperkenalkan oleh Dempster, Laird, dan Rubin (1977), algoritma EM telah menjadi salah satu metode paling berpengaruh dalam statistika komputasional. Metode ini digunakan secara luas dalam berbagai bidang seperti pengelompokan berbasis model (misalnya Gaussian Mixture Models), pemrosesan sinyal, genetika, machine learning, analisis citra, hingga inferensi berbasis model laten. Keunggulan utamanya terletak pada kemampuan menangani struktur data yang rumit dengan pendekatan iteratif yang relatif sederhana dan mudah diimplementasikan.

Dasar Teoretis Algoritma EM

Algoritma Expectation–Maximization (EM) merupakan metode iteratif untuk melakukan pendugaan parameter ketika sebagian data tidak teramati atau bersifat laten. Ide utamanya adalah memaksimalkan expected complete-data log-likelihood alih-alih langsung memaksimalkan observed-data log-likelihood yang sering kali sulit dikerjakan secara analitik.

Pada bagian ini dibahas landasan teoretis yang menjadi dasar algoritma EM.

Observed Data, Missing Data, dan Complete Data

Misalkan:

• \(y\) = observed data
• \(x\) = missing atau latent data
• \(z = (x, y)\) = complete data

Jika distribusi lengkap \(f(z \mid \theta)\) lebih mudah dianalisis dibandingkan \(f(y \mid \theta)\), maka optimasi terhadap likelihood complete-data dapat memberikan pendekatan yang lebih sederhana.

Fungsi Likelihood dan Log-Likelihood

Likelihood data lengkap:

\[ L(\theta \mid z) = f(z \mid \theta). \]

Sementara likelihood data teramati:

\[ L(\theta \mid y) = f(y \mid \theta) = \int f(x, y \mid \theta)\, dx. \]

Integral ini sering sulit dihitung secara eksplisit karena struktur data hilang yang kompleks.

Fungsi Ekspektasi \(Q(\theta \mid \theta^{(i)})\)

Pada iterasi ke-\(i\), algoritma EM membentuk fungsi:

\[ Q(\theta \mid \theta^{(i)}) = \mathbb{E}_{x \mid y, \theta^{(i)}} \left[ \log f(x, y \mid \theta) \right]. \]

Fungsi ini merupakan nilai harapan log-likelihood data lengkap terhadap distribusi bersyarat \(f(x \mid y, \theta^{(i)})\).

Langkah E-step

E-step menghitung:

\[ Q(\theta \mid \theta^{(i)}) = \mathbb{E}\left[ \log f(x, y \mid \theta) \;\middle|\; y, \theta^{(i)} \right]. \]

Proses ini membutuhkan:

1. Bentuk log-likelihood data lengkap
2. Distribusi bersyarat \(f(x \mid y, \theta^{(i)})\).

Output E-step adalah fungsi \(Q(\theta \mid \theta^{(i)})\) yang akan digunakan pada M-step.

Langkah M-step

Pada M-step, parameter diperbarui dengan:

\[ \theta^{(i+1)} = \arg\max_{\theta} Q(\theta \mid \theta^{(i)}). \]

Jika \(Q\) memiliki bentuk sederhana, maka update parameter biasanya dapat diperoleh secara eksplisit.

Sifat Monotonik EM

EM menjamin bahwa:

\[ \ell(\theta^{(i+1)} \mid y) \;\ge\; \ell(\theta^{(i)} \mid y), \]

di mana \(\ell(\theta \mid y) = \log f(y \mid \theta)\).

Sifat ini diperoleh karena:

• E-step membentuk lower bound terhadap log-likelihood
• M-step memaksimalkan lower bound tersebut.

Dengan demikian, nilai log-likelihood tidak menurun di setiap iterasi, meskipun EM tidak selalu dijamin mencapai global maximum.

Interpretasi Probabilistik

EM dapat dipandang sebagai proses yang berulang antara:

• mengestimasi nilai ekspektasi data hilang (E-step), dan
  
• mengoptimasi parameter berbasis data lengkap yang diestimasi (M-step).

Interpretasi ini menjadikan EM sangat relevan untuk model yang mengandung struktur laten seperti mixture models, hidden Markov models, dan pemodelan hierarkis.

Studi Kasus: Pendugaan Parameter Multinomial Menggunakan Algoritma EM

Pada studi ini, Rao (1973) melakukan pembagian secara acak 197 hewan ke dalam empat kategori berdasarkan phenotype. Misalkan vektor jumlah pengamatan dituliskan sebagai:

\[ y = (y_1, y_2, y_3, y_4)^\top = (100, 12, 18, 29)^\top. \]

Peluang untuk setiap kategori diberikan oleh:

\[ \left( p_1 = \frac{1}{2} + \frac{\theta}{4},\;\; p_2 = 1 - \frac{\theta}{4},\;\; p_3 = 1 - \frac{\theta}{4},\;\; p_4 = \frac{\theta}{4} \right), \]

sehingga:

\[ y \sim \text{Multinomial}(p_1, p_2, p_3, p_4). \]

Fungsi kepekatan peluang (pdf) bagi \(y\) adalah:

\[ g(y \mid \theta) = \frac{y_1 + y_2 + y_3 + y_4}{y_1! \, y_2! \, y_3! \, y_4!} \left(\frac{1}{2} + \frac{\theta}{4}\right)^{y_1} \left(1 - \frac{\theta}{4}\right)^{y_2} \left(1 - \frac{\theta}{4}\right)^{y_3} \left(\frac{\theta}{4}\right)^{y_4}. \]

Log-likelihood bagi \(y\) diperoleh dengan:

\[ \log L(\theta \mid y) = c + y_1 \log(2+\theta) + (y_2 + y_3) \log(1-\theta) + y_4 \log \theta, \]

dengan \(c\) merupakan konstanta yang tidak bergantung pada \(\theta\).

Pemecahan kategori dan pembentukan data lengkap

Sekarang peneliti memecah kategori pertama menjadi dua kategori baru sehingga:

\[ x = (x_1, x_2, x_3, x_4, x_5)^\top, \]

dengan relasi:

\[ x_1 + x_2 = y_1,\quad x_3 = y_2,\quad x_4 = y_3,\quad x_5 = y_4. \]

Peluang kategori baru menjadi:

\[ \left( \frac12,\, \frac{\theta}{4},\, 1 - \frac{\theta}{4},\, 1 - \frac{\theta}{4},\, \frac{\theta}{4} \right), \]

sehingga:

\[ x \sim \text{Multinomial}(p_1,p_2,p_3,p_4,p_5). \]

Fungsi peluang bagi \(x\) adalah:

\[ f(x \mid \theta) = \frac{x_1+x_2+x_3+x_4+x_5}{x_1! x_2! x_3! x_4! x_5!} \left(\frac12\right)^{x_1} \left(\frac{\theta}{4}\right)^{x_2} \left(1-\frac{\theta}{4}\right)^{x_3} \left(1-\frac{\theta}{4}\right)^{x_4} \left(\frac{\theta}{4}\right)^{x_5}. \]

Dengan mengambil bagian yang mengandung \(\theta\):

\[ \log L(\theta \mid x) = c + x_1 \log\left(\tfrac12\right) + (x_2 + x_5)\log \theta + (x_3+x_4)\log(1-\theta). \]

Pembentukan complete-data log-likelihood

Karena \(x_1\) dan \(x_2\) tidak teramati, kita nyatakan data laten sebagai:

\[ z = (z_1, z_2) = (x_1, x_2). \]

Dengan substitusi \(x_3 = y_2\), \(x_4 = y_3\), dan \(x_5 = y_4\), diperoleh:

\[ \log L^c(\theta \mid y,z) = c + z_1 \log\left(\frac12\right) + (z_2 + y_4)\log\theta + (y_2+y_3)\log(1-\theta). \]

E-step: Menghitung nilai ekspektasi terhadap data laten

Definisi fungsi \(Q\):

\[ Q(\theta \mid \theta^{(0)}) = \mathbb{E}_{z \mid y,\theta^{(0)}}[\log L^c(\theta \mid y,z)]. \]

Expand:

\[ Q(\theta \mid \theta^{(0)}) = \log\left(\frac12\right)\mathbb{E}[z_1] + (\mathbb{E}[z_2]+y_4)\log\theta + (y_2+y_3)\log(1-\theta). \]

Untuk menghitung ekspektasi, diperlukan pdf dari \(z \mid y\):

\[ f(z \mid \theta, y) = \frac{h(y,z \mid \theta)}{g(y\mid \theta)}, \]

dengan:

\[ h(x,z \mid \theta) = L^c(\theta \mid y,z) = \frac{z_1+z_2+y_2+y_3+y_4}{z_1! z_2! y_2! y_3! y_4!} \left(\frac12\right)^{z_1} \left(\frac{\theta}{4}\right)^{z_2} \left(1-\frac{\theta}{4}\right)^{y_2} \left(1-\frac{\theta}{4}\right)^{y_3} \left(\frac{\theta}{4}\right)^{y_4}, \]

dan

\[ g(y\mid \theta) = \frac{y_1+y_2+y_3+y_4}{y_1! y_2! y_3! y_4!} \left(\frac12+\frac{\theta}{4}\right)^{y_1} \left(1-\frac{\theta}{4}\right)^{y_2} \left(1-\frac{\theta}{4}\right)^{y_3} \left(\frac{\theta}{4}\right)^{y_4}. \]

Menggabungkan, diperoleh bentuk sederhana:

\[ f(z \mid \theta,y) = \frac{y_1!}{z_1! z_2!} \left(\frac12\right)^{z_1} \left(\frac{\theta}{4}\right)^{z_2} \left(\frac{1}{2}+\frac{\theta}{4}\right)^{-y_1}. \]

Karena \(z_1+z_2 = y_1\), distribusinya dapat diturunkan menjadi:

\[ z_1 \sim \text{Binomial}\!\left( y_1,\; \frac{\frac12}{\frac12+\frac{\theta}{4}} \right), \qquad z_2 \sim \text{Binomial}\!\left( y_1,\; \frac{\frac{\theta}{4}}{\frac12+\frac{\theta}{4}} \right). \]

Menggunakan \(\mathbb{E}(X)=np\):

\[ \mathbb{E}[z_1] = y_1 \frac{\frac12}{\frac12+\frac{\theta^{(0)}}{4}}, \quad \mathbb{E}[z_2] = y_1 \frac{\frac{\theta^{(0)}}{4}}{\frac12+\frac{\theta^{(0)}}{4}}. \]

Dengan \(y_1 = 100\):

\[ \hat z_1 = 100 \left( \frac{\frac12}{\frac12+\frac{\theta^{(0)}}{4}} \right), \qquad \hat z_2 = 100 \left( \frac{\frac{\theta^{(0)}}{4}}{\frac12+\frac{\theta^{(0)}}{4}} \right). \]

Sehingga fungsi \(Q\) menjadi:

\[ Q(\theta \mid \theta^{(0)}) = \log\left(\frac12\right)\hat z_1 + (\hat z_2 + y_4)\log\theta + (y_2+y_3)\log(1-\theta). \]

M-step: Memaksimumkan Q untuk memperoleh update parameter

Parameter diperbarui melalui:

\[ \hat\theta^{(m+1)} = \arg\max_{\theta} Q(\theta \mid \theta^{(m)}). \]

Turunan terhadap \(\theta\):

\[ \frac{dQ}{d\theta} = \frac{\hat z_2 + y_4}{\theta} - \frac{y_2+y_3}{1-\theta}. \]

Setara nol:

\[ \frac{\hat z_2+y_4}{\theta} = \frac{y_2+y_3}{1-\theta}. \]

Solusi:

\[ \hat\theta^{(1)} = \frac{\hat z_2 + y_4}{\hat z_2 + y_4 + y_2 + y_3}. \]

Dalam notasi \(x\):

\[ \hat\theta^{(1)} = \frac{x_2 + x_5}{x_2 + x_5 + x_3 + x_4}. \]

Tahapan berikutnya dapat dilanjutkan secara iteratif hingga toleransi konvergensi, misalnya \(10^{-7}\).

## ---------------------------------------------------------
## EM algorithm untuk studi kasus multinomial Rao (1973)
## y = (y1, y2, y3, y4) = (100, 12, 18, 29)
## ---------------------------------------------------------

em_rao <- function(y1 = 100, y2 = 12, y3 = 18, y4 = 29,
                   theta_init = 0.5,   # nilai awal theta
                   tol = 1e-7,         # toleransi konvergensi
                   max_iter = 1000,    # maksimum iterasi
                   verbose = TRUE) {   # TRUE: tampilkan riwayat iterasi
  
  theta_old <- theta_init
  iter <- 0
  diff <- Inf
  
  if (verbose) {
    cat("Iter  theta        z1_hat       z2_hat\n")
  }
  
  while (diff > tol && iter < max_iter) {
    iter <- iter + 1
    
    ## -------------------------
    ## E-step:
    ## z1_hat = E[z1 | y, theta_old]
    ## z2_hat = E[z2 | y, theta_old]
    ## -------------------------
    denom <- 0.5 + theta_old / 4
    
    z1_hat <- y1 * (0.5       / denom)
    z2_hat <- y1 * ((theta_old / 4) / denom)
    
    ## -------------------------
    ## M-step:
    ## theta_new = (z2_hat + y4) / (z2_hat + y4 + y2 + y3)
    ## -------------------------
    theta_new <- (z2_hat + y4) / (z2_hat + y4 + y2 + y3)
    
    diff <- abs(theta_new - theta_old)
    
    if (verbose) {
      cat(sprintf("%4d  % .8f  % .4f  % .4f\n",
                  iter, theta_new, z1_hat, z2_hat))
    }
    
    theta_old <- theta_new
  }
  
  if (diff > tol) {
    warning("Algoritma belum konvergen sampai iterasi maksimum.")
  }
  
  ## Pada saat konvergen, x1_hat = z1_hat, x2_hat = z2_hat
  list(
    theta_hat = theta_old,
    x1_hat    = z1_hat,
    x2_hat    = z2_hat,
    iter      = iter,
    diff      = diff
  )
}

## ---------------------------------------------------------
## Contoh pemanggilan
## ---------------------------------------------------------

hasil_em <- em_rao(theta_init = 0.5, tol = 1e-7, verbose = TRUE)

Iter  theta        z1_hat       z2_hat
   1   0.62025316   80.0000   20.0000
   2   0.63711798   76.3285   23.6715
   3   0.63924804   75.8404   24.1596
   4   0.63951337   75.7792   24.2208
   5   0.63954636   75.7715   24.2285
   6   0.63955047   75.7706   24.2294
   7   0.63955098   75.7705   24.2295
   8   0.63955104   75.7705   24.2295

hasil_em$theta_hat  # estimasi theta

[1] 0.639551

hasil_em$x1_hat     # estimasi x1

[1] 75.77046

hasil_em$x2_hat     # estimasi x2

[1] 24.22954

Muhammad Yusran
Program Studi Magister Statistika dan Sains Data
Sekolah Sains Data, Matematika, dan Informatika
IPB University

Email: muhammadyusran@apps.ipb.ac.id
LinkedIn: Muhammad Yusran
Instagram: @mhammad.yusran_