hsbkw: Hessian Matrix of the Negative Log-Likelihood for the BKw...

In gkwreg: Generalized Kumaraswamy Regression Models for Bounded Data

Hessian Matrix of the Negative Log-Likelihood for the BKw Distribution

Description

Computes the analytic 4x4 Hessian matrix (matrix of second partial derivatives) of the negative log-likelihood function for the Beta-Kumaraswamy (BKw) distribution with parameters alpha (\alpha), beta (\beta), gamma (\gamma), and delta (\delta). This distribution is the special case of the Generalized Kumaraswamy (GKw) distribution where \lambda = 1. The Hessian is useful for estimating standard errors and in optimization algorithms.

Usage

hsbkw(par, data)

Arguments

par	A numeric vector of length 4 containing the distribution parameters in the order: alpha (\alpha > 0), beta (\beta > 0), 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 based on the BKw log-likelihood (\lambda=1 case of GKw, see llbkw):

\ell(\theta | \mathbf{x}) = n[\ln(\alpha) + \ln(\beta) - \ln B(\gamma, \delta+1)] + \sum_{i=1}^{n} [(\alpha-1)\ln(x_i) + (\beta(\delta+1)-1)\ln(v_i) + (\gamma-1)\ln(w_i)]

where \theta = (\alpha, \beta, \gamma, \delta), B(a,b) is the Beta function (beta), and intermediate terms are:

• v_i = 1 - x_i^{\alpha}

• w_i = 1 - v_i^{\beta} = 1 - (1-x_i^{\alpha})^{\beta}

The Hessian matrix returned contains the elements - \frac{\partial^2 \ell(\theta | \mathbf{x})}{\partial \theta_i \partial \theta_j} for \theta_i, \theta_j \in \{\alpha, \beta, \gamma, \delta\}.

Key properties of the returned matrix:

• Dimensions: 4x4.

• Symmetry: The matrix is symmetric.

• Ordering: Rows and columns correspond to the parameters in the order \alpha, \beta, \gamma, \delta.

• Content: Analytic second derivatives of the negative log-likelihood.

This corresponds to the relevant 4x4 submatrix of the 5x5 GKw Hessian (hsgkw) evaluated at \lambda=1. The exact analytical formulas are implemented directly.

Value

Returns a 4x4 numeric matrix representing the Hessian matrix of the negative log-likelihood function, -\partial^2 \ell / (\partial \theta_i \partial \theta_j), where \theta = (\alpha, \beta, \gamma, \delta). Returns a 4x4 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

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

Kumaraswamy, P. (1980). A generalized probability density function for double-bounded random processes. Journal of Hydrology, 46(1-2), 79-88.

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

See Also

hsgkw (parent distribution Hessian), llbkw (negative log-likelihood for BKw), grbkw (gradient for BKw, if available), dbkw (density for BKw), optim, hessian (for numerical Hessian comparison).

Examples

# Assuming existence of rbkw, llbkw, grbkw, hsbkw functions for BKw

# Generate sample data
set.seed(123)
true_par_bkw <- c(alpha = 2, beta = 3, gamma = 1, delta = 0.5)
if (exists("rbkw")) {
  sample_data_bkw <- rbkw(100, alpha = true_par_bkw[1], beta = true_par_bkw[2],
                         gamma = true_par_bkw[3], delta = true_par_bkw[4])
} else {
  sample_data_bkw <- rgkw(100, alpha = true_par_bkw[1], beta = true_par_bkw[2],
                         gamma = true_par_bkw[3], delta = true_par_bkw[4], lambda = 1)
}
hist(sample_data_bkw, breaks = 20, main = "BKw(2, 3, 1, 0.5) Sample")

# --- Find MLE estimates ---
start_par_bkw <- c(1.5, 2.5, 0.8, 0.3)
mle_result_bkw <- stats::optim(par = start_par_bkw,
                               fn = llbkw,
                               gr = if (exists("grbkw")) grbkw else NULL,
                               method = "BFGS",
                               hessian = TRUE, # Ask optim for numerical Hessian
                               data = sample_data_bkw)

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

  mle_par_bkw <- mle_result_bkw$par
  cat("\nComparing Hessians for BKw at MLE estimates:\n")

  # Numerical Hessian of llbkw
  num_hess_bkw <- numDeriv::hessian(func = llbkw, x = mle_par_bkw, data = sample_data_bkw)

  # Analytical Hessian from hsbkw
  ana_hess_bkw <- hsbkw(par = mle_par_bkw, data = sample_data_bkw)

  cat("Numerical Hessian (BKw):\n")
  print(round(num_hess_bkw, 4))
  cat("Analytical Hessian (BKw):\n")
  print(round(ana_hess_bkw, 4))

  # Check differences
  cat("Max absolute difference between BKw Hessians:\n")
  print(max(abs(num_hess_bkw - ana_hess_bkw)))

  # Optional: Use analytical Hessian for Standard Errors
  # tryCatch({
  #   cov_matrix_bkw <- solve(ana_hess_bkw)
  #   std_errors_bkw <- sqrt(diag(cov_matrix_bkw))
  #   cat("Std. Errors from Analytical BKw Hessian:\n")
  #   print(std_errors_bkw)
  # }, error = function(e) {
  #   warning("Could not invert analytical BKw Hessian: ", e$message)
  # })

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

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