Algoritma Newton-Raphson

Membaca data

data <- read.csv("D:/dataset.csv")
Y <- data$satellites
Y

##   [1]  8  0  9  0  4  0  0  0  0  0  0  0 11  0 14  8  1  1  0  5  4  3  1  2  3
##  [26]  0  3  5  0  0  4  0  0  8  5  0  0  6  0  6  3  5  6  5  9  4  6  4  3  3
##  [51]  5  5  6  4  5 15  3  3  0  0  0  5  3  5  1  8 10  0  0  3  7  1  0  6  0
##  [76]  0  3  4  0  5  0  0  0  4  0  3  0  0  0  0  5  0  0  0  0  1  0  1  1  1
## [101]  1  1  1  4  1  1  1  1  2  4  3  6  0  2  2  0 12  0  5  6  6  2  0  2  3
## [126]  0  3  4  2  6  6  0  4 10  7  0  5  5  6  6  7  3  3  0  0  8  4  4 10  9
## [151]  4  0  0  0  0  4  0  2  0  4  4  3  8  0  7  0  0  2  3  4  0  0  0

Fungsi log-likelihood untuk distribusi Poisson

log_likelihood <- function(lambda, Y) {
  n <- length(Y)
  sum(Y * log(lambda) - lambda - lfactorial(Y))  # Menggunakan lfactorial untuk stabilitas numerik
}

Gradient (turunan pertama dari log-likelihood)

gradient <- function(lambda, Y) {
  sum(Y) / lambda - length(Y)  # dL/d(lambda)
}

Hessian (turunan kedua dari log-likelihood)

hessian <- function(lambda, Y) {
  -sum(Y) / lambda^2  # d²L/d(lambda²)
}

Algoritma Newton-Raphson untuk estimasi lambda

newton_raphson <- function(Y, lambda_init = 1, tol = 1e-6, max_iter = 100) {
  lambda <- lambda_init
  iter_results <- data.frame(Iteration = integer(0), Lambda = numeric(0), Gradient = numeric(0), Hessian = numeric(0))
  
  for (k in 1:max_iter) {
    grad <- gradient(lambda, Y)
    hess <- hessian(lambda, Y)
    if (hess == 0) {
      cat("Hessian nol pada iterasi ke-", k, ". Algoritma dihentikan.\n")
      break
    }
    
    lambda_new <- lambda - grad / hess
    
    iter_results <- rbind(iter_results, data.frame(Iteration = k, Lambda = lambda_new, Gradient = grad, Hessian = hess))
    
    if (abs(lambda_new - lambda) < tol) {
      cat("Konvergen pada iterasi ke-", k, "\n")
      break
    }
    
    lambda <- lambda_new
  }
  
  return(list(estimated_lambda = lambda_new, iter_results = iter_results))
}

Jalankan algoritma Newton-Raphson

result <- newton_raphson(Y)

## Konvergen pada iterasi ke- 7

Tampilkan estimasi parameter lambda

cat("Estimasi parameter lambda:", result$estimated_lambda, "\n")

## Estimasi parameter lambda: 2.919075

Tampilkan hasil semua iterasi

cat("\nHasil Iterasi:\n")

## 
## Hasil Iterasi:

print(result$iter_results)

##   Iteration   Lambda     Gradient    Hessian
## 1         1 1.657426 3.320000e+02 -505.00000
## 2         2 2.373779 1.316894e+02 -183.83289
## 3         3 2.817212 3.974091e+01  -89.62118
## 4         4 2.915521 6.255263e+00  -63.62861
## 5         5 2.919071 2.109223e-01  -59.40995
## 6         6 2.919075 2.565319e-04  -59.26552
## 7         7 2.919075 3.803677e-10  -59.26535