Algoritma Fisher Scoring

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 - log(factorial(Y)))  # Log-likelihood Poisson
}

Gradient (turunan pertama)

gradient <- function(lambda, Y) {
  sum(Y) / lambda - length(Y)  # Turunan pertama log-likelihood
}

Fisher Information

fisher_information <- function(lambda, Y) {
  length(Y) / lambda  # Informasi Fisher untuk distribusi Poisson
}

Algoritma Fisher Scoring dengan output semua hasil iterasi

fisher_scoring <- 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), FisherInfo = numeric(0))
  
  for (k in 1:max_iter) {
    grad <- gradient(lambda, Y)
    fisher_info <- fisher_information(lambda, Y)
    
    lambda_new <- lambda + grad / fisher_info
    
    iter_results <- rbind(iter_results, data.frame(Iteration = k, Lambda = lambda_new, Gradient = grad, FisherInfo = fisher_info))
    
    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 Fisher Scoring

result <- fisher_scoring(Y)

## Konvergen pada iterasi ke- 2

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 FisherInfo
## 1         1 2.919075 3.320000e+02  173.00000
## 2         2 2.919075 2.842171e-14   59.26535