Computes the gradient vector (vector of partial derivatives) of the negative log-likelihood function for the five-parameter Generalized Kumaraswamy (GKw) distribution. This provides the analytical gradient, often used for efficient optimization via maximum likelihood estimation.
Arguments
- par
A numeric vector of length 5 containing the distribution parameters in the order:
alpha(\(\alpha > 0\)),beta(\(\beta > 0\)),gamma(\(\gamma > 0\)),delta(\(\delta \ge 0\)),lambda(\(\lambda > 0\)).- data
A numeric vector of observations. All values must be strictly between 0 and 1 (exclusive).
Value
Returns a numeric vector of length 5 containing the partial derivatives
of the negative log-likelihood function \(-\ell(\theta | \mathbf{x})\) with
respect to each parameter:
\((-\partial \ell/\partial \alpha, -\partial \ell/\partial \beta, -\partial \ell/\partial \gamma, -\partial \ell/\partial \delta, -\partial \ell/\partial \lambda)\).
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).
Details
The components of the gradient vector of the negative log-likelihood (\(-\nabla \ell(\theta | \mathbf{x})\)) are:
$$ -\frac{\partial \ell}{\partial \alpha} = -\frac{n}{\alpha} - \sum_{i=1}^{n}\ln(x_i) + \sum_{i=1}^{n}\left[x_i^{\alpha} \ln(x_i) \left(\frac{\beta-1}{v_i} - \frac{(\gamma\lambda-1) \beta v_i^{\beta-1}}{w_i} + \frac{\delta \lambda \beta v_i^{\beta-1} w_i^{\lambda-1}}{z_i}\right)\right] $$ $$ -\frac{\partial \ell}{\partial \beta} = -\frac{n}{\beta} - \sum_{i=1}^{n}\ln(v_i) + \sum_{i=1}^{n}\left[v_i^{\beta} \ln(v_i) \left(\frac{\gamma\lambda-1}{w_i} - \frac{\delta \lambda w_i^{\lambda-1}}{z_i}\right)\right] $$ $$ -\frac{\partial \ell}{\partial \gamma} = n[\psi(\gamma) - \psi(\gamma+\delta+1)] - \lambda\sum_{i=1}^{n}\ln(w_i) $$ $$ -\frac{\partial \ell}{\partial \delta} = n[\psi(\delta+1) - \psi(\gamma+\delta+1)] - \sum_{i=1}^{n}\ln(z_i) $$ $$ -\frac{\partial \ell}{\partial \lambda} = -\frac{n}{\lambda} - \gamma\sum_{i=1}^{n}\ln(w_i) + \delta\sum_{i=1}^{n}\frac{w_i^{\lambda}\ln(w_i)}{z_i} $$
where:
\(v_i = 1 - x_i^{\alpha}\)
\(w_i = 1 - v_i^{\beta} = 1 - (1-x_i^{\alpha})^{\beta}\)
\(z_i = 1 - w_i^{\lambda} = 1 - [1-(1-x_i^{\alpha})^{\beta}]^{\lambda}\)
\(\psi(\cdot)\) is the digamma function (
digamma).
Numerical stability is ensured through careful implementation, including checks for valid inputs and handling of intermediate calculations involving potentially small or large numbers, often leveraging the Armadillo C++ library for efficiency.
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.
Examples
# \donttest{
## Example 1: Basic Gradient Evaluation
# Generate sample data
set.seed(123)
n <- 1000
true_params <- c(alpha = 2.0, beta = 3.0, gamma = 1.5, delta = 2.0, lambda = 1.8)
data <- rgkw(n, alpha = true_params[1], beta = true_params[2],
gamma = true_params[3], delta = true_params[4],
lambda = true_params[5])
# Evaluate gradient at true parameters
grad_true <- grgkw(par = true_params, data = data)
cat("Gradient at true parameters:\n")
#> Gradient at true parameters:
print(grad_true)
#> [1] -34.386342 12.010575 -19.736267 7.392701 -22.078415
cat("Norm:", sqrt(sum(grad_true^2)), "\n")
#> Norm: 47.52161
# Evaluate at different parameter values
test_params <- rbind(
c(1.5, 2.5, 1.2, 1.5, 1.5),
c(2.0, 3.0, 1.5, 2.0, 1.8),
c(2.5, 3.5, 1.8, 2.5, 2.0)
)
grad_norms <- apply(test_params, 1, function(p) {
g <- grgkw(p, data)
sqrt(sum(g^2))
})
results <- data.frame(
Alpha = test_params[, 1],
Beta = test_params[, 2],
Gamma = test_params[, 3],
Delta = test_params[, 4],
Lambda = test_params[, 5],
Grad_Norm = grad_norms
)
print(results, digits = 4)
#> Alpha Beta Gamma Delta Lambda Grad_Norm
#> 1 1.5 2.5 1.2 1.5 1.5 1504.78
#> 2 2.0 3.0 1.5 2.0 1.8 47.52
#> 3 2.5 3.5 1.8 2.5 2.0 1402.07
## Example 2: Gradient in Optimization
# Optimization with analytical gradient
fit_with_grad <- optim(
par = c(1.5, 2.5, 1.2, 1.5, 1.5),
fn = llgkw,
gr = grgkw,
data = data,
method = "BFGS",
hessian = TRUE,
control = list(trace = 0, maxit = 1000)
)
# Optimization without gradient
fit_no_grad <- optim(
par = c(1.5, 2.5, 1.2, 1.5, 1.5),
fn = llgkw,
data = data,
method = "BFGS",
hessian = TRUE,
control = list(trace = 0, maxit = 1000)
)
comparison <- data.frame(
Method = c("With Gradient", "Without Gradient"),
Alpha = c(fit_with_grad$par[1], fit_no_grad$par[1]),
Beta = c(fit_with_grad$par[2], fit_no_grad$par[2]),
Gamma = c(fit_with_grad$par[3], fit_no_grad$par[3]),
Delta = c(fit_with_grad$par[4], fit_no_grad$par[4]),
Lambda = c(fit_with_grad$par[5], fit_no_grad$par[5]),
NegLogLik = c(fit_with_grad$value, fit_no_grad$value),
Iterations = c(fit_with_grad$counts[1], fit_no_grad$counts[1])
)
print(comparison, digits = 4, row.names = FALSE)
#> Method Alpha Beta Gamma Delta Lambda NegLogLik Iterations
#> With Gradient 1.205 3.288 0.3823 1.462 13.88 -704.3 386
#> Without Gradient 1.256 3.386 0.3747 1.403 13.43 -704.3 367
## Example 3: Verifying Gradient at MLE
mle <- fit_with_grad$par
names(mle) <- c("alpha", "beta", "gamma", "delta", "lambda")
# At MLE, gradient should be approximately zero
gradient_at_mle <- grgkw(par = mle, data = data)
cat("\nGradient at MLE:\n")
#>
#> Gradient at MLE:
print(gradient_at_mle)
#> [1] -0.040704813 0.024508274 0.007276465 0.017706753 -0.010821899
cat("Max absolute component:", max(abs(gradient_at_mle)), "\n")
#> Max absolute component: 0.04070481
cat("Gradient norm:", sqrt(sum(gradient_at_mle^2)), "\n")
#> Gradient norm: 0.05235577
## Example 4: Numerical vs Analytical Gradient
# Manual finite difference gradient
numerical_gradient <- function(f, x, data, h = 1e-7) {
grad <- numeric(length(x))
for (i in seq_along(x)) {
x_plus <- x_minus <- x
x_plus[i] <- x[i] + h
x_minus[i] <- x[i] - h
grad[i] <- (f(x_plus, data) - f(x_minus, data)) / (2 * h)
}
return(grad)
}
# Compare at MLE
grad_analytical <- grgkw(par = mle, data = data)
grad_numerical <- numerical_gradient(llgkw, mle, data)
comparison_grad <- data.frame(
Parameter = c("alpha", "beta", "gamma", "delta", "lambda"),
Analytical = grad_analytical,
Numerical = grad_numerical,
Abs_Diff = abs(grad_analytical - grad_numerical),
Rel_Error = abs(grad_analytical - grad_numerical) /
(abs(grad_analytical) + 1e-10)
)
print(comparison_grad, digits = 8)
#> Parameter Analytical Numerical Abs_Diff Rel_Error
#> 1 alpha -0.0407048130 -0.0407061407 1.3276474e-06 3.2616473e-05
#> 2 beta 0.0245082736 0.0245046294 3.6441642e-06 1.4869118e-04
#> 3 gamma 0.0072764653 0.0072782314 1.7660659e-06 2.4270931e-04
#> 4 delta 0.0177067532 0.0177067250 2.8273547e-08 1.5967663e-06
#> 5 lambda -0.0108218991 -0.0108224185 5.1944308e-07 4.7999253e-05
## Example 5: Score Test Statistic
# Score test for H0: theta = theta0
theta0 <- c(1.8, 2.8, 1.3, 1.8, 1.6)
score_theta0 <- grgkw(par = theta0, data = data)
# Fisher information at theta0
fisher_info <- hsgkw(par = theta0, data = data)
# Score test statistic
score_stat <- t(score_theta0) %*% solve(fisher_info) %*% score_theta0
p_value <- pchisq(score_stat, df = 5, lower.tail = FALSE)
cat("\nScore Test:\n")
#>
#> Score Test:
cat("H0: alpha=1.8, beta=2.8, gamma=1.3, delta=1.8, lambda=1.6\n")
#> H0: alpha=1.8, beta=2.8, gamma=1.3, delta=1.8, lambda=1.6
cat("Test statistic:", score_stat, "\n")
#> Test statistic: 258.9207
cat("P-value:", format.pval(p_value, digits = 4), "\n")
#> P-value: < 2.2e-16
## Example 6: Confidence Ellipse (Alpha vs Beta)
# Observed information
obs_info <- hsgkw(par = mle, data = data)
vcov_full <- solve(obs_info)
vcov_2d <- vcov_full[1:2, 1:2]
# Create confidence ellipse
theta <- seq(0, 2 * pi, length.out = round(n/4))
chi2_val <- qchisq(0.95, df = 2)
eig_decomp <- eigen(vcov_2d)
ellipse <- matrix(NA, nrow = round(n/4), ncol = 2)
for (i in 1:round(n/4)) {
v <- c(cos(theta[i]), sin(theta[i]))
ellipse[i, ] <- mle[1:2] + sqrt(chi2_val) *
(eig_decomp$vectors %*% diag(sqrt(eig_decomp$values)) %*% v)
}
# Marginal confidence intervals
se_2d <- sqrt(diag(vcov_2d))
ci_alpha <- mle[1] + c(-1, 1) * 1.96 * se_2d[1]
ci_beta <- mle[2] + c(-1, 1) * 1.96 * se_2d[2]
# Plot
plot(ellipse[, 1], ellipse[, 2], type = "l", lwd = 2, col = "#2E4057",
xlab = expression(alpha), ylab = expression(beta),
main = "95% Confidence Region (Alpha vs Beta)", las = 1)
# Add marginal CIs
abline(v = ci_alpha, col = "#808080", lty = 3, lwd = 1.5)
abline(h = ci_beta, col = "#808080", lty = 3, lwd = 1.5)
points(mle[1], mle[2], pch = 19, col = "#8B0000", cex = 1.5)
points(true_params[1], true_params[2], pch = 17, col = "#006400", cex = 1.5)
legend("topright",
legend = c("MLE", "True", "95% CR", "Marginal 95% CI"),
col = c("#8B0000", "#006400", "#2E4057", "#808080"),
pch = c(19, 17, NA, NA),
lty = c(NA, NA, 1, 3),
lwd = c(NA, NA, 2, 1.5),
bty = "n")
grid(col = "gray90")
## Example 7: Confidence Ellipse (Gamma vs Delta)
# Extract 2x2 submatrix for gamma and delta
vcov_2d_gd <- vcov_full[3:4, 3:4]
# Create confidence ellipse
eig_decomp_gd <- eigen(vcov_2d_gd)
ellipse_gd <- matrix(NA, nrow = round(n/4), ncol = 2)
for (i in 1:round(n/4)) {
v <- c(cos(theta[i]), sin(theta[i]))
ellipse_gd[i, ] <- mle[3:4] + sqrt(chi2_val) *
(eig_decomp_gd$vectors %*% diag(sqrt(eig_decomp_gd$values)) %*% v)
}
# Marginal confidence intervals
se_2d_gd <- sqrt(diag(vcov_2d_gd))
ci_gamma <- mle[3] + c(-1, 1) * 1.96 * se_2d_gd[1]
ci_delta <- mle[4] + c(-1, 1) * 1.96 * se_2d_gd[2]
# Plot
plot(ellipse_gd[, 1], ellipse_gd[, 2], type = "l", lwd = 2, col = "#2E4057",
xlab = expression(gamma), ylab = expression(delta),
main = "95% Confidence Region (Gamma vs Delta)", las = 1)
# Add marginal CIs
abline(v = ci_gamma, col = "#808080", lty = 3, lwd = 1.5)
abline(h = ci_delta, col = "#808080", lty = 3, lwd = 1.5)
points(mle[3], mle[4], pch = 19, col = "#8B0000", cex = 1.5)
points(true_params[3], true_params[4], pch = 17, col = "#006400", cex = 1.5)
legend("topright",
legend = c("MLE", "True", "95% CR", "Marginal 95% CI"),
col = c("#8B0000", "#006400", "#2E4057", "#808080"),
pch = c(19, 17, NA, NA),
lty = c(NA, NA, 1, 3),
lwd = c(NA, NA, 2, 1.5),
bty = "n")
grid(col = "gray90")
# }