# \donttest{
# Generate sample data from a known GKw distribution
set.seed(123)
true_par <- c(alpha = 2, beta = 3, gamma = 1.0, delta = 0.5, lambda = 0.5)
sample_data <- rgkw(100, alpha = true_par[1], beta = true_par[2],
gamma = true_par[3], delta = true_par[4], lambda = true_par[5])
hist(sample_data, breaks = 20, main = "Sample GKw Data")
# --- Maximum Likelihood Estimation using optim ---
# Initial parameter guess (can be crucial)
start_par <- c(1.5, 2.5, 1.2, 0.3, 0.6)
# Perform optimization (minimizing negative log-likelihood)
# Ensure data is passed correctly to llgkw
mle_result <- stats::optim(par = start_par,
fn = llgkw,
method = "BFGS", # Method supporting bounds might be safer
hessian = TRUE,
data = sample_data)
# Check convergence and results
if (mle_result$convergence == 0) {
print("Optimization converged successfully.")
mle_par <- mle_result$par
print("Estimated parameters:")
print(mle_par)
print("True parameters:")
print(true_par)
# Standard errors from Hessian (optional)
# fisher_info <- solve(mle_result$hessian) # Need positive definite Hessian
# std_errors <- sqrt(diag(fisher_info))
# print("Approximate Standard Errors:")
# print(std_errors)
} else {
warning("Optimization did not converge!")
print(mle_result$message)
}
# --- Compare numerical and analytical derivatives (if available) ---
# Requires the 'numDeriv' package and analytical functions 'grgkw', 'hsgkw'
if (requireNamespace("numDeriv", quietly = TRUE) &&
exists("grgkw") && exists("hsgkw") && mle_result$convergence == 0) {
cat("\nComparing Derivatives at MLE estimates:\n")
# Numerical derivatives
num_grad <- numDeriv::grad(func = llgkw, x = mle_par, data = sample_data)
num_hess <- numDeriv::hessian(func = llgkw, x = mle_par, data = sample_data)
# Analytical derivatives (assuming they exist)
# Note: grgkw/hsgkw might compute derivatives of log-likelihood,
# while llgkw is negative log-likelihood. Adjust signs if needed.
# Assuming grgkw/hsgkw compute derivatives of NEGATIVE log-likelihood here:
ana_grad <- grgkw(par = mle_par, data = sample_data)
ana_hess <- hsgkw(par = mle_par, data = sample_data)
# Check differences (should be small if analytical functions are correct)
cat("Difference between numerical and analytical gradient:\n")
print(summary(abs(num_grad - ana_grad)))
cat("Difference between numerical and analytical Hessian:\n")
print(summary(abs(num_hess - ana_hess)))
} else {
cat("\nSkipping derivative comparison.\n")
cat("Requires 'numDeriv' package and functions 'grgkw', 'hsgkw'.\n")
}
# }
Run the code above in your browser using DataLab