grbeta: Gradient of the Negative Log-Likelihood for the Beta Distribution (gamma, delta+1 Parameterization)

Description

Computes the gradient vector (vector of first 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 gradient is useful for optimization algorithms.

Usage

grbeta(par, data)

Value

Returns a numeric vector of length 2 containing the partial derivatives of the negative log-likelihood function $-\ell(\theta | \mathbf{x})$ with respect to each parameter: $(-\partial \ell/\partial \gamma, -\partial \ell/\partial \delta)$. Returns a vector of 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 gradient of the negative log-likelihood for a Beta distribution with parameters shape1 = gamma ($\gamma$) and shape2 = delta + 1 ($\delta+1$). The components of the gradient vector ($-\nabla \ell(\theta | \mathbf{x})$) are:

$$ -\frac{\partial \ell}{\partial \gamma} = n[\psi(\gamma) - \psi(\gamma+\delta+1)] - \sum_{i=1}^{n}\ln(x_i) $$ $$ -\frac{\partial \ell}{\partial \delta} = n[\psi(\delta+1) - \psi(\gamma+\delta+1)] - \sum_{i=1}^{n}\ln(1-x_i) $$

where $\psi(\cdot)$ is the digamma function (digamma). These formulas represent the derivatives of $-\ell(\theta)$, consistent with minimizing the negative log-likelihood. They correspond to the relevant components of the general GKw gradient (grgkw) evaluated at $\alpha=1, \beta=1, \lambda=1$. Note the parameterization: the standard Beta shape parameters are $\gamma$ and $\delta+1$.

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 gradient 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 = grbeta, # Use analytical gradient
                               method = "L-BFGS-B",
                               lower = c(1e-6, 1e-6), # Bounds: gamma>0, delta>=0
                               hessian = TRUE,
                               data = sample_data_beta)

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

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

  # Numerical gradient of llbeta
  num_grad_beta <- numDeriv::grad(func = llbeta, x = mle_par_beta, data = sample_data_beta)

  # Analytical gradient from grbeta
  ana_grad_beta <- grbeta(par = mle_par_beta, data = sample_data_beta)

  cat("Numerical Gradient (Beta):\n")
  print(num_grad_beta)
  cat("Analytical Gradient (Beta):\n")
  print(ana_grad_beta)

  # Check differences
  cat("Max absolute difference between Beta gradients:\n")
  print(max(abs(num_grad_beta - ana_grad_beta)))

} else {
  cat("\nSkipping Beta gradient comparison.\n")
}

# Example with Hessian comparison (if hsbeta exists)
if (mle_result_beta$convergence == 0 &&
    requireNamespace("numDeriv", quietly = TRUE) && exists("hsbeta")) {

  num_hess_beta <- numDeriv::hessian(func = llbeta, x = mle_par_beta, data = sample_data_beta)
  ana_hess_beta <- hsbeta(par = mle_par_beta, data = sample_data_beta)
  cat("\nMax absolute difference between Beta Hessians:\n")
  print(max(abs(num_hess_beta - ana_hess_beta)))

}

# }

Run the code above in your browser using DataLab