Dữ liệu sau đây là mẫu i.i.d. (independent, identical distributed – độc lập cùng phân phối) từ phân phối Cauchy(θ, 1): 1.77, -0.23, 2.76, 3.80, 3.47, 56.75, -1.34, 4.24, -2.44, 3.29, 3.71, -2.40, 4.53, -0.07, -1.05, -13.87, -2.53, -1.75, 0.27, 43.21.
Vẽ biểu đồ hàm log-likelihood của θ dựa trên dữ liệu.
Áp dụng thuật toán Newton-Raphson để tìm ước lượng hợp lý cực đại cho θ, với mỗi điểm bắt đầu sau: -11, -1, 0, 1.5, 4, 4.7, 7, 8, và 38. Hãy nhận xét kết quả nghiệm thu được (hội tụ, ổn định). Giá trị trung bình của dữ liệu có phải là điểm khởi đầu tốt không?
Áp dụng phương pháp bisection với điểm bắt đầu là [−1, 1]. Sử dụng các lần chạy bổ sung để minh họa cách thức mà phương pháp bisection có thể không tìm được giá trị cực đại toàn cục.
Từ các giá trị khởi đầu của (θ (0), θ(1)) = (−2, −1), áp dụng phương pháp secant để tìm ước lượng hợp lý cực đại cho θ. Điều gì xảy ra khi (θ (0), θ(1)) = (−3, 3), và đối với các lựa chọn khởi đầu khác?
Sử dụng ví dụ này để so sánh tốc độ và tính ổn định của phương pháp Newton–Raphson, bisection, fixed-point và phương pháp secant. Kết luận của bạn có thay đổi khi bạn áp dụng các phương pháp này cho một mẫu ngẫu nhiên có kích thước 20 từ phân phối chuẩn (θ, 1) không?
Với hàm mật độ xác xuất cauchy f(x;θ,γ) = \[ f(x; \theta, \gamma) = \frac{1}{\pi \gamma \left[1 + \left(\frac{x - \theta}{\gamma}\right)^2 \right]} \]
Lấy lograrit hàm mật độ xác xuất cauchy: \[ \log f(x; \theta, \gamma) = -\log(\pi \gamma) - \log\left(1 + \left(\frac{x - \theta}{\gamma}\right)^2 \right) \] **Khởi tạo hàm tính log-likelihood
# step 1: Create log function
cauchy_data <- c( 1.77, -0.23, 2.76, 3.80, 3.47, 56.75, -1.34, 4.24, -2.44, 3.29, 3.71, -2.40, 4.53, -0.07, -1.05, -13.87, -2.53, -1.75, 0.27, 43.21)
logL_cauchy <- function(data, theta){
sum(dcauchy(x = data, location = theta, scale = 1, log = TRUE))
}
# step 2: sum all vecter theta
logL_cauchy_vec <- Vectorize(FUN = logL_cauchy, vectorize.args = "theta")sum(dcauchy(x = cauchy_data, location = -10, scale = 1, log = TRUE))## [1] -122.5858với dcauchy là hàm lấy log của cauchy: -log(π) - log(1+(x-θ)^2). Thử tạo hàm đối chứng kết quả hàm dcauchy
test_dcauchy <- function(data, theta) {
log_cauchy = -log(pi) - log(1+(data-theta)^2)
return(log_cauchy)
}
cat("Tính tổng log cauchy bằng hàm test_dcauchy trên tập cauchy_data", sum(test_dcauchy(cauchy_data, 2)))## Tính tổng log cauchy bằng hàm test_dcauchy trên tập cauchy_data -74.6098cat("Tính tổng log cauchy bằng hàm logL_cauchy trên tập cauchy_data", logL_cauchy(cauchy_data, 2))## Tính tổng log cauchy bằng hàm logL_cauchy trên tập cauchy_data -74.6098Với θ ∈ (−20, 60), ta tính các giá trị của hàm log-likelihood dựa trên dữ liệu.
# step 3: Create sequence theta from -20 to 60
theta_x <- seq(from = -20, to = 60, by = 0.1)
# step 4: Return vector log_cauchy
logL_cauchy_data <- logL_cauchy_vec(data = cauchy_data, theta = theta_x)Vẽ biểu đồ
# Biểu đồ với Plotly
p <- plot_ly(
x = ~theta_x,
y = ~logL_cauchy_data,
type = 'scatter',
mode = 'lines',
line = list(dash = 'solid') # Tương ứng với lty = 1
) %>%
layout(
xaxis = list(title = "theta"),
yaxis = list(title = "logL Cauchy"),
title = "Log-Likelihood for Cauchy Distribution"
)
# Hiển thị biểu đồ
p- Cách 1
## [1] "Giá trị trung bình là: 5.106"Đạo hàm bậc 1 của log-likelihood: \[ \ell'(\theta) = \sum_{i=1}^n \frac{2(x_i - \theta)}{1 + (x_i - \theta)^2} \]
Đạo hàm bậc 2 của log-likelihood: \[ \ell''(\theta) = \sum_{i=1}^n \frac{-2\left[1 - (x_i - \theta)^2\right]}{\left[1 + (x_i - \theta)^2\right]^2} \]
Khởi tạo hàm tính tổng đạo hàm bậc 1, bậc 2
# f <- function(x) {
# log(x)/(1 + x)
# }
f_prime_log <- function(x, theta){
sum(2*(x-theta)/(1+(x-theta)^2))
}
f_2prime_log <- function(x, theta){
sum(-2*(1-(x-theta)^2)/(1+(x-theta)^2)^2)
}Khi này, thuật toán Newton sẽ tiến hành, với điểm bắt đầu θ(0) = (-11, -1, 0, 1.5, 4, 4.7, 7, 8, và 38), như sau:
# Nhập danh sách khởi đầu theta, cập nhập điểm trung bình
start_theta <- c(-11, -1, 0, 1.5, 4, 4.7, 5.1, 7, 8, 38)
itr <- 20 # vòng lặp tối thiểu 20 lần để dễ quan sát hội tụ
#Tạo biến results dạng list để lưu trữ kết quả chạy thuật toán Newton với mỗi điểm start_theta
results <- list()
# MAIN
for (theta in start_theta) { #Vòng lặp chạy từng điểm start_theta
# Lưu điểm khởi đầu
theta_iter = theta
# Với mỗi điểm start_theta thuật toán Newton tiến hành, với số vòng lặp xác định ở trên
for(i in 1:itr){
# Với mỗi vòng lặp cập nhập lại giá trị theta
theta_iter[i+1] <- theta_iter[i] - f_prime_log(cauchy_data, theta_iter[i] )/f_2prime_log(cauchy_data, theta_iter[i])
}
# Lưu kết quả vào results dạng list
results[[paste0("θ_", theta)]] <- theta_iter
}## # A tibble: 21 × 10
## `θ_-11` `θ_-1` θ_0 θ_1.5 θ_4 θ_4.7 θ_5.1 θ_7 θ_8 θ_38
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -11 -1 0 1.5 4 4.7 5.1 7 8 38
## 2 -16.5 -0.241 -0.196 1.68 2.94 -0.747 17.7 11.8 13.8 43.0
## 3 -25.3 -0.192 -0.192 1.71 2.83 -0.0852 32.6 21.5 25.4 42.8
## 4 -48.4 -0.192 -0.192 1.71 2.82 -0.192 47.8 38.9 44.2 42.8
## 5 -97.3 -0.192 -0.192 1.71 2.82 -0.192 55.6 42.0 -67.6 42.8
## 6 -197. -0.192 -0.192 1.71 2.82 -0.192 54.4 40.6 -136. 42.8
## 7 -397. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.1 -275. 42.8
## 8 -797. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.0 -555. 42.8
## 9 -1599. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.0 -1114. 42.8
## 10 -3203. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.0 -2232. 42.8
## 11 -6410. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.0 -4468. 42.8
## 12 -12826. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.0 -8942. 42.8
## 13 -25657. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.0 -17888. 42.8
## 14 -51319. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.0 -35782. 42.8
## 15 -102643. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.0 -71569. 42.8
## 16 -205291. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.0 -143143. 42.8
## 17 -410587. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.0 -286291. 42.8
## 18 -821179. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.0 -572588. 42.8
## 19 -1642363. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.0 -1145180. 42.8
## 20 -3284731. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.0 -2290365. 42.8
## 21 -6569467. -0.192 -0.192 1.71 2.82 -0.192 54.9 41.0 -4580736. 42.8- Cách 2: Sử dụng vòng lặp dừng while thu gọn kết quả trên. Với mỗi điểm bắt đầu nào có hội tụ tại vòng lặp nào, nghiệm gần chính xác là bao nhiêu
# Nhập danh sách khởi đầu theta
start_theta <- c(-11, -1, 0, 1.5, 4, 4.7, 5.1, 7, 8, 38)
itr <- 200
tolerance <- 1e-6 # Ngưỡng hội tụ
# MAIN
for (theta in start_theta) { #Vòng lặp chạy từng điểm start_theta
# Lưu điểm khởi đầu
theta_iter = theta
i <- 1
# Với mỗi điểm start_theta thuật toán Newton tiến hành, với số vòng lặp xác định ở trên
while(i <= itr){
# Với mỗi vòng lặp cập nhập lại giá trị theta
theta_iter[i+1] <- theta_iter[i] - f_prime_log(cauchy_data, theta_iter[i] )/f_2prime_log(cauchy_data, theta_iter[i])
# Kiểm tra điều kiện hội tụ
if (abs(theta_iter[i+1] - theta_iter[i]) < tolerance) {
break # Dừng vòng lặp nếu hội tụ
}
i <- i + 1
}
is_converge <- ifelse(i <= itr, "Hội tụ", "Không hội tụ")
theta_tail <- ifelse(i <= itr, tail(theta_iter, 1), "Không có")
# In kết quả
combine_string <- paste0("Với điểm theta bắt đầu:", theta, ", có nghiệm gần đúng:", theta_tail, ", Số vòng lặp:", i-1, "\n")
cat(combine_string)
}## Với điểm theta bắt đầu:-11, có nghiệm gần đúng:Không có, Số vòng lặp:200
## Với điểm theta bắt đầu:-1, có nghiệm gần đúng:-0.192286613229646, Số vòng lặp:3
## Với điểm theta bắt đầu:0, có nghiệm gần đúng:-0.192286613229651, Số vòng lặp:3
## Với điểm theta bắt đầu:1.5, có nghiệm gần đúng:1.71358683480959, Số vòng lặp:4
## Với điểm theta bắt đầu:4, có nghiệm gần đúng:2.81747216557324, Số vòng lặp:4
## Với điểm theta bắt đầu:4.7, có nghiệm gần đúng:-0.192286613229649, Số vòng lặp:4
## Với điểm theta bắt đầu:5.1, có nghiệm gần đúng:54.8766179072917, Số vòng lặp:8
## Với điểm theta bắt đầu:7, có nghiệm gần đúng:41.0408478179443, Số vòng lặp:8
## Với điểm theta bắt đầu:8, có nghiệm gần đúng:Không có, Số vòng lặp:200
## Với điểm theta bắt đầu:38, có nghiệm gần đúng:42.7953774720175, Số vòng lặp:4## khoi tao khoang bat dau va diem bat dau
theta_a <- -1
theta_b <- 1
theta_avg <- theta_a + (theta_b - theta_a)/2
## xac dinh so vong lap
itr <- 40
## MAIN
for (i in 1:itr){
if (f_prime_log(cauchy_data,theta_a[i])*f_prime_log(cauchy_data,theta_avg[i]) <= 0) {
# nghiệm nằm trong khoảng [theta_a, theta_avg], chỉ cập nhập lại giá trị theta_b
theta_a[i+1] <- theta_a[i]
theta_b[i + 1] <- theta_avg[i]
} else {
# nghiệm nằm trong khoảng [theta_b, theta_avg], chỉ cập nhập lại giá trị theta_a
theta_a[i + 1] <- theta_avg[i]
theta_b[i + 1] <- theta_b[i]
}
theta_avg[i+1] <- theta_a[i+1] + (theta_b[i+1] - theta_a[i+1])/2
}
head(cbind(theta_a, theta_b, theta_avg), 40)## theta_a theta_b theta_avg
## [1,] -1.0000000 1.0000000 0.0000000
## [2,] -1.0000000 0.0000000 -0.5000000
## [3,] -0.5000000 0.0000000 -0.2500000
## [4,] -0.2500000 0.0000000 -0.1250000
## [5,] -0.2500000 -0.1250000 -0.1875000
## [6,] -0.2500000 -0.1875000 -0.2187500
## [7,] -0.2187500 -0.1875000 -0.2031250
## [8,] -0.2031250 -0.1875000 -0.1953125
## [9,] -0.1953125 -0.1875000 -0.1914062
## [10,] -0.1953125 -0.1914062 -0.1933594
## [11,] -0.1933594 -0.1914062 -0.1923828
## [12,] -0.1923828 -0.1914062 -0.1918945
## [13,] -0.1923828 -0.1918945 -0.1921387
## [14,] -0.1923828 -0.1921387 -0.1922607
## [15,] -0.1923828 -0.1922607 -0.1923218
## [16,] -0.1923218 -0.1922607 -0.1922913
## [17,] -0.1922913 -0.1922607 -0.1922760
## [18,] -0.1922913 -0.1922760 -0.1922836
## [19,] -0.1922913 -0.1922836 -0.1922874
## [20,] -0.1922874 -0.1922836 -0.1922855
## [21,] -0.1922874 -0.1922855 -0.1922865
## [22,] -0.1922874 -0.1922865 -0.1922870
## [23,] -0.1922870 -0.1922865 -0.1922867
## [24,] -0.1922867 -0.1922865 -0.1922866
## [25,] -0.1922867 -0.1922866 -0.1922867
## [26,] -0.1922867 -0.1922866 -0.1922866
## [27,] -0.1922866 -0.1922866 -0.1922866
## [28,] -0.1922866 -0.1922866 -0.1922866
## [29,] -0.1922866 -0.1922866 -0.1922866
## [30,] -0.1922866 -0.1922866 -0.1922866
## [31,] -0.1922866 -0.1922866 -0.1922866
## [32,] -0.1922866 -0.1922866 -0.1922866
## [33,] -0.1922866 -0.1922866 -0.1922866
## [34,] -0.1922866 -0.1922866 -0.1922866
## [35,] -0.1922866 -0.1922866 -0.1922866
## [36,] -0.1922866 -0.1922866 -0.1922866
## [37,] -0.1922866 -0.1922866 -0.1922866
## [38,] -0.1922866 -0.1922866 -0.1922866
## [39,] -0.1922866 -0.1922866 -0.1922866
## [40,] -0.1922866 -0.1922866 -0.1922866Công thức lặp nghiệm \[ x^{(t+1)} = x^{(t)} - f'(x^{(t)}) \cdot \frac{x^{(t)} - x^{(t-1)}}{f'(x^{(t)}) - f'(x^{(t-1)})} \]
theta_0 <- -2
theta_1 <- -1
itr <- 20
theta <- theta_1
## MAIN
for(i in 1:itr){
# Cập nhập lại theta sau khi mỗi lần lặp nghiệm
theta[i + 1] <- theta[i] - f_prime_log(cauchy_data,theta[i])*((theta[i]-theta_0[i])/(f_prime_log(cauchy_data,theta[i]) - f_prime_log(cauchy_data,theta_0[i])))
# cập nhập lại điểm khởi đầu
theta_0[i+1] = theta[i]
}
print(theta)## [1] -1.0000000 -0.4412058 -0.1245642 -0.1984022 -0.1923655 -0.1922865
## [7] -0.1922866 -0.1922866 -0.1922866 -0.1922866 -0.1922866 -0.1922866
## [13] NaN NaN NaN NaN NaN NaN
## [19] NaN NaN NaNSử dụng vòng lặp điều kiện để dừng lại khi hàm hội tụ
theta_0 <- -2
theta_1 <- -1
theta <- theta_1
# Ngưỡng hội tụ
tol <- 1e-6
# Vòng lặp secant
i <- 1
while (abs(theta[i] - theta_0[i]) >= tol) {
theta[i + 1] <- theta[i] - f_prime_log(cauchy_data,theta[i])*((theta[i]-theta_0[i])/(f_prime_log(cauchy_data,theta[i]) - f_prime_log(cauchy_data,theta_0[i])))
# cập nhập lại điểm khởi đầu
theta_0[i+1] = theta[i]
i <- i + 1
}
root <- tail(theta, 1)
# In kết quả
cat("Nghiệm gần đúng:", root, "\n")## Nghiệm gần đúng: -0.1922866cat("Số vòng lặp:", i-1, "\n")## Số vòng lặp: 6# Đoạn mã 1
time1 <- system.time({
result1 <- sum(runif(1e6))
})
# Đoạn mã 2
time2 <- system.time({
result2 <- prod(1:1000)
})
# In kết quả
cat("Thời gian chạy đoạn mã 1:", time1["elapsed"], "giây\n")## Thời gian chạy đoạn mã 1: 0.01 giâycat("Thời gian chạy đoạn mã 2:", time2["elapsed"], "giây\n")## Thời gian chạy đoạn mã 2: 0 giâySo sánh các phương pháp
Tiêu chí so sánh:
Tốc độ hội tụ: Số lần lặp để đạt giá trị θ gần cực đại.
Tính ổn định: Khả năng hội tụ đúng điểm cực đại từ các điểm khởi đầu khác nhau.
So sánh:
Xét mật độ \(f(x) = \frac{1 - \cos(x - \theta)}{2\pi}\) trên 0 ≤ x ≤ 2π, trong đó θ là tham số và θ ∈ [−π, π]. Dữ liệu i.i.d. sau phát sinh từ mật độ này: 3.91, 4.85, 2.28, 4.06, 3.70, 4.04, 5.46, 3.53, 2.28, 1.96, 2.53, 3.88, 2.22, 3.47, 4.82, 2.46, 2.99, 2.54, 0.52, 2.50. Chúng ta muốn ước tính θ.
Vẽ đồ thị hàm log-likelihood với θ ∈ [−π, π].
Tìm ước lượng method-of-moment cho θ.
Tìm MLE cho θ bằng phương pháp Newton–Raphson, sử dụng kết quả từ (b) làm giá trị bắt đầu. Bạn tìm thấy những nghiệm nào khi bắt đầu ở -2.7 và 2.7?
Lặp lại phần (c) bằng cách sử dụng 200 giá trị bắt đầu cách đều nhau giữa −π và π. Nói cách khác, chia tập hợp các giá trị bắt đầu thành các nhóm riêng biệt, với mỗi nhóm tương ứng với một kết quả duy nhất riêng biệt của quá trình tối ưu hóa (một chế độ cục bộ). Thảo luận về kết quả của bạn.
Tìm hai giá trị bắt đầu, gần bằng nhau nhất có thể, mà phương pháp Newton–Raphson hội tụ thành hai nghiệm khác nhau.