Package 'gkwdist'

Density of the Beta Distribution (gamma, delta+1 Parameterization)

Description

Computes the probability density function (PDF) 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 distribution is defined on the interval (0, 1).

Usage

dbeta_(x, gamma = 1, delta = 0, log = FALSE)

Arguments

x	Vector of quantiles (values between 0 and 1).
gamma	First shape parameter (shape1), $\gamma > 0$. Can be a scalar or a vector. Default: 1.0.
delta	Second shape parameter is delta + 1 (shape2), requires $\delta \ge 0$ so that shape2 >= 1. Can be a scalar or a vector. Default: 0.0 (leading to shape2 = 1).
log	Logical; if TRUE, the logarithm of the density is returned ($\log(f(x))$). Default: FALSE.

Details

The probability density function (PDF) calculated by this function corresponds to a standard Beta distribution $Beta(\gamma, \delta+1)$:

$f(x; \gamma, \delta) = \frac{x^{\gamma-1} (1-x)^{(\delta+1)-1}}{B(\gamma, \delta+1)} = \frac{x^{\gamma-1} (1-x)^{\delta}}{B(\gamma, \delta+1)}$

for $0 < x < 1$, where $B(a,b)$ is the Beta function (beta).

This specific parameterization arises as a special case of the five-parameter Generalized Kumaraswamy (GKw) distribution (dgkw) obtained by setting the parameters $\alpha = 1$, $\beta = 1$, and $\lambda = 1$. It is therefore equivalent to the McDonald (Mc)/Beta Power distribution (dmc) with $\lambda = 1$.

Note the difference in the second parameter compared to dbeta, where dbeta(x, shape1, shape2) uses shape2 directly. Here, shape1 = gamma and shape2 = delta + 1.

Value

A vector of density values ($f(x)$) or log-density values ($\log(f(x))$). The length of the result is determined by the recycling rule applied to the arguments (x, gamma, delta). Returns 0 (or -Inf if log = TRUE) for x outside the interval (0, 1), or NaN if parameters are invalid (e.g., gamma <= 0, delta < 0).

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

dbeta (standard R implementation), dgkw (parent distribution density), dmc (McDonald/Beta Power density), pbeta_, qbeta_, rbeta_ (other functions for this parameterization, if they exist).

Examples

# Example values
x_vals <- c(0.2, 0.5, 0.8)
gamma_par <- 2.0 # Corresponds to shape1
delta_par <- 3.0 # Corresponds to shape2 - 1
shape1 <- gamma_par
shape2 <- delta_par + 1

# Calculate density using dbeta_
densities <- dbeta_(x_vals, gamma_par, delta_par)
print(densities)

# Compare with stats::dbeta
densities_stats <- stats::dbeta(x_vals, shape1 = shape1, shape2 = shape2)
print(paste("Max difference vs stats::dbeta:", max(abs(densities - densities_stats))))

# Compare with dgkw setting alpha=1, beta=1, lambda=1
densities_gkw <- dgkw(x_vals,
  alpha = 1.0, beta = 1.0, gamma = gamma_par,
  delta = delta_par, lambda = 1.0
)
print(paste("Max difference vs dgkw:", max(abs(densities - densities_gkw))))

# Compare with dmc setting lambda=1
densities_mc <- dmc(x_vals, gamma = gamma_par, delta = delta_par, lambda = 1.0)
print(paste("Max difference vs dmc:", max(abs(densities - densities_mc))))

# Calculate log-density
log_densities <- dbeta_(x_vals, gamma_par, delta_par, log = TRUE)
print(log_densities)
print(stats::dbeta(x_vals, shape1 = shape1, shape2 = shape2, log = TRUE))

# Plot the density
curve_x <- seq(0.001, 0.999, length.out = 200)
curve_y <- dbeta_(curve_x, gamma = 2, delta = 3) # Beta(2, 4)
plot(curve_x, curve_y,
  type = "l", main = "Beta(2, 4) Density via dbeta_",
  xlab = "x", ylab = "f(x)", col = "blue"
)
curve(stats::dbeta(x, 2, 4), add = TRUE, col = "red", lty = 2)
legend("topright",
  legend = c("dbeta_(gamma=2, delta=3)", "stats::dbeta(shape1=2, shape2=4)"),
  col = c("blue", "red"), lty = c(1, 2), bty = "n"
)

---

Density of the Beta-Kumaraswamy (BKw) Distribution

Description

Computes the probability density function (PDF) for the Beta-Kumaraswamy (BKw) distribution with parameters alpha ($\alpha$), beta ($\beta$), gamma ($\gamma$), and delta ($\delta$). This distribution is defined on the interval (0, 1).

Usage

dbkw(x, alpha = 1, beta = 1, gamma = 1, delta = 0, log = FALSE)

Arguments

x	Vector of quantiles (values between 0 and 1).
alpha	Shape parameter alpha > 0. Can be a scalar or a vector. Default: 1.0.
beta	Shape parameter beta > 0. Can be a scalar or a vector. Default: 1.0.
gamma	Shape parameter gamma > 0. Can be a scalar or a vector. Default: 1.0.
delta	Shape parameter delta >= 0. Can be a scalar or a vector. Default: 0.0.
log	Logical; if TRUE, the logarithm of the density is returned ($\log(f(x))$). Default: FALSE.

Details

The probability density function (PDF) of the Beta-Kumaraswamy (BKw) distribution is given by:

$f(x; \alpha, \beta, \gamma, \delta) = \frac{\alpha \beta}{B(\gamma, \delta+1)} x^{\alpha - 1} \bigl(1 - x^\alpha\bigr)^{\beta(\delta+1) - 1} \bigl[1 - \bigl(1 - x^\alpha\bigr)^\beta\bigr]^{\gamma - 1}$

for $0 < x < 1$, where $B(a,b)$ is the Beta function (beta).

The BKw distribution is a special case of the five-parameter Generalized Kumaraswamy (GKw) distribution (dgkw) obtained by setting the parameter $\lambda = 1$. Numerical evaluation is performed using algorithms similar to those for dgkw, ensuring stability.

Value

A vector of density values ($f(x)$) or log-density values ($\log(f(x))$). The length of the result is determined by the recycling rule applied to the arguments (x, alpha, beta, gamma, delta). Returns 0 (or -Inf if log = TRUE) for x outside the interval (0, 1), or NaN if parameters are invalid (e.g., alpha <= 0, beta <= 0, gamma <= 0, delta < 0).

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.

See Also

dgkw (parent distribution density), pbkw, qbkw, rbkw (other BKw functions),

Examples

# Example values
x_vals <- c(0.2, 0.5, 0.8)
alpha_par <- 2.0
beta_par <- 1.5
gamma_par <- 1.0 # Equivalent to Kw when gamma=1
delta_par <- 0.5

# Calculate density
densities <- dbkw(x_vals, alpha_par, beta_par, gamma_par, delta_par)
print(densities)

# Calculate log-density
log_densities <- dbkw(x_vals, alpha_par, beta_par, gamma_par, delta_par,
  log = TRUE
)
print(log_densities)
# Check: should match log(densities)
print(log(densities))

# Compare with dgkw setting lambda = 1
densities_gkw <- dgkw(x_vals, alpha_par, beta_par,
  gamma = gamma_par,
  delta = delta_par, lambda = 1.0
)
print(paste("Max difference:", max(abs(densities - densities_gkw)))) # Should be near zero

# Plot the density for different gamma values
curve_x <- seq(0.01, 0.99, length.out = 200)
curve_y1 <- dbkw(curve_x, alpha = 2, beta = 3, gamma = 0.5, delta = 1)
curve_y2 <- dbkw(curve_x, alpha = 2, beta = 3, gamma = 1.0, delta = 1)
curve_y3 <- dbkw(curve_x, alpha = 2, beta = 3, gamma = 2.0, delta = 1)

plot(curve_x, curve_y1,
  type = "l", main = "BKw Density Examples (alpha=2, beta=3, delta=1)",
  xlab = "x", ylab = "f(x)", col = "blue", ylim = range(0, curve_y1, curve_y2, curve_y3)
)
lines(curve_x, curve_y2, col = "red")
lines(curve_x, curve_y3, col = "green")
legend("topright",
  legend = c("gamma=0.5", "gamma=1.0", "gamma=2.0"),
  col = c("blue", "red", "green"), lty = 1, bty = "n"
)

---

Density of the Exponentiated Kumaraswamy (EKw) Distribution

Description

Computes the probability density function (PDF) for the Exponentiated Kumaraswamy (EKw) distribution with parameters alpha ($\alpha$), beta ($\beta$), and lambda ($\lambda$). This distribution is defined on the interval (0, 1).

Usage

dekw(x, alpha = 1, beta = 1, lambda = 1, log = FALSE)

Arguments

x	Vector of quantiles (values between 0 and 1).
alpha	Shape parameter alpha > 0. Can be a scalar or a vector. Default: 1.0.
beta	Shape parameter beta > 0. Can be a scalar or a vector. Default: 1.0.
lambda	Shape parameter lambda > 0 (exponent parameter). Can be a scalar or a vector. Default: 1.0.
log	Logical; if TRUE, the logarithm of the density is returned ($\log(f(x))$). Default: FALSE.

Details

The probability density function (PDF) of the Exponentiated Kumaraswamy (EKw) distribution is given by:

$f(x; \alpha, \beta, \lambda) = \lambda \alpha \beta x^{\alpha-1} (1 - x^\alpha)^{\beta-1} \bigl[1 - (1 - x^\alpha)^\beta \bigr]^{\lambda - 1}$

for $0 < x < 1$.

The EKw distribution is a special case of the five-parameter Generalized Kumaraswamy (GKw) distribution (dgkw) obtained by setting the parameters $\gamma = 1$ and $\delta = 0$. When $\lambda = 1$, the EKw distribution reduces to the standard Kumaraswamy distribution.

Value

A vector of density values ($f(x)$) or log-density values ($\log(f(x))$). The length of the result is determined by the recycling rule applied to the arguments (x, alpha, beta, lambda). Returns 0 (or -Inf if log = TRUE) for x outside the interval (0, 1), or NaN if parameters are invalid (e.g., alpha <= 0, beta <= 0, lambda <= 0).

Author(s)

Lopes, J. E.

References

Nadarajah, S., Cordeiro, G. M., & Ortega, E. M. (2012). The exponentiated Kumaraswamy distribution. Journal of the Franklin Institute, 349(3),

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.

See Also

dgkw (parent distribution density), pekw, qekw, rekw (other EKw functions),

Examples

# Example values
x_vals <- c(0.2, 0.5, 0.8)
alpha_par <- 2.0
beta_par <- 3.0
lambda_par <- 1.5 # Exponent parameter

# Calculate density
densities <- dekw(x_vals, alpha_par, beta_par, lambda_par)
print(densities)

# Calculate log-density
log_densities <- dekw(x_vals, alpha_par, beta_par, lambda_par, log = TRUE)
print(log_densities)
# Check: should match log(densities)
print(log(densities))

# Compare with dgkw setting gamma = 1, delta = 0
densities_gkw <- dgkw(x_vals, alpha_par, beta_par,
  gamma = 1.0, delta = 0.0,
  lambda = lambda_par
)
print(paste("Max difference:", max(abs(densities - densities_gkw)))) # Should be near zero

# Plot the density for different lambda values
curve_x <- seq(0.01, 0.99, length.out = 200)
curve_y1 <- dekw(curve_x, alpha = 2, beta = 3, lambda = 0.5) # less peaked
curve_y2 <- dekw(curve_x, alpha = 2, beta = 3, lambda = 1.0) # standard Kw
curve_y3 <- dekw(curve_x, alpha = 2, beta = 3, lambda = 2.0) # more peaked

plot(curve_x, curve_y2,
  type = "l", main = "EKw Density Examples (alpha=2, beta=3)",
  xlab = "x", ylab = "f(x)", col = "red", ylim = range(0, curve_y1, curve_y2, curve_y3)
)
lines(curve_x, curve_y1, col = "blue")
lines(curve_x, curve_y3, col = "green")
legend("topright",
  legend = c("lambda=0.5", "lambda=1.0 (Kw)", "lambda=2.0"),
  col = c("blue", "red", "green"), lty = 1, bty = "n"
)

---

Density of the Generalized Kumaraswamy Distribution

Description

Computes the probability density function (PDF) for the five-parameter Generalized Kumaraswamy (GKw) distribution, defined on the interval (0, 1).

Usage

dgkw(x, alpha = 1, beta = 1, gamma = 1, delta = 0, lambda = 1, log = FALSE)

Arguments

x	Vector of quantiles (values between 0 and 1).
alpha	Shape parameter alpha > 0. Can be a scalar or a vector. Default: 1.0.
beta	Shape parameter beta > 0. Can be a scalar or a vector. Default: 1.0.
gamma	Shape parameter gamma > 0. Can be a scalar or a vector. Default: 1.0.
delta	Shape parameter delta >= 0. Can be a scalar or a vector. Default: 0.0.
lambda	Shape parameter lambda > 0. Can be a scalar or a vector. Default: 1.0.
log	Logical; if TRUE, the logarithm of the density is returned. Default: FALSE.

Details

The probability density function of the Generalized Kumaraswamy (GKw) distribution with parameters alpha ($\alpha$), beta ($\beta$), gamma ($\gamma$), delta ($\delta$), and lambda ($\lambda$) is given by:

$f(x; \alpha, \beta, \gamma, \delta, \lambda) = \frac{\lambda \alpha \beta x^{\alpha-1}(1-x^{\alpha})^{\beta-1}} {B(\gamma, \delta+1)} [1-(1-x^{\alpha})^{\beta}]^{\gamma\lambda-1} [1-[1-(1-x^{\alpha})^{\beta}]^{\lambda}]^{\delta}$

for $x \in (0,1)$, where $B(a, b)$ is the Beta function beta.

This distribution was proposed by Cordeiro & de Castro (2011) and includes several other distributions as special cases:

• Kumaraswamy (Kw): gamma = 1, delta = 0, lambda = 1

• Exponentiated Kumaraswamy (EKw): gamma = 1, delta = 0

• Beta-Kumaraswamy (BKw): lambda = 1

• Generalized Beta type 1 (GB1 - implies McDonald): alpha = 1, beta = 1

• Beta distribution: alpha = 1, beta = 1, lambda = 1

The function includes checks for valid parameters and input values x. It uses numerical stabilization for x close to 0 or 1.

Value

A vector of density values ($f(x)$) or log-density values ($\log(f(x))$). The length of the result is determined by the recycling rule applied to the arguments (x, alpha, beta, gamma, delta, lambda). Returns 0 (or -Inf if log = TRUE) for x outside the interval (0, 1), or NaN if parameters are invalid.

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, 81(7), 883-898.

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

See Also

pgkw, qgkw, rgkw (if these exist), dbeta, integrate

Examples

# Simple density evaluation at a point
dgkw(0.5, alpha = 2, beta = 3, gamma = 1, delta = 0, lambda = 1) # Kw case

# Plot the PDF for various parameter sets
x_vals <- seq(0.01, 0.99, by = 0.01)

# Standard Kumaraswamy (gamma=1, delta=0, lambda=1)
pdf_kw <- dgkw(x_vals, alpha = 2, beta = 3, gamma = 1, delta = 0, lambda = 1)

# Beta equivalent (alpha=1, beta=1, lambda=1) - Beta(gamma, delta+1)
pdf_beta <- dgkw(x_vals, alpha = 1, beta = 1, gamma = 2, delta = 3, lambda = 1)
# Compare with stats::dbeta
pdf_beta_check <- stats::dbeta(x_vals, shape1 = 2, shape2 = 3 + 1)
# max(abs(pdf_beta - pdf_beta_check)) # Should be close to zero

# Exponentiated Kumaraswamy (gamma=1, delta=0)
pdf_ekw <- dgkw(x_vals, alpha = 2, beta = 3, gamma = 1, delta = 0, lambda = 2)

plot(x_vals, pdf_kw,
  type = "l", ylim = range(c(pdf_kw, pdf_beta, pdf_ekw)),
  main = "GKw Densities Examples", ylab = "f(x)", xlab = "x", col = "blue"
)
lines(x_vals, pdf_beta, col = "red")
lines(x_vals, pdf_ekw, col = "green")
legend("topright",
  legend = c("Kw(2,3)", "Beta(2,4) equivalent", "EKw(2,3, lambda=2)"),
  col = c("blue", "red", "green"), lty = 1, bty = "n"
)

# Log-density
log.pdf_val <- dgkw(0.5, 2, 3, 1, 0, 1, log = TRUE)
print(log.pdf_val)
print(log(dgkw(0.5, 2, 3, 1, 0, 1))) # Should match

---

Density of the Kumaraswamy-Kumaraswamy (kkw) Distribution

Description

Computes the probability density function (PDF) for the Kumaraswamy-Kumaraswamy (kkw) distribution with parameters alpha ($\alpha$), beta ($\beta$), delta ($\delta$), and lambda ($\lambda$). This distribution is defined on the interval (0, 1).

Usage

dkkw(x, alpha = 1, beta = 1, delta = 0, lambda = 1, log = FALSE)

Arguments

x	Vector of quantiles (values between 0 and 1).
alpha	Shape parameter alpha > 0. Can be a scalar or a vector. Default: 1.0.
beta	Shape parameter beta > 0. Can be a scalar or a vector. Default: 1.0.
delta	Shape parameter delta >= 0. Can be a scalar or a vector. Default: 0.0.
lambda	Shape parameter lambda > 0. Can be a scalar or a vector. Default: 1.0.
log	Logical; if TRUE, the logarithm of the density is returned ($\log(f(x))$). Default: FALSE.

Details

The Kumaraswamy-Kumaraswamy (kkw) distribution is a special case of the five-parameter Generalized Kumaraswamy distribution (dgkw) obtained by setting the parameter $\gamma = 1$.

The probability density function is given by:

$f(x; \alpha, \beta, \delta, \lambda) = (\delta + 1) \lambda \alpha \beta x^{\alpha - 1} (1 - x^\alpha)^{\beta - 1} \bigl[1 - (1 - x^\alpha)^\beta\bigr]^{\lambda - 1} \bigl\{1 - \bigl[1 - (1 - x^\alpha)^\beta\bigr]^\lambda\bigr\}^{\delta}$

for $0 < x < 1$. Note that $1/(\delta+1)$ corresponds to the Beta function term $B(1, \delta+1)$ when $\gamma=1$.

Numerical evaluation follows similar stability considerations as dgkw.

Value

A vector of density values ($f(x)$) or log-density values ($\log(f(x))$). The length of the result is determined by the recycling rule applied to the arguments (x, alpha, beta, delta, lambda). Returns 0 (or -Inf if log = TRUE) for x outside the interval (0, 1), or NaN if parameters are invalid (e.g., alpha <= 0, beta <= 0, delta < 0, lambda <= 0).

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.

See Also

dgkw (parent distribution density), pkkw, qkkw, rkkw (if they exist), dbeta

Examples

# Example values
x_vals <- c(0.2, 0.5, 0.8)
alpha_par <- 2.0
beta_par <- 3.0
delta_par <- 0.5
lambda_par <- 1.5

# Calculate density
densities <- dkkw(x_vals, alpha_par, beta_par, delta_par, lambda_par)
print(densities)

# Calculate log-density
log_densities <- dkkw(x_vals, alpha_par, beta_par, delta_par, lambda_par,
  log = TRUE
)
print(log_densities)
# Check: should match log(densities)
print(log(densities))

# Compare with dgkw setting gamma = 1
densities_gkw <- dgkw(x_vals, alpha_par, beta_par,
  gamma = 1.0,
  delta_par, lambda_par
)
print(paste("Max difference:", max(abs(densities - densities_gkw)))) # Should be near zero

# Plot the density
curve_x <- seq(0.01, 0.99, length.out = 200)
curve_y <- dkkw(curve_x, alpha_par, beta_par, delta_par, lambda_par)
plot(curve_x, curve_y,
  type = "l", main = "kkw Density Example",
  xlab = "x", ylab = "f(x)", col = "blue"
)

---

Density of the Kumaraswamy (Kw) Distribution

Description

Computes the probability density function (PDF) for the two-parameter Kumaraswamy (Kw) distribution with shape parameters alpha ($\alpha$) and beta ($\beta$). This distribution is defined on the interval (0, 1).

Usage

dkw(x, alpha = 1, beta = 1, log = FALSE)

Arguments

x	Vector of quantiles (values between 0 and 1).
alpha	Shape parameter alpha > 0. Can be a scalar or a vector. Default: 1.0.
beta	Shape parameter beta > 0. Can be a scalar or a vector. Default: 1.0.
log	Logical; if TRUE, the logarithm of the density is returned ($\log(f(x))$). Default: FALSE.

Details

The probability density function (PDF) of the Kumaraswamy (Kw) distribution is given by:

$f(x; \alpha, \beta) = \alpha \beta x^{\alpha-1} (1 - x^\alpha)^{\beta-1}$

for $0 < x < 1$, $\alpha > 0$, and $\beta > 0$.

The Kumaraswamy distribution is identical to the Generalized Kumaraswamy (GKw) distribution (dgkw) with parameters $\gamma = 1$, $\delta = 0$, and $\lambda = 1$. It is also a special case of the Exponentiated Kumaraswamy (dekw) with $\lambda = 1$, and the Kumaraswamy-Kumaraswamy (dkkw) with $\delta = 0$ and $\lambda = 1$.

Value

A vector of density values ($f(x)$) or log-density values ($\log(f(x))$). The length of the result is determined by the recycling rule applied to the arguments (x, alpha, beta). Returns 0 (or -Inf if log = TRUE) for x outside the interval (0, 1), or NaN if parameters are invalid (e.g., alpha <= 0, beta <= 0).

Author(s)

Lopes, J. E.

References

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

Jones, M. C. (2009). Kumaraswamy's distribution: A beta-type distribution with some tractability advantages. Statistical Methodology, 6(1), 70-81.

See Also

dgkw (parent distribution density), dekw, dkkw, pkw, qkw, rkw (other Kw functions), dbeta

Examples

# Example values
x_vals <- c(0.2, 0.5, 0.8)
alpha_par <- 2.0
beta_par <- 3.0

# Calculate density using dkw
densities <- dkw(x_vals, alpha_par, beta_par)
print(densities)

# Calculate log-density
log_densities <- dkw(x_vals, alpha_par, beta_par, log = TRUE)
print(log_densities)
# Check: should match log(densities)
print(log(densities))

# Compare with dgkw setting gamma = 1, delta = 0, lambda = 1
densities_gkw <- dgkw(x_vals, alpha_par, beta_par,
  gamma = 1.0, delta = 0.0,
  lambda = 1.0
)
print(paste("Max difference:", max(abs(densities - densities_gkw)))) # Should be near zero

# Plot the density for different shape parameter combinations
curve_x <- seq(0.001, 0.999, length.out = 200)
plot(curve_x, dkw(curve_x, alpha = 2, beta = 3),
  type = "l",
  main = "Kumaraswamy Density Examples", xlab = "x", ylab = "f(x)",
  col = "blue", ylim = c(0, 4)
)
lines(curve_x, dkw(curve_x, alpha = 3, beta = 2), col = "red")
lines(curve_x, dkw(curve_x, alpha = 0.5, beta = 0.5), col = "green") # U-shaped
lines(curve_x, dkw(curve_x, alpha = 5, beta = 1), col = "purple") # J-shaped
lines(curve_x, dkw(curve_x, alpha = 1, beta = 3), col = "orange") # J-shaped (reversed)
legend("top",
  legend = c("a=2, b=3", "a=3, b=2", "a=0.5, b=0.5", "a=5, b=1", "a=1, b=3"),
  col = c("blue", "red", "green", "purple", "orange"), lty = 1, bty = "n", ncol = 2
)

---

Density of the McDonald (Mc)/Beta Power Distribution Distribution

Description

Computes the probability density function (PDF) for the McDonald (Mc) distribution (also previously referred to as Beta Power) with parameters gamma ($\gamma$), delta ($\delta$), and lambda ($\lambda$). This distribution is defined on the interval (0, 1).

Usage

dmc(x, gamma = 1, delta = 0, lambda = 1, log = FALSE)

Arguments

x	Vector of quantiles (values between 0 and 1).
gamma	Shape parameter gamma > 0. Can be a scalar or a vector. Default: 1.0.
delta	Shape parameter delta >= 0. Can be a scalar or a vector. Default: 0.0.
lambda	Shape parameter lambda > 0. Can be a scalar or a vector. Default: 1.0.
log	Logical; if TRUE, the logarithm of the density is returned ($\log(f(x))$). Default: FALSE.

Details

The probability density function (PDF) of the McDonald (Mc) distribution is given by:

$f(x; \gamma, \delta, \lambda) = \frac{\lambda}{B(\gamma,\delta+1)} x^{\gamma \lambda - 1} (1 - x^\lambda)^\delta$

for $0 < x < 1$, where $B(a,b)$ is the Beta function (beta).

The Mc distribution is a special case of the five-parameter Generalized Kumaraswamy (GKw) distribution (dgkw) obtained by setting the parameters $\alpha = 1$ and $\beta = 1$. It was introduced by McDonald (1984) and is related to the Generalized Beta distribution of the first kind (GB1). When $\lambda=1$, it simplifies to the standard Beta distribution with parameters $\gamma$ and $\delta+1$.

Value

A vector of density values ($f(x)$) or log-density values ($\log(f(x))$). The length of the result is determined by the recycling rule applied to the arguments (x, gamma, delta, lambda). Returns 0 (or -Inf if log = TRUE) for x outside the interval (0, 1), or NaN if parameters are invalid (e.g., gamma <= 0, delta < 0, lambda <= 0).

Author(s)

Lopes, J. E.

References

McDonald, J. B. (1984). Some generalized functions for the size distribution of income. Econometrica, 52(3), 647-663.

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.

See Also

dgkw (parent distribution density), pmc, qmc, rmc (other Mc functions), dbeta

Examples

# Example values
x_vals <- c(0.2, 0.5, 0.8)
gamma_par <- 2.0
delta_par <- 1.5
lambda_par <- 1.0 # Equivalent to Beta(gamma, delta+1)

# Calculate density using dmc
densities <- dmc(x_vals, gamma_par, delta_par, lambda_par)
print(densities)
# Compare with Beta density
print(stats::dbeta(x_vals, shape1 = gamma_par, shape2 = delta_par + 1))

# Calculate log-density
log_densities <- dmc(x_vals, gamma_par, delta_par, lambda_par, log = TRUE)
print(log_densities)

# Compare with dgkw setting alpha = 1, beta = 1
densities_gkw <- dgkw(x_vals,
  alpha = 1.0, beta = 1.0, gamma = gamma_par,
  delta = delta_par, lambda = lambda_par
)
print(paste("Max difference:", max(abs(densities - densities_gkw)))) # Should be near zero

# Plot the density for different lambda values
curve_x <- seq(0.01, 0.99, length.out = 200)
curve_y1 <- dmc(curve_x, gamma = 2, delta = 3, lambda = 0.5)
curve_y2 <- dmc(curve_x, gamma = 2, delta = 3, lambda = 1.0) # Beta(2, 4)
curve_y3 <- dmc(curve_x, gamma = 2, delta = 3, lambda = 2.0)

plot(curve_x, curve_y2,
  type = "l", main = "McDonald (Mc) Density (gamma=2, delta=3)",
  xlab = "x", ylab = "f(x)", col = "red", ylim = range(0, curve_y1, curve_y2, curve_y3)
)
lines(curve_x, curve_y1, col = "blue")
lines(curve_x, curve_y3, col = "green")
legend("topright",
  legend = c("lambda=0.5", "lambda=1.0 (Beta)", "lambda=2.0"),
  col = c("blue", "red", "green"), lty = 1, bty = "n"
)

---

Estimate Distribution Parameters Using Method of Moments

Description

Estimates parameters for various distribution families from the Generalized Kumaraswamy family using the method of moments. The implementation is optimized for numerical stability and computational efficiency through Nelder-Mead optimization and adaptive numerical integration.

Usage

gkwgetstartvalues(x, family = "gkw", n_starts = 5L)

Arguments

x	Numeric vector of observations. All values must be in the open interval (0,1). Values outside this range will be automatically truncated to avoid numerical issues.
family	Character string specifying the distribution family. Valid options are: "gkw" (Generalized Kumaraswamy - 5 parameters), "bkw" (Beta-Kumaraswamy - 4 parameters), "kkw" (Kumaraswamy-Kumaraswamy - 4 parameters), "ekw" (Exponentiated Kumaraswamy - 3 parameters), "mc" (McDonald - 3 parameters), "kw" (Kumaraswamy - 2 parameters), "beta" (Beta - 2 parameters). The string is case-insensitive. Default is "gkw".
n_starts	Integer specifying the number of different initial parameter values to try during optimization. More starting points increase the probability of finding the global optimum at the cost of longer computation time. Default is 5.

Details

The function uses the method of moments to estimate distribution parameters by minimizing the weighted sum of squared relative errors between theoretical and sample moments (orders 1 through 5). The optimization employs the Nelder-Mead simplex algorithm, which is derivative-free and particularly robust for this problem.

Key implementation features: logarithmic calculations for numerical stability, adaptive numerical integration using Simpson's rule with fallback to trapezoidal rule, multiple random starting points to avoid local minima, decreasing weights for higher-order moments (1.0, 0.8, 0.6, 0.4, 0.2), and automatic parameter constraint enforcement.

Parameter Constraints: All parameters are constrained to positive values. Additionally, family-specific constraints are enforced: alpha and beta in (0.1, 50.0), gamma in (0.1, 10.0) for GKw-related families or (0.1, 50.0) for Beta, delta in (0.01, 10.0), and lambda in (0.1, 20.0).

The function will issue warnings for empty input vectors, sample sizes less than 10 (unreliable estimation), or failure to find valid parameter estimates (returns defaults).

Value

Named numeric vector containing the estimated parameters for the specified distribution family. Parameter names correspond to the distribution specification. If estimation fails, returns a vector of NA values with appropriate parameter names.

References

Jones, M. C. (2009). Kumaraswamy's distribution: A beta-type distribution with some tractability advantages. Statistical Methodology, 6(1), 70-81.

Examples

# Generate sample data from Beta distribution
set.seed(123)
x <- rbeta(100, shape1 = 2, shape2 = 3)

# Estimate Beta parameters
params_beta <- gkwgetstartvalues(x, family = "beta")
print(params_beta)

# Estimate Kumaraswamy parameters
params_kw <- gkwgetstartvalues(x, family = "kw")
print(params_kw)

# Estimate GKw parameters with more starting points
params_gkw <- gkwgetstartvalues(x, family = "gkw", n_starts = 10)
print(params_gkw)

---

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)

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

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

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,

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

See Also

grgkw, grmc (related gradients), llbeta (negative log-likelihood function), hsbeta (Hessian, if available), dbeta_, pbeta_, qbeta_, rbeta_, optim, grad (for numerical gradient comparison), digamma.

Examples

## Example 1: Basic Gradient Evaluation

# Generate sample data
set.seed(123)
n <- 1000
true_params <- c(gamma = 2.0, delta = 3.0)
data <- rbeta_(n, gamma = true_params[1], delta = true_params[2])

# Evaluate gradient at true parameters
grad_true <- grbeta(par = true_params, data = data)
cat("Gradient at true parameters:\n")
print(grad_true)
cat("Norm:", sqrt(sum(grad_true^2)), "\n")

# Evaluate at different parameter values
test_params <- rbind(
  c(1.5, 2.5),
  c(2.0, 3.0),
  c(2.5, 3.5)
)

grad_norms <- apply(test_params, 1, function(p) {
  g <- grbeta(p, data)
  sqrt(sum(g^2))
})

results <- data.frame(
  Gamma = test_params[, 1],
  Delta = test_params[, 2],
  Grad_Norm = grad_norms
)
print(results, digits = 4)


## Example 2: Gradient in Optimization

# Optimization with analytical gradient
fit_with_grad <- optim(
  par = c(1.5, 2.5),
  fn = llbeta,
  gr = grbeta,
  data = data,
  method = "L-BFGS-B",
  lower = c(0.01, 0.01),
  upper = c(100, 100),
  hessian = TRUE,
  control = list(trace = 0)
)

# Optimization without gradient
fit_no_grad <- optim(
  par = c(1.5, 2.5),
  fn = llbeta,
  data = data,
  method = "L-BFGS-B",
  lower = c(0.01, 0.01),
  upper = c(100, 100),
  hessian = TRUE,
  control = list(trace = 0)
)

comparison <- data.frame(
  Method = c("With Gradient", "Without Gradient"),
  Gamma = c(fit_with_grad$par[1], fit_no_grad$par[1]),
  Delta = c(fit_with_grad$par[2], fit_no_grad$par[2]),
  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)


## Example 3: Verifying Gradient at MLE

mle <- fit_with_grad$par
names(mle) <- c("gamma", "delta")

# At MLE, gradient should be approximately zero
gradient_at_mle <- grbeta(par = mle, data = data)
cat("\nGradient at MLE:\n")
print(gradient_at_mle)
cat("Max absolute component:", max(abs(gradient_at_mle)), "\n")
cat("Gradient norm:", sqrt(sum(gradient_at_mle^2)), "\n")


## 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 <- grbeta(par = mle, data = data)
grad_numerical <- numerical_gradient(llbeta, mle, data)

comparison_grad <- data.frame(
  Parameter = c("gamma", "delta"),
  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)


## Example 5: Score Test Statistic

# Score test for H0: theta = theta0
theta0 <- c(1.8, 2.8)
score_theta0 <- -grbeta(par = theta0, data = data)

# Fisher information at theta0
fisher_info <- hsbeta(par = theta0, data = data)

# Score test statistic
score_stat <- t(score_theta0) %*% solve(fisher_info) %*% score_theta0
p_value <- pchisq(score_stat, df = 2, lower.tail = FALSE)

cat("\nScore Test:\n")
cat("H0: gamma=1.8, delta=2.8\n")
cat("Test statistic:", score_stat, "\n")
cat("P-value:", format.pval(p_value, digits = 4), "\n")


## Example 6: Confidence Ellipse (Gamma vs Delta)

# Observed information
obs_info <- hsbeta(par = mle, data = data)
vcov_full <- solve(obs_info)

# Create confidence ellipse
theta <- seq(0, 2 * pi, length.out = 100)
chi2_val <- qchisq(0.95, df = 2)

eig_decomp <- eigen(vcov_full)
ellipse <- matrix(NA, nrow = 100, ncol = 2)
for (i in 1:100) {
  v <- c(cos(theta[i]), sin(theta[i]))
  ellipse[i, ] <- mle + sqrt(chi2_val) *
    (eig_decomp$vectors %*% diag(sqrt(eig_decomp$values)) %*% v)
}

# Marginal confidence intervals
se_2d <- sqrt(diag(vcov_full))
ci_gamma <- mle[1] + c(-1, 1) * 1.96 * se_2d[1]
ci_delta <- mle[2] + c(-1, 1) * 1.96 * se_2d[2]

# Plot

plot(ellipse[, 1], ellipse[, 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[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")

---

Gradient of the Negative Log-Likelihood for the BKw Distribution

Description

Computes the gradient vector (vector of first 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 gradient is typically used in optimization algorithms for maximum likelihood estimation.

Usage

grbkw(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

The components of the gradient vector of the negative log-likelihood ($-\nabla \ell(\theta | \mathbf{x})$) for the BKw ($\lambda=1$) model 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(\delta+1)-1}{v_i} - \frac{(\gamma-1) \beta v_i^{\beta-1}}{w_i}\right)\right]$

$-\frac{\partial \ell}{\partial \beta} = -\frac{n}{\beta} - (\delta+1)\sum_{i=1}^{n}\ln(v_i) + \sum_{i=1}^{n}\left[\frac{(\gamma-1) v_i^{\beta} \ln(v_i)}{w_i}\right]$

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

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

where:

• $v_i = 1 - x_i^{\alpha}$

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

• $\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 general GKw gradient (grgkw) components for $\alpha, \beta, \gamma, \delta$ evaluated at $\lambda=1$. Note that the component for $\lambda$ is omitted. Numerical stability is maintained through careful implementation.

Value

Returns a numeric vector of length 4 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)$. 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).

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 gradient formulas might be derived or sourced from additional references).

See Also

grgkw (parent distribution gradient), llbkw (negative log-likelihood for BKw), hsbkw (Hessian for BKw, if available), dbkw (density for BKw), optim, grad (for numerical gradient comparison), digamma.

Examples

## Example 1: Basic Gradient Evaluation
# Generate sample data
set.seed(2203)
n <- 1000
true_params <- c(alpha = 2.0, beta = 1.5, gamma = 1.5, delta = 0.5)
data <- rbkw(n,
  alpha = true_params[1], beta = true_params[2],
  gamma = true_params[3], delta = true_params[4]
)

# Evaluate gradient at true parameters
grad_true <- grbkw(par = true_params, data = data)
cat("Gradient at true parameters:\n")
print(grad_true)
cat("Norm:", sqrt(sum(grad_true^2)), "\n")

# Evaluate at different parameter values
test_params <- rbind(
  c(1.5, 1.0, 1.0, 0.3),
  c(2.0, 1.5, 1.5, 0.5),
  c(2.5, 2.0, 2.0, 0.7)
)

grad_norms <- apply(test_params, 1, function(p) {
  g <- grbkw(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],
  Grad_Norm = grad_norms
)
print(results, digits = 4)


## Example 2: Gradient in Optimization

# Optimization with analytical gradient
fit_with_grad <- optim(
  par = c(1.8, 1.2, 1.1, 0.3),
  fn = llbkw,
  gr = grbkw,
  data = data,
  method = "Nelder-Mead",
  hessian = TRUE,
  control = list(trace = 0)
)

# Optimization without gradient
fit_no_grad <- optim(
  par = c(1.8, 1.2, 1.1, 0.3),
  fn = llbkw,
  data = data,
  method = "Nelder-Mead",
  hessian = TRUE,
  control = list(trace = 0)
)

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]),
  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)


## Example 3: Verifying Gradient at MLE

mle <- fit_with_grad$par
names(mle) <- c("alpha", "beta", "gamma", "delta")

# At MLE, gradient should be approximately zero
gradient_at_mle <- grbkw(par = mle, data = data)
cat("\nGradient at MLE:\n")
print(gradient_at_mle)
cat("Max absolute component:", max(abs(gradient_at_mle)), "\n")
cat("Gradient norm:", sqrt(sum(gradient_at_mle^2)), "\n")


## 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 <- grbkw(par = mle, data = data)
grad_numerical <- numerical_gradient(llbkw, mle, data)

comparison_grad <- data.frame(
  Parameter = c("alpha", "beta", "gamma", "delta"),
  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)


## Example 5: Score Test Statistic

# Score test for H0: theta = theta0
theta0 <- c(1.8, 1.3, 1.2, 0.4)
score_theta0 <- -grbkw(par = theta0, data = data)

# Fisher information at theta0
fisher_info <- hsbkw(par = theta0, data = data)

# Score test statistic
score_stat <- t(score_theta0) %*% solve(fisher_info) %*% score_theta0
p_value <- pchisq(score_stat, df = 4, lower.tail = FALSE)

cat("\nScore Test:\n")
cat("H0: alpha=1.8, beta=1.3, gamma=1.2, delta=0.4\n")
cat("Test statistic:", score_stat, "\n")
cat("P-value:", format.pval(p_value, digits = 4), "\n")


## Example 6: Confidence Ellipses (Selected pairs)

# Observed information
obs_info <- hsbkw(par = mle, data = data)
vcov_full <- solve(obs_info)

# Create confidence ellipses
theta <- seq(0, 2 * pi, length.out = 100)
chi2_val <- qchisq(0.95, df = 2)

# Alpha vs Beta ellipse
vcov_ab <- vcov_full[1:2, 1:2]
eig_decomp_ab <- eigen(vcov_ab)
ellipse_ab <- matrix(NA, nrow = 100, ncol = 2)
for (i in 1:100) {
  v <- c(cos(theta[i]), sin(theta[i]))
  ellipse_ab[i, ] <- mle[1:2] + sqrt(chi2_val) *
    (eig_decomp_ab$vectors %*% diag(sqrt(eig_decomp_ab$values)) %*% v)
}

# Alpha vs Gamma ellipse
vcov_ag <- vcov_full[c(1, 3), c(1, 3)]
eig_decomp_ag <- eigen(vcov_ag)
ellipse_ag <- matrix(NA, nrow = 100, ncol = 2)
for (i in 1:100) {
  v <- c(cos(theta[i]), sin(theta[i]))
  ellipse_ag[i, ] <- mle[c(1, 3)] + sqrt(chi2_val) *
    (eig_decomp_ag$vectors %*% diag(sqrt(eig_decomp_ag$values)) %*% v)
}

# Beta vs Delta ellipse
vcov_bd <- vcov_full[c(2, 4), c(2, 4)]
eig_decomp_bd <- eigen(vcov_bd)
ellipse_bd <- matrix(NA, nrow = 100, ncol = 2)
for (i in 1:100) {
  v <- c(cos(theta[i]), sin(theta[i]))
  ellipse_bd[i, ] <- mle[c(2, 4)] + sqrt(chi2_val) *
    (eig_decomp_bd$vectors %*% diag(sqrt(eig_decomp_bd$values)) %*% v)
}

# Marginal confidence intervals
se_ab <- sqrt(diag(vcov_ab))
ci_alpha_ab <- mle[1] + c(-1, 1) * 1.96 * se_ab[1]
ci_beta_ab <- mle[2] + c(-1, 1) * 1.96 * se_ab[2]

se_ag <- sqrt(diag(vcov_ag))
ci_alpha_ag <- mle[1] + c(-1, 1) * 1.96 * se_ag[1]
ci_gamma_ag <- mle[3] + c(-1, 1) * 1.96 * se_ag[2]

se_bd <- sqrt(diag(vcov_bd))
ci_beta_bd <- mle[2] + c(-1, 1) * 1.96 * se_bd[1]
ci_delta_bd <- mle[4] + c(-1, 1) * 1.96 * se_bd[2]

# Plot selected ellipses

# Alpha vs Beta
plot(ellipse_ab[, 1], ellipse_ab[, 2],
  type = "l", lwd = 2, col = "#2E4057",
  xlab = expression(alpha), ylab = expression(beta),
  main = "Alpha vs Beta", las = 1, xlim = range(ellipse_ab[, 1], ci_alpha_ab),
  ylim = range(ellipse_ab[, 2], ci_beta_ab)
)
abline(v = ci_alpha_ab, col = "#808080", lty = 3, lwd = 1.5)
abline(h = ci_beta_ab, 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)
grid(col = "gray90")

# Alpha vs Gamma
plot(ellipse_ag[, 1], ellipse_ag[, 2],
  type = "l", lwd = 2, col = "#2E4057",
  xlab = expression(alpha), ylab = expression(gamma),
  main = "Alpha vs Gamma", las = 1, xlim = range(ellipse_ag[, 1], ci_alpha_ag),
  ylim = range(ellipse_ag[, 2], ci_gamma_ag)
)
abline(v = ci_alpha_ag, col = "#808080", lty = 3, lwd = 1.5)
abline(h = ci_gamma_ag, col = "#808080", lty = 3, lwd = 1.5)
points(mle[1], mle[3], pch = 19, col = "#8B0000", cex = 1.5)
points(true_params[1], true_params[3], pch = 17, col = "#006400", cex = 1.5)
grid(col = "gray90")

# Beta vs Delta
plot(ellipse_bd[, 1], ellipse_bd[, 2],
  type = "l", lwd = 2, col = "#2E4057",
  xlab = expression(beta), ylab = expression(delta),
  main = "Beta vs Delta", las = 1, xlim = range(ellipse_bd[, 1], ci_beta_bd),
  ylim = range(ellipse_bd[, 2], ci_delta_bd)
)
abline(v = ci_beta_bd, col = "#808080", lty = 3, lwd = 1.5)
abline(h = ci_delta_bd, col = "#808080", lty = 3, lwd = 1.5)
points(mle[2], mle[4], pch = 19, col = "#8B0000", cex = 1.5)
points(true_params[2], true_params[4], pch = 17, col = "#006400", cex = 1.5)
grid(col = "gray90")

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", cex = 0.8
)

---

Gradient of the Negative Log-Likelihood for the EKw Distribution

Description

Computes the gradient vector (vector of first partial derivatives) of the negative log-likelihood function for the Exponentiated Kumaraswamy (EKw) distribution with parameters alpha ($\alpha$), beta ($\beta$), and lambda ($\lambda$). This distribution is the special case of the Generalized Kumaraswamy (GKw) distribution where $\gamma = 1$ and $\delta = 0$. The gradient is useful for optimization.

Usage

grekw(par, data)

Arguments

par	A numeric vector of length 3 containing the distribution parameters in the order: alpha ($\alpha > 0$), beta ($\beta > 0$), lambda ($\lambda > 0$).
data	A numeric vector of observations. All values must be strictly between 0 and 1 (exclusive).

Details

The components of the gradient vector of the negative log-likelihood ($-\nabla \ell(\theta | \mathbf{x})$) for the EKw ($\gamma=1, \delta=0$) model 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{(\lambda-1) \beta v_i^{\beta-1}}{w_i}\right)\right]$

$-\frac{\partial \ell}{\partial \beta} = -\frac{n}{\beta} - \sum_{i=1}^{n}\ln(v_i) + \sum_{i=1}^{n}\left[\frac{(\lambda-1) v_i^{\beta} \ln(v_i)}{w_i}\right]$

$-\frac{\partial \ell}{\partial \lambda} = -\frac{n}{\lambda} - \sum_{i=1}^{n}\ln(w_i)$

where:

• $v_i = 1 - x_i^{\alpha}$

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

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 $\gamma=1, \delta=0$.

Value

Returns a numeric vector of length 3 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 \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).

Author(s)

Lopes, J. E.

References

Nadarajah, S., Cordeiro, G. M., & Ortega, E. M. (2012). The exponentiated Kumaraswamy distribution. Journal of the Franklin Institute, 349(3),

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 gradient formulas might be derived or sourced from additional references).

See Also

grgkw (parent distribution gradient), llekw (negative log-likelihood for EKw), hsekw (Hessian for EKw, if available), dekw (density for EKw), optim, grad (for numerical gradient comparison).

Examples

## Example 1: Basic Gradient Evaluation

# Generate sample data
set.seed(123)
n <- 1000
true_params <- c(alpha = 2.5, beta = 3.5, lambda = 2.0)
data <- rekw(n,
  alpha = true_params[1], beta = true_params[2],
  lambda = true_params[3]
)

# Evaluate gradient at true parameters
grad_true <- grekw(par = true_params, data = data)
cat("Gradient at true parameters:\n")
print(grad_true)
cat("Norm:", sqrt(sum(grad_true^2)), "\n")

# Evaluate at different parameter values
test_params <- rbind(
  c(2.0, 3.0, 1.5),
  c(2.5, 3.5, 2.0),
  c(3.0, 4.0, 2.5)
)

grad_norms <- apply(test_params, 1, function(p) {
  g <- grekw(p, data)
  sqrt(sum(g^2))
})

results <- data.frame(
  Alpha = test_params[, 1],
  Beta = test_params[, 2],
  Lambda = test_params[, 3],
  Grad_Norm = grad_norms
)
print(results, digits = 4)


## Example 2: Gradient in Optimization

# Optimization with analytical gradient
fit_with_grad <- optim(
  par = c(2, 3, 1.5),
  fn = llekw,
  gr = grekw,
  data = data,
  method = "BFGS",
  hessian = TRUE,
  control = list(trace = 0)
)

# Optimization without gradient
fit_no_grad <- optim(
  par = c(2, 3, 1.5),
  fn = llekw,
  data = data,
  method = "BFGS",
  hessian = TRUE,
  control = list(trace = 0)
)

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]),
  Lambda = c(fit_with_grad$par[3], fit_no_grad$par[3]),
  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)


## Example 3: Verifying Gradient at MLE

mle <- fit_with_grad$par
names(mle) <- c("alpha", "beta", "lambda")

# At MLE, gradient should be approximately zero
gradient_at_mle <- grekw(par = mle, data = data)
cat("\nGradient at MLE:\n")
print(gradient_at_mle)
cat("Max absolute component:", max(abs(gradient_at_mle)), "\n")
cat("Gradient norm:", sqrt(sum(gradient_at_mle^2)), "\n")


## 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 <- grekw(par = mle, data = data)
grad_numerical <- numerical_gradient(llekw, mle, data)

comparison_grad <- data.frame(
  Parameter = c("alpha", "beta", "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)


## Example 5: Score Test Statistic

# Score test for H0: theta = theta0
theta0 <- c(2.2, 3.2, 1.8)
score_theta0 <- -grekw(par = theta0, data = data)

# Fisher information at theta0
fisher_info <- hsekw(par = theta0, data = data)

# Score test statistic
score_stat <- t(score_theta0) %*% solve(fisher_info) %*% score_theta0
p_value <- pchisq(score_stat, df = 3, lower.tail = FALSE)

cat("\nScore Test:\n")
cat("H0: alpha=2.2, beta=3.2, lambda=1.8\n")
cat("Test statistic:", score_stat, "\n")
cat("P-value:", format.pval(p_value, digits = 4), "\n")


## Example 6: Confidence Ellipse (Alpha vs Beta)

# Observed information
obs_info <- hsekw(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 = 100)
chi2_val <- qchisq(0.95, df = 2)

eig_decomp <- eigen(vcov_2d)
ellipse <- matrix(NA, nrow = 100, ncol = 2)
for (i in 1:100) {
  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 (Alpha vs Lambda)

# Extract 2x2 submatrix for alpha and lambda
vcov_2d_al <- vcov_full[c(1, 3), c(1, 3)]

# Create confidence ellipse
eig_decomp_al <- eigen(vcov_2d_al)
ellipse_al <- matrix(NA, nrow = 100, ncol = 2)
for (i in 1:100) {
  v <- c(cos(theta[i]), sin(theta[i]))
  ellipse_al[i, ] <- mle[c(1, 3)] + sqrt(chi2_val) *
    (eig_decomp_al$vectors %*% diag(sqrt(eig_decomp_al$values)) %*% v)
}

# Marginal confidence intervals
se_2d_al <- sqrt(diag(vcov_2d_al))
ci_alpha_2 <- mle[1] + c(-1, 1) * 1.96 * se_2d_al[1]
ci_lambda <- mle[3] + c(-1, 1) * 1.96 * se_2d_al[2]

# Plot

plot(ellipse_al[, 1], ellipse_al[, 2],
  type = "l", lwd = 2, col = "#2E4057",
  xlab = expression(alpha), ylab = expression(lambda),
  main = "95% Confidence Region (Alpha vs Lambda)", las = 1
)

# Add marginal CIs
abline(v = ci_alpha_2, col = "#808080", lty = 3, lwd = 1.5)
abline(h = ci_lambda, col = "#808080", lty = 3, lwd = 1.5)

points(mle[1], mle[3], pch = 19, col = "#8B0000", cex = 1.5)
points(true_params[1], true_params[3], 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 8: Confidence Ellipse (Beta vs Lambda)

# Extract 2x2 submatrix for beta and lambda
vcov_2d_bl <- vcov_full[2:3, 2:3]

# Create confidence ellipse
eig_decomp_bl <- eigen(vcov_2d_bl)
ellipse_bl <- matrix(NA, nrow = 100, ncol = 2)
for (i in 1:100) {
  v <- c(cos(theta[i]), sin(theta[i]))
  ellipse_bl[i, ] <- mle[2:3] + sqrt(chi2_val) *
    (eig_decomp_bl$vectors %*% diag(sqrt(eig_decomp_bl$values)) %*% v)
}

# Marginal confidence intervals
se_2d_bl <- sqrt(diag(vcov_2d_bl))
ci_beta_2 <- mle[2] + c(-1, 1) * 1.96 * se_2d_bl[1]
ci_lambda_2 <- mle[3] + c(-1, 1) * 1.96 * se_2d_bl[2]

# Plot

plot(ellipse_bl[, 1], ellipse_bl[, 2],
  type = "l", lwd = 2, col = "#2E4057",
  xlab = expression(beta), ylab = expression(lambda),
  main = "95% Confidence Region (Beta vs Lambda)", las = 1
)

# Add marginal CIs
abline(v = ci_beta_2, col = "#808080", lty = 3, lwd = 1.5)
abline(h = ci_lambda_2, col = "#808080", lty = 3, lwd = 1.5)

points(mle[2], mle[3], pch = 19, col = "#8B0000", cex = 1.5)
points(true_params[2], true_params[3], 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")

---