# \donttest{
# Assuming existence of rekw, llekw, grekw, hsekw functions for EKw
# Generate sample data
set.seed(123)
true_par_ekw <- c(alpha = 2, beta = 3, lambda = 0.5)
if (exists("rekw")) {
sample_data_ekw <- rekw(100, alpha = true_par_ekw[1], beta = true_par_ekw[2],
lambda = true_par_ekw[3])
} else {
sample_data_ekw <- rgkw(100, alpha = true_par_ekw[1], beta = true_par_ekw[2],
gamma = 1, delta = 0, lambda = true_par_ekw[3])
}
hist(sample_data_ekw, breaks = 20, main = "EKw(2, 3, 0.5) Sample")
# --- Find MLE estimates ---
start_par_ekw <- c(1.5, 2.5, 0.8)
mle_result_ekw <- stats::optim(par = start_par_ekw,
fn = llekw,
gr = if (exists("grekw")) grekw else NULL,
method = "BFGS",
hessian = TRUE, # Ask optim for numerical Hessian
data = sample_data_ekw)
# --- Compare analytical Hessian to numerical Hessian ---
if (mle_result_ekw$convergence == 0 &&
requireNamespace("numDeriv", quietly = TRUE) &&
exists("hsekw")) {
mle_par_ekw <- mle_result_ekw$par
cat("\nComparing Hessians for EKw at MLE estimates:\n")
# Numerical Hessian of llekw
num_hess_ekw <- numDeriv::hessian(func = llekw, x = mle_par_ekw, data = sample_data_ekw)
# Analytical Hessian from hsekw
ana_hess_ekw <- hsekw(par = mle_par_ekw, data = sample_data_ekw)
cat("Numerical Hessian (EKw):\n")
print(round(num_hess_ekw, 4))
cat("Analytical Hessian (EKw):\n")
print(round(ana_hess_ekw, 4))
# Check differences
cat("Max absolute difference between EKw Hessians:\n")
print(max(abs(num_hess_ekw - ana_hess_ekw)))
# Optional: Use analytical Hessian for Standard Errors
# tryCatch({
# cov_matrix_ekw <- solve(ana_hess_ekw)
# std_errors_ekw <- sqrt(diag(cov_matrix_ekw))
# cat("Std. Errors from Analytical EKw Hessian:\n")
# print(std_errors_ekw)
# }, error = function(e) {
# warning("Could not invert analytical EKw Hessian: ", e$message)
# })
} else {
cat("\nSkipping EKw Hessian comparison.\n")
cat("Requires convergence, 'numDeriv' package, and function 'hsekw'.\n")
}
# }
Run the code above in your browser using DataLab