hsbeta: Hessian Matrix of the Negative Log-Likelihood for the Beta Distribution (gamma, delta+1 Parameterization)

Description

Computes the analytic 2x2 Hessian matrix (matrix of second partial derivatives) of the negative log-likelihood function for the standard Beta distribution, using a parameterization common in generalized distribution families. The distribution is parameterized by gamma ($\gamma$) and delta ($\delta$), corresponding to the standard Beta distribution with shape parameters shape1 = gamma and shape2 = delta + 1. The Hessian is useful for estimating standard errors and in optimization algorithms.

Usage

hsbeta(par, data)

Value

Returns a 2x2 numeric matrix representing the Hessian matrix of the negative log-likelihood function, $-\partial^2 \ell / (\partial \theta_i \partial \theta_j)$, where $\theta = (\gamma, \delta)$. Returns a 2x2 matrix populated with NaN if any parameter values are invalid according to their constraints, or if any value in data is not in the interval (0, 1).

Arguments

par: A numeric vector of length 2 containing the distribution parameters in the order: gamma ($\gamma > 0$), delta ($\delta \ge 0$).
data: A numeric vector of observations. All values must be strictly between 0 and 1 (exclusive).

Author

Lopes, J. E.

Details

This function calculates the analytic second partial derivatives of the negative log-likelihood function ($-\ell(\theta|\mathbf{x})$) for a Beta distribution with parameters shape1 = gamma ($\gamma$) and shape2 = delta + 1 ($\delta+1$). The components of the Hessian matrix ($-\mathbf{H}(\theta)$) are:

$$ -\frac{\partial^2 \ell}{\partial \gamma^2} = n[\psi'(\gamma) - \psi'(\gamma+\delta+1)] $$ $$ -\frac{\partial^2 \ell}{\partial \gamma \partial \delta} = -n\psi'(\gamma+\delta+1) $$ $$ -\frac{\partial^2 \ell}{\partial \delta^2} = n[\psi'(\delta+1) - \psi'(\gamma+\delta+1)] $$

where $\psi'(\cdot)$ is the trigamma function (trigamma). These formulas represent the second derivatives of $-\ell(\theta)$, consistent with minimizing the negative log-likelihood. They correspond to the relevant 2x2 submatrix of the general GKw Hessian (hsgkw) evaluated at $\alpha=1, \beta=1, \lambda=1$. Note the parameterization difference from the standard Beta distribution (shape2 = delta + 1).

The returned matrix is symmetric.

References

Johnson, N. L., Kotz, S., & Balakrishnan, N. (1995). Continuous Univariate Distributions, Volume 2 (2nd ed.). Wiley.

Cordeiro, G. M., & de Castro, M. (2011). A new family of generalized distributions. Journal of Statistical Computation and Simulation,

(Note: Specific Hessian formulas might be derived or sourced from additional references).

Examples

Run this code

# \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")

# --- Find MLE estimates ---
start_par_beta <- c(1.5, 2.5)
mle_result_beta <- stats::optim(par = start_par_beta,
                               fn = llbeta,
                               gr = if (exists("grbeta")) grbeta else NULL,
                               method = "L-BFGS-B",
                               lower = c(1e-6, 1e-6), # Bounds: gamma>0, delta>=0
                               hessian = TRUE, # Ask optim for numerical Hessian
                               data = sample_data_beta)

# --- Compare analytical Hessian to numerical Hessian ---
if (mle_result_beta$convergence == 0 &&
    requireNamespace("numDeriv", quietly = TRUE) &&
    exists("hsbeta")) {

  mle_par_beta <- mle_result_beta$par
  cat("\nComparing Hessians for Beta at MLE estimates:\n")

  # Numerical Hessian of llbeta
  num_hess_beta <- numDeriv::hessian(func = llbeta, x = mle_par_beta, data = sample_data_beta)

  # Analytical Hessian from hsbeta
  ana_hess_beta <- hsbeta(par = mle_par_beta, data = sample_data_beta)

  cat("Numerical Hessian (Beta):\n")
  print(round(num_hess_beta, 4))
  cat("Analytical Hessian (Beta):\n")
  print(round(ana_hess_beta, 4))

  # Check differences
  cat("Max absolute difference between Beta Hessians:\n")
  print(max(abs(num_hess_beta - ana_hess_beta)))

  # Optional: Use analytical Hessian for Standard Errors
  # tryCatch({
  #   cov_matrix_beta <- solve(ana_hess_beta) # ana_hess_beta is already Hessian of negative LL
  #   std_errors_beta <- sqrt(diag(cov_matrix_beta))
  #   cat("Std. Errors from Analytical Beta Hessian:\n")
  #   print(std_errors_beta)
  # }, error = function(e) {
  #   warning("Could not invert analytical Beta Hessian: ", e$message)
  # })

} else {
  cat("\nSkipping Beta Hessian comparison.\n")
  cat("Requires convergence, 'numDeriv' package, and function 'hsbeta'.\n")
}

# }

Run the code above in your browser using DataLab