hsbeta: Hessian Matrix of the Negative Log-Likelihood for the Beta...

In gkwreg: Generalized Kumaraswamy Regression Models for Bounded Data

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)

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).

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.

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).

Author(s)

Lopes, J. E.

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).

See Also

hsgkw, hsmc (related Hessians), llbeta (negative log-likelihood function), grbeta (gradient, if available), dbeta_, pbeta_, qbeta_, rbeta_, optim, hessian (for numerical Hessian comparison), trigamma.

Examples

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

gkwreg documentation built on April 16, 2025, 1:10 a.m.