Learn R Programming

gkwreg (version 1.0.7)

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

Description

Computes 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. This function is suitable for maximum likelihood estimation.

Usage

llbeta(par, data)

Value

Returns a single double value representing the negative log-likelihood (\(-\ell(\theta|\mathbf{x})\)). Returns Inf

if any parameter values in par 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 negative log-likelihood for a Beta distribution with parameters shape1 = gamma (\(\gamma\)) and shape2 = delta + 1 (\(\delta+1\)). The probability density function (PDF) is: $$ f(x | \gamma, \delta) = \frac{x^{\gamma-1} (1-x)^{\delta}}{B(\gamma, \delta+1)} $$ for \(0 < x < 1\), where \(B(a,b)\) is the Beta function (beta). The log-likelihood function \(\ell(\theta | \mathbf{x})\) for a sample \(\mathbf{x} = (x_1, \dots, x_n)\) is \(\sum_{i=1}^n \ln f(x_i | \theta)\): $$ \ell(\theta | \mathbf{x}) = \sum_{i=1}^{n} [(\gamma-1)\ln(x_i) + \delta\ln(1-x_i)] - n \ln B(\gamma, \delta+1) $$ where \(\theta = (\gamma, \delta)\). This function computes and returns the negative log-likelihood, \(-\ell(\theta|\mathbf{x})\), suitable for minimization using optimization routines like optim. It is equivalent to the negative log-likelihood of the GKw distribution (llgkw) evaluated at \(\alpha=1, \beta=1, \lambda=1\), and also to the negative log-likelihood of the McDonald distribution (llmc) evaluated at \(\lambda=1\). The term \(\ln B(\gamma, \delta+1)\) is typically computed using log-gamma functions (lgamma) for numerical stability.

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,

See Also

llgkw, llmc (related negative log-likelihoods), dbeta_, pbeta_, qbeta_, rbeta_, grbeta (gradient, if available), hsbeta (Hessian, if available), optim, lbeta.

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

# --- 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