Penaksiran Parameter MLE dengan Metode Numerik Iteratif

1 Pendahuluan

Maximum Likelihood Estimation (MLE) merupakan metode penaksiran parameter yang bertujuan memaksimumkan fungsi log-likelihood \(\ell(\boldsymbol{\theta})\) terhadap parameter \(\boldsymbol{\theta}\). Secara umum, penaksiran MLE dilakukan dengan menyelesaikan persamaan score, yaitu:

\[ U(\boldsymbol{\theta}) = \frac{\partial \ell(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = \boldsymbol{0}. \]

Namun, persamaan tersebut sering kali tidak dapat diselesaikan secara eksplisit karena melibatkan fungsi khusus, model campuran, atau variabel laten. Oleh karena itu, diperlukan metode numerik iteratif untuk memperoleh taksiran parameter yang optimal. Tiga metode yang dibahas dalam materi ini adalah Newton-Raphson, BFGS, dan Expectation-Maximization (EM).

1.1 Metode Newton-Raphson

Newton-Raphson merupakan metode iteratif yang menggunakan turunan pertama (gradien atau score) dan turunan kedua (Hessian) dari fungsi log-likelihood. Pembaruan parameter dilakukan berdasarkan pendekatan ekspansi Taylor orde kedua, dengan rumus:

\[ \boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} - \left[H(\boldsymbol{\theta}^{(t)})\right]^{-1} U(\boldsymbol{\theta}^{(t)}). \]

dengan \(U(\boldsymbol{\theta})\) merupakan vektor score dan \(H(\boldsymbol{\theta})\) merupakan matriks Hessian. Metode ini memiliki konvergensi kuadratik apabila kondisi konvergensinya terpenuhi, tetapi memerlukan perhitungan Hessian pada setiap iterasi. Metode ini juga sensitif terhadap pemilihan nilai awal.

1.2 Metode BFGS

BFGS (Broyden-Fletcher-Goldfarb-Shanno) merupakan metode kuasi-Newton yang mengurangi kebutuhan menghitung dan menginvers matriks Hessian secara langsung pada setiap iterasi. Metode ini membangun aproksimasi Hessian atau invers Hessian secara iteratif menggunakan informasi gradien dan perubahan parameter.

Pada iterasi ke-\(t\), perubahan parameter dan gradien didefinisikan sebagai:

\[ \mathbf{s}^{(t)} = \boldsymbol{\theta}^{(t+1)} - \boldsymbol{\theta}^{(t)}, \]

\[ \mathbf{y}^{(t)} = U(\boldsymbol{\theta}^{(t+1)}) - U(\boldsymbol{\theta}^{(t)}). \]

Aproksimasi invers Hessian diperbarui menggunakan informasi tersebut. Metode BFGS umumnya lebih efisien ketika perhitungan Hessian mahal dan dapat mencapai konvergensi superlinear dalam kondisi yang sesuai.

1.3 Metode Expectation-Maximization (EM)

EM (Expectation-Maximization) merupakan algoritma iteratif yang digunakan untuk penaksiran parameter pada model yang melibatkan data hilang atau variabel laten, seperti model campuran. Algoritma ini terdiri atas dua tahap utama.

Tahap Expectation (E-step) menghitung ekspektasi log-likelihood data lengkap berdasarkan data teramati dan taksiran parameter saat ini, yang secara umum dinyatakan sebagai:

\[ Q(\boldsymbol{\theta}\mid\boldsymbol{\theta}^{(t)}) = E\left[ \ell_c(\boldsymbol{\theta}) \mid \mathbf{y},\boldsymbol{\theta}^{(t)} \right]. \]

Selanjutnya, tahap Maximization (M-step) memperbarui parameter dengan memaksimumkan fungsi tersebut:

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

Kedua tahap dilakukan secara bergantian hingga memenuhi kriteria konvergensi. Salah satu sifat utama EM adalah nilai likelihood data teramati tidak menurun pada setiap iterasi, meskipun algoritma dapat berakhir pada maksimum lokal.

Berdasarkan ketiga metode tersebut, latihan ini menerapkan Newton-Raphson untuk memperoleh MLE distribusi Beta, BFGS untuk penaksiran parameter regresi Poisson, serta EM untuk menaksir parameter model campuran dua distribusi Poisson. Penerapan ini bertujuan menunjukkan bagaimana metode numerik iteratif digunakan untuk menyelesaikan permasalahan penaksiran parameter yang tidak selalu memiliki solusi analitik sederhana.

2 Menurunkan Iterasi Newton-Raphson untuk MLE Distribusi Beta

2.1 Model dan fungsi likelihood

Misalkan \(X_1,\ldots,X_n\) merupakan sampel acak independen dari distribusi Beta dengan parameter \(\alpha>0\) dan \(\beta>0\). Fungsi kepadatan peluangnya adalah:

\[ f(x_i\mid\alpha,\beta) =\frac{x_i^{\alpha-1}(1-x_i)^{\beta-1}}{B(\alpha,\beta)}, \qquad 0<x_i<1, \]

dengan fungsi beta:

\[ B(\alpha,\beta)= \frac{\Gamma(\alpha)\Gamma(\beta)} {\Gamma(\alpha+\beta)}. \]

Likelihood sampel adalah:

\[ L(\alpha,\beta) =\prod_{i=1}^{n} \frac{x_i^{\alpha-1}(1-x_i)^{\beta-1}} {B(\alpha,\beta)}. \]

Dengan mengambil logaritma, diperoleh log-likelihood:

\[ \begin{aligned} \ell(\alpha,\beta) &=\log L(\alpha,\beta)\\ &=(\alpha-1)\sum_{i=1}^{n}\log x_i +(\beta-1)\sum_{i=1}^{n}\log(1-x_i)\\ &\quad-n\log\Gamma(\alpha)-n\log\Gamma(\beta) +n\log\Gamma(\alpha+\beta). \end{aligned} \]

MLE \((\widehat\alpha,\widehat\beta)\) diperoleh dengan memaksimumkan \(\ell(\alpha,\beta)\). Secara umum, persamaan score tidak dapat diselesaikan secara eksplisit karena melibatkan fungsi digamma.

2.2 Turunan pertama: Score

Gunakan identitas:

\[ \frac{d}{dx}\log\Gamma(x)=\psi(x), \]

dengan \(\psi(x)\) adalah fungsi digamma. Turunan terhadap \(\alpha\) adalah:

\[ \frac{\partial\ell}{\partial\alpha} =\sum_{i=1}^{n}\log x_i -n\psi(\alpha)+n\psi(\alpha+\beta). \]

Turunan terhadap \(\beta\) adalah:

\[ \frac{\partial\ell}{\partial\beta} =\sum_{i=1}^{n}\log(1-x_i) -n\psi(\beta)+n\psi(\alpha+\beta). \]

Dengan demikian, vektor score dapat ditulis:

\[ U(\boldsymbol\theta)= \begin{pmatrix} U_\alpha\\ U_\beta \end{pmatrix} = \begin{pmatrix} \sum_i\log x_i-n\psi(\alpha)+n\psi(\alpha+\beta)\\ \sum_i\log(1-x_i)-n\psi(\beta)+n\psi(\alpha+\beta) \end{pmatrix}, \]

dengan \(\boldsymbol\theta=(\alpha,\beta)^\mathsf T\). Persamaan score yang hendak diselesaikan adalah \(U(\boldsymbol\theta)=\boldsymbol 0\).

2.3 Turunan kedua: Hessian

Turunan fungsi digamma adalah fungsi trigamma, ditulis \(\psi_1(x)\). Elemen-elemen Hessian log-likelihood adalah:

\[ \frac{\partial^2\ell}{\partial\alpha^2} =-n\psi_1(\alpha)+n\psi_1(\alpha+\beta), \]

\[ \frac{\partial^2\ell}{\partial\beta^2} =-n\psi_1(\beta)+n\psi_1(\alpha+\beta), \]

dan turunan silang:

\[ \frac{\partial^2\ell}{\partial\alpha\,\partial\beta} =\frac{\partial^2\ell}{\partial\beta\,\partial\alpha} =n\psi_1(\alpha+\beta). \]

Jadi Hessian-nya adalah:

\[ H(\boldsymbol\theta)= n\begin{pmatrix} \psi_1(\alpha+\beta)-\psi_1(\alpha)& \psi_1(\alpha+\beta)\\ \psi_1(\alpha+\beta)& \psi_1(\alpha+\beta)-\psi_1(\beta) \end{pmatrix}. \]

2.4 Penurunan iterasi Newton-Raphson

Untuk memaksimumkan log-likelihood, pendekatan Taylor pada score di sekitar nilai saat ini \(\boldsymbol\theta^{(t)}\) memberikan:

\[ U(\boldsymbol\theta^{(t+1)}) \approx U(\boldsymbol\theta^{(t)}) +H(\boldsymbol\theta^{(t)}) (\boldsymbol\theta^{(t+1)}-\boldsymbol\theta^{(t)}). \]

Dengan menyamakan pendekatan tersebut dengan nol dan menyelesaikan terhadap parameter berikutnya, diperoleh:

\[ \boxed{ \boldsymbol\theta^{(t+1)} =\boldsymbol\theta^{(t)} -\left[H(\boldsymbol\theta^{(t)})\right]^{-1} U(\boldsymbol\theta^{(t)}) }. \]

Secara komputasi, lebih baik menyelesaikan sistem linear \(H(\boldsymbol\theta^{(t)})\boldsymbol d^{(t)} =-U(\boldsymbol\theta^{(t)})\), lalu menetapkan \(\boldsymbol\theta^{(t+1)}=\boldsymbol\theta^{(t)}+\boldsymbol d^{(t)}\), daripada menghitung invers matriks secara langsung.

Kriteria berhenti dapat menggunakan \(\|\boldsymbol\theta^{(t+1)}-\boldsymbol\theta^{(t)}\|<\varepsilon\) atau \(\|U(\boldsymbol\theta^{(t+1)})\|<\varepsilon\). Karena parameter Beta harus positif, nilai awal harus positif dan langkah iterasi perlu diperiksa agar tidak menghasilkan parameter nonpositif.

2.5 Implementasi R

Fungsi digamma() dan trigamma() tersedia di R. Kode berikut menjalankan Newton-Raphson dengan backtracking sederhana agar parameter tetap positif dan log-likelihood tidak menurun.

# Contoh Data Simulasi
set.seed(123)
x_beta <- rbeta(200, shape1 = 2.5, shape2 = 5)

beta_loglik <- function(par, x) {
  alpha <- par[1]
  beta  <- par[2]
  if (alpha <= 0 || beta <= 0) return(-Inf)
  n <- length(x)
  (alpha - 1) * sum(log(x)) +
    (beta - 1) * sum(log1p(-x)) -
    n * lbeta(alpha, beta)
}

beta_score <- function(par, x) {
  alpha <- par[1]
  beta  <- par[2]
  n <- length(x)
  c(
    sum(log(x)) - n * digamma(alpha) + n * digamma(alpha + beta),
    sum(log1p(-x)) - n * digamma(beta) + n * digamma(alpha + beta)
  )
}

beta_hessian <- function(par, x) {
  alpha <- par[1]
  beta  <- par[2]
  n <- length(x)
  common <- n * trigamma(alpha + beta)
  matrix(c(
    n * (trigamma(alpha + beta) - trigamma(alpha)), common,
    common, n * (trigamma(alpha + beta) - trigamma(beta))
  ), nrow = 2, byrow = TRUE)
}

newton_beta <- function(x, start = c(1, 1),
                        tol = 1e-8, maxit = 1000) {
  par <- start
  
  if (any(par <= 0))
    stop("Nilai awal alpha dan beta harus positif.")
  
  history <- data.frame(
    iter = 0,
    alpha = par[1],
    beta = par[2],
    logLik = beta_loglik(par, x),
    score_norm = sqrt(sum(beta_score(par, x)^2))
  )
  
  converged <- FALSE
  reason <- "Mencapai batas maksimum iterasi"
  
  for (it in 1:maxit) {
    U <- beta_score(par, x)
    H <- beta_hessian(par, x)
    
    # Arah Newton
    direction <- tryCatch(
      as.vector(solve(H, -U)),
      error = function(e) rep(NA_real_, length(par))
    )
    
    if (any(!is.finite(direction))) {
      stop("Sistem Newton gagal diselesaikan; coba nilai awal lain.")
    }
    
    # Backtracking
    step <- 1
    old_ll <- beta_loglik(par, x)
    
    repeat {
      candidate <- par + step * direction
      
      if (all(candidate > 0) &&
          beta_loglik(candidate, x) >= old_ll) break
      
      step <- step / 2
      
      if (step < 1e-10) {
        stop("Langkah terlalu kecil; coba nilai awal lain.")
      }
    }
    
    change <- max(abs(candidate - par))
    par <- candidate
    
    score_norm <- sqrt(sum(beta_score(par, x)^2))
    
    history <- rbind(history, data.frame(
      iter = it,
      alpha = par[1],
      beta = par[2],
      logLik = beta_loglik(par, x),
      score_norm = score_norm
    ))
    
    # Kriteria konvergensi
    if (change < tol || score_norm < tol) {
      converged <- TRUE
      
      if (score_norm < tol) {
        reason <- "Norma score memenuhi toleransi"
      } else {
        reason <- "Perubahan parameter memenuhi toleransi"
      }
      
      break
    }
  }
  
  list(
    estimate = setNames(par, c("alpha", "beta")),
    logLik = beta_loglik(par, x),
    iterations = nrow(history) - 1,
    converged = converged,
    reason = reason,
    grad_norm = sqrt(sum(beta_score(par, x)^2)),
    history = history
  )
}

# Hasil
hasil_beta <- newton_beta(
  x_beta,
  start = c(1, 1),
  tol = 1e-7,
  maxit = 1000
)

# Estimasi parameter
hasil_beta$estimate

##   alpha    beta 
## 2.80589 5.71917

# Nilai log-likelihood
hasil_beta$logLik

## [1] 99.6724

# Jumlah iterasi
hasil_beta$iterations

## [1] 7

# Status konvergensi
hasil_beta$converged

## [1] TRUE

# Keterangan konvergensi
hasil_beta$reason

## [1] "Norma score memenuhi toleransi"

# Norma gradien akhir (jika disimpan oleh fungsi)
hasil_beta$grad_norm

## [1] 9.68868e-09

# Riwayat iterasi
head(hasil_beta$history)

tail(hasil_beta$history)

2.6 Interpretasi

Berdasarkan hasil estimasi parameter distribusi Beta menggunakan metode Newton-Raphson, diperoleh estimasi parameter pertama sebesar \(\hat{\alpha}=2.80589\) dan parameter kedua sebesar \(\hat{\beta}=5.71917\). Nilai tersebut menunjukkan parameter distribusi Beta yang diperoleh dari data simulasi melalui proses maksimum likelihood. Nilai log-likelihood yang dihasilkan adalah \(99.6724\), dengan algoritma mencapai konvergensi pada iterasi ke-7.

Konvergensi algoritma ditunjukkan oleh norma skor akhir sebesar \(9.68868\times10^{-9}\), yang telah memenuhi toleransi sebesar \(10^{-7}\). Selama iterasi, nilai log-likelihood meningkat dan norma skor menurun hingga mendekati nol. Hal ini menunjukkan bahwa proses iterasi berhasil menemukan estimasi parameter yang memenuhi kriteria konvergensi metode Newton-Raphson.

3 Mengimplementasikan BFGS untuk Regresi Poisson

3.1 Model regresi Poisson

Misalkan \(Y_i\) adalah respons cacah dan \(x_i\) adalah vektor prediktor. Model regresi Poisson menggunakan:

\[ Y_i\mid x_i\sim\operatorname{Poisson}(\mu_i), \qquad \log(\mu_i)=x_i^\mathsf T\boldsymbol\beta, \]

sehingga \(\mu_i=\exp(x_i^\mathsf T\boldsymbol\beta)\). Untuk \(n\) pengamatan independen, likelihood-nya adalah:

\[ L(\boldsymbol\beta)= \prod_{i=1}^{n}\frac{\exp(-\mu_i)\mu_i^{y_i}}{y_i!}. \]

Log-likelihood adalah:

\[ \ell(\boldsymbol\beta) =\sum_{i=1}^{n} \left[y_i x_i^\mathsf T\boldsymbol\beta -\exp(x_i^\mathsf T\boldsymbol\beta)-\log(y_i!)\right]. \]

Maksimisasi \(\ell(\boldsymbol\beta)\) ekuivalen dengan minimisasi negatif log-likelihood \(f(\boldsymbol\beta)=-\ell(\boldsymbol\beta)\).

3.2 Gradien dan Hessian

Tuliskan \(X\) sebagai matriks desain, \(\boldsymbol y=(y_1,\ldots,y_n)^\mathsf T\), dan \(\boldsymbol\mu=(\mu_1,\ldots,\mu_n)^\mathsf T\). Karena \(\partial\mu_i/\partial\boldsymbol\beta=\mu_i x_i\), gradien log-likelihood adalah:

\[ \nabla\ell(\boldsymbol\beta) =X^\mathsf T(\boldsymbol y-\boldsymbol\mu). \]

Dengan demikian, gradien negatif log-likelihood adalah:

\[ \boxed{\nabla f(\boldsymbol\beta) =X^\mathsf T(\boldsymbol\mu-\boldsymbol y)}. \]

Hessian log-likelihood adalah:

\[ \nabla^2\ell(\boldsymbol\beta) =-X^\mathsf T W X, \qquad W=\operatorname{diag}(\mu_1,\ldots,\mu_n), \]

sehingga Hessian negatif log-likelihood adalah \(\nabla^2 f(\boldsymbol\beta)=X^\mathsf T W X\). BFGS tidak perlu menghitung Hessian ini pada setiap iterasi; metode tersebut membangun aproksimasi invers Hessian dari perubahan gradien.

3.3 Pembaruan BFGS

Misalkan \(f\) adalah fungsi yang diminimalkan dan \(g_t=\nabla f(\boldsymbol\beta_t)\). Arah pencarian ditentukan oleh:

\[ \boldsymbol d_t=-M_tg_t, \]

dengan \(M_t\) sebagai aproksimasi invers Hessian. Dengan ukuran langkah \(\eta_t>0\),

\[ \boldsymbol\beta_{t+1} =\boldsymbol\beta_t+\eta_t\boldsymbol d_t. \]

Definisikan:

\[ s_t=\boldsymbol\beta_{t+1}-\boldsymbol\beta_t, \qquad y_t=g_{t+1}-g_t, \qquad \rho_t=\frac{1}{y_t^\mathsf T s_t}. \]

Pembaruan aproksimasi invers Hessian adalah:

\[ \boxed{ M_{t+1}= (I-\rho_t s_ty_t^\mathsf T)M_t (I-\rho_t y_ts_t^\mathsf T) +\rho_t s_ts_t^\mathsf T }. \]

Pembaruan digunakan ketika \(y_t^\mathsf T s_t>0\), yang membantu menjaga sifat positif-definit aproksimasi invers Hessian. Ukuran langkah dipilih melalui backtracking line search agar fungsi objektif menurun. Iterasi berhenti jika norma gradien atau perubahan parameter cukup kecil.

3.4 Implementasi R

Kode berikut mengimplementasikan pembaruan invers Hessian BFGS dan backtracking line search tanpa memanggil optim(method = "BFGS"). Untuk kestabilan numerik, fungsi objektif yang diminimalkan adalah negatif log-likelihood.

# Contoh Data Simulasi
set.seed(456)
n <- 250
x1 <- rnorm(n)
x2 <- rbinom(n, size = 1, prob = 0.4)
X <- model.matrix(~ x1 + x2)
beta_true <- c(0.3, 0.5, -0.4)
mu_true <- exp(as.vector(X %*% beta_true))
y <- rpois(n, lambda = mu_true)

# Negatif log-likelihood Poisson
poisson_nll <- function(beta, X, y) {
  eta <- as.vector(X %*% beta)
  mu <- exp(pmin(eta, 700))  # mengurangi risiko overflow
  -sum(dpois(y, lambda = mu, log = TRUE))
}

# Gradien negatif log-likelihood
poisson_grad <- function(beta, X, y) {
  eta <- as.vector(X %*% beta)
  mu <- exp(pmin(eta, 700))
  as.vector(t(X) %*% (mu - y))
}


bfgs_poisson <- function(X, y, start = rep(0, ncol(X)),
                         tol = 1e-7, maxit = 1000) {
  p <- ncol(X)
  beta <- as.vector(start)
  M <- diag(p)
  
  f <- poisson_nll(beta, X, y)
  g <- poisson_grad(beta, X, y)
  grad_norm <- sqrt(sum(g^2))
  
  history <- data.frame(
    iter = 0,
    nll = f,
    grad_norm = grad_norm
  )
  
  converged <- grad_norm < tol
  reason <- if (converged) "Konvergen pada nilai awal" else ""
  
  it <- 0
  
  # Iterasi berjalan selama belum konvergen
  while (!converged && it < maxit) {
    
    # Arah pencarian BFGS
    direction <- -as.vector(M %*% g)
    
    # Jika bukan arah turun, gunakan steepest descent
    if (!all(is.finite(direction)) ||
        sum(g * direction) >= 0) {
      M <- diag(p)
      direction <- -g
    }
    
    slope <- sum(g * direction)
    step <- 1
    
    # Backtracking line search
    repeat {
      beta_new <- beta + step * direction
      f_new <- poisson_nll(beta_new, X, y)
      
      if (is.finite(f_new) &&
          f_new <= f + 1e-4 * step * slope) {
        break
      }
      
      step <- step / 2
      
      if (step < 1e-12) break
    }
    
    # Jika line search gagal, reset arah dan coba lagi
    if (step < 1e-12) {
      M <- diag(p)
      direction <- -g
      slope <- sum(g * direction)
      step <- 1
      
      repeat {
        beta_new <- beta + step * direction
        f_new <- poisson_nll(beta_new, X, y)
        
        if (is.finite(f_new) &&
            f_new <= f + 1e-4 * step * slope) {
          break
        }
        
        step <- step / 2
        
        if (step < 1e-12) break
      }
    }
    
    # Tidak dapat melanjutkan dengan langkah yang valid
    if (step < 1e-12) {
      reason <- "Line search gagal; iterasi tidak dapat dilanjutkan"
      break
    }
    
    # Gradien pada parameter baru
    g_new <- poisson_grad(beta_new, X, y)
    
    # Perubahan parameter dan gradien
    s <- beta_new - beta
    yy <- g_new - g
    ys <- sum(yy * s)
    
    # Pembaruan aproksimasi invers Hessian
    if (is.finite(ys) && ys > 1e-12) {
      rho <- 1 / ys
      I <- diag(p)
      A <- I - rho * outer(s, yy)
      M <- A %*% M %*% t(A) + rho * outer(s, s)
    } else {
      M <- diag(p)
    }
    
    # Memperbarui parameter dan fungsi objektif
    beta <- beta_new
    f <- f_new
    g <- g_new
    grad_norm <- sqrt(sum(g^2))
    it <- it + 1
    
    history <- rbind(history, data.frame(
      iter = it,
      nll = f,
      grad_norm = grad_norm
    ))
    
    # Periksa konvergensi setelah pembaruan
    converged <- is.finite(grad_norm) &&
      grad_norm < tol
    
    if (converged) {
      reason <- "Norma gradien memenuhi toleransi"
    }
  }
  
  if (!converged && it >= maxit) {
    reason <- "Mencapai iterasi maksimum sebelum konvergen"
  }
  
  list(
    coefficients = setNames(beta, colnames(X)),
    logLik = -f,
    iterations = it,
    converged = converged,
    reason = reason,
    history = history,
    grad_norm = grad_norm
  )
}

# Hasil
hasil_bfgs <- bfgs_poisson(
  X, y,
  tol = 1e-7,
  maxit = 1000
)

# Estimasi parameter
hasil_bfgs$coefficients

## (Intercept)          x1          x2 
##    0.261051    0.541544   -0.280723

# Nilai log-likelihood
hasil_bfgs$logLik

## [1] -346.25

# Jumlah iterasi
hasil_bfgs$iterations

## [1] 11

# Status konvergensi
hasil_bfgs$converged

## [1] TRUE

# Keterangan konvergensi
hasil_bfgs$reason

## [1] "Norma gradien memenuhi toleransi"

# Norma gradien akhir
hasil_bfgs$grad_norm

## [1] 3.39367e-10

# Riwayat iterasi
head(hasil_bfgs$history)

tail(hasil_bfgs$history)

Pemeriksaan dengan fungsi bawaan R

Sebagai pemeriksaan tambahan, hasil implementasi manual dapat dibandingkan dengan hasil glm() yang menggunakan model linear tergeneralisasi dengan keluarga Poisson. Perbandingan ini bukan bagian dari algoritma BFGS manual, tetapi berguna untuk memeriksa apakah estimasi koefisien berada pada nilai yang serupa.

fit_glm <- glm(y ~ x1 + x2, family = poisson(link = "log"))
cbind(
  BFGS_manual = hasil_bfgs$coefficients,
  GLM_R = coef(fit_glm)
)

##             BFGS_manual     GLM_R
## (Intercept)    0.261051  0.261051
## x1             0.541544  0.541544
## x2            -0.280723 -0.280723

3.5 Interpretasi

Berdasarkan hasil estimasi regresi Poisson menggunakan metode BFGS, diperoleh estimasi koefisien intersep sebesar \(\hat{\beta}_0=0.261051\), koefisien variabel prediktor pertama sebesar \(\hat{\beta}_1=0.541544\), dan koefisien variabel prediktor kedua sebesar \(\hat{\beta}2=-0.280723\). Dengan demikian, model regresi Poisson yang diperoleh dapat dituliskan sebagai \(\log(\hat{\mu}i)=0.261051+0.541544x{1i}-0.280723x{2i}\), dengan \(\hat{\mu}_i\) menyatakan estimasi nilai harapan variabel respons untuk pengamatan ke-\(i\).

Koefisien variabel prediktor pertama bernilai positif, sehingga kenaikan satu unit pada variabel tersebut berkaitan dengan kenaikan nilai harapan respons, dengan variabel lain tetap. Sebaliknya, koefisien variabel prediktor kedua bernilai negatif, sehingga nilai satu pada variabel tersebut berkaitan dengan nilai harapan respons yang lebih rendah dibandingkan nilai nol, dengan variabel lain tetap. Model konvergen pada iterasi ke-11 dengan nilai log-likelihood sebesar \(-346.25\) dan norma gradien akhir sebesar \(3.39367\times10^{-10}\). Koefisien yang dihasilkan juga sama dengan koefisien dari fungsi glm() dengan distribusi Poisson dan fungsi penghubung log, sehingga hasil implementasi BFGS sesuai dengan hasil pembanding tersebut.

4 Menerapkan EM untuk Model Campuran Dua Poisson

4.1 Definisi model campuran

Misalkan data cacah \(Y_1,\ldots,Y_n\) berasal dari campuran dua distribusi Poisson. Setiap pengamatan berasal dari komponen pertama dengan probabilitas \(\pi\), dan dari komponen kedua dengan probabilitas \(1-\pi\). Modelnya adalah:

\[ P(Y_i=y_i) =\pi\frac{e^{-\lambda_1}\lambda_1^{y_i}}{y_i!} +(1-\pi)\frac{e^{-\lambda_2}\lambda_2^{y_i}}{y_i!}, \]

dengan \(0<\pi<1\), \(\lambda_1>0\), dan \(\lambda_2>0\).

Misalkan \(Z_i\in\{0,1\}\) adalah indikator laten, dengan \(Z_i=1\) jika pengamatan berasal dari komponen pertama dan \(Z_i=0\) jika berasal dari komponen kedua. Maka:

\[ P(Z_i=1)=\pi,\qquad P(Z_i=0)=1-\pi, \]

dan,

\[ Y_i\mid Z_i=1\sim\operatorname{Poisson}(\lambda_1), \qquad Y_i\mid Z_i=0\sim\operatorname{Poisson}(\lambda_2). \]

Log-likelihood data teramati adalah:

\[ \ell_{\mathrm{obs}}(\pi,\lambda_1,\lambda_2) =\sum_{i=1}^{n} \log\left[ \pi f_1(y_i)+(1-\pi)f_2(y_i) \right], \]

dengan \(f_k(y_i)=P(Y_i=y_i\mid Z_i=k)\). Karena terdapat penjumlahan di dalam logaritma, optimasi langsung tidak menghasilkan pembaruan parameter sederhana. EM mengatasi hal ini dengan memperlakukan \(Z_i\) sebagai data laten.

4.2 Complete-data log-likelihood

Jika \(Z_i\) diketahui, log-likelihood lengkap adalah:

\[ \begin{aligned} \ell_c =&\sum_{i=1}^{n}\Big[ Z_i\{\log\pi+y_i\log\lambda_1-\lambda_1-\log(y_i!)\}\\ &+(1-Z_i)\{\log(1-\pi)+y_i\log\lambda_2 -\lambda_2-\log(y_i!)\}\Big]. \end{aligned} \]

EM mengganti \(Z_i\) yang tidak teramati dengan ekspektasi posteriornya pada E-step, kemudian memaksimumkan ekspektasi log-likelihood lengkap pada M-step.

4.3 E-step: menghitung responsibility

Pada iterasi ke-\(t\), hitung probabilitas posterior bahwa pengamatan ke-\(i\) berasal dari komponen pertama:

\[ \boxed{ \tau_i^{(t)} =E(Z_i\mid y_i,\boldsymbol\theta^{(t)}) =\frac{\pi^{(t)}f_1(y_i;\lambda_1^{(t)})} {\pi^{(t)}f_1(y_i;\lambda_1^{(t)}) +(1-\pi^{(t)})f_2(y_i;\lambda_2^{(t)})} }. \]

Nilai \(\tau_i^{(t)}\) disebut responsibility. Posterior untuk komponen kedua adalah \(1-\tau_i^{(t)}\). Dalam implementasi, probabilitas Poisson dapat dihitung pada skala log untuk mengurangi risiko underflow.

4.4 M-step: memperbarui parameter

Ekspektasi log-likelihood lengkap, dengan mengabaikan suku yang tidak bergantung pada parameter, adalah:

\[ \begin{aligned} Q(\boldsymbol\theta\mid\boldsymbol\theta^{(t)}) =&\sum_{i=1}^{n}\tau_i^{(t)}\log\pi +\sum_{i=1}^{n}(1-\tau_i^{(t)})\log(1-\pi)\\ &+\sum_{i=1}^{n}\tau_i^{(t)} (y_i\log\lambda_1-\lambda_1)\\ &+\sum_{i=1}^{n}(1-\tau_i^{(t)}) (y_i\log\lambda_2-\lambda_2). \end{aligned} \]

Turunan terhadap \(\pi\) adalah:

\[ \frac{\partial Q}{\partial\pi} =\frac{\sum_i\tau_i^{(t)}}{\pi} -\frac{\sum_i(1-\tau_i^{(t)})}{1-\pi}. \]

Menyamakan turunan dengan nol menghasilkan:

\[ \boxed{\pi^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n}\tau_i^{(t)}}. \]

Turunan terhadap \(\lambda_1\) adalah \(\sum_i\tau_i^{(t)}y_i/\lambda_1-\sum_i\tau_i^{(t)}\). Dengan menyamakan turunan dengan nol, diperoleh:

\[ \boxed{ \lambda_1^{(t+1)} =\frac{\sum_i\tau_i^{(t)}y_i}{\sum_i\tau_i^{(t)}} }. \]

Dengan cara yang sama,

\[ \boxed{ \lambda_2^{(t+1)} =\frac{\sum_i(1-\tau_i^{(t)})y_i} {\sum_i(1-\tau_i^{(t)})} }. \]

Jadi, pada setiap iterasi, E-step menghitung responsibility dengan parameter saat ini, lalu M-step menghitung rata-rata tertimbang untuk memperbarui bobot campuran dan kedua parameter Poisson.

4.5 Kriteria konvergensi

Log-likelihood data teramati dihitung dengan:

\[ \ell_{\mathrm{obs}}^{(t)} =\sum_{i=1}^{n}\log\left[ \pi^{(t)}f_1(y_i;\lambda_1^{(t)}) +(1-\pi^{(t)})f_2(y_i;\lambda_2^{(t)}) \right]. \]

Sesuai sifat EM, likelihood teramati tidak menurun pada setiap iterasi ideal. Algoritma dapat dihentikan ketika:

\[ \left|\ell_{\mathrm{obs}}^{(t+1)} -\ell_{\mathrm{obs}}^{(t)}\right|<\varepsilon, \]

atau ketika perubahan parameter cukup kecil. EM dapat konvergen ke maksimum lokal, sehingga beberapa nilai awal dapat dicoba.

4.6 Implementasi R

# Fungsi Log-likelihood Teramati untuk Campuran Dua Poisson
mix2_loglik <- function(y, pi, lambda1, lambda2) {
  if (pi <= 0 || pi >= 1 || lambda1 <= 0 || lambda2 <= 0) {
    return(-Inf)
  }
  log_f1 <- log(pi) + dpois(y, lambda1, log = TRUE)
  log_f2 <- log1p(-pi) + dpois(y, lambda2, log = TRUE)

  # log(exp(a) + exp(b)) yang lebih stabil secara numerik
  m <- pmax(log_f1, log_f2)
  sum(m + log(exp(log_f1 - m) + exp(log_f2 - m)))
}

em_poisson2 <- function(y, pi = 0.5, lambda1 = 2, lambda2 = 8,
                        tol = 1e-8, maxit = 1000) {
  if (any(y < 0) || any(y != floor(y))) {
    stop("Semua pengamatan harus berupa bilangan cacah nonnegatif.")
  }
  
  if (pi <= 0 || pi >= 1 || lambda1 <= 0 || lambda2 <= 0) {
    stop("Gunakan 0 < pi < 1 dan lambda1, lambda2 > 0.")
  }
  
  # Menata label awal
  if (lambda1 > lambda2) {
    temp <- lambda1
    lambda1 <- lambda2
    lambda2 <- temp
    pi <- 1 - pi
  }
  
  ll <- mix2_loglik(y, pi, lambda1, lambda2)
  
  history <- data.frame(
    iter = 0,
    pi = pi,
    lambda1 = lambda1,
    lambda2 = lambda2,
    logLik = ll
  )
  
  converged <- FALSE
  reason <- "Mencapai batas maksimum iterasi"
  
  for (it in 1:maxit) {
    
    # E-step
    log_num1 <- log(pi) + dpois(y, lambda1, log = TRUE)
    log_num2 <- log1p(-pi) + dpois(y, lambda2, log = TRUE)
    
    m <- pmax(log_num1, log_num2)
    den <- m + log(exp(log_num1 - m) +
                     exp(log_num2 - m))
    
    tau <- exp(log_num1 - den)
    
    # M-step
    sum_tau <- sum(tau)
    sum_one_minus_tau <- sum(1 - tau)
    
    if (sum_tau <= 0 || sum_one_minus_tau <= 0) {
      stop("Satu komponen kehilangan seluruh bobot; coba nilai awal lain.")
    }
    
    pi_new <- mean(tau)
    lambda1_new <- sum(tau * y) / sum_tau
    lambda2_new <- sum((1 - tau) * y) / sum_one_minus_tau
    
    lambda1_new <- max(lambda1_new, 1e-10)
    lambda2_new <- max(lambda2_new, 1e-10)
    
    # Menata label komponen
    if (lambda1_new > lambda2_new) {
      temp <- lambda1_new
      lambda1_new <- lambda2_new
      lambda2_new <- temp
      pi_new <- 1 - pi_new
    }
    
    ll_new <- mix2_loglik(
      y, pi_new, lambda1_new, lambda2_new
    )
    
    history <- rbind(history, data.frame(
      iter = it,
      pi = pi_new,
      lambda1 = lambda1_new,
      lambda2 = lambda2_new,
      logLik = ll_new
    ))
    
    change <- max(abs(c(
      pi_new - pi,
      lambda1_new - lambda1,
      lambda2_new - lambda2
    )))
    
    ll_change <- abs(ll_new - ll)
    
    pi <- pi_new
    lambda1 <- lambda1_new
    lambda2 <- lambda2_new
    ll <- ll_new
    
    # Kriteria konvergensi
    if (change < tol || ll_change < tol) {
      converged <- TRUE
      
      if (change < tol && ll_change < tol) {
        reason <- "Perubahan parameter dan log-likelihood memenuhi toleransi"
      } else if (change < tol) {
        reason <- "Perubahan parameter memenuhi toleransi"
      } else {
        reason <- "Perubahan log-likelihood memenuhi toleransi"
      }
      
      break
    }
  }
  
  list(
    estimate = c(
      pi = pi,
      lambda1 = lambda1,
      lambda2 = lambda2
    ),
    logLik = ll,
    iterations = nrow(history) - 1,
    converged = converged,
    reason = reason,
    grad_norm = NA_real_,
    history = history
  )
}

# Contoh data simulasi dari campuran dua Poisson
set.seed(789)
n_mix <- 500
z <- rbinom(n_mix, size = 1, prob = 0.4)
y_mix <- ifelse(z == 1,
                rpois(n_mix, lambda = 2),
                rpois(n_mix, lambda = 8))

# Hasil
hasil_em <- em_poisson2(
  y_mix,
  pi = 0.5,
  lambda1 = 1.5,
  lambda2 = 7
)

# Estimasi parameter
hasil_em <- em_poisson2(
  y_mix,
  pi = 0.5,
  lambda1 = 1.5,
  lambda2 = 7,
  tol = 1e-7,
  maxit = 1000
)

# Estimasi parameter
hasil_em$estimate

##       pi  lambda1  lambda2 
## 0.381047 1.836703 7.619526

# Nilai log-likelihood
hasil_em$logLik

## [1] -1315.27

# Jumlah iterasi
hasil_em$iterations

## [1] 29

# Status konvergensi
hasil_em$converged

## [1] TRUE

# Keterangan konvergensi
hasil_em$reason

## [1] "Perubahan log-likelihood memenuhi toleransi"

# Riwayat iterasi
head(hasil_em$history)

tail(hasil_em$history)

Memeriksa kenaikan log-likelihood dan visualisasi

# Selisih log-likelihood antar iterasi
diff(hasil_em$history$logLik)

##  [1] 2.80186e+01 7.18046e-02 3.59529e-02 2.16578e-02 1.31390e-02 7.98434e-03
##  [7] 4.85791e-03 2.95855e-03 1.80316e-03 1.09962e-03 6.70894e-04 4.09467e-04
## [13] 2.49980e-04 1.52647e-04 9.32273e-05 5.69452e-05 3.47869e-05 2.12525e-05
## [19] 1.29847e-05 7.93371e-06 4.84771e-06 2.96218e-06 1.81007e-06 1.10608e-06
## [25] 6.75905e-07 4.13038e-07 2.52404e-07 1.54244e-07 9.42584e-08

# Plot Konvergensi EM untuk Campuran Dua Poisson
plot(hasil_em$history$iter,
     hasil_em$history$logLik,
     type = "o",
     col = "#7B2CBF",
     lwd = 2.5,
     pch = 16,
     cex = 0.7,
     xlab = "Iterasi",
     ylab = "Log-likelihood Teramati",
     main = "Konvergensi EM untuk Campuran Dua Poisson",
     col.main = "#343A40",
     col.lab = "#343A40",
     col.axis = "#495057",
     cex.main = 1.2,
     cex.lab = 1.05,
     las = 1)

# Garis bantu horizontal pada nilai akhir
abline(h = tail(hasil_em$history$logLik, 1),
       col = "#E76F51",
       lty = 2,
       lwd = 1.5)

# Anotasi log-likelihood akhir
text(x = max(hasil_em$history$iter) * 0.55,
     y = tail(hasil_em$history$logLik, 1) + 1,
     labels = paste0("Log-likelihood akhir = ",
                     round(tail(hasil_em$history$logLik, 1), 2)),
     col = "#E76F51",
     pos = 3,
     cex = 0.85)

# Grid tipis agar plot mudah dibaca
grid(col = "gray85", lty = 3)

4.7 Interpretasi

Berdasarkan hasil estimasi model campuran dua distribusi Poisson menggunakan algoritma Expectation-Maximization (EM), diperoleh estimasi proporsi komponen pertama sebesar \(\hat{\pi}=0.381047\), parameter rata-rata komponen pertama sebesar \(\hat{\lambda}_1=1.836703\), dan parameter rata-rata komponen kedua sebesar \(\hat{\lambda}_2=7.619526\). Estimasi \(\hat{\pi}=0.381047\) menunjukkan bahwa proporsi campuran yang diestimasi untuk komponen pertama sekitar 38.10%, sedangkan proporsi komponen kedua adalah \(1-\hat{\pi}=0.618953\), atau sekitar 61.90%.

Hasil estimasi menunjukkan bahwa komponen pertama memiliki rata-rata yang lebih rendah daripada komponen kedua, yaitu \(\hat{\lambda}_1=1.836703\) dan \(\hat{\lambda}_2=7.619526\). Algoritma EM mencapai konvergensi pada iterasi ke-29 dengan nilai log-likelihood sebesar \(-1315.27\). Selama proses iterasi, nilai log-likelihood meningkat hingga perubahan antara iterasi berturut-turut memenuhi toleransi sebesar \(10^{-7}\). Dengan demikian, kriteria konvergensi yang digunakan telah terpenuhi dan parameter campuran dua distribusi Poisson berhasil diestimasi dari data simulasi.

Berdasarkan plot konvergensi algoritma Expectation-Maximization (EM) pada model campuran dua distribusi Poisson, terlihat bahwa nilai log-likelihood meningkat dari sekitar \(-1343.45\) pada iterasi awal menjadi \(-1315.27\) pada iterasi ke-29. Peningkatan ini menunjukkan bahwa algoritma EM secara bertahap menemukan estimasi parameter yang memberikan kecocokan model yang lebih baik terhadap data. Pada iterasi-iterasi awal, peningkatan log-likelihood terlihat cukup besar, kemudian perubahan nilainya semakin kecil seiring bertambahnya iterasi.

Pada iterasi selanjutnya, kurva log-likelihood cenderung mendatar karena perubahan nilai log-likelihood ant ariterasi semakin kecil hingga memenuhi toleransi sebesar \(10^{-7}\). Kondisi tersebut menunjukkan bahwa algoritma EM telah mencapai konvergensi pada iterasi ke-29. Dengan demikian, plot mendukung hasil perhitungan numerik bahwa proses estimasi telah mencapai kondisi konvergen berdasarkan kriteria yang digunakan.

5 Summary

1. Newton–Raphson untuk Distribusi Beta

   Metode Newton–Raphson digunakan untuk mengestimasi parameter distribusi Beta melalui pemaksimuman log-likelihood. Prosesnya meliputi penurunan log-likelihood, pembentukan vektor score menggunakan fungsi digamma dan matriks Hessian menggunakan fungsi trigamma, kemudian pembaruan parameter dengan rumus:

   \[ \boldsymbol{\theta}^{(t+1)} = \boldsymbol{\theta}^{(t)} - [H(\boldsymbol{\theta}^{(t)})]^{-1} U(\boldsymbol{\theta}^{(t)}). \]

   Berdasarkan simulasi, diperoleh \(\hat{\alpha}=2.80589\) dan \(\hat{\beta}=5.71917\). Algoritma konvergen pada iterasi ke-7 dengan log-likelihood sebesar \(99.6724\) dan norma skor akhir sebesar \(9.68868\times10^{-9}\). Hasil ini menunjukkan bahwa estimasi parameter telah memenuhi kriteria konvergensi yang ditetapkan.

2. BFGS untuk Regresi Poisson

   Metode BFGS diterapkan untuk mengestimasi koefisien regresi Poisson dengan meminimumkan fungsi negatif log-likelihood. Gradien fungsi objektif dirumuskan sebagai:

   \[ \nabla f(\boldsymbol{\beta}) = X^\mathsf{T}(\boldsymbol{\mu}-\boldsymbol{y}), \]

   dengan \(\boldsymbol{\mu}=\exp(X\boldsymbol{\beta})\). Selanjutnya, aproksimasi invers Hessian diperbarui menggunakan formula BFGS dan line search untuk menentukan langkah iterasi. Hasil simulasi menghasilkan \(\hat{\beta}_0=0.261051\), \(\hat{\beta}_1=0.541544\), dan \(\hat{\beta}_2=-0.280723\). Algoritma konvergen pada iterasi ke-11 dengan log-likelihood sebesar \(-346.25\). Koefisien yang diperoleh sama dengan hasil fungsi glm() pada R, sehingga implementasi BFGS memberikan hasil yang konsisten dengan metode pembanding tersebut.

3. EM untuk Campuran Dua Distribusi Poisson

   Algoritma Expectation-Maximization (EM) digunakan untuk mengestimasi parameter model campuran dua distribusi Poisson ketika keanggotaan setiap observasi pada komponen campuran tidak diketahui. Pada tahap E-step, dihitung probabilitas tanggung jawab (responsibility) observasi ke-\(i\) terhadap komponen pertama:

   \[ \tau_i = \frac{\pi f(y_i;\lambda_1)} {\pi f(y_i;\lambda_1)+(1-\pi)f(y_i;\lambda_2)}. \]

   Pada tahap M-step, parameter diperbarui menggunakan rata-rata tertimbang:

   \[ \hat{\pi} = \frac{1}{n}\sum_{i=1}^{n}\tau_i, \]

   \[ \hat{\lambda}_1 = \frac{\sum_{i=1}^{n}\tau_i y_i} {\sum_{i=1}^{n}\tau_i}, \qquad \hat{\lambda}_2 = \frac{\sum_{i=1}^{n}(1-\tau_i)y_i} {\sum_{i=1}^{n}(1-\tau_i)}. \]

   Berdasarkan simulasi, diperoleh \(\hat{\pi}=0.381047\), \(\hat{\lambda}_1=1.836703\), dan \(\hat{\lambda}_2=7.619526\). Algoritma konvergen pada iterasi ke-29 dengan log-likelihood sebesar \(-1315.27\). Plot konvergensi menunjukkan bahwa log-likelihood meningkat dari sekitar \(-1343.45\) menjadi \(-1315.27\), kemudian kurva semakin mendatar hingga perubahan log-likelihood memenuhi toleransi \(10^{-7}\). Hal ini menunjukkan bahwa algoritma EM telah mencapai konvergensi berdasarkan kriteria yang digunakan.