# \donttest{
# Assuming existence of rmc, grmc, hsmc functions for Mc distribution
# Generate sample data from a known Mc distribution
set.seed(123)
true_par_mc <- c(gamma = 2, delta = 3, lambda = 0.5)
# Use rmc for data generation
sample_data_mc <- rmc(100, gamma = true_par_mc[1], delta = true_par_mc[2],
lambda = true_par_mc[3])
hist(sample_data_mc, breaks = 20, main = "Mc(2, 3, 0.5) Sample")
# --- Maximum Likelihood Estimation using optim ---
# Initial parameter guess
start_par_mc <- c(1.5, 2.5, 0.8)
# Perform optimization (minimizing negative log-likelihood)
mle_result_mc <- stats::optim(par = start_par_mc,
fn = llmc, # Use the Mc neg-log-likelihood
method = "BFGS", # Or "L-BFGS-B" with lower=1e-6
hessian = TRUE,
data = sample_data_mc)
# Check convergence and results
if (mle_result_mc$convergence == 0) {
print("Optimization converged successfully.")
mle_par_mc <- mle_result_mc$par
print("Estimated Mc parameters:")
print(mle_par_mc)
print("True Mc parameters:")
print(true_par_mc)
} else {
warning("Optimization did not converge!")
print(mle_result_mc$message)
}
# --- Compare numerical and analytical derivatives (if available) ---
# Requires 'numDeriv' package and analytical functions 'grmc', 'hsmc'
if (mle_result_mc$convergence == 0 &&
requireNamespace("numDeriv", quietly = TRUE) &&
exists("grmc") && exists("hsmc")) {
cat("\nComparing Derivatives at Mc MLE estimates:\n")
# Numerical derivatives of llmc
num_grad_mc <- numDeriv::grad(func = llmc, x = mle_par_mc, data = sample_data_mc)
num_hess_mc <- numDeriv::hessian(func = llmc, x = mle_par_mc, data = sample_data_mc)
# Analytical derivatives (assuming they return derivatives of negative LL)
ana_grad_mc <- grmc(par = mle_par_mc, data = sample_data_mc)
ana_hess_mc <- hsmc(par = mle_par_mc, data = sample_data_mc)
# Check differences
cat("Max absolute difference between gradients:\n")
print(max(abs(num_grad_mc - ana_grad_mc)))
cat("Max absolute difference between Hessians:\n")
print(max(abs(num_hess_mc - ana_hess_mc)))
} else {
cat("\nSkipping derivative comparison for Mc.\n")
cat("Requires convergence, 'numDeriv' package and functions 'grmc', 'hsmc'.\n")
}
# }
Run the code above in your browser using DataLab