Dempster, Laird và Rubin (1977) [1] đã xem xét phân phối đa thức với bốn nhóm danh mục, tức là multinomial có hàm xác suất,
\[ p(x_1, x_2, x_3, x_4) = \frac{n!}{x_1! x_2! x_3! x_4!} p_1^{x_1} p_2^{x_2} p_3^{x_3} p_4^{x_4} \]
trong đó \(n = x_1 + x_2 + x_3 + x_4\) và \(p_1 + p_2 + p_3 + p_4 = 1\). Họ cho rằng xác suất có liên quan đến một tham số duy nhất, \(\theta\):
\[ p_1 = \frac{1}{2} + \frac{1}{4}\theta \] \[ p_2 = \frac{1}{4} - \frac{1}{4}\theta \] \[ p_3 = \frac{1}{4} - \frac{1}{4}\theta \] \[ p_4 = \frac{1}{4}\theta, \]
với \(0 \leq \theta \leq 1\). (Mô hình này quay trở lại ví dụ được Fisher thảo luận năm 1925 trong Statistical Methods for Research Workers.) Với một quan sát \((x_1, x_2, x_3, x_4)\), hàm log-likelihood là:
\[ \ell(\theta) = x_1 \log(2 + \theta) + (x_2 + x_3) \log(1 - \theta) + x_4 \log(\theta) + c \]
và
\[ \frac{d\ell(\theta)}{d\theta} = \frac{x_1}{2 + \theta} - \frac{x_2 + x_3}{1 - \theta} + \frac{x_4}{\theta}. \]
Mục tiêu là ước lượng \(\theta\).
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 = (125, 18, 20, 34)\).
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 việc sử dụng phương pháp của Newton là rất hữu ích. Viết các bước giải thuật để xác định ước lượng theo phương pháp của Newton. Viết một chương trình thi hành giải thuật này, bắt đầu với \(\theta^{(0)} = 0.5\).
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^{(0)} = 0.5\).
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^{(0)} = 0.5\).
Công thức Gradient
\[ \frac{d\ell(\theta)}{d\theta} = \frac{x_1}{2 + \theta} - \frac{x_2 + x_3}{1 - \theta} + \frac{x_4}{\theta}. \]
Công thức Hessian
\[ \frac{d^2\ell(\theta)}{d\theta^2} = -\frac{x_1}{(2 + \theta)^2} - \frac{x_2 + x_3}{(1 - \theta)^2} - \frac{x_4}{\theta^2}. \]
# Dữ liệu
x1 <- 125
x2 <- 18
x3 <- 20
x4 <- 34
#Hàm log-likelihood
logL <- function(theta) {
x1 * log(2 + theta) + (x2 + x3) * log(1 - theta) + x4 * log(theta)
}
# Hàm tính gradient
gradient <- function(theta) {
x1 / (2 + theta) - (x2 + x3) / (1 - theta) + x4 / theta
}
# Hàm tính Hessian
hessian <- function(theta) {
-x1 / (2 + theta)^2 - (x2 + x3) / (1 - theta)^2 - x4 / theta^2
}Giải phương trình sử dụng uniroot
theta_est <- uniroot(gradient, lower = 0.01, upper = 0.99)$root
#In kết quả
cat("Ước lượng theta (MLE) bằng phương pháp giải phương trình đa thức:", theta_est)## Ước lượng theta (MLE) bằng phương pháp giải phương trình đa thức: 0.6268217☑ Giải thuật:
Khởi tạo giá trị ban đầu: Chọn giá trị ban đầu cho \(\theta_0\)
Tính Gradient và Hessian của hàm log-likelihood
Cập nhật tham số: Sử dụng phương pháp Newton để cập nhật giá trị của \(\theta\)
Kiểm tra điều kiện dừng: Lặp lại các bước 2 và 3 cho đến khi đạt được điều kiện dừng
☑ Chương trình thi hành:
newton_method <- function(theta0, tol = 1e-8, max_iter = 100) {
theta <- theta0
for (i in 1:max_iter) {
grad <- gradient(theta)
hess <- hessian(theta)
theta_new <- theta - grad / hess
if (abs(theta_new - theta) < tol) {
break
}
theta <- theta_new
}
return(list(theta = theta, iterations = i))
}
# Ước lượng theta
theta_est <- newton_method(0.5)
# In ra kết quả
cat(glue("Ước lượng theta (MLE) bằng phương pháp Newton: {theta_est$theta}, sau {theta_est$iterations} vòng lặp"))## Ước lượng theta (MLE) bằng phương pháp Newton: 0.626821497870984, sau 5 vòng lặp☑ Giải thuật:
Khởi tạo giá trị ban đầu: Chọn giá trị ban đầu cho \(\theta_0\)
Tính gradient và ma trận thông tin Fisher(kỳ vọng của Hessian âm)
Cập nhật tham số: Sử dụng phương pháp Newton để cập nhật giá trị của \(\theta\)
Kiểm tra điều kiện dừng: Lặp lại các bước 2 và 3 cho đến khi đạt được điều kiện dừng.
Ma trận Fisher Information được xác định là:
\[ I(\theta) = -E \left[ \frac{d^2 \ell (\theta)}{d\theta^2} \right] \]
Do \(E[x_i] = np_i\), ta thay kỳ vọng của \(x_1, x_2, x_3\), và \(x_4\):
\[ E[x_1] = np_1 = n \left( \frac{1}{2} + \frac{1}{4} \theta \right) \]
\[ E[x_2 + x_3] = n(p_2 + p_3) = n \left( \frac{1}{4} - \frac{1}{4} \theta + \frac{1}{4} - \frac{1}{4} \theta \right) = n \left( \frac{1}{2} - \frac{1}{2} \theta \right) \]
\[ E[x_4] = np_4 = n \frac{1}{4} \theta \]
Do đó, kỳ vọng của hessian:
\[ E \left[ \frac{d^2 \ell (\theta)}{d\theta^2} \right] = -\frac{n \left( \frac{1}{2} + \frac{1}{4} \theta \right)}{(2 + \theta)^2} - \frac{n \left( \frac{1}{2} - \frac{1}{2} \theta \right)}{(1 - \theta)^2} - \frac{n \frac{1}{4} \theta}{\theta^2} \]
=> Ma trận Fisher:
\[ I(\theta) = \frac{n \left( \frac{1}{2} + \frac{1}{4} \theta \right)}{(2 + \theta)^2} + \frac{n \left( \frac{1}{2} - \frac{1}{2} \theta \right)}{(1 - \theta)^2} + \frac{n \frac{1}{4} \theta}{\theta^2} \]
☑ Chương trình thi hành:
# Hàm ma trận thông tin Fisher
fisher_information <- function(theta, n) {
p1 <- 0.5 + 0.25 * theta
p2 <- 0.25 - 0.25 * theta
p3 <- 0.25 - 0.25 * theta
p4 <- 0.25 * theta
return((n * p1 / (2 + theta)^2) + (n * (p2 + p3) / (1 - theta)^2) + (n * p4 / theta^2))
}
# Phương pháp Fisher Scoring
fisher_scoring <- function(theta0, n = 197, tol = 1e-6, max_iter = 100) {
theta <- theta0
for (i in 1:max_iter) {
grad <- gradient(theta)
info <- fisher_information(theta,n)
theta_new <- theta + grad / info
if (abs(theta_new - theta) < tol) {
break
}
theta <- theta_new
}
return(list(theta = theta, iterations = i))
}
# Ước lượng theta bằng PP Fisher Scoring
theta_est <- fisher_scoring(0.5)
# In ra kết quả
cat(glue("Ước lượng theta (MLE) bằng phương pháp Fisher Scoring: {theta_est$theta}, sau {theta_est$iterations} vòng lặp"))## Ước lượng theta (MLE) bằng phương pháp Fisher Scoring: 0.626820916319874, sau 5 vòng lặp☑ Giải thuật:
Khởi tạo giá trị ban đầu: Chọn giá trị ban đầu cho \(\theta_0\)
Ma trận khởi tạo: Ma trận đơn vị
Xấp xỉ ma trận Hessian: Sử dụng phương pháp quasi-Newton (BFGS) để xấp xỉ Hessian
Cập nhật θ: \[ \theta^{(k+1)} = \theta^{(k)} - H^{(k)} \cdot \frac{d\ell(\theta^{(k)})}{d\theta} \]
Kiểm tra điều kiện dừng: Lặp lại các bước 2 và 3 cho đến khi đạt được điều kiện dừng.
☑ Chương trình thi hành:
# Quasi-Newton method (BFGS)
quasi_newton <- function(theta0 = 0.5, tol = 1e-7, epsilon = 1e-6, max_iter = 1000) {
theta <- theta0
M <- diag(1) # Khởi tạo ma trận Hessian xấp xỉ là ma trận đơn vị
theta_path <- matrix(NA, nrow = max_iter, ncol = length(theta0)) # Lưu trữ các bước của tham số
iter <- 0
for (i in 1:max_iter) {
grad <- gradient(theta) # Tính gradient tại theta hiện tại
grad
# Cập nhật theta mới bằng công thức Quasi-Newton
delta <- solve(M, -grad) # Tính delta
theta_new <- theta + delta # Cập nhật theta mới
# Ràng buộc theta trong khoảng (0, 1)
theta_new <- max(min(theta_new, 1 - epsilon), epsilon)
# Kiểm tra điều kiện dừng (khi sự thay đổi của tham số nhỏ hơn tol)
if (abs(theta_new - theta) < tol) {
break
}
# Tính toán delta_theta và delta_grad
delta_theta <- theta_new - theta
delta_grad <- gradient(theta_new) - grad
# Cập nhật ma trận Hessian xấp xỉ theo công thức BFGS
M <- M + (delta_grad %*% t(delta_grad)) / as.numeric(t(delta_grad) %*% delta_theta) -
(M %*% delta_theta %*% t(delta_theta) %*% M) / as.numeric(t(delta_theta) %*% M %*% delta_theta)
# Cập nhật giá trị tham số
theta <- theta_new
iter <- iter +1
}
# Trả về kết quả với giá trị tham số tối ưu và số lần lặp
return(list(theta = theta, iterations = iter, path = theta_path))
}
# Ước lượng theta bằng PP Fisher Scoring
theta_est <- quasi_newton(0.5)
# In ra kết quả
cat(glue("Ước lượng theta (MLE) bằng phương pháp quasi-Newton: {theta_est$theta}, sau {theta_est$iterations} vòng lặp"))logL_optimization <- function(theta) {
if (theta <= 0 || theta >= 1) {
return(Inf) # Trả về giá trị rất lớn nếu theta nằm ngoài khoảng [0, 1]
}
- (x1 * log(2 + theta) + (x2 + x3) * log(1 - theta) + x4 * log(theta))
}
# Sử dụng phương pháp Quasi - Newton
theta_est <- optim(0.5, logL_optimization, method = "BFGS")$par
cat("Ước lượng theta (MLE) bằng phương pháp Quasi- Newton:", theta_est)## Ước lượng theta (MLE) bằng phương pháp Quasi- Newton: 0.6268212Phương pháp giải phương trình đa thức
Ưu điểm: Chính xác, dễ sử dụng cho phương trình 1 biến.
Nhược điểm: Cần cung cấp khoảng nghiệm phù hợp.
Sử dụng khi: Khi có phương trình bất kỳ cần tìm nghiệm trong khoảng.
Phương pháp Newton:
Ưu điểm: Hội tụ nhanh nếu nghiệm khởi tạo tốt.
Nhược điểm: Cần tính đạo hàm bậc hai, dễ bị sai số.
Sử dụng khi: Khi có thể tính đạo hàm bậc hai dễ dàng.
Phương pháp Fisher Scoring
Ưu điểm: Ổn định hơn Newton, không cần đạo hàm bậc hai trực tiếp.
Nhược điểm: Phụ thuộc vào ma trận Fisher, có thể phức tạp.
Sử dụng khi: Khi cần một phương pháp ổn định hơn Newton.
Phương pháp Quasi - Newton
Ưu điểm: Không cần tính đạo hàm bậc hai, hội tụ nhanh.
Nhược điểm: Có thể chậm hơn Newton trong một số trường hợp.
Sử dụng khi: Khi không muốn tính đạo hàm bậc hai trực tiếp.
Nhận xét chung:
Phương pháp giải phương trình đa thức là lựa chọn tốt nhất vì bài toán này chỉ có một phương trình một ẩn, không cần đạo hàm bậc hai.
Phương pháp Newton hội tụ nhanh nhưng dễ gặp sai số nếu chọn điểm khởi tạo xấu.
Phương pháp Fisher Scoring ổn định hơn Newton nhưng yêu cầu tính toán kỳ vọng hàm Fisher.
Phương pháp Quasi-Newton phù hợp nếu mở rộng bài toán sang nhiều biến nhưng không cần thiết trong bài toán này.
Giả sử thời gian sống sót \(t\) cho các cá thể trong quần thể có hàm mật độ \(f\) và hàm phân phối tích lũy \(F\). Khi đó hàm sống sót (survival function) là \(S(t) = 1 - F(t)\). Hàm nguy cơ (hazard function) là \(h(t) = \frac{f(t)}{1 - F(t)}\), hàm này đo lường rủi ro tử vong tức thời tại thời điểm \(t\) khi sống sót đến thời điểm \(t\). Một mô hình nguy cơ tỷ lệ đặt ra rằng hàm nguy cơ phụ thuộc vào cả thời gian và một vectơ gồm các hiệp biến, \(x\), thông qua mô hình
\[ h(t; x) = \lambda(t) \exp \left( x^T \beta\ \right), \]
trong đó \(\beta\) là vectơ tham số. Đặt \(\Lambda(t) = \int_{-\infty}^t \lambda(u) \, du\), ta chứng minh được rằng
\[ S(t) = \exp \{-\Lambda(t) \exp \left( x^T \beta\ \right)\} \]
và \(f(t) = \lambda(t) \exp \{x^T \beta - \Lambda(t) \exp \left( x^T \beta\ \right)\}\).
\[ \ell(\beta) = \sum_{i=1}^n (w_i \log(\mu_i) - \mu_i) + \sum_{i=1}^n w_i \log \left( \frac{\lambda(t_i)}{\Lambda(t_i)} \right), \]
trong đó, \(\mu_i = \Lambda(t_i) \exp \left( x_i^T \beta\ \right)\).
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 [2]. 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 cho mỗi bệnh nhân đã được ghi lại (xem file
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(t) =
t^\alpha\) với \(\alpha >
0\), khi đó, ta xác định được \(\lambda(t) = \alpha t^{\alpha-1}\) và \(f(t) = \alpha t^{\alpha-1} \exp \{x^T \beta -
t^\alpha \exp \{x^T \beta\}\}\) (hàm mật độ Weibull). Áp dụng
tham số hóa biến phụ thuộc: \(x_i^T \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) 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 bước giải thuật để xác định ước lượng MLE của \(\alpha\), \(\beta_0\) và \(\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.
Viết các bước giải thuật để xác định ước lượng MLE của \(\alpha\), \(\beta_0\) và \(\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 này.
Giả thiết bài toán:
Giả sử dữ liệu là thời gian sống sót đã bị che khuất (censored survival time) \(t_i\) với \(i = 1, \dots, n\).
Đặt \(w_i = 1\) nếu \(t_i\) là thời gian không bị che khuất, \(w_i = 0\) nếu \(t_i\) là thời gian bị che khuất.
Hàm nguy cơ (hazard function) Hàm nguy cơ \(h(t; x)\) được định nghĩa: \[ h(t; x) = \lambda(t) \exp(x^\top \beta) \]
Hàm sống sót (survival function) Hàm sống sót \(S(t)\) mà ta xác suất sống sót đến thời điểm \(t\): \[ S(t) = \exp(-\Lambda(t) \exp(x^\top \beta)) \]
Hàm mật độ (Density Function) Hàm mật độ \(f(t)\) là đạo hàm của hàm phân phối tích lũy \(F(t) = 1 - S(t)\): \[ f(t) = \frac{d}{dt} F(t) = \frac{d}{dt} [1 - S(t)] = -\frac{d}{dt} S(t) \]
\[ f(t) = -\frac{d}{dt} [\exp(-\Lambda(t) \exp(x^\top \beta))] \]
\[ f(t) = \exp(-\Lambda(t) \exp(x^\top \beta)) \cdot (-\Lambda(t) \exp(x^\top \beta)) \]
\[ \Rightarrow f(t) = \lambda(t) \exp \left\{ x^\top \beta - \Lambda(t) \exp(x^\top \beta) \right\} \]
Likelihood cho quan sát không bị che khuất Đối với một quan sát không bị che khuất (\(w_i = 1\)), chúng ta biết chính xác thời gian sống sót \(t_i\). Do đó, likelihood cho quan sát này chính là giá trị của hàm mật độ tại thời điểm \(t_i\): \[ L_i = f(t_i) = \lambda(t_i) \exp(x_i^\top \beta - \Lambda(t_i) \exp(x_i^\top \beta)) \]
Likelihood cho quan sát bị che khuất Một quan sát bị che khuất (\(w_i = 0\)) có nghĩa là chúng ta không biết chính xác thời gian sống sót của đối tượng, nhưng biết rằng nó đã sống sót ít nhất đến thời điểm \(t_i\). Trong trường hợp này, thông tin duy nhất chúng ta có là đối tượng vẫn còn sống tại thời điểm \(t_i\).
\(\Rightarrow\) Do đó, likelihood cho quan sát này chính là xác suất sống sót đến thời điểm \(t_i\), tức là giá trị của hàm sống sót tại thời điểm \(t_i\): \[ L_i = S(t_i) = \exp(-\Lambda(t_i) \exp(x_i^\top \beta)) \]
Hàm log-likelihood tổng hợp cho cả quan sát bị che khuất và không bị che khuất \[ L = \prod_{i=1}^{n} \left[ (f(t_i))^{w_i} S(t_i)^{1-w_i} \right] \]
Log-likelihood: \[ l(\beta) = \sum_{i=1}^{n} \left[ w_i \log f(t_i) + (1-w_i) \log S(t_i) \right] \]
Thay \(f(t_i)\) và \(S(t_i)\) vào công thức trên: \[ l(\beta) = \sum_{i=1}^{n} \left[ w_i (\log \lambda(t_i) + x_i^\top \beta - \Lambda(t_i) \exp(x_i^\top \beta)) + (1-w_i) (-\Lambda(t_i) \exp(x_i^\top \beta)) \right] \]
Thay \(\mu_i = \Lambda(t_i) \exp(x_i^\top \beta)\), ta được Hàm log-likelihood cần chứng minh: \[ l(\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) \]
Hàm mật độ Weibull: \[ f(t; \alpha, \beta) = \alpha t^{\alpha-1} \exp \left\{ x^\top \beta - t^\alpha \exp(x^\top \beta) \right\} \]
Hàm log-likelihood đã được xác định như sau: \[ l(\alpha, \beta_0, \beta_1) = \sum_{i=1}^{n} w_i \left[ \log \alpha + (\alpha - 1) \log t_i + x_i^\top \beta - t_i^\alpha \exp(x_i^\top \beta) \right] \]
Trong đó:
\(x_i^\top \beta = \beta_0 + \delta_i \beta_1\)
\(\delta_i = 1\) nếu bệnh nhân thứ \(i\) nằm trong nhóm điều trị (treatment) và \(\delta_i = 0\) nếu bệnh nhân đó nằm trong nhóm đối chứng (control).
Đạo hàm theo \(\alpha\): \[ \frac{\partial l}{\partial \alpha} = \sum_{i=1}^{n} w_i \left( \frac{1}{\alpha} + \log t_i - t_i^\alpha \log t_i \exp(x_i^\top \beta) \right) \]
Đạo hàm theo \(\beta_0\): \[ \frac{\partial l}{\partial \beta_0} = \sum_{i=1}^{n} w_i \left( 1 - t_i^\alpha \exp(x_i^\top \beta) \right) \]
Đạo hàm theo \(\beta_1\): \[ \frac{\partial l}{\partial \beta_1} = \sum_{i=1}^{n} w_i \delta_i \left( 1 - t_i^\alpha \exp(x_i^\top \beta) \right) \]
Đạo hàm bậc 2
Đạo hàm bậc 2 theo \(\alpha\): \[ \frac{\partial^2 l}{\partial \alpha^2} = \sum_{i=1}^{n} w_i \left( -\frac{1}{\alpha^2} - t_i^\alpha \log t_i^2 \exp(x_i^\top \beta) \right) \]
Đạo hàm bậc 2 theo \(\beta_0\): \[ \frac{\partial^2 l}{\partial \beta_0^2} = - \sum_{i=1}^{n} w_i t_i^\alpha \exp(x_i^\top \beta) \]
Đạo hàm bậc 2 theo \(\beta_1\): \[ \frac{\partial^2 l}{\partial \beta_1^2} = - \sum_{i=1}^{n} w_i \delta_i t_i^\alpha \exp(x_i^\top \beta) \]
Đạo hàm hỗn hợp:
\[ \frac{\partial^2 l}{\partial \alpha \partial \beta_0} = - \sum_{i=1}^{n} w_i t_i^\alpha \log t_i \exp(x_i^\top \beta) \]
\[ \frac{\partial^2 l}{\partial \alpha \partial \beta_1} = - \sum_{i=1}^{n} w_i \delta_i t_i^\alpha \log t_i \exp(x_i^\top \beta) \]
\[ \frac{\partial^2 l}{\partial \beta_0 \partial \beta_1} = - \sum_{i=1}^{n} w_i \delta_i t_i^\alpha \exp(x_i^\top \beta) \]
☑ Giải thuật Newton: 1. Khởi tạo giá trị ban đầu: Chọn giá trị ban đầu cho \(\alpha_0\), \(\beta_0\), \(\beta_1\)
Tính Gradient và Hessian của hàm log-likelihood
Cập nhật tham số: Sử dụng phương pháp Newton để cập nhật giá trị của (_0), \(\beta_0\), \(\beta_1\)
Kiểm tra điều kiện dừng: Lặp lại các bước 2 và 3 cho đến khi đạt được điều kiện dừng
☑ Giải thuật F: 1. Khởi tạo giá trị ban đầu: Chọn giá trị ban đầu cho \(\theta_0\)
Tính gradient và ma trận thông tin Fisher(kỳ vọng của Hessian âm)
Cập nhật tham số: Sử dụng phương pháp Newton để cập nhật giá trị của (_0), \(\beta_0\), \(\beta_1\)
Kiểm tra điều kiện dừng: Lặp lại các bước 2 và 3 cho đến khi đạt được điều kiện dừng.
☑ Chương trình thi hành:
# Đọc dữ liệu
data <- read.csv("mercaptopurine.csv")
# Tạo biến nhị phân cho nhóm điều trị
data$group_treatment <- ifelse(data$group == "treatment", 1, 0)
# Hàm log-likelihood
log_likelihood <- function(params, data) {
alpha <- params[1]
beta0 <- params[2]
beta1 <- params[3]
mu <- data$time^alpha * exp(beta0 + data$group_treatment * beta1)
sum(data$censored * (log(alpha) + (alpha - 1) *ifelse(data$time > 0, log(data$time), 0) + beta0 + data$group_treatment * beta1 - mu))
}
# Hàm gradient
gradient <- function(params, data) {
alpha <- params[1]
beta0 <- params[2]
beta1 <- params[3]
mu <- data$time^alpha * exp(beta0 + data$group_treatment * beta1)
log_time <- ifelse(data$time > 0, log(data$time), 0)
grad_alpha <- sum(data$censored * (1/alpha + log(data$time) - mu * log(data$time)))
grad_beta0 <- sum(data$censored * (1 - mu))
grad_beta1 <- sum(data$censored * data$group_treatment * (1 - mu))
return(c(grad_alpha, grad_beta0, grad_beta1))
}
#Hàm Hessian (với điều chỉnh để tránh ma trận suy biến)
hessian <- function(params, data) {
alpha <- params[1]
beta0 <- params[2]
beta1 <- params[3]
mu <- data$time^alpha * exp(beta0 + data$group_treatment * beta1)
log_time <- ifelse(data$time > 0, log(data$time), 0)
h11 <- -sum(data$censored * (1/alpha^2 + mu * log_time^2))
h12 <- -sum(data$censored * mu * log_time)
h13 <- -sum(data$censored * data$group_treatment * mu * log_time)
h22 <- -sum(data$censored * mu)
h23 <- -sum(data$censored * data$group_treatment * mu)
h33 <- -sum(data$censored * data$group_treatment^2 * mu)
# Thêm epsilon vào đường chéo chính để tránh ma trận suy biến
epsilon <- 1e-6
return(matrix(c(h11, h12, h13, h12, h22, h23, h13, h23, h33), nrow = 3) + diag(epsilon, 3))
}
# Thuật toán Newton-Raphson
newton_raphson <- function(params, data, tol = 1e-6, max_iter = 100) {
for (i in 1:max_iter) {
g <- gradient(params, data)
H <- hessian(params, data)
params_new <- params - solve(H, g)
if (max(abs(params_new - params)) < tol) break
params <- params_new
}
return(params)
}
# Khởi tạo giá trị tham số ban đầu
initial_params <- c(1, 0, 0)
# Áp dụng thuật toán Newton-Raphson
result_new <- newton_raphson(initial_params, data)
cat("Ước lượng tham số alpha, beta0, beta1:\n")
print(result_new)## Ước lượng tham số alpha, beta0, beta1:## [1] 2.217536 -3.505163 -3.505163# Hàm tính ma trận thông tin Fisher
fisher_information_matrix <- function(params, data) {
H <- hessian(params, data)
return(-H)
}
# Tính ma trận thông tin Fisher
I <- fisher_information_matrix(result_new, data)
cat("Ma trận thông tin Fisher:\n")
print(I)
# Thuật toán Fisher Scoring
fisher_scoring <- function(params, data, tol=1e-6, max_iter=100) {
for (i in 1:max_iter) {
g <- gradient(params, data)
H <- hessian(params, data)
# Cập nhật tham số
params_new <- params + solve(I, g)
# Kiểm tra điều kiện hội tụ
if (max(abs(params_new - params)) < tol) break
params <- params_new
}
return(params)
}
# Áp dụng thuật toán Fisher Scoring
result_fs <- fisher_scoring(initial_params, data)
print("Kết quả từ thuật toán Fisher Scoring:")
print(result_fs)## [1] "Kết quả từ thuật toán Fisher Scoring:"## [1] 2.217536 -3.505163 -3.505163# Phương pháp Quasi-Newton BFGS
quasi_newton <- function(init_params, data, tol = 1e-6, max_iter = 100) {
optim_result <- optim(
init_params,
log_likelihood,
data = data,
method = "BFGS",
)
return(optim_result$par)
}
#Khởi tạo giá trị ban đầu cho α, β0 và β1
init_params <- c(1, 0, 0)
# Chạy các phương pháp
result_quasi <- quasi_newton(init_params, data)
cat("Kết quả phương pháp Quasi-Newton:", result_quasi, "\n")## Kết quả phương pháp Quasi-Newton: 149.7742 47.60001 47.60001# Ước lượng sai số chuẩn
standard_errors <- sqrt(diag(solve(I)))
cat("Sai số chuẩn:\n")
cat("alpha:", standard_errors[1], "\n")
cat("beta0:", standard_errors[2], "\n")
cat("beta1:", standard_errors[3], "\n")## Sai số chuẩn:## alpha: 0.5267836## beta0: NaN## beta1: NaN# Tính ma trận tương quan
calculate_correlation_matrix <- function(I) {
covariance_matrix <- solve(I)
correlation_matrix <- covariance_matrix / sqrt(outer(diag(covariance_matrix), diag(covariance_matrix)))
return(correlation_matrix)
}
# Tính ma trận tương quan
calculate_correlation_matrix <- function(I) {
covariance_matrix <- solve(I)
correlation_matrix <- covariance_matrix / sqrt(outer(diag(covariance_matrix), diag(covariance_matrix)))
return(correlation_matrix)
}
# Hiển thị và kiểm tra tương quan mạnh
correlation_matrix <- calculate_correlation_matrix(I)
cat("Ma trận tương quan:\n")
print(correlation_matrix)## Ma trận tương quan:## [,1] [,2] [,3]
## [1,] 1 NaN NaN
## [2,] NaN -1.000000 1.000003
## [3,] NaN 1.000003 -1.000000# Kiểm tra tương quan mạnh (|r| > 0.8)
strong_correlations <- which(abs(correlation_matrix) > 0.8 & abs(correlation_matrix) < 1, arr.ind = TRUE)
if (nrow(strong_correlations) > 0) {
apply(strong_correlations, 1, function(idx) {
cat("Tham số", idx[1], "và", idx[2], "có tương quan mạnh:", correlation_matrix[idx[1], idx[2]], "\n")
})
} else {
cat("Không có tham số nào có tương quan mạnh.\n")
}## Không có tham số nào có tương quan mạnh.Có 46 vụ tràn dầu thỏ với khối lượng ít nhất là 1000 thùng từ tàu chở dầu trên vùng biển Hoa Kỳ trong giai đoạn 1974-1999. Tập tin oilspills.dat chứa dữ liệu sau:
Dữ liệu được điều chỉnh từ [3]. Lượng dầu được vận chuyển được đo bằng tỷ thùng (Bbbl).
Khối lượng dầu được vận chuyển là thước đo mức độ tiếp xúc với rủi ro tràn dầu. Giả sử chúng ta sử dụng giả định quá trình Poisson được đưa ra bởi \(N_i | b_{i1}, b_{i2} \sim P(\lambda_i)\) trong đó \(\lambda_i = \alpha_1 b_{i1} + \alpha_2 b_{i2}\). Các tham số của mô hình này là \(\alpha_1\) và \(\alpha_2\), biểu thị tỷ lệ xảy ra tràn dầu trên mỗi Bbbl dầu được vận chuyển trong quá trình xuất nhập khẩu và trong nước.
Xác định hàm log-likelihood cho các hệ số \(\alpha_1\) và \(\alpha_2\).
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.
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 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 lượng sai số chuẩn (standard errors) của ước lượng MLE.
Áp dụng phương pháp quasi-Newton với hai cách chọn \(M^{(0)}\): (1) ma trận đơn vị âm và (2) ma trận thông tin Fisher, \(-I(\alpha^{(0)})\). So sánh tính ổn định, tốc độ của hai cách chọn này.
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.
Hàm likelihood cho mô hình Poisson là:
\[ L(\alpha_1, \alpha_2) = \prod_{i=1}^{n} \frac{e^{-\lambda_i} \lambda_i^{N_i}}{N_i!} \]
Lấy logarit hàm likelihood ta được Hàm log-likelihood:
\[ \ell(\alpha_1, \alpha_2) = \sum_{i=1}^{n} \left( -\lambda_i + N_i \log(\lambda_i) - \log(N_i!) \right) \]
Thay \(\lambda_i = \alpha_1 b_{i1} + \alpha_2 b_{i2}\) vào, ta có:
\[ \ell(\alpha_1, \alpha_2) = \sum_{i=1}^{n} \left( -(\alpha_1 b_{i1} + \alpha_2 b_{i2}) + N_i \log(\alpha_1 b_{i1} + \alpha_2 b_{i2}) - \log(N_i!) \right) \]
Vì \(\log(N_i!)\) không phụ thuộc vào \(\alpha_1\) và \(\alpha_2\), nó có thể được bỏ qua khi tối ưu hóa hàm log-likelihood. Do đó, hàm log-likelihood rút gọn là:
\[ \ell(\alpha_1, \alpha_2) = \sum_{i=1}^{n} \left( -(\alpha_1 b_{i1} + \alpha_2 b_{i2}) + N_i \log(\alpha_1 b_{i1} + \alpha_2 b_{i2}) \right) \]
☑ Giải thuật:
Khởi tạo giá trị ban đầu: Chọn giá trị ban đầu cho \(\alpha_1\) và \(\alpha_2\).
Tính Gradient và Hessian của hàm log-likelihood
Cập nhật tham số: Sử dụng phương pháp Newton để cập nhật giá trị của \(\alpha_1\) và \(\alpha_2\).
Kiểm tra điều kiện dừng: Lặp lại các bước 2 và 3 cho đến khi đạt được điều kiện dừng
Công thức Gradient:
\[ \nabla \ell(\alpha_1, \alpha_2) = \begin{bmatrix} \frac{\partial \ell}{\partial \alpha_1} \\ \frac{\partial \ell}{\partial \alpha_2} \end{bmatrix} \]
Trong đó:
\[ \frac{\partial \ell}{\partial \alpha_1} = \sum_{i=1}^{n} \left( \frac{N_i}{\lambda_i} b_{i1} - b_{i1} \right) \]
\[ \frac{\partial \ell}{\partial \alpha_2} = \sum_{i=1}^{n} \left( \frac{N_i}{\lambda_i} b_{i2} - b_{i2} \right) \]
Với \(\lambda_i = \alpha_1 b_{i1} + \alpha_2 b_{i2}\).
Công thức Hessian:
\[ H(\alpha_1, \alpha_2) = \begin{bmatrix} \frac{\partial^2 \ell}{\partial \alpha_1^2} & \frac{\partial^2 \ell}{\partial \alpha_1 \partial \alpha_2} \\ \frac{\partial^2 \ell}{\partial \alpha_2 \partial \alpha_1} & \frac{\partial^2 \ell}{\partial \alpha_2^2} \end{bmatrix} \]
Trong đó:
\[ H_{11} = \frac{\partial^2 \ell}{\partial \alpha_1^2} = -\sum_{i=1}^{n} \left( \frac{N_i}{\lambda_i^2} b_{i1}^2 \right) \]
\[ H_{12} = \frac{\partial^2 \ell}{\partial \alpha_1 \partial \alpha_2} = -\sum_{i=1}^{n} \left( \frac{N_i}{\lambda_i^2} b_{i1} b_{i2} \right) \]
\[ H_{21} = \frac{\partial^2 \ell}{\partial \alpha_2 \partial \alpha_1} = H_{12} \]
\[ H_{22}= \frac{\partial^2 \ell}{\partial \alpha_2^2} = -\sum_{i=1}^{n} \left( \frac{N_i}{\lambda_i^2} b_{i2}^2 \right) \]
Với \(\lambda_i = \alpha_1 b_{i1} + \alpha_2 b_{i2}\).
☑ Chương trình thi hành:
# đọc dữ liệu
data_bt3 <- read.table("oilspills.dat", header = TRUE)
# Khởi tạo giá trị ban đầu
alpha0 <- c(0.1, 0.1) # Giá trị ban đầu cho alpha1 và alpha2
tol <- 1e-6 # Ngưỡng hội tụ
max_iter <- 100 # Số lần lặp tối đa
# Hàm log-likelihood
log_likelihood <- function(alpha, data) {
lambda <- alpha[1] * data$importexport + alpha[2] * data$domestic
sum(-lambda + data$spills * log(lambda)) # trả kết quả dưới dạng vector
}
# Gradient của hàm log-likelihood
gradient <- function(alpha, data) {
lambda <- alpha[1] * data$importexport + alpha[2] * data$domestic
c(sum(data$spills / lambda * data$importexport - data$importexport), # trả kết quả dưới dạng vector
sum(data$spills / lambda * data$domestic - data$domestic))
}
# Hessian của hàm log-likelihood
hessian <- function(alpha, data) {
lambda <- alpha[1] * data$importexport + alpha[2] * data$domestic
H11 <- -sum(data$spills / lambda^2 * data$importexport^2)
H12 <- -sum(data$spills / lambda^2 * data$importexport * data$domestic)
H21 <- H12
H22 <- -sum(data$spills / lambda^2 * data$domestic^2)
matrix(c(H11, H12, H21, H22), nrow = 2) # trả kết quả dưới dạng matrix
}
# Phương pháp Newton
newton_method <- function(data, alpha0, tol = 1e-6, max_iter = 100) {
alpha <- alpha0
alpha_path <- matrix(NA, nrow = max_iter, ncol = 2)
for (i in 1:max_iter) {
grad <- gradient(alpha, data)
hess <- hessian(alpha, data)
# Cập nhập alpha
delta <- solve(hess, -grad)
alpha_new <- alpha + delta
alpha_path[i, ] <- alpha_new
# Kiểm tra điều kiện dừng
if (sqrt(sum((alpha_new - alpha)^2)) < tol) {
alpha_path <- alpha_path[1:i,]
break
}
alpha <- alpha_new
}
return(list(alpha = alpha, path = alpha_path))
}
result_newton_method <- newton_method(data_bt3, alpha0, tol, max_iter)
# Kết quả
cat("Ước lượng MLE của alpha1 và alpha2:\n")
cat("alpha1 =", result_newton_method$alpha[1], "\n")
cat("alpha2 =", result_newton_method$alpha[2], "\n")## Ước lượng MLE của alpha1 và alpha2:## Ước lượng alpha1: 1.097153## Ước lượng alpha2: 0.9375546# Cách 2: Sử dụng hàm optim sẵn có để tối đa hóa log-likelihood
# Đặt giá trị khởi tạo cho alpha1 và alpha2
init_params <- c(1, 1)
# Tối đa hóa hàm log-likelihood bằng cách sử dụng hàm optim trong R
result <- optim(init_params, log_likelihood, data = data_bt3, method = "BFGS", control = list(fnscale = -1))
# In ra kết quả ước lượng cho alpha1 và alpha2
alpha1_hat <- result$par[1]
alpha2_hat <- result$par[2]
cat("Ước lượng alpha1:", alpha1_hat, "\n")
cat("Ước lượng alpha2:", alpha2_hat, "\n")## Ước lượng MLE của alpha1 và alpha2 theo cách 2:## Ước lượng alpha1: 1.097153## Ước lượng alpha2: 0.9375544☑ Giải thuật:
Khởi tạo giá trị ban đầu: Chọn giá trị ban đầu cho \(\alpha_1\) và \(\alpha_2\).
Tính gradient và ma trận thông tin Fisher(kỳ vọng của Hessian âm):
Cập nhật tham số: Sử dụng phương pháp Newton để cập nhật giá trị của \(\alpha_1\) và \(\alpha_2\).
Kiểm tra điều kiện dừng: Lặp lại các bước 2 và 3 cho đến khi đạt được điều kiện dừng.
Phương pháp Fisher Scoring sử dụng công thức sau để cập nhật các tham số \(\alpha_1\) và \(\alpha_2\):
\[ θ^{(k+1)} = θ^{(k)} + I(θ^{(k)})^{-1} \nabla L(θ^{(k)}) \]
Trong đó:
\[I(\alpha) = -E[H(\alpha)]\]
Trong mô hình Poisson, vì \(N_i\) có phân phối Poisson với tham số \(\lambda_i\), ta có \(E[N_i] = \lambda_i\). Do đó, ma trận Fisher Information được tính như sau:
\[I_{11} = \sum_{i=1}^{n} \frac{b_{i1}^2}{\lambda_i}, \quad I_{12} = I_{21} = \sum_{i=1}^{n} \frac{b_{i1} b_{i2}}{\lambda_i}, \quad I_{22} = \sum_{i=1}^{n} \frac{b_{i2}^2}{\lambda_i} \]
☑ Chương trình thi hành:
# Khởi tạo giá trị ban đầu
alpha0 <- c(0.1, 0.1) # Giá trị ban đầu cho alpha1 và alpha2
tol <- 1e-6 # Ngưỡng hội tụ
max_iter <- 100 # Số lần lặp tối đa
# Ma trận thông tin Fisher
fisher_information <- function(alpha, data) {
lambda <- alpha[1] * data$importexport + alpha[2] * data$domestic
I11 <- sum(data$importexport^2 / lambda)
I12 <- sum(data$importexport * data$domestic / lambda)
I22 <- sum(data$domestic^2 / lambda)
matrix(c(I11, I12, I12, I22), nrow = 2)
}
# Phương pháp Fisher Scoring
fisher_scoring_method <- function(data, alpha0, tol = 1e-6, max_iter = 100) {
alpha <- alpha0
alpha_path <- matrix(NA, nrow = max_iter, ncol = 2)
for (i in 1:max_iter) {
grad <- gradient(alpha, data)
fisher <- fisher_information(alpha, data)
delta <- solve(fisher, grad)
alpha_new <- alpha + delta
alpha_path[i, ] <- alpha_new
# Kiểm tra điều kiện dừng
if (sqrt(sum((alpha_new - alpha)^2)) < tol) {
alpha_path <- alpha_path[1:i, ]
break
}
alpha <- alpha_new
}
return(list(alpha = alpha, path = alpha_path))
}
result_fisher_scoring_method <- fisher_scoring_method(data_bt3, alpha0, tol, max_iter)
# Kết quả
cat("Ước lượng MLE của alpha1 và alpha2 bằng phương pháp Fisher Scoring:\n")
cat("alpha1 =", result_fisher_scoring_method$alpha[1], "\n")
cat("alpha2 =", result_fisher_scoring_method$alpha[2], "\n")## Ước lượng MLE của alpha1 và alpha2 bằng phương pháp Fisher Scoring:## alpha1 = 1.097153## alpha2 = 0.9375543☑ Giải thuật:
Tính ma trận thông tin Fisher: Sử dụng các ước lượng MLE của \(\alpha_1\) và \(\alpha_2\)
Tính ma trận hiệp phương sai: Ma trận hiệp phương sai của các ước lượng MLE là nghịch đảo của ma trận thông tin Fisher.
Tính sai số chuẩn: Sai số chuẩn là căn bậc hai của các phần tử trên đường chéo chính của ma trận hiệp phương sai.
☑ Chương trình thi hành:
# Giả sử chúng ta đã có ước lượng MLE của alpha1 và alpha2 từ câu 1b hoặc 1c
alpha <- c(1.097153, 0.9375544) # Thay thế bằng ước lượng MLE thực tế
# Tính ma trận thông tin Fisher tại ước lượng MLE
fisher <- fisher_information(alpha, data_bt3)
# Tính ma trận hiệp phương sai
cov_matrix <- solve(fisher)
# Tính sai số chuẩn
standard_errors <- sqrt(diag(cov_matrix))
# Kết quả
cat("Sai số chuẩn của ước lượng MLE:\n")
cat("Standard Error of alpha1 =", standard_errors[1], "\n")
cat("Standard Error of alpha2 =", standard_errors[2], "\n")## Sai số chuẩn của ước lượng MLE:## Standard Error of alpha1 = 0.437556## Standard Error of alpha2 = 0.6314688☑ Giải thuật:
Khởi tạo giá trị ban đầu: Chọn giá trị ban đầu cho \(\alpha_1\) và \(\alpha_2\).
Chọn ma trận khởi tạo: Sử dụng hai cách chọn ma trận khởi tạo \(M^{(0)}\):
Cập nhật tham số: Sử dụng phương pháp quasi-Newton (BFGS) để cập nhật giá trị của \(\alpha_1\) và \(\alpha_2\).
So sánh tính ổn định và tốc độ hội tụ: Theo dõi số lần lặp và sự thay đổi trong giá trị của \(\alpha_1\) và \(\alpha_2\) để so sánh hiệu suất của hai cách chọn ma trận khởi tạo.
☑ Chương trình thi hành:
# Khởi tạo giá trị ban đầu
alpha <- c(0.1, 0.1) # Giá trị ban đầu cho alpha1 và alpha2
tol <- 1e-6 # Ngưỡng hội tụ
max_iter <- 1000 # Số lần lặp tối đa
# Hàm log-likelihood
log_likelihood <- function(alpha, data) {
lambda <- alpha[1] * data$importexport + alpha[2] * data$domestic
sum(-lambda + data$spills * log(lambda))
}
# Gradient của hàm log-likelihood
gradient <- function(alpha, data) {
lambda <- alpha[1] * data$importexport + alpha[2] * data$domestic
c(sum(data$spills / lambda * data$importexport - data$importexport),
sum(data$spills / lambda * data$domestic - data$domestic))
}
# Ma trận thông tin Fisher
fisher_information <- function(alpha, data) {
lambda <- alpha[1] * data$importexport + alpha[2] * data$domestic
I11 <- sum(data$importexport^2 / lambda)
I12 <- sum(data$importexport * data$domestic / lambda)
I21 <- I12
I22 <- sum(data$domestic^2 / lambda)
matrix(c(I11, I12, I21, I22), nrow = 2)
}
# Phương pháp quasi-Newton với ma trận khởi tạo khác nhau
quasi_newton_method <- function(data, alpha0 , M0, tol = 1e-6, max_iter = 100) {
alpha <- alpha0
M <- M0
alpha_path <- matrix(NA, nrow = max_iter, ncol = 2)
for (i in 1:max_iter) {
M <- M0
grad <- gradient(alpha, data)
delta <- solve(M0, -grad)
alpha_new <- alpha + delta
alpha_path[i, ] <- alpha_new
# Kiểm tra điều kiện dừng
if (sqrt(sum((alpha_new - alpha)^2)) < tol) {
alpha_path <- alpha_path[1:i, ]
break
}
# Cập nhật ma trận M sử dụng BFGS
delta_alpha <- alpha_new - alpha
delta_grad <- gradient(alpha_new, data) - grad
M <- M + (delta_grad %*% t(delta_grad)) / as.numeric(t(delta_grad) %*% delta_alpha) -
(M %*% delta_alpha %*% t(delta_alpha) %*% M) / as.numeric(t(delta_alpha) %*% M %*% delta_alpha)
# Cập nhật alpha
alpha <- alpha_new
}
return(list(alpha = alpha, iterations = i, path = alpha_path))
}
# Cách 1: Ma trận đơn vị âm
M0_identity <- -diag(2)
result_identity <- quasi_newton_method(data_bt3, alpha, M0_identity, tol, max_iter)
# Cách 2: Ma trận thông tin Fisher
fisher <- fisher_information(alpha, data_bt3)
M0_fisher <- -fisher
result_fisher <- quasi_newton_method(data_bt3, alpha, M0_fisher, tol, max_iter)
# Kết quả
cat("Kết quả với ma trận đơn vị âm:\n")
cat("alpha1 =", result_identity$alpha[1], "\n")
cat("alpha2 =", result_identity$alpha[2], "\n")
cat("Số lần lặp:", result_identity$iterations, "\n\n")
cat("Kết quả với ma trận thông tin Fisher:\n")
cat("alpha1 =", result_fisher$alpha[1], "\n")
cat("alpha2 =", result_fisher$alpha[2], "\n")
cat("Số lần lặp:", result_fisher$iterations, "\n")## Kết quả với ma trận đơn vị âm:## alpha1 = -26536.16## alpha2 = -17512.67## Số lần lặp: 1000## Kết quả với ma trận thông tin Fisher:## alpha1 = 1.097158## alpha2 = 0.9375495## Số lần lặp: 65Với cách chọn ma trận \(M^{(0)}\) là ma trận thông tin Fisher \(-I(\alpha^{(0)})\) hội tụ sau 65 lần lặp, với ma trận đơn vị âm với 100 lần lặp vẫn chưa hội tụ. => Ma trận Fisher hội tụ nhanh hơn và ổn định hơn ma trận đơn vị
☑ Chương trình thi hành:
# Tạo lưới giá trị cho alpha1 và alpha2
alpha1_seq <- seq(-0.5, 3, length = 100)
alpha2_seq <- seq(-0.5, 3, length = 100)
# Tính toán log-likelihood cho từng giá trị alpha1, alpha2
log_likelihood_grid <- matrix(0, nrow = length(alpha1_seq), ncol = length(alpha2_seq))
for (i in 1:length(alpha1_seq)) {
for (j in 1:length(alpha2_seq)) {
log_likelihood_grid[i, j] <- log_likelihood(c(alpha1_seq[i], alpha2_seq[j]), data_bt3)
}
}
# Lưu đường đi của các phương pháp
result_newton <- newton_method(data_bt3, c(0.1, 0.1))
result_fisher <- fisher_scoring_method(data_bt3, c(0.1, 0.1))
result_quasi_newton <- quasi_newton_method(data_bt3, c(0.1, 0.1), -fisher_information(c(0.1, 0.1), data_bt3))
# Tạo đồ thị contour với plotly
p <- plot_ly(x = alpha1_seq, y = alpha2_seq, z = log_likelihood_grid, type = "contour",
contours = list(start = -200, end = -10, size = 5, coloring = "none", showlabels = TRUE,
labelfont = list(size = 10, color = "black")),
line = list(color = 'black'),
showlegend = TRUE) %>% # Remove the contour trace from the legend
layout(title = "Log-Likelihood",
xaxis = list(title = "alpha1"),
yaxis = list(title = "alpha2"))
# Thêm đường đi vào biểu đồ cho các phương pháp
p <- p %>%
add_trace(x = result_newton$path[, 1], y = result_newton$path[, 2], type = "scatter", mode = "lines+markers",
line = list(color = "red", width = 2), marker = list(color = "red", size = 6), name = "Newton") %>%
add_trace(x = result_fisher$path[, 1], y = result_fisher$path[, 2], type = "scatter", mode = "lines+markers",
line = list(color = "blue", width = 2), marker = list(color = "blue", size = 6), name = "Fisher Scoring") %>%
add_trace(x = result_quasi_newton$path[, 1], y = result_quasi_newton$path[, 2], type = "scatter", mode = "lines+markers",
line = list(color = "green", width = 2), marker = list(color = "green", size = 6), name = "Quasi-Newton")
# Hiển thị đồ thị
p