llbeta: Negative Log-Likelihood for the Beta Distribution (gamma,...

View source: R/RcppExports.R

llbetaR Documentation

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)

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

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

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,

See Also

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

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

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




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