Bài tập 3

a. Xác hàm log-likelihood cho các hệ số α1 và α2.

Giả sử số vụ tràn dầu \(N_i\) trong năm thứ \(i\) tuân theo phân phối Poisson với tham số \(\lambda_i\):

\[ N_i \mid b_{i1}, b_{i2} \sim \text{Poisson}(\lambda_i) \]

Trong đó:

\[ \lambda_i = \alpha_1 b_{i1} + \alpha_2 b_{i2} \]

  • \(b_{i1}\): Lượng dầu vận chuyển qua hoạt động xuất nhập khẩu.
  • \(b_{i2}\): Lượng dầu vận chuyển qua hoạt động trong nước.
Hàm log-likelihood

Hàm likelihood của mô hình Poisson là:

\[ L(\alpha_1, \alpha_2) = \prod_{i=1}^{n} \frac{e^{-\lambda_i} \lambda_i^{N_i}}{N_i!} \]

Hàm log-likelihood là:

\[ \ell(\alpha_1, \alpha_2) = \sum_{i=1}^{n} \left[ N_i \log(\lambda_i) - \lambda_i - \log(N_i!) \right] \]

Thay \(\lambda_i = \alpha_1 b_{i1} + \alpha_2 b_{i2}\), ta có:

\[ \ell(\alpha_1, \alpha_2) = \sum_{i=1}^{n} \left[ N_i \log\left(\alpha_1 b_{i1} + \alpha_2 b_{i2}\right) - \left(\alpha_1 b_{i1} + \alpha_2 b_{i2}\right) - \log(N_i!) \right] \]

b. Viết các bước giải thuật để xác định ước lượng MLE của α1 và α2 bằng phương pháp Newton

Các bước giải thuật 1. Khởi tạo: - Chọn giá trị ban đầu: \[ \alpha^{(0)} = \bigl(\alpha_1^{(0)},\, \alpha_2^{(0)}\bigr) \]

  1. Tính gradient và Hessian:

    • Tính gradient:

      • Đạo hàm theo \(\alpha_1\): \[ \frac{\partial \ell}{\partial \alpha_1} = \sum_{i=1}^{n} b_{i1}\left(\frac{N_i}{\lambda_i} - 1\right). \]

      • Đạo hàm theo \(\alpha_2\): \[ \frac{\partial \ell}{\partial \alpha_2} = \sum_{i=1}^{n} b_{i2}\left(\frac{N_i}{\lambda_i} - 1\right). \]

      • Vậy vector gradient là: \[ \nabla \ell(\alpha_1, \alpha_2) = \begin{pmatrix} \displaystyle \sum_{i=1}^{n} b_{i1}\left(\frac{N_i}{\lambda_i} - 1\right) \\ \\ \displaystyle \sum_{i=1}^{n} b_{i2}\left(\frac{N_i}{\lambda_i} - 1\right) \end{pmatrix}. \]

    • Tính Hessian:

      • Đạo hàm riêng cấp 2 theo \(\alpha_1\): \[ \frac{\partial^2 \ell}{\partial \alpha_1^2} = -\sum_{i=1}^{n} \frac{N_i \, b_{i1}^2}{\lambda_i^2}. \]

      • Đạo hàm riêng cấp 2 theo \(\alpha_2\): \[ \frac{\partial^2 \ell}{\partial \alpha_2^2} = -\sum_{i=1}^{n} \frac{N_i \, b_{i2}^2}{\lambda_i^2}. \]

      • Đạo hàm hỗn hợp theo \(\alpha_1\) và \(\alpha_2\): \[ \frac{\partial^2 \ell}{\partial \alpha_1 \partial \alpha_2} = -\sum_{i=1}^{n} \frac{N_i \, b_{i1} b_{i2}}{\lambda_i^2}. \]
        => Do đó, ma trận Hessian của hàm log-likelihood là: \[ H(\alpha_1, \alpha_2) = \begin{pmatrix} -\displaystyle \sum_{i=1}^{n} \frac{N_i \, b_{i1}^2}{\lambda_i^2} & -\displaystyle \sum_{i=1}^{n} \frac{N_i \, b_{i1} b_{i2}}{\lambda_i^2} \\[2mm] -\displaystyle \sum_{i=1}^{n} \frac{N_i \, b_{i1} b_{i2}}{\lambda_i^2} & -\displaystyle \sum_{i=1}^{n} \frac{N_i \, b_{i2}^2}{\lambda_i^2} \end{pmatrix}. \]

  2. Cập nhật tham số: \[ \alpha^{(k+1)} = \alpha^{(k)} - H^{-1}\bigl(\alpha^{(k)}\bigr) \, \nabla \ell\bigl(\alpha^{(k)}\bigr) \]

  3. Kiểm tra điều kiện dừng:

    • Nếu \[ \|\alpha^{(k+1)} - \alpha^{(k)}\| < \epsilon \] (với \(\epsilon\) là ngưỡng nhỏ), dừng lại.
    • Ngược lại, quay lại bước 2.

Chương trình thi hành giải thuật

#Đọc dữ liệu
data <- read.table("/Users/kieuhoa.vo/Downloads/0. Study/Master/MTT018 - Thong ke tinh toan/Assignment/Bai tap 1/oilspills.dat", header = TRUE)
spills <- data$spills
importexport <- data$importexport
domestic <- data$domestic

#Hàm log-likelihood
log_likelihood <- function(alpha, spills, importexport, domestic)
{
    lambda <- alpha[1] * importexport + alpha[2] * domestic
    l_likelihood <- sum(-lambda + spills * log(lambda) - log(factorial(spills)))
    return (l_likelihood)
}

#Tính đạo hàm 
gradient <- function(alpha, splills, importexport, domestic)
{
  lambda <- alpha[1] * importexport + alpha[2] * domestic
  d_alpha1 <- sum(-importexport + (splills*importexport)/lambda)
  d_alpha2 <- sum(-domestic + (splills*domestic)/lambda)
  return(c(d_alpha1, d_alpha2))
}

#Tính Hessian matrix
hessian <- function(alpha, splills, importexport, domestic)
{
  lambda <- alpha[1] * importexport + alpha[2] * domestic
  h11 <- sum( (-importexport^2 *splills) / lambda^2)
  h22 <- sum( (-domestic^2 *splills) / lambda^2)
  h12 <- sum( (-spills *importexport *domestic)/lambda^2 )  #đạo hàm hỗn hợp 
  return(matrix(c(h11, h12, h12, h22), nrow = 2, byrow = TRUE))
 }

newton <- function(alpha_init, spills, importexport, domestic, tolerance = 1e-6, max_iter = 150)
{
  alpha <- alpha_init
  path <- matrix(alpha, nrow = 1)  # Khởi tạo đường đi với điểm ban đầu
  for (i in 1:max_iter)
  {
    grad <- gradient(alpha, spills, importexport, domestic)
    hess <- hessian(alpha, spills, importexport, domestic)
    if (det(hess) != 0)
    {
      alpha_new <- alpha - solve(hess) %*% grad
      
      path <- rbind(path, t(alpha_new))  # Lưu giá trị mới vào đường đi
      
      if (max(abs(alpha_new - alpha )) < tolerance)
      {
        #print(cat("Break at irr: ", i ))
        break
      }
      alpha <- alpha_new
    }
  }
  return(list(alpha = alpha, path = path))
}

#Ước lượng MLE bằng pp  newton
alpha_init_nt <- c(0.7, 0.7)

result <- newton(alpha_init_nt, spills, importexport, domestic)
alpha_estimated_nt <- result$alpha
newton_path <- result$path

print(cat("Kết quả Ước lượng MLE bằng pp  newton","alpha1 = ",alpha_estimated_nt[1], ",alpha2 = ", alpha_estimated_nt[2], "\n"))
## Kết quả Ước lượng MLE bằng pp  newton alpha1 =  1.097152 ,alpha2 =  0.9375546 
## NULL

c. Viết các bước giải thuật để xác định ước lượng MLE của α1 và α2 bằng phương pháp Fisher Scoring. Viết một chương trình thi hành giải thuật này. So sánh với kết quả câu (b)

Các bước giải thuật:

1.Khởi tạo Chọn giá trị ban đầu: \[ \alpha^{(0)} = \bigl(\alpha_1^{(0)},\, \alpha_2^{(0)}\bigr). \]

2: Tính gradient và Fisher Information

  • Tính gradient: \[ \nabla \ell\bigl(\alpha^{(k)}\bigr) = \begin{pmatrix} \displaystyle \frac{\partial \ell}{\partial \alpha_1} \\[1mm] \displaystyle \frac{\partial \ell}{\partial \alpha_2} \end{pmatrix}. \]

  • Tính Fisher Information: \[ I\bigl(\alpha^{(k)}\bigr) = \begin{pmatrix} \displaystyle \sum_{i=1}^{n} \frac{b_{i1}^2}{\lambda_i} & \displaystyle \sum_{i=1}^{n} \frac{b_{i1}b_{i2}}{\lambda_i} \\[2mm] \displaystyle \sum_{i=1}^{n} \frac{b_{i1}b_{i2}}{\lambda_i} & \displaystyle \sum_{i=1}^{n} \frac{b_{i2}^2}{\lambda_i} \end{pmatrix}. \]

3: Cập nhật tham số
Cập nhật tham số theo công thức: \[ \alpha^{(k+1)} = \alpha^{(k)} + I^{-1}\bigl(\alpha^{(k)}\bigr)\, \nabla \ell\bigl(\alpha^{(k)}\bigr). \]

4: Kiểm tra điều kiện dừng
Nếu \[ \|\alpha^{(k+1)} - \alpha^{(k)}\| < \epsilon, \]
với \(\epsilon\) là ngưỡng nhỏ ( \(10^{-6}\)), dừng lại. Ngược lại, quay lại Bước 2.

Chương trình thi hành giải thuật

#Tính fisher information matrix 
fisher_information <- function(alpha, importexport, domestic) 
{
  lambda <- alpha[1] * importexport + alpha[2] * domestic
  i11 <- sum((importexport^2)/lambda)
  i22 <- sum((domestic^2)/lambda)
  i12 <- sum((importexport*domestic)/lambda)
  return(matrix(c(i11, i12, i12, i22), nrow = 2, byrow = TRUE))
}

#Phương pháp Fisher Scoring
fisher_scoring <- function (alpha_init, spills, importexport, domestic,  tol = 1e-6, max_iter = 150)
{
  alpha <- alpha_init
  path <- matrix(alpha, nrow = 1)
  for (i in 1:max_iter) 
    {
    # Tính gradient
    grad <- gradient(alpha, spills, importexport, domestic)
    
    # Tính Fisher Information
    fisher_info <- fisher_information(alpha, importexport, domestic)
    
    # Cập nhật tham số
    alpha_new <- alpha + solve(fisher_info) %*% grad
    path <- rbind(path, t(alpha_new))  # Lưu giá trị mới vào đường đi
    
    # Kiểm tra điều kiện dừng
    if (max(abs(alpha_new - alpha)) < tol) {
      break
    }
    alpha <- alpha_new
  }
  return(list(alpha = alpha, path = path))
}
# Ước lượng MLE bằng Fisher Scoring
alpha_init_Fisher <- c(0.7, 0.7)
result <- fisher_scoring(alpha_init_Fisher, spills, importexport, domestic)
alpha_mle_fisher <- result$alpha
fisher_path  <- result$path

print(cat("Kết quả ướclượng MLE bằng Fisher Scoring","alpha1 = ",alpha_mle_fisher[1], ",alpha2 = ", alpha_mle_fisher[2], "\n"))
## Kết quả ướclượng MLE bằng Fisher Scoring alpha1 =  1.097153 ,alpha2 =  0.9375543 
## NULL

d. Ước lượng sai số chuẩn (standard errors) của ước lượng MLE

Sai số chuẩn của các ước lượng MLE \(\alpha_1\) và \(\alpha_2\) được tính bằng căn bậc hai của các phần tử trên đường chéo chính của ma trận nghịch đảo Fisher:

\[ SE(\alpha_1) = \sqrt{[I^{-1}(\alpha)]_{11}}, \]

\[ SE(\alpha_2) = \sqrt{[I^{-1}(\alpha)]_{22}}. \]

# Ước lượng sai số chuẩn
standard_errors <- function(alpha, importexport, domestic) 
  {
    # Tính ma trận thông tin Fisher
    fisher_info <- fisher_information(alpha, importexport, domestic)
    
    # Tính ma trận nghịch đảo của Fisher Information
    inv_fisher_info <- solve(fisher_info)
    
    # Tính sai số chuẩn
    se_alpha1 <- sqrt(inv_fisher_info[1, 1])
    se_alpha2 <- sqrt(inv_fisher_info[2, 2])
    
    return(c(se_alpha1, se_alpha2))
}

#Sử dụng ước lượng MLE từ câu c
standard_err <- standard_errors(alpha_mle_fisher, importexport, domestic)
print(cat("Ước ượng sai số chuẩn (standard errors) của ước lượng MLE","se_alpha1 = ",standard_err[1], ",se_alpha2 = ", standard_err[2], "\n"))
## Ước ượng sai số chuẩn (standard errors) của ước lượng MLE se_alpha1 =  0.437556 ,se_alpha2 =  0.6314687 
## NULL

e. Áp dụng phương pháp quasi-Newton với hai cách chọn M(0): (1) ma trận đơn vị âmvà (2) ma trận thông tin Fisher, −Iα(0). So sánh tính ổn định, tốc độ của hai cách chọn này.

Phương pháp quasi-Newton cập nhật tham số:

\[ \alpha^{(k+1)} = \alpha^{(k)} - \eta_k \, M^{(k)} \, \nabla \ell\bigl(\alpha^{(k)}\bigr) \]

Trong đó: - \(\eta_k\) : Bước học (learning rate) tại lần lặp thứ \(k\). - \(M^{(k)}\) : Xấp xỉ của ma trận nghịch đảo Hessian tại lần lặp thứ \(k\). - \(\nabla \ell\bigl(\alpha^{(k)}\bigr)\) : Gradient của hàm log-likelihood tại \(\alpha^{(k)}\).

Cập nhật ma trận \(M^{(k)}\) theo công thức BFGS:

\[ M^{(k+1)} = M^{(k)} + \frac{\left(s^{(k)} - M^{(k)} y^{(k)}\right)\left(s^{(k)} - M^{(k)} y^{(k)}\right)^T}{\left(s^{(k)} - M^{(k)} y^{(k)}\right)^T y^{(k)}} \]

với:

\[ s^{(k)} = \alpha^{(k+1)} - \alpha^{(k)} \]

# Phương pháp quasi-Newton
quasi_newton <- function(alpha_init, spills, importexport, domestic, M_init, tol = 1e-6, max_iter = 100) 
{
  alpha <- alpha_init
  path <- matrix(alpha, nrow = 1)
  M <- M_init
  for (i in 1:max_iter) {
    grad <- gradient(alpha, spills, importexport, domestic)
    
    alpha_new <- alpha - M %*% grad
    path <- rbind(path, t(alpha_new))  # Lưu giá trị mới vào đường đi

    s <- alpha_new - alpha
    y <- gradient(alpha_new, spills, importexport, domestic) - grad
    M <- M + (s - M %*% y) %*% t(s - M %*% y) / as.numeric(t(s - M %*% y) %*% y)
    if (max(abs(alpha_new - alpha)) < tol)
    {
      break
    }
    alpha <- alpha_new
  }
  return(list(alpha = alpha,path = path ) )
}

alpha_init_quasi <- c(1, 1)

# Cách 1: Ma trận đơn vị âm
M_init_1 <- -diag(2)
cat("Thời gian chạy với ma trận đơn vị âm:\n")
## Thời gian chạy với ma trận đơn vị âm:
time_1 <- round(system.time({
  alpha_mle_1 <- quasi_newton(alpha_init_quasi, spills, importexport, domestic, M_init_1)
}), 12)
cat("Ước lượng MLE với ma trận đơn vị âm:")
## Ước lượng MLE với ma trận đơn vị âm:
print(alpha_mle_1$alpha)
##           [,1]
## [1,] 1.0971525
## [2,] 0.9375546
print(time_1)
##    user  system elapsed 
##   0.005   0.000   0.005
# Cách 2: Ma trận thông tin Fisher
fisher_info <- fisher_information(alpha_init_quasi, importexport, domestic)
M_init_2 <- -solve(fisher_info)
cat("Thời gian chạy với ma trận thông tin Fisher:\n")
## Thời gian chạy với ma trận thông tin Fisher:
time_2 <- round(system.time({
  alpha_mle_2 <- quasi_newton(alpha_init_quasi, spills, importexport, domestic, M_init_2)
}), 12)
cat("Ước lượng MLE với ma trận thông tin Fisher:")
## Ước lượng MLE với ma trận thông tin Fisher:
print(alpha_mle_2$alpha)
##           [,1]
## [1,] 1.0971525
## [2,] 0.9375547
print(time_2)
##    user  system elapsed 
##       0       0       0
# So sánh thời gian chạy
cat("So sánh thời gian chạy:\n")
## So sánh thời gian chạy:
cat("Ma trận đơn vị âm:", time_1["elapsed"], "giây\n")
## Ma trận đơn vị âm: 0.005 giây
cat("Ma trận thông tin Fisher:", time_2["elapsed"], "giây\n")
## Ma trận thông tin Fisher: 0 giây
cat("So sánh số vòng lặp:\n")
## So sánh số vòng lặp:
cat("Ma trận đơn vị âm:", length(alpha_mle_1$path[,1]), "vòng lặp\n")
## Ma trận đơn vị âm: 10 vòng lặp
cat("Ma trận thông tin Fisher:", length(alpha_mle_2$path[,2]), "vòng lặp\n")
## Ma trận thông tin Fisher: 6 vòng lặp

Kết luận so sánh:

  • Ma trận đơn vị âm hội tụ chậm hơn, đặc biệt khi giá trị ban đầu xa giá trị tối ưu. Cách này có thể ổn định hơn khi giá trị ban đầu gần với giá trị tối ưu.

  • Ma trận thông tin Fisher: Hội tụ nhanh hơn vì nó cung cấp một xấp xỉ tốt hơn cho Hessian, tuy nhiên có thể không ổn định nếu Fisher information matrix không được tính chính xác.

f. Xây dựng một đồ thị biểu diễn vùng tối ưu của hàm log-likelihood và để so sánh các đường đi được thực hiện bởi các phương pháp được sử dụng trong (b), (c) và (e). Chọn vùng vẽ đồ thị và điểm bắt đầu để minh họa tốt nhất các tính năng về hiệu suất của thuật toán.

  • Chọn khoảng cho alpha1, alpha2 [0.5,2]
# Tạo lưới giá trị cho alpha1 và alpha2
alpha1 <- seq(0.5, 2, length.out = 100)
alpha2 <- seq(0.5, 2, length.out = 100)
loglik_mat <- matrix(NA, nrow = length(alpha1), ncol = length(alpha2))

# Tính giá trị log-likelihood cho mỗi điểm trên lưới
for (i in 1:length(alpha1)) {
  for (j in 1:length(alpha2)) 
    {
    loglik_mat[i, j] <- log_likelihood(c(alpha1[i], alpha2[j]), spills, importexport, domestic)
    }
}

# Vẽ đồ thị contour
contour(
  x = alpha1, 
  y = alpha2, 
  z = loglik_mat,
  xlab = "α₁",
  ylab = "α₂",
  main = "Contour của Log-Likelihood Function",
  nlevels = 100 # Số đường đồng mức
)

# Thêm đường đi của Newton Method
lines(newton_path[, 1], newton_path[, 2], col = "red", lwd = 1, lty = 1)

  # Thêm điểm bắt đầu và kết thúc
points(alpha_init_nt[1], alpha_init_nt[2], col = "black", pch = 16, cex = 1)  # Điểm bắt đầu
points(alpha_estimated_nt[1], alpha_estimated_nt[2], col = "blue", pch = 16, cex = 0.7)   # Điểm kết thúc

#-----------------
# Thêm đường đi của Fisher scoring
lines(fisher_path[, 1], fisher_path[, 2], col = "green", lwd = 1, lty = 1)

# Thêm điểm bắt đầu và kết thúc
points(alpha_init_Fisher[1], alpha_init_Fisher[2], col = "black", pch = 16, cex = 0.7)  # Điểm bắt đầu
points(alpha_mle_fisher[1], alpha_mle_fisher[2], col = "blue", pch = 4, cex = 1)   # Điểm kết thúc

#-----------------
# Thêm đường đi của quasi-Newton
path_quasi <- alpha_mle_2$path
lines(path_quasi[, 1], path_quasi[, 2], col = "blue", lwd = 1, lty = 1)

# Thêm điểm bắt đầu và kết thúc
points(alpha_init_quasi[1], alpha_init_quasi[2], col = "black", pch = 16, cex = 0.7)  # Điểm bắt đầu
points(alpha_mle_2$alpha[1], alpha_mle_2$alpha[2], col = "red", pch = 16, cex = 0.7)   # Điểm kết thúc

Trong đó:

  • Đường màu đỏ là đường đi của phương pháp Newton
  • Đường màu xanh lá là đường đi của phương pháp Fisher Scoring.
  • Đường màu xanh dương là đường đi của phương pháp quasi-Newton.
  • Màu đen là điểm bắt đầu, điểm màu đỏ là điểm kết thúc (trong trường hợp này các phương pháp gần như cho kết quả xấp xỉ nhau)