Bài tập nhóm số 1

(a) Xác định MLE cho \(\theta\). (Chỉ cần giải một phương trình đa thức đơn giản). Đánh giá ước lượng bằng cách sử dụng dữ liệu mà Dempster, Laird và Rubin đã sử dụng: \(n=197\) và \(x=\left(125, 18, 20 , 34\right)\).

Để xác định MLE cho \(\theta\) ta cần giải phương trình \(\dfrac{dl\left(\theta\right)}{d\theta}=0\) với điều kiện \(0 \le \theta \le 1\).

Bằng các biến đổi sơ cấp ta dễ dàng biến đổi được:

\[\begin{array}{l} \dfrac{x_{1}}{2+\theta } -\dfrac{x_{2} +x_{3}}{1-\theta } +\dfrac{x_{4}}{\theta } =0\\ \Leftrightarrow x_{1}( 1-\theta ) \theta -( x_{2} +x_{3})( 2+\theta ) \theta +x_{4}( 2+\theta )( 1-\theta ) =0\\ \Leftrightarrow -( x_{1} +x_{2} +x_{3} +x_{4}) \theta ^{2} +( x_{1} -2( x_{2} +x_{3}) -x_{4}) \theta +2x_{4} \theta =0 \end{array}\]

Đặt \(a=-\left(x_1+x_2+x_3+x_4\right); b = x_1 -2 \left(x_2+x_3\right) -x_4; c=2x_4.\)

Giải phương trình bậc \(2\) với biến \(\theta\) và kết hợp với điều kiện dễ dàng tính được \(\theta=\)0.627.

Code R:

# Định nghĩa hàm tính MLE cho theta
mle_theta <- function(x1, x2, x3, x4) {
  # Định nghĩa các hệ số của phương trình bậc hai
  a <- -(x1 + x2 + x3 + x4)
  b <- x1 - 2 * (x2 + x3) - x4
  c <- 2 * x4

  # Tính delta (b^2 - 4ac)
  delta <- b^2 - 4 * a * c

  # Kiểm tra nếu delta < 0 thì phương trình vô nghiệm
  if (delta < 0) {
    return("Vô nghiệm")  # Không có nghiệm thực hợp lệ
  }

  # Tính hai nghiệm của phương trình bậc hai
  theta1 <- (-b + sqrt(delta)) / (2 * a)
  theta2 <- (-b - sqrt(delta)) / (2 * a)

  # Kiểm tra nghiệm nằm trong khoảng [0,1]
  theta_valid <- c(theta1, theta2)[c(theta1, theta2) >= 0 & c(theta1, theta2) <= 1]

  # Trả về nghiệm hợp lệ
  if (length(theta_valid) > 0) {
    return(theta_valid)
  } else {
    return("Vô nghiệm")  # Trả về NA nếu không có nghiệm hợp lệ
  }
}

# Ví dụ sử dụng hàm với bộ dữ liệu (125, 18, 20, 34)
theta_mle <- mle_theta(125, 18, 20, 34)
print(theta_mle)

## [1] 0.6268215

(b) Mặc dù dễ dàng tìm thấy nghiệm tối ưu như trong phần trước của bài tập này, nhưng chúng ta cũng có thể sử dụng phương pháp của Newton. Viết chương trình thi hành giải thuật này, bắt đầu với \(\theta^{\left(0\right)}=0.5.\)

newton_method <- function(x1, x2, x3, x4, theta_begin, tol = 1e-10, max_iter = 100){
  theta <- theta_begin #Giá trị khởi tạo
  
  results <- data.frame(Iteration = integer(), Theta = numeric()) #Khởi tạo bảng kết quả
  
  for (i in 1:max_iter){
  #Đạo hàm bậc nhất
  l_prime <- (x1 / (2 + theta)) - ((x2 + x3) / (1 - theta)) + (x4 / theta)
  
  #Đạo hàm bậc hai
  l_double_prime <- -(x1 / (2 + theta)^2) - (x2 + x3) / (1 - theta)^2 - x4 / theta^2
  
  #Cặp nhật theta theo phương pháp Newton
  theta_new <- theta - l_prime / l_double_prime
  
  #Lưu kết quả vào bảng
  results <- rbind(results, data.frame(Iteration = i, Theta = theta_new))
  
  
  #Kiểm tra điều kiện dừng
  if (abs(theta_new - theta) < tol){
    return(results) #Trả về kết quả nếu thoả điều kiện dừng
  }
  
  #Cập nhật theta cho vòng lặp tiếp theo
  theta <- theta_new
  }
  
  #Trả về kết quả nếu không hội tụ
  return("Không hội tụ")
}

# Áp dụng phương pháp Newton với dữ liệu (125, 18, 20, 34) và theta khởi tạo là 0.5
results_b <- newton_method(125, 18, 20, 34, 0.5)
print(format(results_b, digits = 10))

##   Iteration        Theta
## 1         1 0.6363636364
## 2         2 0.6269686722
## 3         3 0.6268215315
## 4         4 0.6268214979
## 5         5 0.6268214979

(c) Viết các bước giải thuật để xác định ước lượng theo phương pháp Fisher Scoring. Viết một chương trình thi hành giải thuật này, bắt đầu với \(\theta^{\left(0\right)}=0.5.\)

fisher_scoring <- function(x1, x2, x3, x4, theta_begin, tol = 1e-10, max_iter = 100) {
  theta <- theta_begin  # Giá trị khởi tạo
  iter_results <- data.frame(Iteration = integer(), Theta = numeric())  # Tạo bảng kết quả

  for (i in 1:max_iter) {
    # Tính đạo hàm bậc nhất l'(theta)
    l_prime <- (x1 / (2 + theta)) - ((x2 + x3) / (1 - theta)) + (x4 / theta)

    # Tính ma trận thông tin Fisher I(theta)
    I_theta <- (x1 / (2 + theta)^2) + ((x2 + x3) / (1 - theta)^2) + (x4 / theta^2)

    # Cập nhật theta theo phương pháp Fisher Scoring
    theta_new <- theta + l_prime / I_theta

    # Lưu kết quả vào bảng
    iter_results <- rbind(iter_results, data.frame(Iteration = i, Theta = theta_new))

    # Kiểm tra điều kiện hội tụ
    if (abs(theta_new - theta) < tol) {
      print("Thuật toán hội tụ!")
      return(iter_results)  # Trả về bảng kết quả
    }

    # Cập nhật giá trị theta
    theta <- theta_new
  }

  print("Thuật toán không hội tụ sau số lần lặp tối đa")
  return(iter_results)  # Trả về bảng kết quả ngay khi không hội tụ
}

# Chạy thử nghiệm với dữ liệu trong đề bài
result_c <- fisher_scoring(125, 18, 20, 34, 0.5)

## [1] "Thuật toán hội tụ!"

print(format(result_c, digits = 10))

##   Iteration        Theta
## 1         1 0.6363636364
## 2         2 0.6269686722
## 3         3 0.6268215315
## 4         4 0.6268214979
## 5         5 0.6268214979

(d) Viết các bước giải thuật để xác định ước lượng theo phương pháp quasi-Newton. Viết một chương trình thi hành giải thuật này, bắt đầu với \(\theta^{\left(0\right)}=0.5.\)

# Khởi tạo dữ liệu
x1 <- 125
x2 <- 18
x3 <- 20
x4 <- 34

# Định nghĩa hàm log-likelihood
logL <- function(theta) {
  if (theta <= 0 || theta >= 1) {
    return(-Inf)  # Tránh log của số âm hoặc 0
  }
  return(x1 * log(2 + theta) + (x2 + x3) * log(1 - theta) + x4 * log(theta))
}

# Định nghĩa đạo hàm của log-likelihood
logL_prime <- function(theta) {
  if (theta <= 0 || theta >= 1) {
    return(NA)  # Tránh giá trị không xác định
  }
  return(x1 / (2 + theta) - (x2 + x3) / (1 - theta) + x4 / theta)
}

# BFGS Algorithm
alpha_default <- 1
alpha <- alpha_default
M <- matrix(-1, 1, 1)  # Hessian xấp xỉ nghịch đảo
iter <- 20  # Số vòng lặp tối đa
theta_0 <- 0.5  # Giá trị khởi tạo
theta_vals <- numeric(iter + 1)
theta_vals[1] <- theta_0
epsilon <- 1e-10  # Ngưỡng hội tụ

for (i in 1:iter) {
  hessian_inv <- solve(M)  # Nghịch đảo của ma trận Hessian xấp xỉ
  score <- logL_prime(theta_0)  # Gradient (đạo hàm cấp 1)
  
  theta_t <- theta_0 - alpha * hessian_inv * score  # Cập nhật theta

  # Điều chỉnh alpha nếu cần
  f_t <- logL(theta_t)
  f_t0 <- logL(theta_0)
  
  while (f_t < f_t0) {  # Nếu log-likelihood giảm thì giảm alpha
    alpha <- alpha / 2
    theta_t <- theta_0 - alpha * hessian_inv * score
    f_t <- logL(theta_t)
  }

  theta_vals[i + 1] <- theta_t  # Lưu giá trị theta mới

  # Cập nhật ma trận M theo BFGS
  ht <- theta_t - theta_0
  score_t <- logL_prime(theta_t)
  ut <- score_t - score
  v <- ut - M * ht
  M_old <- M
  
  if (ht != 0 && ut != 0) {
    M <- M - ((M * ht * t(M * ht)) / (t(ht) * M * ht)) + ((ut * t(ut)) / (t(ht) * ut))
  }

  # Kiểm tra điều kiện cập nhật M
  if (abs(t(v) * ht) < epsilon) {
    M <- M_old
  }

  # Reset alpha và cập nhật theta_0
  alpha <- alpha_default
  theta_0 <- theta_t
}

# In kết quả
print(theta_vals)

##  [1] 0.5000000 0.6640625 0.6208978 0.6264557 0.6268248 0.6268215 0.6268215
##  [8] 0.6268215 0.6268215 0.6268215 0.6268215 0.6268215 0.6268215 0.6268215
## [15] 0.6268215 0.6268215 0.6268215 0.6268215 0.6268215 0.6268215 0.6268215

# Vẽ đồ thị hội tụ của theta
plot(theta_vals, type = "b", pch = 16, col = "blue",
     xlab = "Số vòng lặp", ylab = expression(theta),
     main = "Quá trình hội tụ của theta",
     ylim = c(min(theta_vals), max(theta_vals)))

# Đánh dấu các điểm trên đường đi hội tụ
points(theta_vals, pch = 16, col = "red")

(e) So sánh các phương pháp này như thế nào

3️⃣ Yêu cầu tính toán

• Newton-Raphson: Cần tính đạo hàm cấp 1 và cấp 2 (Hessian).
  
• Fisher Scoring: Cần tính đạo hàm cấp 1 và ma trận thông tin Fisher \(I(\theta)\).
  
• Quasi-Newton: Chỉ cần đạo hàm cấp 1, nhưng phải lưu trữ và cập nhật ma trận Hessian xấp xỉ.

Giả sử dữ liệu của chúng ta là thời gian sống sót đã bị che khuất (censored survival time) \(t_i\) đối với \(i=1,...,n\). Vào cuối nghiên cứu, một bệnh nhân hoặc đã chết (thời gian sống sót đã biết) hoặc vẫn còn sống (thời gian bị che khuất; biết là sống sót ít nhất cho đến khi kết thúc nghiên cứu). Đặt \(w_i\) là \(1\) nếu \(t_i\) là thời gian không bị che khuất và \(0\) nếu \(t_i\) là thời gian bị che khuất. Chứng minh rằng log-likelihood có dạng

\[ \ell(\beta) = \sum_{i=1}^{n} \left( w_i \log(\mu_i) - \mu_i \right) + \sum_{i=1}^{n} w_i \log \left( \frac{\lambda(t_i)}{\Lambda(t_i)} \right) \]

với \(\mu_i = \Lambda(t_i) \exp(x_i^T \beta)\).

Ta có:
Đối với một quan sát không bị che khuất \(t_i\), hàm Likelihood là:

\[ f(t_i) = \lambda(t_i) \exp \left\{ x_i^T \beta - \Lambda(t_i) \exp(x_i^T \beta) \right\} \]

Đối với một quan sát bị che khuất \(t_i\), hàm Likelihood là:

\[ S(t_i) = \exp \left\{ -\Lambda(t_i) \exp(x_i^T \beta) \right\} \]

Chúng ta sử dụng \(w_i\) như một biến chỉ báo, trong đó \(w_i = 1\) nếu \(t_i\) không bị che khuất và \(w_i = 0\) nếu \(t_i\) bị che khuất.

Do đó, hàm log-likelihood kết hợp là:

\[ \ell(\beta) = \sum_{i=1}^{n} w_i \log(f(t_i)) + \sum_{i=1}^{n} (1 - w_i) \log(S(t_i)) \]

Thay thế các biểu thức cho \(\log(f(t_i))\) và \(\log(S(t_i))\):

\[ \begin{align*} \ell(\beta) &= \sum_{i=1}^{n} w_i \left[ \log(\lambda(t_i)) + x_i^T \beta - \Lambda(t_i) \exp(x_i^T \beta) \right] \\ &\quad + \sum_{i=1}^{n} (1 - w_i) \left[ -\Lambda(t_i) \exp(x_i^T \beta) \right] \\ &= \sum_{i=1}^{n} w_i \log(\lambda(t_i)) + \sum_{i=1}^{n} w_i x_i^T \beta - \sum_{i=1}^{n} w_i\Lambda(t_i) \exp(x_i^T \beta) \\ &\quad -\sum_{i=1}^{n}\Lambda(t_i) \exp(x_i^T \beta)+\sum_{i=1}^{n} w_i\Lambda(t_i) \exp(x_i^T \beta) \\ &= \sum_{i=1}^{n} w_i \log(\lambda(t_i)) + \sum_{i=1}^{n} w_i x_i^T \beta - \sum_{i=1}^{n} \Lambda(t_i) \exp(x_i^T \beta) \end{align*} \]

Sử dụng \(\mu_i = \Lambda(t_i) \exp(x_i^T \beta)\):

\[ \ell(\beta) = \sum_{i=1}^{n} w_i \log(\lambda(t_i)) + \sum_{i=1}^{n} w_i x_i^T \beta - \sum_{i=1}^{n} \mu_i \]

Ta có: \(\mu_i = \Lambda(t_i) \exp(x_i^T \beta)\), nên \(\log(\mu_i) = \log(\Lambda(t_i)) + x_i^T \beta\) và \(x_i^T \beta= \log(\mu_i)-\log(\Lambda(t_i))\)

Suy ra:

\[\begin{align*} \ell(\beta) &= \sum_{i=1}^{n} w_i (\log(\lambda(t_i)) + x_i^T \beta) - \sum_{i=1}^{n} \mu_i \\ &= \sum_{i=1}^{n} w_i (\log(\lambda(t_i)) + \log(\mu_i) - \log(\Lambda(t_i))) - \sum_{i=1}^{n} \mu_i \\ &= \sum_{i=1}^{n} \left( w_i \log(\mu_i) - \mu_i \right) + \sum_{i=1}^{n} w_i (\log (\lambda(t_i)) - \log(\Lambda(t_i))) \\ &= \sum_{i=1}^{n} \left( w_i \log(\mu_i) - \mu_i \right) + \sum_{i=1}^{n} w_i \log \left( \frac{\lambda(t_i)}{\Lambda(t_i)} \right) \end{align*} \] Hoàn tất chứng minh.

Hãy xem xét một mô hình về thời gian thuyên giảm cho bệnh nhân mắc bệnh bạch cầu cấp tính trong một thử nghiệm lâm sàng. Bệnh nhân được điều trị bằng thuốc 6-mercaptopurine. (6-MP) hoặc giả dược (E. J. Freireich, E. Gehan, E. Frei III, L. R. Schroeder, I. J. Wolman, R. Anabari, E. O. Burgert, S. D. Mills, D. Pinkel, O. S. Selawry, J. H. Moon, B. R. Gendel, C. L. Spurr, R. Storrs, F. Haurani, B. Hoogstraten, and S. Lee, n.d.). Một năm sau khi bắt đầu nghiên cứu, thời gian (đơn vị: tuần) của thời gian thuyên giảm giảm cho mỗi bệnh nhân đã được ghi lại mercaptopurine.csv. Một số kết quả đã bị che khuất vì thời gian thuyên giảm kéo dài sau thời gian nghiên cứu. Mục tiêu là xác định liệu phương pháp điều trị có kéo dài thời gian thuyên giảm hay không. Giả sử chúng ta đặt \(\Lambda\left(t\right)=t^\alpha\) với \(\alpha >0\), khi đó, ta xác định được \(\lambda\left(t\right)=\alpha t^{\alpha-1}\) và \(f\left(t\right)=\alpha t^{\alpha-1}exp\{x^T\beta-t^{\alpha}exp\left(x^T\beta\right)\}\) (hàm mật độ Weibull). Áp dụng tham số hoá biến phụ thuộc: \(x^T_i\beta=\beta_0+\delta_i\beta_i\), trong đó \(\delta_i\) là \(1\) nếu bệnh nhân thứ \(i\) nằm trong nhóm điều trị (treatment) và bằng \(0\) nếu bệnh nhân đó nằm trong nhóm đối chứng (nhóm control). Viết các giải thuật để xác định ước lượng MLE của \(\alpha, \beta_0, \beta_1\) bằng phương pháp Newton và Fisher Scoring. Viết chương trình thi hành các giải thuật này.

Xét mô hình thời gian thuyên giảm với hàm nguy cơ tỷ lệ như sau:

\[\Lambda(t) = t^\alpha \quad \text{với} \quad \alpha > 0\]

Hàm nguy cơ tỷ lệ:

\[\lambda(t) = \alpha t^{\alpha-1}\]

Hàm mật độ Weibull:

\[f(t) = \alpha t^{\alpha-1} \exp \left( x^\top \beta - t^\alpha \exp \left( x^\top \beta \right) \right)\]

Tham số hóa biến phụ thuộc được đưa ra bởi \(x_i^\top \beta = \beta_0 + \delta_i \beta_1\), trong đó:

• \(\delta_i\) là 1 nếu bệnh nhân thứ i nằm trong nhóm điều trị (treatment).

• \(\delta_i\) là 0 nếu bệnh nhân thứ i nằm trong nhóm kiểm soát (control).

\(\mu_i = \Lambda(t_i) \exp(x_i^T \beta) = t_i^\alpha\exp(\beta_0 + \delta_i \beta_1)\)

Giải thuật bằng Newton-Raphson và Fisher Scoring

1. Phương pháp Newton-Raphson:

Ta cần tìm ước lượng cực đại của hàm likelihood

Phương pháp Newton-Raphson cập nhật các tham số bằng cách:

\[ \theta^{(k+1)} = \theta^{(k)} - \left( \nabla^2 \ell(\theta^{(k)}) \right)^{-1} \nabla \ell(\theta^{(k)}) \]

Trong đó:

• \(\theta = (\alpha, \beta_0, \beta_1)\)

• \(\nabla \ell(\theta)\) là gradient của hàm log-likelihood.

• \(\nabla^2 \ell(\theta)\) là ma trận Hessian của hàm log-likelihood.

2. Phương pháp Fisher Scoring:

Tương tự như giải thuật Newton, nhưng thay ma trận Hessian bằng kỳ vọng của ma trận Hessian:

\[\theta^{(k+1)} = \theta^{(k)} - \left( \mathbb{E}\left[ \nabla^2 \ell(\theta^{(k)}) \right] \right)^{-1} \nabla \ell(\theta^{(k)})\] Code R:

# Load necessary libraries
  library(readr)

# Đọc dữ liệu từ file CSV

data <- read.csv("C:/Users/AZ/Desktop/Study/Master_Data/Thongke/CodeR/BT1/mercaptopurine.csv")  # Đường dẫn file CSV

treatment <- ifelse(data$group == "treatment", 1, 0)
status <- as.numeric(data$censored)
time <- as.numeric(data$time)

# Định nghĩa hàm log-likelihood Weibull
logLik_Weibull <- function(params, time, status, treatment) {
  alpha <- params[1]
  beta0 <- params[2]
  beta1 <- params[3]
  
  x_beta <- beta0 + beta1 * treatment
  lambda_t <- alpha * time^(alpha - 1)  
  Lambda_t <- time^alpha
  mu_i <- Lambda_t * exp(x_beta)
  
  log_lik <- sum(status * log(mu_i) - mu_i) + sum(status * log(alpha / time))
  
  return(-log_lik)  # Đảo dấu vì optim() tối thiểu hóa
}

# Hàm tính gradient (đạo hàm bậc nhất)
gradient <- function(params, time, status, treatment) {
  alpha <- params[1]
  beta0 <- params[2]
  beta1 <- params[3]
  
  x_beta <- beta0 + beta1 * treatment
  lambda_t <- alpha * time^(alpha - 1)
  Lambda_t <- time^alpha
  mu_i <- Lambda_t * exp(x_beta)
  
  d_alpha <- sum(status/alpha + status * log(time) - mu_i * log(time))
  d_beta0 <- sum(status - mu_i)
  d_beta1 <- sum(treatment * (status - mu_i))
  
  return(c(d_alpha, d_beta0, d_beta1))
}

# Khởi tạo giá trị ban đầu cho các tham số
initial_params <- c(1, 0, 0)

# Hàm tính Hessian (đạo hàm bậc hai, dùng cho Newton-Raphson)
hessian <- function(params, time, status, treatment) {
  alpha <- params[1]
  beta0 <- params[2]
  beta1 <- params[3]
  
  x_beta <- beta0 + beta1 * treatment
  lambda_t <- alpha * time^(alpha - 1)
  Lambda_t <- time^alpha
  mu_i <- Lambda_t * exp(x_beta)
  
  d2_alpha <- sum(-status/alpha^2 - mu_i * log(time)^2)
  d2_beta0 <- -sum(mu_i)
  d2_beta1 <- -sum(treatment^2 * mu_i)
  d_alpha_beta0 <- -sum(mu_i * log(time))
  d_alpha_beta1 <- -sum(treatment * mu_i * log(time))
  d2_beta0_beta1 <- -sum(treatment * mu_i)
  
  H <- matrix(c(d2_alpha, d_alpha_beta0, d_alpha_beta1,
                d_alpha_beta0, d2_beta0, d2_beta0_beta1,
                d_alpha_beta1, d2_beta0_beta1, d2_beta1), 3, 3, byrow = TRUE)
  
  return(-H)
}

# Newton-Raphson
newton_raphson <- function(gradient, initial_params, time, status, treatment, tol = 1e-6, max_iter = 100) {
  params <- initial_params
  for (i in 1:max_iter) {
    grad <- gradient(params, time, status, treatment)
    H <- -hessian(params, time, status, treatment)  
    H <- H + diag(1e-6, nrow(H))# Thêm số nhỏ vào đường chéo để tính được nghịch đảo khi ma trận gần suy biến.
    step <- solve(H) %*% grad
    params <- params - step
    # In thông tin hội tụ
    cat(sprintf("Iter %d | log-likelihood: %.6f\n", i, -logLik_Weibull(params, time, status, treatment)))
    if (sqrt(sum(step^2)) < tol) {
  
      break
    }
  }
  return(params)
}

# Fisher Scoring
fisher_scoring <- function(initial_params, time, status, treatment, tol = 1e-6, max_iter = 100) {
  params <- initial_params
  for (i in 1:max_iter) {
    grad <- gradient(params, time, status, treatment)
    I_fisher <- -hessian(params, time, status, treatment)  
    I_fisher <- I_fisher + diag(1e-6, nrow(I_fisher))  # Thêm số nhỏ vào đường chéo để tính được nghịch đảo khi ma trận gần suy biến.
    inv_I_fisher <- tryCatch(solve(I_fisher), error = function(e) NA)
    update <- inv_I_fisher %*% grad  
    params <- params - update
    # In thông tin hội tụ
    cat(sprintf("Iter %d | log-likelihood: %.6f\n", i, -logLik_Weibull(params, time, status, treatment)))
    if (max(abs(update)) < tol) {
      
      break
    }
  }
  return(params)
}

# Chạy phương pháp Newton-Raphson
cat("Ước lượng bằng phương pháp Newton-Raphson:\n")

## Ước lượng bằng phương pháp Newton-Raphson:

params_newton <- newton_raphson(gradient, initial_params, time, status, treatment)

## Iter 1 | log-likelihood: -214.626911
## Iter 2 | log-likelihood: -100.472229
## Iter 3 | log-likelihood: -62.830467
## Iter 4 | log-likelihood: -51.421739
## Iter 5 | log-likelihood: -48.146876
## Iter 6 | log-likelihood: -47.211070
## Iter 7 | log-likelihood: -46.906381
## Iter 8 | log-likelihood: -46.795223
## Iter 9 | log-likelihood: -46.754331
## Iter 10 | log-likelihood: -46.739288
## Iter 11 | log-likelihood: -46.733754
## Iter 12 | log-likelihood: -46.731718
## Iter 13 | log-likelihood: -46.730969
## Iter 14 | log-likelihood: -46.730694
## Iter 15 | log-likelihood: -46.730592
## Iter 16 | log-likelihood: -46.730555
## Iter 17 | log-likelihood: -46.730541
## Iter 18 | log-likelihood: -46.730536
## Iter 19 | log-likelihood: -46.730691
## Iter 20 | log-likelihood: -46.730591
## Iter 21 | log-likelihood: -46.730555
## Iter 22 | log-likelihood: -46.730541
## Iter 23 | log-likelihood: -46.730536
## Iter 24 | log-likelihood: -46.730613
## Iter 25 | log-likelihood: -46.730563
## Iter 26 | log-likelihood: -46.730544
## Iter 27 | log-likelihood: -46.730537
## Iter 28 | log-likelihood: -46.730535
## Iter 29 | log-likelihood: -46.730535
## Iter 30 | log-likelihood: -46.730535
## Iter 31 | log-likelihood: -46.730535
## Iter 32 | log-likelihood: -46.730535
## Iter 33 | log-likelihood: -46.730535
## Iter 34 | log-likelihood: -46.730535
## Iter 35 | log-likelihood: -46.730535
## Iter 36 | log-likelihood: -46.730535
## Iter 37 | log-likelihood: -46.730535
## Iter 38 | log-likelihood: -46.730535
## Iter 39 | log-likelihood: -46.730535
## Iter 40 | log-likelihood: -46.730535
## Iter 41 | log-likelihood: -46.730535
## Iter 42 | log-likelihood: -46.730535
## Iter 43 | log-likelihood: -46.730535
## Iter 44 | log-likelihood: -46.730535
## Iter 45 | log-likelihood: -46.730535
## Iter 46 | log-likelihood: -46.730535
## Iter 47 | log-likelihood: -46.730535
## Iter 48 | log-likelihood: -46.730535
## Iter 49 | log-likelihood: -46.730535
## Iter 50 | log-likelihood: -46.730535
## Iter 51 | log-likelihood: -46.730535
## Iter 52 | log-likelihood: -46.730535
## Iter 53 | log-likelihood: -46.730535
## Iter 54 | log-likelihood: -46.730535
## Iter 55 | log-likelihood: -46.730535
## Iter 56 | log-likelihood: -46.730535
## Iter 57 | log-likelihood: -46.730535
## Iter 58 | log-likelihood: -46.730536
## Iter 59 | log-likelihood: -46.730537
## Iter 60 | log-likelihood: -46.730535
## Iter 61 | log-likelihood: -46.730535
## Iter 62 | log-likelihood: -46.730535
## Iter 63 | log-likelihood: -46.730535
## Iter 64 | log-likelihood: -46.730535
## Iter 65 | log-likelihood: -46.730535
## Iter 66 | log-likelihood: -46.730535
## Iter 67 | log-likelihood: -46.730535
## Iter 68 | log-likelihood: -46.730535
## Iter 69 | log-likelihood: -46.730535
## Iter 70 | log-likelihood: -46.730535
## Iter 71 | log-likelihood: -46.730535
## Iter 72 | log-likelihood: -46.730535
## Iter 73 | log-likelihood: -46.730535
## Iter 74 | log-likelihood: -46.730535
## Iter 75 | log-likelihood: -46.730535
## Iter 76 | log-likelihood: -46.730535
## Iter 77 | log-likelihood: -46.730535
## Iter 78 | log-likelihood: -46.730535
## Iter 79 | log-likelihood: -46.730535
## Iter 80 | log-likelihood: -46.730535
## Iter 81 | log-likelihood: -46.730535
## Iter 82 | log-likelihood: -46.730535
## Iter 83 | log-likelihood: -46.730535
## Iter 84 | log-likelihood: -46.730535
## Iter 85 | log-likelihood: -46.730535
## Iter 86 | log-likelihood: -46.730535
## Iter 87 | log-likelihood: -46.730535
## Iter 88 | log-likelihood: -46.730535
## Iter 89 | log-likelihood: -46.730535
## Iter 90 | log-likelihood: -46.730535
## Iter 91 | log-likelihood: -46.730535
## Iter 92 | log-likelihood: -46.730535
## Iter 93 | log-likelihood: -46.730535
## Iter 94 | log-likelihood: -46.730535
## Iter 95 | log-likelihood: -46.730535
## Iter 96 | log-likelihood: -46.730535
## Iter 97 | log-likelihood: -46.730535
## Iter 98 | log-likelihood: -46.730535
## Iter 99 | log-likelihood: -46.730535
## Iter 100 | log-likelihood: -46.730535

print(params_newton)

##           [,1]
## [1,]   2.50892
## [2,] -27.30059
## [3,]  19.12388

# Chạy phương pháp Fisher Scoring
cat("\nƯớc lượng bằng phương pháp Fisher Scoring:\n")

## 
## Ước lượng bằng phương pháp Fisher Scoring:

params_fisher <- fisher_scoring(initial_params, time, status, treatment)

## Iter 1 | log-likelihood: -214.626911
## Iter 2 | log-likelihood: -100.472229
## Iter 3 | log-likelihood: -62.830467
## Iter 4 | log-likelihood: -51.421739
## Iter 5 | log-likelihood: -48.146876
## Iter 6 | log-likelihood: -47.211070
## Iter 7 | log-likelihood: -46.906381
## Iter 8 | log-likelihood: -46.795223
## Iter 9 | log-likelihood: -46.754331
## Iter 10 | log-likelihood: -46.739288
## Iter 11 | log-likelihood: -46.733754
## Iter 12 | log-likelihood: -46.731718
## Iter 13 | log-likelihood: -46.730969
## Iter 14 | log-likelihood: -46.730694
## Iter 15 | log-likelihood: -46.730592
## Iter 16 | log-likelihood: -46.730555
## Iter 17 | log-likelihood: -46.730541
## Iter 18 | log-likelihood: -46.730536
## Iter 19 | log-likelihood: -46.730691
## Iter 20 | log-likelihood: -46.730591
## Iter 21 | log-likelihood: -46.730555
## Iter 22 | log-likelihood: -46.730541
## Iter 23 | log-likelihood: -46.730536
## Iter 24 | log-likelihood: -46.730613
## Iter 25 | log-likelihood: -46.730563
## Iter 26 | log-likelihood: -46.730544
## Iter 27 | log-likelihood: -46.730537
## Iter 28 | log-likelihood: -46.730535
## Iter 29 | log-likelihood: -46.730535
## Iter 30 | log-likelihood: -46.730535
## Iter 31 | log-likelihood: -46.730535
## Iter 32 | log-likelihood: -46.730535
## Iter 33 | log-likelihood: -46.730535
## Iter 34 | log-likelihood: -46.730535
## Iter 35 | log-likelihood: -46.730535
## Iter 36 | log-likelihood: -46.730535
## Iter 37 | log-likelihood: -46.730535
## Iter 38 | log-likelihood: -46.730535
## Iter 39 | log-likelihood: -46.730535
## Iter 40 | log-likelihood: -46.730535
## Iter 41 | log-likelihood: -46.730535
## Iter 42 | log-likelihood: -46.730535
## Iter 43 | log-likelihood: -46.730535
## Iter 44 | log-likelihood: -46.730535
## Iter 45 | log-likelihood: -46.730535
## Iter 46 | log-likelihood: -46.730535
## Iter 47 | log-likelihood: -46.730535
## Iter 48 | log-likelihood: -46.730535
## Iter 49 | log-likelihood: -46.730535
## Iter 50 | log-likelihood: -46.730535
## Iter 51 | log-likelihood: -46.730535
## Iter 52 | log-likelihood: -46.730535
## Iter 53 | log-likelihood: -46.730535
## Iter 54 | log-likelihood: -46.730535
## Iter 55 | log-likelihood: -46.730535
## Iter 56 | log-likelihood: -46.730535
## Iter 57 | log-likelihood: -46.730535
## Iter 58 | log-likelihood: -46.730536
## Iter 59 | log-likelihood: -46.730537
## Iter 60 | log-likelihood: -46.730535
## Iter 61 | log-likelihood: -46.730535
## Iter 62 | log-likelihood: -46.730535
## Iter 63 | log-likelihood: -46.730535
## Iter 64 | log-likelihood: -46.730535
## Iter 65 | log-likelihood: -46.730535
## Iter 66 | log-likelihood: -46.730535
## Iter 67 | log-likelihood: -46.730535
## Iter 68 | log-likelihood: -46.730535
## Iter 69 | log-likelihood: -46.730535
## Iter 70 | log-likelihood: -46.730535
## Iter 71 | log-likelihood: -46.730535
## Iter 72 | log-likelihood: -46.730535
## Iter 73 | log-likelihood: -46.730535
## Iter 74 | log-likelihood: -46.730535
## Iter 75 | log-likelihood: -46.730535
## Iter 76 | log-likelihood: -46.730535
## Iter 77 | log-likelihood: -46.730535
## Iter 78 | log-likelihood: -46.730535
## Iter 79 | log-likelihood: -46.730535
## Iter 80 | log-likelihood: -46.730535
## Iter 81 | log-likelihood: -46.730535
## Iter 82 | log-likelihood: -46.730535
## Iter 83 | log-likelihood: -46.730535
## Iter 84 | log-likelihood: -46.730535
## Iter 85 | log-likelihood: -46.730535
## Iter 86 | log-likelihood: -46.730535
## Iter 87 | log-likelihood: -46.730535
## Iter 88 | log-likelihood: -46.730535
## Iter 89 | log-likelihood: -46.730535
## Iter 90 | log-likelihood: -46.730535
## Iter 91 | log-likelihood: -46.730535
## Iter 92 | log-likelihood: -46.730535
## Iter 93 | log-likelihood: -46.730535
## Iter 94 | log-likelihood: -46.730535
## Iter 95 | log-likelihood: -46.730535
## Iter 96 | log-likelihood: -46.730535
## Iter 97 | log-likelihood: -46.730535
## Iter 98 | log-likelihood: -46.730535
## Iter 99 | log-likelihood: -46.730535
## Iter 100 | log-likelihood: -46.730535

print(params_fisher)

##           [,1]
## [1,]   2.50892
## [2,] -27.30059
## [3,]  19.12388

(c) Viết các bước giải thuật để xác định ước lượng MLE của \(\alpha, \beta_0, \beta_1\) bằng phương pháp quasi-Newton. Viết một chương trình thi hành giải thuật.

Code R:

# Đọc file CSV

library(readr)
data <- read.csv("C:/Users/AZ/Desktop/Study/Master_Data/Thongke/CodeR/BT1/mercaptopurine.csv")  # Đường dẫn file CSV

treatment <- ifelse(data$group == "treatment", 1, 0)
status <- as.numeric(data$censored)
time <- as.numeric(data$time)

# Định nghĩa log-likelihood Weibull
f <- function(x, time, status, treatment) {
  alpha <- x[1]
  beta0 <- x[2]
  beta1 <- x[3]
  
  if (alpha <= 0) {
    return(Inf) # Tránh trường hợp alpha <= 0
  }
  
  x_beta <- beta0 + beta1 * treatment
  Lambda_t <- time^alpha
  mu_i <- Lambda_t * exp(x_beta)
  
  if (any(mu_i <= 0) || any(time <= 0)) {
    return(Inf)  # Tránh các giá trị log không hợp lệ
  }
  
  log_lik <- sum(status * log(mu_i) - mu_i) + sum(status * log(alpha / time))
  
  return(-log_lik)  # Đảo dấu vì optim() tối thiểu hóa
}

f_prime <- function(x, time, status, treatment) {
  alpha <- x[1]
  beta0 <- x[2]
  beta1 <- x[3]
  
  if (alpha <= 0) {
    return(c(NA, NA, NA))  # Tránh trường hợp alpha <= 0
  }
  
  x_beta <- beta0 + beta1 * treatment
  Lambda_t <- time^alpha
  mu_i <- Lambda_t * exp(x_beta)
  
  if (any(mu_i <= 0) || any(time <= 0)) {
    return(c(NA, NA, NA))  # Tránh các giá trị log không hợp lệ
  }
  
  f_t <- alpha * (time^(alpha-1)) * exp((time^alpha) * exp(beta0 + treatment * beta1))
  
  d_alpha <- sum(status/alpha + status * log(time) - mu_i * log(time))
  d_beta0 <- sum(status - mu_i)
  d_beta1 <- sum(treatment * (status - mu_i))
  
  return(c(d_alpha, d_beta0, d_beta1))
}


# BFGS Quasi-Newton
M <- diag(c(-1,-1,-1), 3, 3)
k_default <- 1
k <- k_default
x0 <- c(1,0,0)
iter <- 42
x_vals_bfgs <- matrix(NA, nrow=iter+1, ncol=3)
x_vals_bfgs[1,] <- x0
epsilon <- 1e-10

for (i in 1:iter) {
  hessian_inv <- solve(M)
  score <- f_prime(x0, time, status, treatment)
  xt <- x0 - k * hessian_inv %*% score
  
  # Kia;m tra gia:#m log-likelihood
  f_t <- f(xt, time, status, treatment)
  f_t0 <- f(x0, time, status, treatment)
  
  while (f_t > f_t0) {  # Nếu log-likelihood tăng thì giảm k
    k <- k / 2
    xt <- x0 - k * hessian_inv %*% score
    f_t <- f(xt, time, status, treatment)
  }
  
  x_vals_bfgs[i+1,] <- xt
  
  ## Cập nhật ma trận M theo BFGS
  ht <- xt - x0
  score_t <- f_prime(xt, time, status, treatment)
  ut <- score_t - score
  
  v <- ut - M %*% ht
  M_old <- M
  
  M <- M - ((M %*% ht %*% t(M %*% ht)) / as.numeric(t(ht) %*% M %*% ht)) +
    ((ut %*% t(ut)) / as.numeric(t(ht) %*% ut))
  
  # Kiểm tra điều kiện cập nhật M
  if (abs(t(ht) %*% ut) < epsilon) {
    M <- M_old
  }
  
  ## Reset alpha và  cập nhật x_0
  k <- k_default
  x0 <- xt
}

# In giá trị của nghiệm lập
head(x_vals_bfgs, n = 40)

##            [,1]        [,2]       [,3]
##  [1,] 1.0000000   0.0000000  0.0000000
##  [2,] 0.2814039  -0.2583008 -0.1694336
##  [3,] 4.1582689 -11.1667462 -0.8451386
##  [4,] 3.1096409  -8.5220573 -0.6953237
##  [5,] 2.2006902  -7.6508261  1.6644693
##  [6,] 1.8782574  -6.4730475  0.9458558
##  [7,] 1.8782939  -6.4325883  0.6771135
##  [8,] 1.9309632  -6.7158016  0.5363962
##  [9,] 2.0407991  -7.3065544  0.5954402
## [10,] 2.2114175  -8.2071530  0.8406507
## [11,] 2.4114451  -9.3010503  1.2784471
## [12,] 2.5880960 -10.4091982  1.8865781
## [13,] 2.6556742 -11.1022185  2.4341223
## [14,] 2.6445246 -11.4769415  2.8807247
## [15,] 2.5794038 -11.8985936  3.5320519
## [16,] 2.5046246 -12.3323591  4.1984237
## [17,] 2.4740589 -12.7730744  4.7205880
## [18,] 2.4739449 -13.2749694  5.2093697
## [19,] 2.4900680 -13.9195983  5.7944325
## [20,] 2.5079524 -14.6073201  6.4251363
## [21,] 2.5167438 -15.2487213  7.0429060
## [22,] 2.5169243 -15.8363542  7.6340572
## [23,] 2.5130041 -16.4421786  8.2545307
## [24,] 2.5088772 -17.0809888  8.9065172
## [25,] 2.5069858 -17.7203059  9.5507692
## [26,] 2.5071317 -18.3527094 10.1816750
## [27,] 2.5081365 -18.9984955 10.8237651
## [28,] 2.5090362 -19.6614727 11.4839017
## [29,] 2.5093918 -20.3246733 12.1462378
## [30,] 2.5093031 -20.9804524 12.8025378
## [31,] 2.5090576 -21.6394037 13.4623723
## [32,] 2.5088679 -22.3086421 14.1321980
## [33,] 2.5088080 -22.9811909 14.8048729
## [34,] 2.5088408 -23.6504861 15.4740116
## [35,] 2.5088984 -24.3208091 16.1441320
## [36,] 2.5089369 -24.9976326 16.8208393
## [37,] 2.5089455 -25.6772309 17.5004253
## [38,] 2.5089356 -26.3541355 18.1773726
## [39,] 2.5089226 -27.0310145 18.8542965
## [40,] 2.5089150 -27.7127451 19.5360490

(d) Ước lượng sai số chuẩn (standard errors) của ước lượng MLE. Có bất kỳ ước lượng MLE nào có tương quan mạnh hay không?

# Ước lượng sai số chuẩn của các ước lượng MLE
standard_errors <- function(params, time, status, treatment) {
  H <- -hessian(params, time, status, treatment)
  H <-H + diag(1e-6, nrow(H)) 
  H_inv<-solve(H)
  se <- sqrt(diag(H_inv))
  return(se)
}
# Chạy phương pháp Newton-Raphson
cat("Ước lượng bằng phương pháp Newton-Raphson:\n")

## Ước lượng bằng phương pháp Newton-Raphson:

params_newton <- newton_raphson(gradient, initial_params, time, status, treatment)

## Iter 1 | log-likelihood: -214.626911
## Iter 2 | log-likelihood: -100.472229
## Iter 3 | log-likelihood: -62.830467
## Iter 4 | log-likelihood: -51.421739
## Iter 5 | log-likelihood: -48.146876
## Iter 6 | log-likelihood: -47.211070
## Iter 7 | log-likelihood: -46.906381
## Iter 8 | log-likelihood: -46.795223
## Iter 9 | log-likelihood: -46.754331
## Iter 10 | log-likelihood: -46.739288
## Iter 11 | log-likelihood: -46.733754
## Iter 12 | log-likelihood: -46.731718
## Iter 13 | log-likelihood: -46.730969
## Iter 14 | log-likelihood: -46.730694
## Iter 15 | log-likelihood: -46.730592
## Iter 16 | log-likelihood: -46.730555
## Iter 17 | log-likelihood: -46.730541
## Iter 18 | log-likelihood: -46.730536
## Iter 19 | log-likelihood: -46.730691
## Iter 20 | log-likelihood: -46.730591
## Iter 21 | log-likelihood: -46.730555
## Iter 22 | log-likelihood: -46.730541
## Iter 23 | log-likelihood: -46.730536
## Iter 24 | log-likelihood: -46.730613
## Iter 25 | log-likelihood: -46.730563
## Iter 26 | log-likelihood: -46.730544
## Iter 27 | log-likelihood: -46.730537
## Iter 28 | log-likelihood: -46.730535
## Iter 29 | log-likelihood: -46.730535
## Iter 30 | log-likelihood: -46.730535
## Iter 31 | log-likelihood: -46.730535
## Iter 32 | log-likelihood: -46.730535
## Iter 33 | log-likelihood: -46.730535
## Iter 34 | log-likelihood: -46.730535
## Iter 35 | log-likelihood: -46.730535
## Iter 36 | log-likelihood: -46.730535
## Iter 37 | log-likelihood: -46.730535
## Iter 38 | log-likelihood: -46.730535
## Iter 39 | log-likelihood: -46.730535
## Iter 40 | log-likelihood: -46.730535
## Iter 41 | log-likelihood: -46.730535
## Iter 42 | log-likelihood: -46.730535
## Iter 43 | log-likelihood: -46.730535
## Iter 44 | log-likelihood: -46.730535
## Iter 45 | log-likelihood: -46.730535
## Iter 46 | log-likelihood: -46.730535
## Iter 47 | log-likelihood: -46.730535
## Iter 48 | log-likelihood: -46.730535
## Iter 49 | log-likelihood: -46.730535
## Iter 50 | log-likelihood: -46.730535
## Iter 51 | log-likelihood: -46.730535
## Iter 52 | log-likelihood: -46.730535
## Iter 53 | log-likelihood: -46.730535
## Iter 54 | log-likelihood: -46.730535
## Iter 55 | log-likelihood: -46.730535
## Iter 56 | log-likelihood: -46.730535
## Iter 57 | log-likelihood: -46.730535
## Iter 58 | log-likelihood: -46.730536
## Iter 59 | log-likelihood: -46.730537
## Iter 60 | log-likelihood: -46.730535
## Iter 61 | log-likelihood: -46.730535
## Iter 62 | log-likelihood: -46.730535
## Iter 63 | log-likelihood: -46.730535
## Iter 64 | log-likelihood: -46.730535
## Iter 65 | log-likelihood: -46.730535
## Iter 66 | log-likelihood: -46.730535
## Iter 67 | log-likelihood: -46.730535
## Iter 68 | log-likelihood: -46.730535
## Iter 69 | log-likelihood: -46.730535
## Iter 70 | log-likelihood: -46.730535
## Iter 71 | log-likelihood: -46.730535
## Iter 72 | log-likelihood: -46.730535
## Iter 73 | log-likelihood: -46.730535
## Iter 74 | log-likelihood: -46.730535
## Iter 75 | log-likelihood: -46.730535
## Iter 76 | log-likelihood: -46.730535
## Iter 77 | log-likelihood: -46.730535
## Iter 78 | log-likelihood: -46.730535
## Iter 79 | log-likelihood: -46.730535
## Iter 80 | log-likelihood: -46.730535
## Iter 81 | log-likelihood: -46.730535
## Iter 82 | log-likelihood: -46.730535
## Iter 83 | log-likelihood: -46.730535
## Iter 84 | log-likelihood: -46.730535
## Iter 85 | log-likelihood: -46.730535
## Iter 86 | log-likelihood: -46.730535
## Iter 87 | log-likelihood: -46.730535
## Iter 88 | log-likelihood: -46.730535
## Iter 89 | log-likelihood: -46.730535
## Iter 90 | log-likelihood: -46.730535
## Iter 91 | log-likelihood: -46.730535
## Iter 92 | log-likelihood: -46.730535
## Iter 93 | log-likelihood: -46.730535
## Iter 94 | log-likelihood: -46.730535
## Iter 95 | log-likelihood: -46.730535
## Iter 96 | log-likelihood: -46.730535
## Iter 97 | log-likelihood: -46.730535
## Iter 98 | log-likelihood: -46.730535
## Iter 99 | log-likelihood: -46.730535
## Iter 100 | log-likelihood: -46.730535

print(params_newton)

##           [,1]
## [1,]   2.50892
## [2,] -27.30059
## [3,]  19.12388

se_newton <- standard_errors(params_newton, time, status, treatment)

## Warning in sqrt(diag(H_inv)): NaNs produced

cat("Sai số chuẩn (standard errors) cho phương pháp Newton-Raphson:\n")

## Sai số chuẩn (standard errors) cho phương pháp Newton-Raphson:

print(se_newton)

## [1]      NaN 709.5399 709.5399

# Chạy phương pháp Fisher Scoring
cat("\nƯớc lượng bằng phương pháp Fisher Scoring:\n")

## 
## Ước lượng bằng phương pháp Fisher Scoring:

params_fisher <- fisher_scoring(initial_params, time, status, treatment)

## Iter 1 | log-likelihood: -214.626911
## Iter 2 | log-likelihood: -100.472229
## Iter 3 | log-likelihood: -62.830467
## Iter 4 | log-likelihood: -51.421739
## Iter 5 | log-likelihood: -48.146876
## Iter 6 | log-likelihood: -47.211070
## Iter 7 | log-likelihood: -46.906381
## Iter 8 | log-likelihood: -46.795223
## Iter 9 | log-likelihood: -46.754331
## Iter 10 | log-likelihood: -46.739288
## Iter 11 | log-likelihood: -46.733754
## Iter 12 | log-likelihood: -46.731718
## Iter 13 | log-likelihood: -46.730969
## Iter 14 | log-likelihood: -46.730694
## Iter 15 | log-likelihood: -46.730592
## Iter 16 | log-likelihood: -46.730555
## Iter 17 | log-likelihood: -46.730541
## Iter 18 | log-likelihood: -46.730536
## Iter 19 | log-likelihood: -46.730691
## Iter 20 | log-likelihood: -46.730591
## Iter 21 | log-likelihood: -46.730555
## Iter 22 | log-likelihood: -46.730541
## Iter 23 | log-likelihood: -46.730536
## Iter 24 | log-likelihood: -46.730613
## Iter 25 | log-likelihood: -46.730563
## Iter 26 | log-likelihood: -46.730544
## Iter 27 | log-likelihood: -46.730537
## Iter 28 | log-likelihood: -46.730535
## Iter 29 | log-likelihood: -46.730535
## Iter 30 | log-likelihood: -46.730535
## Iter 31 | log-likelihood: -46.730535
## Iter 32 | log-likelihood: -46.730535
## Iter 33 | log-likelihood: -46.730535
## Iter 34 | log-likelihood: -46.730535
## Iter 35 | log-likelihood: -46.730535
## Iter 36 | log-likelihood: -46.730535
## Iter 37 | log-likelihood: -46.730535
## Iter 38 | log-likelihood: -46.730535
## Iter 39 | log-likelihood: -46.730535
## Iter 40 | log-likelihood: -46.730535
## Iter 41 | log-likelihood: -46.730535
## Iter 42 | log-likelihood: -46.730535
## Iter 43 | log-likelihood: -46.730535
## Iter 44 | log-likelihood: -46.730535
## Iter 45 | log-likelihood: -46.730535
## Iter 46 | log-likelihood: -46.730535
## Iter 47 | log-likelihood: -46.730535
## Iter 48 | log-likelihood: -46.730535
## Iter 49 | log-likelihood: -46.730535
## Iter 50 | log-likelihood: -46.730535
## Iter 51 | log-likelihood: -46.730535
## Iter 52 | log-likelihood: -46.730535
## Iter 53 | log-likelihood: -46.730535
## Iter 54 | log-likelihood: -46.730535
## Iter 55 | log-likelihood: -46.730535
## Iter 56 | log-likelihood: -46.730535
## Iter 57 | log-likelihood: -46.730535
## Iter 58 | log-likelihood: -46.730536
## Iter 59 | log-likelihood: -46.730537
## Iter 60 | log-likelihood: -46.730535
## Iter 61 | log-likelihood: -46.730535
## Iter 62 | log-likelihood: -46.730535
## Iter 63 | log-likelihood: -46.730535
## Iter 64 | log-likelihood: -46.730535
## Iter 65 | log-likelihood: -46.730535
## Iter 66 | log-likelihood: -46.730535
## Iter 67 | log-likelihood: -46.730535
## Iter 68 | log-likelihood: -46.730535
## Iter 69 | log-likelihood: -46.730535
## Iter 70 | log-likelihood: -46.730535
## Iter 71 | log-likelihood: -46.730535
## Iter 72 | log-likelihood: -46.730535
## Iter 73 | log-likelihood: -46.730535
## Iter 74 | log-likelihood: -46.730535
## Iter 75 | log-likelihood: -46.730535
## Iter 76 | log-likelihood: -46.730535
## Iter 77 | log-likelihood: -46.730535
## Iter 78 | log-likelihood: -46.730535
## Iter 79 | log-likelihood: -46.730535
## Iter 80 | log-likelihood: -46.730535
## Iter 81 | log-likelihood: -46.730535
## Iter 82 | log-likelihood: -46.730535
## Iter 83 | log-likelihood: -46.730535
## Iter 84 | log-likelihood: -46.730535
## Iter 85 | log-likelihood: -46.730535
## Iter 86 | log-likelihood: -46.730535
## Iter 87 | log-likelihood: -46.730535
## Iter 88 | log-likelihood: -46.730535
## Iter 89 | log-likelihood: -46.730535
## Iter 90 | log-likelihood: -46.730535
## Iter 91 | log-likelihood: -46.730535
## Iter 92 | log-likelihood: -46.730535
## Iter 93 | log-likelihood: -46.730535
## Iter 94 | log-likelihood: -46.730535
## Iter 95 | log-likelihood: -46.730535
## Iter 96 | log-likelihood: -46.730535
## Iter 97 | log-likelihood: -46.730535
## Iter 98 | log-likelihood: -46.730535
## Iter 99 | log-likelihood: -46.730535
## Iter 100 | log-likelihood: -46.730535

print(params_fisher)

##           [,1]
## [1,]   2.50892
## [2,] -27.30059
## [3,]  19.12388

se_fisher <- standard_errors(params_fisher, time, status, treatment)

## Warning in sqrt(diag(H_inv)): NaNs produced

cat("Sai số chuẩn (standard errors) cho phương pháp Fisher Scoring:\n")

## Sai số chuẩn (standard errors) cho phương pháp Fisher Scoring:

print(se_fisher)

## [1]      NaN 709.5399 709.5399

cat("\nƯớc lượng bằng phương pháp Quasi-Newton:\n")

## 
## Ước lượng bằng phương pháp Quasi-Newton:

print(print(x_vals_bfgs[40,]))

## [1]   2.508915 -27.712745  19.536049
## [1]   2.508915 -27.712745  19.536049

se_quasi_newton <- standard_errors(print(x_vals_bfgs[40,]), time, status, treatment)

## [1]   2.508915 -27.712745  19.536049

## Warning in sqrt(diag(H_inv)): NaNs produced

cat("Sai số chuẩn (standard errors) cho phương pháp Quasi-Newton:\n")

## Sai số chuẩn (standard errors) cho phương pháp Quasi-Newton:

print(se_quasi_newton)

## [1]     NaN 708.715 708.715

(a) Xác định log-likelihood cho các hệ số \(\alpha_1\) và \(\alpha_2\).

(b) Viết các bước giải thuật để xác định ước lượng MLE của \(\alpha_1\) và \(\alpha_2\) bằng phương pháp Newton. Viết một chương trình thi hành giải thuật này.

## (b) Newton-Raphson Algorithm

newton_raphson <- function(alpha0, tol = 1e-8, max_iter = 1000) {
  path <- matrix(alpha0, ncol = 2)
  for (i in 1:max_iter) {
    g <- gradient(path[i, ])
    H <- hessian(path[i, ])
    delta <- solve(H, -g)
    new_point <- path[i, ] + delta
    path <- rbind(path, new_point)
    cat("Iteration:", i, "| alpha: (", new_point, ")\n")
    if (max(abs(delta)) < tol) break
  }
  return(path)
}

nr_path <- newton_raphson(alpha0)

## Iteration: 1 | alpha: ( 1.090831 0.9426757 )
## Iteration: 2 | alpha: ( 1.097128 0.9375748 )
## Iteration: 3 | alpha: ( 1.097153 0.9375546 )
## Iteration: 4 | alpha: ( 1.097153 0.9375546 )

print(nr_path)

##               [,1]      [,2]
##           1.000000 1.0000000
## new_point 1.090831 0.9426757
## new_point 1.097128 0.9375748
## new_point 1.097153 0.9375546
## new_point 1.097153 0.9375546

(c) Viết các bước giải thuật để xác định MLE của \(\alpha_1\) và \(\alpha_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) Fisher Scoring Algorithm

fisher_scoring <- function(alpha0, tol = 1e-6, max_iter = 100) {
  path <- matrix(alpha0, ncol = 2)
  for (i in 1:max_iter) {
    g <- gradient(path[i, ])
    I <- fisher_info(path[i, ])
    if (any(is.na(c(g, I)))) break
    delta <- solve(I, g)
    new_point <- path[i, ] + delta
    path <- rbind(path, new_point)
    cat("Iteration:", i, "| alpha: (", new_point, ")\n")
    if (max(abs(delta)) < tol) break
  }
  return(path)
}

fs_path <- fisher_scoring(alpha0)

## Iteration: 1 | alpha: ( 1.117785 0.9062845 )
## Iteration: 2 | alpha: ( 1.090702 0.9473314 )
## Iteration: 3 | alpha: ( 1.099195 0.9344591 )
## Iteration: 4 | alpha: ( 1.096508 0.9385308 )
## Iteration: 5 | alpha: ( 1.097356 0.9372463 )
## Iteration: 6 | alpha: ( 1.097088 0.9376519 )
## Iteration: 7 | alpha: ( 1.097173 0.9375239 )
## Iteration: 8 | alpha: ( 1.097146 0.9375643 )
## Iteration: 9 | alpha: ( 1.097155 0.9375515 )
## Iteration: 10 | alpha: ( 1.097152 0.9375556 )
## Iteration: 11 | alpha: ( 1.097153 0.9375543 )
## Iteration: 12 | alpha: ( 1.097152 0.9375547 )

print(fs_path)

##               [,1]      [,2]
##           1.000000 1.0000000
## new_point 1.117785 0.9062845
## new_point 1.090702 0.9473314
## new_point 1.099195 0.9344591
## new_point 1.096508 0.9385308
## new_point 1.097356 0.9372463
## new_point 1.097088 0.9376519
## new_point 1.097173 0.9375239
## new_point 1.097146 0.9375643
## new_point 1.097155 0.9375515
## new_point 1.097152 0.9375556
## new_point 1.097153 0.9375543
## new_point 1.097152 0.9375547

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

\[\begin{aligned} \hat{\lambda}_i = \hat{\alpha}_1 b_{i1} + \hat{\alpha}_2 b_{i2} \\ \text{SE}(\hat{\alpha}_j) = \sqrt{\text{diag} (I^{-1}(\hat{\alpha}))_{j}} \\ \text{SE}(\hat{\alpha}_1) = \sqrt{\text{diag} (I^{-1}(\hat{\alpha}))_{1}} \\ \text{SE}(\hat{\alpha}_1) = \sqrt{\left[ I(\hat{\alpha})^{-1} \right]_{11}} \\ \text{SE}(\hat{\alpha}_2) = \sqrt{\left[ I(\hat{\alpha})^{-1} \right]_{22}} \\ \end{aligned}\]

Code R

## (d) Standard Errors of MLE Estimates

alpha_newton <- as.vector(tail(newton_raphson(alpha0), 1))

## Iteration: 1 | alpha: ( 1.090831 0.9426757 )
## Iteration: 2 | alpha: ( 1.097128 0.9375748 )
## Iteration: 3 | alpha: ( 1.097153 0.9375546 )
## Iteration: 4 | alpha: ( 1.097153 0.9375546 )

se_newton <- sqrt(diag(solve(-hessian(alpha_newton))))
alpha_fisher <- as.vector(tail(fisher_scoring(alpha0), 1))

## Iteration: 1 | alpha: ( 1.117785 0.9062845 )
## Iteration: 2 | alpha: ( 1.090702 0.9473314 )
## Iteration: 3 | alpha: ( 1.099195 0.9344591 )
## Iteration: 4 | alpha: ( 1.096508 0.9385308 )
## Iteration: 5 | alpha: ( 1.097356 0.9372463 )
## Iteration: 6 | alpha: ( 1.097088 0.9376519 )
## Iteration: 7 | alpha: ( 1.097173 0.9375239 )
## Iteration: 8 | alpha: ( 1.097146 0.9375643 )
## Iteration: 9 | alpha: ( 1.097155 0.9375515 )
## Iteration: 10 | alpha: ( 1.097152 0.9375556 )
## Iteration: 11 | alpha: ( 1.097153 0.9375543 )
## Iteration: 12 | alpha: ( 1.097152 0.9375547 )

se_fisher <- sqrt(diag(solve(fisher_info(alpha_fisher))))

cat("MLE Estimates (Newton):", round(alpha_newton, 5), "\n")

## MLE Estimates (Newton): 1.09715 0.93755

cat("MLE Estimates (Fisher):", round(as.numeric(alpha_fisher), 5), "\n")

## MLE Estimates (Fisher): 1.09715 0.93755

cat("Standard Errors (Newton):", round(se_newton, 5), "\n")

## Standard Errors (Newton): 0.38961 0.55468

cat("Standard Errors (Fisher):", round(se_fisher, 5), "\n")

## Standard Errors (Fisher): 0.43756 0.63147

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

## (e) Quasi-Newton Method with Two Initial Matrix Choices

quasi_newton <- function(alpha0, M0, tol = 1e-6, max_iter = 100) {
  path <- matrix(alpha0, nrow = 1)
  M <- M0
  g_current <- gradient(alpha0)
  for (i in 1:max_iter) {
    step <- as.vector(-M %*% g_current)
    new_alpha <- path[i, ] + step
    # Ensure alpha > 0 to keep lambda valid
    if (any(new_alpha <= 0)) {
      warning("Negative alpha values: ", new_alpha, ". Step rejected.")
      break
    }
    s <- as.vector(new_alpha - path[i, ])
    g_new <- gradient(new_alpha)
    y <- as.vector(g_new - g_current)

    # Compute rho and update M
    sy <- sum(y * s)
    if (sy == 0) break # Avoid division by zero
    rho <- 1 / sy

    # Update M using outer products
    I <- diag(2)
    term1 <- I - rho * outer(s, y)
    term2 <- I - rho * outer(y, s)
    M <- term1 %*% M %*% term2 + rho * outer(s, s)

    path <- rbind(path, new_alpha)
    g_current <- g_new
    if (max(abs(step)) < tol) break
  }
  path
}

qn_negI_path <- quasi_newton(alpha0, M0 = -diag(2))
qn_negFisher_path <- quasi_newton(alpha0, M0 = -solve(fisher_info(alpha0)))
print(qn_negI_path)

##                [,1]      [,2]
##           1.0000000 1.0000000
## new_alpha 2.1128229 1.3871771
## new_alpha 1.2246525 0.8085567
## new_alpha 0.9207694 1.1448675
## new_alpha 1.1084478 0.9222804
## new_alpha 1.0981072 0.9387168
## new_alpha 1.0963802 0.9371610
## new_alpha 1.0971540 0.9375539
## new_alpha 1.0971525 0.9375546
## new_alpha 1.0971525 0.9375546

print(qn_negFisher_path)

##               [,1]      [,2]
##           1.000000 1.0000000
## new_alpha 1.117785 0.9062845
## new_alpha 1.097394 0.9395147
## new_alpha 1.097107 0.9375701
## new_alpha 1.097153 0.9375540
## new_alpha 1.097153 0.9375546

(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 hoạ tốt nhất các tính năng về hiệu suất của thuật toán.

Code R

## (f) Plot Optimization Paths

library(ggplot2)
set.seed(123)

print(nr_path)

##               [,1]      [,2]
##           1.000000 1.0000000
## new_point 1.090831 0.9426757
## new_point 1.097128 0.9375748
## new_point 1.097153 0.9375546
## new_point 1.097153 0.9375546

print(fs_path)

##               [,1]      [,2]
##           1.000000 1.0000000
## new_point 1.117785 0.9062845
## new_point 1.090702 0.9473314
## new_point 1.099195 0.9344591
## new_point 1.096508 0.9385308
## new_point 1.097356 0.9372463
## new_point 1.097088 0.9376519
## new_point 1.097173 0.9375239
## new_point 1.097146 0.9375643
## new_point 1.097155 0.9375515
## new_point 1.097152 0.9375556
## new_point 1.097153 0.9375543
## new_point 1.097152 0.9375547

print(qn_negI_path)

##                [,1]      [,2]
##           1.0000000 1.0000000
## new_alpha 2.1128229 1.3871771
## new_alpha 1.2246525 0.8085567
## new_alpha 0.9207694 1.1448675
## new_alpha 1.1084478 0.9222804
## new_alpha 1.0981072 0.9387168
## new_alpha 1.0963802 0.9371610
## new_alpha 1.0971540 0.9375539
## new_alpha 1.0971525 0.9375546
## new_alpha 1.0971525 0.9375546

print(qn_negFisher_path)

##               [,1]      [,2]
##           1.000000 1.0000000
## new_alpha 1.117785 0.9062845
## new_alpha 1.097394 0.9395147
## new_alpha 1.097107 0.9375701
## new_alpha 1.097153 0.9375540
## new_alpha 1.097153 0.9375546

# Generate grid coordinates
alpha1_grid <- seq(0.8, 1.8, length.out = 200)
alpha2_grid <- seq(0.7, 1.7, length.out = 200)

# Calculate log-likelihood surface
logL_grid <- outer(
  alpha1_grid, alpha2_grid,
  Vectorize(function(a1, a2) log_likelihood(c(a1, a2)))
)

# Create contour data using base R
contour_data <- expand.grid(alpha1 = alpha1_grid, alpha2 = alpha2_grid)
contour_data$logL <- as.vector(logL_grid)

# Convert path matrices to data frames
nr_df <- data.frame(nr_path, Method = "Newton-Raphson")
fs_df <- data.frame(fs_path, Method = "Fisher-Scoring")
qn_negI_df <- data.frame(qn_negI_path, Method = "Quasi-Newton (-I)")
qn_negFisher_df <- data.frame(qn_negFisher_path, Method = "Quasi-Newton (-Fisher)")

# Combine all paths
all_paths <- rbind(nr_df, fs_df, qn_negI_df, qn_negFisher_df)
colnames(all_paths) <- c("alpha1", "alpha2", "Method")

# Create plot
ggplot(contour_data, aes(x = alpha1, y = alpha2)) +
  geom_contour(aes(z = logL), bins = 30, color = "gray70") +
  geom_path(
    data = all_paths, aes(color = Method, group = Method),
    linewidth = 0.8, alpha = 0.8
  ) +
  geom_point(data = all_paths, aes(color = Method), size = 1.5) +
  geom_point(aes(x = 1.097153, y = 0.9375547),
    color = "red", size = 4, shape = 18
  ) +
  scale_color_manual(values = c(
    "Newton-Raphson" = "#1f77b4",
    "Fisher-Scoring" = "#ff7f0e",
    "Quasi-Newton (-I)" = "#2ca02c",
    "Quasi-Newton (-Fisher)" = "#d62728"
  )) +
  labs(
    title = "Optimization Paths Comparison",
    x = expression(alpha[1]),
    y = expression(alpha[2]),
    subtitle = "Red diamond shows true MLE values"
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    panel.grid = element_blank(),
    panel.border = element_rect(fill = NA, color = "gray80")
  )

## Warning in geom_point(aes(x = 1.097153, y = 0.9375547), color = "red", size = 4, : All aesthetics have length 1, but the data has 40000 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
##   a single row.

a) Xây dựng hàm mục tiêu để ước lượng cho các tham số r và K của mô hình, sao cho sai số của mô hình là nhỏ nhất

Ta có biểu thức cho giá trị của quần thể tại thời điểm \(t\) như sau: \[ N_t = \frac{K\,N_0}{N_0 + \bigl(K - N_0\bigr)\exp(-\,r\,t)} \]

Giả sử ta có \(n\) dữ liệu quan sát \((t_i, N_{t_i})\), trong đó

• \(t_i\) (ngày) thời gian đo

• \(N_{t_i}\) là số lượng bọ cánh cứng thực tế quan sát được tại thời điểm \(t_i\)

Để ước lượng cho các tham số \((r,K)\) sao cho sai số của mô hình là nhỏ nhất, hàm mục tiêu được xây dựng theo tổng bình phương sai số (SSE - Sum of Squared Errors) \[ f(r,K) =\text{SSE}(r, K) = \sum_{i=1}^{n} \biggl[ N_{t_i} - \frac{K\,N_0}{ N_0 + (K - N_0)e^{-r\,t_i} } \biggr]^2 \] Trong đó:

• \(n\) : là số điểm dữ liệu

• \(t_i\) : (ngày) thời gian đo

• \(N_{t_i}\) : là số lượng bọ cánh cứng thực tế quan sát được tại thời điểm \(t_i\)

• \(N_0\) = \(N(0)\) : kích thước quần thể ban đầu

• \(r\) : là tốc độ tăng trưởng nội tại

• \(K\) : là sức chứa của môi trường

b) Viết các bước giải thuật để xác định ước lượng của r và K bằng phương pháp Newton. Viết một chương trình thi hành giải thuật này.

\[ f(r,K) =\text{SSE}(r, K) = \sum_{i=1}^{n} \biggl[ N_{t_i} - \frac{K\,N_0}{ N_0 + (K - N_0)e^{-r\,t_i} } \biggr]^2 = \sum_{i=1}^{n} \biggl[ N_{t_i} - \frac{K\,N_0\,e^{rt_i}}{ \,N_0\,e^{rt_i} + K - N_0 } \biggr]^2 \]

Đạo hàm bậc 1 theo K

\[ \frac{\partial }{\partial K}f(r,K) = 2\sum_{i=1}^{n} \biggl[ N_{t_i} - \frac{K\,N_0\,e^{rt_i}}{ \,N_0\,e^{rt_i} + K - N_0 } \biggr] \frac{\partial }{\partial K} \biggl[ N_{t_i} - \frac{K\,N_0\,e^{rt_i}}{ \,N_0\,e^{rt_i} + K - N_0 } \biggr] \] Với

Suy ra

\[ \frac{\partial}{\partial K}f(r,K) = -2\sum_{i=1}^{n} \biggl[ N_{t_i} - \frac{K\,N_0\,e^{rt_i}}{ \,N_0\,e^{rt_i} + K - N_0 } \biggr] \biggl[\frac{\,N_0^2\,e^{rt_i}(\,e^{rt_i} -1)}{ (\,N_0\,e^{rt_i} + K - N_0)^2 } \biggr] \]

Đạo hàm bậc 1 theo r

\[ \frac{\partial f}{\partial r}(r,K) = 2\sum_{i=1}^{n} \biggl[ N_{t_i} - \frac{K\,N_0\,e^{rt_i}}{ \,N_0\,e^{rt_i} + K - N_0 } \biggr] \frac{\partial }{\partial r} \biggl[ N_{t_i} - \frac{K\,N_0\,e^{rt_i}}{ \,N_0\,e^{rt_i} + K - N_0 } \biggr] \]

Với

Suy ra

\[ \frac{\partial f}{\partial r}(r,K) = -2\sum_{i=1}^{n} \biggl[ N_{t_i} - \frac{K\,N_0\,e^{rt_i}}{ \,N_0\,e^{rt_i} + K - N_0 } \biggr] \biggl[ \frac{K\,N_0\,t_i\,e^{rt_i}(K -N_0)}{ (\,N_0\,e^{rt_i} + K - N_0)^2 } \biggr] \]

Đạo hàm bậc 2 theo K

\[ \frac{\partial^2}{\partial K^2}f(r,K) = -2\sum_{i=1}^{n} \biggl[ \frac{\partial}{\partial K}\biggl[ N_{t_i} - \frac{K\,N_0\,e^{rt_i}}{ \,N_0\,e^{rt_i} + K - N_0 } \biggr] \biggl[\frac{\,N_0^2\,e^{rt_i}(\,e^{rt_i} -1)}{ (\,N_0\,e^{rt_i} + K - N_0)^2 } \biggr] + \biggl[ N_{t_i} - \frac{K\,N_0\,e^{rt_i}}{ \,N_0\,e^{rt_i} + K - N_0 } \biggr] \frac{\partial}{\partial K} \biggl[\frac{\,N_0^2\,e^{rt_i}(\,e^{rt_i} -1)}{ (\,N_0\,e^{rt_i} + K - N_0)^2 } \biggr] \biggr] \] Với

\[ \frac{\partial}{\partial K} \biggl[ N_{t_i} - \frac{K\,N_0\,e^{rt_i}}{ \,N_0\,e^{rt_i} + K - N_0 } \biggr] = - \frac{\,N_0^2\,e^{rt_i}(\,e^{rt_i} -1)}{ (\,N_0\,e^{rt_i} + K - N_0)^2 } \]

\[ \frac{\partial}{\partial K} \biggl[\frac{\,N_0\,e^{2rt_i}(\,e^{rt_i} -1)}{ (\,N_0\,e^{rt_i} + K - N_0)^2 } \biggr] = -\frac{2\,N_0^2\,e^{rt_i}(\,e^{rt_i} -1)}{ (\,N_0\,e^{rt_i} + K - N_0)^3 } \]

Suy ra

Đạo hàm bậc 2 theo r

\[ \frac{\partial^2}{\partial r^2}f(r,K) = -2\sum_{i=1}^{n} \biggl[ \frac{\partial}{\partial r} \biggl[ N_{t_i} - \frac{K\,N_0\,e^{rt_i}}{ \,N_0\,e^{rt_i} + K - N_0 } \biggr] \biggl[ \frac{K\,N_0\,t_i\,e^{rt_i}(K -N_0)}{ (\,N_0\,e^{rt_i} + K - N_0)^2 } \biggr] + \biggl[ N_{t_i} - \frac{K\,N_0\,e^{rt_i}}{ \,N_0\,e^{rt_i} + K - N_0 } \biggr] \frac{\partial}{\partial r} \biggl[ \frac{K\,N_0\,t_i\,e^{rt_i}(K -N_0)}{ (\,N_0\,e^{rt_i} + K - N_0)^2 } \biggr] \biggr] \]

Với

\[ \frac{\partial}{\partial r} \biggl[ N_{t_i} - \frac{K\,N_0\,e^{rt_i}}{ \,N_0\,e^{rt_i} + K - N_0 } \biggr] = -\frac{K\,N_0\,t_i\,e^{rt_i}(K -N_0)}{ (\,N_0\,e^{rt_i} + K - N_0)^2 } \]

Suy ra

Đạo hàm bậc 2 theo r và K

Với \[ \frac{\partial}{\partial K} \biggl[ N_{t_i} - \frac{K\,N_0\,e^{rt_i}}{ \,N_0\,e^{rt_i} + K - N_0 } \biggr] = - \frac{\,N_0^2\,e^{rt_i}(\,e^{rt_i} -1)}{ (\,N_0\,e^{rt_i} + K - N_0)^2 } \]

Suy ra

Với dữ liệu quan sát được tại nhiều thời điểm khác nhau sau 154 ngày như sau

# Chuỗi CSV dưới dạng text
txt <- "days,beetles
0,2
8,47
28,192
41,256
63,768
79,896
97,1120
117,896
135,1184
154,1024
"

# Đọc vào R như là data.frame:
beetles_data <- read.csv(text = txt, header = TRUE)

# Hiển thị bảng:
knitr::kable(beetles_data)

Các bước giải thuật theo phương pháp Newton

Bước 1. Khởi tạo

• Chọn giá trị ban đầu cho vector tham số \(\mathbf{x_0} = (r_0,K_0)\)

• Chọn số vòng lặp tối đa (hoặc ngưỡng dừng).

Bước 2. Tính vector đạo hàm cấp 1 và ma trận đạo hàm cấp 2 (ma trận Hessian) của hàm \(f(r,K)\)

• Ở bước lặp \(k\), lấy \(\mathbf{x} = \mathbf{x}^{(k)}\)

• Tính đạo hàm cấp 1 \(\nabla f(x)\)

• Tính ma trận đạo hàm cấp 2 Hessian

Bước 3. Cập nhật vòng lặp

• Xác định bước Newton \(\Delta\mathbf{x} = \mathbf{H}(\mathbf{x})^{-1} \nabla f(x)\)

• Cập nhật \(\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} - \Delta\mathbf{x}\)

Bước 4. Kiểm tra điều kiện dừng

• Nếu \(\frac{\| \mathbf{x}^{(k+1)} - \mathbf{x}^{(k)}\|}{\|\mathbf{x}^{(k)} + \epsilon \|} \leq \epsilon\) thì dừng.

• Ngược lại, lặp lại bước 2 cho \(k \leftarrow k +1\)

fx_prime <- function(r, K, N0, t, N){
  # ham Nt
  func_Nt <- (K * N0 * exp(r * t)) / (N0 * exp(r * t) + K - N0)
  # mau so
  denom <- (N0 * exp(r * t) + K - N0)
  #dao ham theo K
  df_wrtK <- N0^2 * exp(r * t) * (exp(r * t) - 1) / denom^2
  #dao ham cua SSE theo K
  prime_K <- sum( -2 * (N - func_Nt) * df_wrtK )
  
  #dao ham theo r
  df_wrtr <- K*N0*t * exp(r * t)*(K-N0) / denom^2
  #dao ham cua SSE theo r
  prime_r <- sum( -2 * (N - func_Nt) * df_wrtr )
  res <- c(prime_r,prime_K)
  return(res)
} # tạo đạo hàm cấp 1


f2_prime <- function(r, K, N0, t, N){
  # ham Nt
  func_Nt <- (K * N0 * exp(r * t)) / (N0 * exp(r * t) + K - N0)
  A1 <- N0^2*exp(r*t)*(exp(r * t) - 1)
  A2 <- K*N0*t*exp(r*t)*(K-N0)
  # mau so
  denom <- (N0 * exp(r * t) + K - N0) 
  prime_kk <- sum(2*(A1/denom^3)*(A1/denom + 2*(N-func_Nt)))
  prime_rr <- sum(2* A2/denom^3 * (-t*(K-N0-N0*exp(r*t)) * (N-func_Nt) +A2/denom) )
  prime_rk <- sum(2* (N0^2*t*exp(r*t)/denom^4)* (K*N0*exp(r*t)*(exp(r * t) - 1)*(K - N0) - denom*(2*K*exp(r*t) - denom) * (N-func_Nt)) )
  res <- matrix(c(prime_rr,prime_rk, prime_rk, prime_kk),
                nrow = 2, byrow = TRUE)
  return(res)
} # tạo đạo hàm cấp 2

# Multi Newton 

mul_newton <- function(r, K, N0, t, N,iter,tol = 1e-6){
  
  x <- c(r,K)
  x_vals <- matrix(NA, nrow = iter+1, ncol = 2,
                   dimnames = list(NULL, c("r","K")))
  x_vals[1, ] <- x
  
  for (i in 1:iter){
    r_current <- x[1]
    K_current <- x[2]
    # gradient và Hessian
    grad <- fx_prime(r=r_current, K=K_current, N0, t, N)
    H    <- f2_prime(r=r_current, K=K_current, N0, t, N)
    
    # cập nhật x <- x - H^{-1} * grad
    h <- solve(H, grad)
    
    x_new <- x - h
    
    x_vals[i+1, ] <- x_new
    
    # 4) Tính tỉ lệ sai số tương đối
    #    ratio = ||x_new - x|| / ||x||
    ratio <- norm(x_new - x, type = "2") / (norm(x, type = "2") + tol)
    
    # 5) Kiểm tra điều kiện dừng
    if (ratio < tol) {
      x_vals <- x_vals[1:(i+1), , drop=FALSE]
      break
    }
    x <- x_new
  }
  return(x_vals)
}
t <- c(0,8,28,41,63,79,97,117,135,154)
Num <- c(2,47,192,256,768,896,1120,896,1184,1024)

res <- mul_newton(r=0.1, K=500, N0=2, t=t, N=Num, iter=20)
res

##              r        K
## [1,] 0.1000000  500.000
## [2,] 0.1407367 1092.142
## [3,] 0.1139033 1037.352
## [4,] 0.1173286 1034.866
## [5,] 0.1179416 1033.546
## [6,] 0.1179585 1033.515
## [7,] 0.1179586 1033.515

Kết luận: Ở vòng lặp thứ 7, kết quả là \(r = 0.1179586\) và \(K = 1033.515\)

c) Viết các bước giải thuật để xác định ước lượng của r và K bằng phương pháp quasi-Newton. Viết một chương trình thi hành giải thuật này.

Các bước giải thuật theo phương pháp quasi-Newton

Bước 1. Khởi tạo

• Chọn \(\alpha^{(t)} = 1\) (bước nhảy)

• Chọn ma trận xấp xỉ Hessian \(M^{(t)}\) (ví dụ, lấy \(M^{(0)} = I\))

• Đặt \(\mathbf{x}^{(0)}\) là điểm xuất phát.

Bước 2. Cập nhật \(\mathbf{x}^{(t+1)}\) \[ \mathbf{x}^{(t+1)} \;=\; \mathbf{x}^{(t)} \;-\; \alpha^{(t)}\,\bigl(M^{(t)}\bigr)^{-1}\,\frac{\partial f(\mathbf{x}^{(t)})}{\partial \mathbf{x}} \]

Bước 3. Kiểm tra \(f(\mathbf{x}^{(t+1)})\)
- Nếu \(f(\mathbf{x}^{(t+1)}) < f(\mathbf{x}^{(t)})\) thì tốt, không giảm \(\alpha^{(t)}\).
- Nếu \(f(\mathbf{x}^{(t+1)}) \ge f(\mathbf{x}^{(t)})\) thì đặt \(\alpha^{(t)} := \alpha^{(t)}/2\) rồi tính lại \(\mathbf{x}^{(t+1)}\). Lặp quy trình này cho đến khi \(f(\mathbf{x}^{(t+1)}) - f(\mathbf{x}^{(t)}) > 0\) hoặc \(\alpha^{(t)}\) quá nhỏ.

Bước 4. Cập nhật ma trận \(M^{(t+1)}\) theo công thức BFGS
\[ \begin{aligned} h^{(t)} &= \mathbf{x}^{(t+1)} - \mathbf{x}^{(t)},\\ u^{(t)} &= \frac{\partial f(\mathbf{x}^{(t+1)})}{\partial \mathbf{x}} \;-\; \frac{\partial f(\mathbf{x}^{(t)})}{\partial \mathbf{x}},\\ M^{(t+1)} &= M^{(t)} \;-\; \frac{ M^{(t)} h^{(t)} \,\bigl[M^{(t)} h^{(t)}\bigr]^\top }{ \bigl(h^{(t)}\bigr)^\top\,M^{(t)}\,h^{(t)} } \;+\; \frac{ u^{(t)}\,\bigl[u^{(t)}\bigr]^\top }{ \bigl(h^{(t)}\bigr)^\top\,u^{(t)} }. \end{aligned} \]

Bước 5. Lặp lại cho \(t \leftarrow t+1\) đến khi hội tụ hoặc quá số vòng lặp.

f <- function(r, K, N0, t,N){
  pred <- (K * N0 * exp(r * t)) / (N0 * exp(r * t) + K - N0)
  # Tổng bình phương sai số
  sum((N - pred)^2)
}
quasi_newton <- function(r0, K0, N0, t, N,M,iter,tol = 1e-7) {
  x0 <- c(r0,K0)
  x_val <- matrix(NA,nrow = iter+1, ncol = 2)
  alpha_default <- 1
  alpha <- alpha_default
  x_val[1,] <- x0
  for(i in 1:iter){
    score <- fx_prime(x0[1],x0[2],N0, t, N)
    hessian_inv <- solve(M)
    xt <- x0 - alpha*hessian_inv %*% score
    ## REDUCE ALPHA
    f_xt <- f(xt[1],xt[2],N0, t, N)
    fx <- f(x0[1],x0[2],N0, t, N)
    while(!is.na(f_xt) && !is.na(fx) && f_xt >= fx){
      alpha <- alpha / 2
      xt <- x0 - alpha * (hessian_inv %*% score)
      f_xt <- f(xt[1], xt[2], N0, t, N)
      if(alpha < 1e-12) break
    }
    ##RESET alpha and x0
    x_val[i+1,] <- xt
    ### UPDATE M(t) - BFGS Method
    ht <- xt - x0
    ut <- fx_prime(xt[1],xt[2],N0, t, N) - score
    vt <- ut - M %*% ht
    M_old <- M
    M <- M - (M %*% ht %*% t(M %*% ht))/ (as.numeric(t(ht) %*% M %*% ht)) + (ut%*% t(ut))/(as.numeric(t(ht)%*%ut))
    ### CHECK CONDITION TO UPDATE M
    if(abs(t(vt) %*% vt)[1] <= tol){
      M <- -M_old
    }
    ## RESET alpha and x0
    alpha <- alpha_default
    x0 <- xt
  }
  return(x_val)
}

# Giả sử có dữ liệu t và N, bọ ban đầu N0
t <- c(0,8,28,41,63,79,97,117,135,154)
Num <- c(2,47,192,256,768,896,1120,896,1184,1024)
N0    <- 2

# Gọi thử
r_init <- 0.1
K_init <- 500
M <- diag(c(-1,-1),2,2)
alpha_default <- 1
res <- quasi_newton(r=0.1, K=500, N0=2, t=t, N=Num, 
                    M = M, iter = 40,tol = 1e-7)
res

##             [,1]     [,2]
##  [1,] 0.10000000  500.000
##  [2,] 0.09999021  500.000
##  [3,] 0.09999021  500.000
##  [4,] 0.12215012 1241.266
##  [5,] 0.11882640 1014.959
##  [6,] 0.11962785 1028.823
##  [7,] 0.11962785 1028.823
##  [8,] 0.11907625 1031.612
##  [9,] 0.11883641 1034.580
## [10,] 0.11883641 1034.580
## [11,] 0.11845635 1035.572
## [12,] 0.11845635 1035.572
## [13,] 0.11867225 1033.863
## [14,] 0.11867225 1033.863
## [15,] 0.11867158 1032.815
## [16,] 0.11867158 1032.815
## [17,] 0.11843585 1034.257
## [18,] 0.11843585 1034.257
## [19,] 0.11836987 1032.856
## [20,] 0.11836987 1032.856
## [21,] 0.11812341 1034.349
## [22,] 0.11812341 1034.349
## [23,] 0.11798737 1033.940
## [24,] 0.11789064 1033.983
## [25,] 0.11789064 1033.983
## [26,] 0.11786780 1033.734
## [27,] 0.11786780 1033.734
## [28,] 0.11788314 1033.513
## [29,] 0.11788314 1033.513
## [30,] 0.11790698 1033.488
## [31,] 0.11796578 1033.350
## [32,] 0.11796578 1033.350
## [33,] 0.11797408 1033.455
## [34,] 0.11797408 1033.455
## [35,] 0.11796114 1033.506
## [36,] 0.11795866 1033.517
## [37,] 0.11795866 1033.517
## [38,] 0.11795852 1033.515
## [39,] 0.11795852 1033.515
## [40,] 0.11795856 1033.515
## [41,] 0.11795856 1033.515

Kết luận: Ở vòng lặp thứ 39, kết quả là \(r = 0.1179586\) và \(K = 1033.515\). Kết quả này khớp với kết quả ở câu b) nhưng cần tới 39 vòng lặp mới có thể tìm được tối ưu.

d)Trong nhiều mô hình hóa dân số, người ta áp dụng giả định về log-normality. Giả định đơn giản nhất là \(log\left(N_t\right)\) độc lập và tuân theo phân phối chuẩn với trung bình \(log\{g(t)\}\) và phương sai \(\sigma\). Tìm MLE cho r, K và \(\sigma\) theo giả định này, sử dụng cả 3 phương pháp Newton, Fisher Scoring và quasi-Newton. Cung cấp các sai số chuẩn của ước lượng MLE và ước tính hệ số tương quan giữa chúng. Nhận xét về kết quả.

\[ \log(N_t)\;\sim\;\mathcal{N}\!\Bigl(\log(g(t)),\,\sigma^2\Bigr), \] với \(\displaystyle g(t) = \frac{K\,N_0\,e^{r t}}{N_0\,e^{r t} + (K - N_0)}\).

Khi đó, bạn cần:

Cho các quan sát \(\{(t_i, N_{t_i})\}_{i=1}^n\), giả sử chúng độc lập và

\[ \log(N_{t_i}) \;\sim\; \mathcal{N}\bigl(\,\mu_i,\;\sigma^2\bigr), \quad \text{với } \mu_i = \log\!\bigl(g(t_i;r,K)\bigr). \]

Vậy likelihood (hoặc log-likelihood) cho bộ tham số \((r,\,K,\,\sigma^2)\) là:

\[ \ell(r,K,\sigma^2) \;=\; \sum_{i=1}^n \log\!\Bigl[\, \frac{1}{\sqrt{2\pi\,\sigma^2}} \exp\!\Bigl(\!-\frac{1}{2\sigma^2}\bigl[\log(N_{t_i}) - \log(g(t_i;r,K))\bigr]^2\Bigr) \Bigr]. \]

Rút gọn, ta được:

\[ \ell(r,K,\sigma^2) \;=\; -\,\frac{n}{2}\,\log(2\pi\,\sigma^2) \;-\; \frac{1}{2\sigma^2} \sum_{i=1}^n \bigl[\log(N_{t_i}) - \log(g(t_i))\bigr]^2. \]