# \donttest{
# Assuming existence of rbeta_, llbeta, grbeta, hsbeta functions
# Generate sample data from a Beta(2, 4) distribution
# (gamma=2, delta=3 in this parameterization)
set.seed(123)
true_par_beta <- c(gamma = 2, delta = 3)
sample_data_beta <- rbeta_(100, gamma = true_par_beta[1], delta = true_par_beta[2])
hist(sample_data_beta, breaks = 20, main = "Beta(2, 4) Sample")
# --- Maximum Likelihood Estimation using optim ---
# Initial parameter guess
start_par_beta <- c(1.5, 2.5)
# Perform optimization (minimizing negative log-likelihood)
# Use method="L-BFGS-B" for box constraints (params > 0 / >= 0)
mle_result_beta <- stats::optim(par = start_par_beta,
fn = llbeta, # Use the custom Beta neg-log-likelihood
method = "L-BFGS-B",
lower = c(1e-6, 1e-6), # Bounds: gamma>0, delta>=0
hessian = TRUE,
data = sample_data_beta)
# Check convergence and results
if (mle_result_beta$convergence == 0) {
print("Optimization converged successfully.")
mle_par_beta <- mle_result_beta$par
print("Estimated Beta parameters (gamma, delta):")
print(mle_par_beta)
print("True Beta parameters (gamma, delta):")
print(true_par_beta)
cat(sprintf("MLE corresponds approx to Beta(%.2f, %.2f)\n",
mle_par_beta[1], mle_par_beta[2] + 1))
} else {
warning("Optimization did not converge!")
print(mle_result_beta$message)
}
# --- Compare numerical and analytical derivatives (if available) ---
# Requires 'numDeriv' package and analytical functions 'grbeta', 'hsbeta'
if (mle_result_beta$convergence == 0 &&
requireNamespace("numDeriv", quietly = TRUE) &&
exists("grbeta") && exists("hsbeta")) {
cat("\nComparing Derivatives at Beta MLE estimates:\n")
# Numerical derivatives of llbeta
num_grad_beta <- numDeriv::grad(func = llbeta, x = mle_par_beta, data = sample_data_beta)
num_hess_beta <- numDeriv::hessian(func = llbeta, x = mle_par_beta, data = sample_data_beta)
# Analytical derivatives (assuming they return derivatives of negative LL)
ana_grad_beta <- grbeta(par = mle_par_beta, data = sample_data_beta)
ana_hess_beta <- hsbeta(par = mle_par_beta, data = sample_data_beta)
# Check differences
cat("Max absolute difference between gradients:\n")
print(max(abs(num_grad_beta - ana_grad_beta)))
cat("Max absolute difference between Hessians:\n")
print(max(abs(num_hess_beta - ana_hess_beta)))
} else {
cat("\nSkipping derivative comparison for Beta.\n")
cat("Requires convergence, 'numDeriv' pkg & functions 'grbeta', 'hsbeta'.\n")
}
# }
Run the code above in your browser using DataLab