On the Statistical Properties and Computational Inference of the Generalized Kumaraswamy Distribution Family

1.1 Motivation and Background

The analysis of continuous random variables constrained to the unit interval \((0,1)\) arises naturally in numerous statistical applications, including proportions, rates, percentages, and index measurements. The classical Beta distribution (Johnson et al., 1995) has long served as the canonical model for such data, offering analytical tractability and well-understood properties. However, its cumulative distribution function (CDF) involves the incomplete beta function, requiring numerical evaluation of special functions for quantile computation and simulation.

Kumaraswamy (1980) introduced an alternative two-parameter family with closed-form CDF and quantile function, facilitating computational efficiency while maintaining comparable flexibility to the Beta distribution. Jones (2009) demonstrated that the Kumaraswamy distribution exhibits similar shape characteristics to the Beta family while offering superior computational advantages.

Building upon these foundations, Carrasco, Ferrari and Cordeiro (2010) introduced the five-parameter Generalized Kumaraswamy (GKw) distribution by combining the Kumaraswamy CDF with a generalized-beta generator. Cordeiro and de Castro (2011) developed the related Kumaraswamy-generated (Kw-G) family, which is conceptually connected but is not the five-parameter GKw model studied here. The GKw distribution encompasses a rich hierarchy of submodels and provides substantial flexibility for modeling diverse patterns in bounded data. Identifiability issues for this class were revisited explicitly by Carrasco and Cordeiro (2017), so likelihood asymptotics must be stated conditionally on a regular identifiable parameter region.

Despite its theoretical appeal, a fully explicit and internally consistent analytical treatment of the GKw family—particularly for likelihood-based inference—has remained incomplete in the literature. This vignette fills this gap by providing a rigorous development, including validated expressions for all first and second derivatives of the log-likelihood function, written in a form convenient for implementation in the gkwdist R package.

1.2 Mathematical Preliminaries

We establish notation and fundamental results required for subsequent development.

Notation 1.1. Throughout, we denote

• \(\Gamma(\cdot)\): the gamma function
• \(B(a,b) = \dfrac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\): the beta function
• \(I_z(a,b) = \dfrac{B_z(a,b)}{B(a,b)}\): the regularized incomplete beta function
• \(\psi(x) = \Gamma'(x)/\Gamma(x) = (\ln\Gamma(x))'\): the digamma function
• \(\psi_1(x) = \psi'(x) = (\ln\Gamma(x))''\): the trigamma function
• \(\mathbf{1}_A\): the indicator function of a set \(A\)

We recall basic derivatives of the beta function.

Lemma 1.1 (Derivatives of the beta function). For \(a,b>0\),

\[ \begin{align} \frac{\partial}{\partial a}\ln B(a,b) &= \psi(a) - \psi(a+b), \tag{1.1}\\[3pt] \frac{\partial^2}{\partial a^2}\ln B(a,b) &= \psi_1(a) - \psi_1(a+b), \tag{1.2}\\[3pt] \frac{\partial^2}{\partial a\,\partial b}\ln B(a,b) &= -\psi_1(a+b). \tag{1.3} \end{align} \]

Proof. Since \[ \ln B(a,b) = \ln\Gamma(a) + \ln\Gamma(b) - \ln\Gamma(a+b), \] the identities follow immediately from the definitions of \(\psi\) and \(\psi_1\) and the chain rule. \(\square\)

We will also repeatedly use the following cascade of transformations.

Lemma 1.2 (Cascade transformations). Define, for \(x\in(0,1)\),

\[ \begin{align} v(x; \alpha) &= 1 - x^\alpha, \tag{1.4}\\ w(x; \alpha, \beta) &= 1 - v(x;\alpha)^\beta = 1 - (1-x^\alpha)^\beta, \tag{1.5}\\ z(x; \alpha, \beta, \lambda) &= 1 - w(x;\alpha,\beta)^\lambda = 1 - [1-(1-x^\alpha)^\beta]^\lambda. \tag{1.6} \end{align} \]

Then, for \(\alpha,\beta,\lambda>0\),

\[ \begin{align} \frac{\partial v}{\partial x} &= -\alpha x^{\alpha-1}, \tag{1.7}\\[3pt] \frac{\partial w}{\partial x} &= \alpha\beta x^{\alpha-1}(1-x^\alpha)^{\beta-1}, \tag{1.8}\\[3pt] \frac{\partial z}{\partial x} &= -\alpha\beta\lambda\,x^{\alpha-1} (1-x^\alpha)^{\beta-1}\bigl[1-(1-x^\alpha)^\beta\bigr]^{\lambda-1}. \tag{1.9} \end{align} \]

Proof. Direct differentiation and repeated application of the chain rule. \(\square\)

For brevity we will often write \(v(x)\), \(w(x)\) and \(z(x)\) when the dependence on \((\alpha,\beta,\lambda)\) is clear from the context.

---

2.1 Definition and Fundamental Properties

We start from the five-parameter Generalized Kumaraswamy family.

Definition 2.1 (Generalized Kumaraswamy distribution).
A random variable \(X\) has a Generalized Kumaraswamy distribution with parameter vector \[ \boldsymbol{\theta} = (\alpha,\beta,\gamma,\delta,\lambda)^\top, \] denoted \(X \sim \mathrm{GKw}(\alpha,\beta,\gamma,\delta,\lambda)\), if its probability density function (pdf) is \[ \boxed{ f(x; \boldsymbol{\theta}) = \frac{\lambda\alpha\beta}{B(\gamma,\delta+1)}\; x^{\alpha-1} v(x)^{\beta-1} w(x)^{\gamma\lambda-1} z(x)^{\delta} \,\mathbf{1}_{(0,1)}(x), } \tag{2.1} \] where \[ v(x) = 1-x^\alpha,\qquad w(x) = 1-(1-x^\alpha)^\beta,\qquad z(x) = 1-w(x)^\lambda, \] and the extended parameter space used by the package is \[ \overline{\Theta} = \Bigl\{(\alpha,\beta,\gamma,\delta,\lambda)^\top : \alpha,\beta,\gamma,\lambda>0,\ \delta\ge 0\Bigr\}. \tag{2.2} \] Its interior with respect to the five-dimensional parameterization is \[ \Theta^\circ=\{\boldsymbol\theta:\alpha,\beta,\gamma,\delta,\lambda>0\}. \] Although \(B(\gamma,\delta+1)\) is well-defined for \(\delta>-1\), the package and the submodel hierarchy considered here use \(\delta\ge0\). Boundary points and certain interior manifolds require special care because the full five-parameter representation is not everywhere identifiable.

Proposition 2.1 (Exact non-identifiability manifolds).
For the full five-parameter GKw representation, at least the following parameter manifolds are non-identifiable:

1. If \(\delta=0\), then the density depends on \((\gamma,\lambda)\) only through \(\tau=\gamma\lambda\): \[ f(x)=\tau\alpha\beta\,x^{\alpha-1}v(x)^{\beta-1}w(x)^{\tau-1}. \tag{2.2a} \]
2. If \(\beta=1\), then the density depends on \((\alpha,\lambda)\) only through \(\kappa=\alpha\lambda\): \[ f(x)=\frac{\kappa}{B(\gamma,\delta+1)} x^{\kappa\gamma-1}(1-x^{\kappa})^{\delta}. \tag{2.2b} \]
3. If \(\gamma=1\) and \(\lambda=1\), then the density depends on \((\beta,\delta)\) only through \(\xi=\beta(\delta+1)\): \[ f(x)=\alpha\xi\,x^{\alpha-1}(1-x^\alpha)^{\xi-1}. \tag{2.2c} \] Thus the full five-parameter embedding is not globally identifiable, even though reduced submodels can be identifiable when expressed in their own minimal parameterizations.

Proof. If \(\delta=0\), then \(B(\gamma,1)=1/\gamma\) and \(z(x)^0=1\), so (2.1) reduces immediately to (2.2a). Hence any two pairs \((\gamma,\lambda)\) with the same product \(\gamma\lambda\) generate the same density. If \(\beta=1\), then \(w(x)=x^\alpha\), and (2.1) becomes (2.2b), which depends on \(\alpha\) and \(\lambda\) only through \(\alpha\lambda\). Finally, if \(\gamma=\lambda=1\), then \(B(1,\delta+1)=1/(\delta+1)\), \(w^{\gamma\lambda-1}=1\), and \(z=1-w=v^\beta\). Therefore \[ f(x)=\alpha\beta(\delta+1)x^{\alpha-1}v(x)^{\beta-1+\beta\delta} =\alpha\xi x^{\alpha-1}v(x)^{\xi-1}, \] which is (2.2c). Hence distinct parameter vectors can induce exactly the same distribution on each of these manifolds. Whenever the usual differentiability conditions hold, the Fisher information for the full parameterization is necessarily singular there. \(\square\)

We now verify that (2.1) defines a proper density.

Theorem 2.1 (Validity of the pdf).
For any \(\boldsymbol{\theta} \in \overline{\Theta}\), the function \(f(\cdot; \boldsymbol{\theta})\) in (2.1) is a valid probability density on \((0,1)\).

Proof. Non-negativity is immediate from the definition. To prove normalization, consider the change of variable \[ u = w(x)^\lambda,\qquad 0<u<1. \] From Lemma 1.2 and the chain rule, \[ \frac{du}{dx} = \lambda w(x)^{\lambda-1}\frac{\partial w(x)}{\partial x} = \lambda \alpha\beta\,x^{\alpha-1} v(x)^{\beta-1} w(x)^{\lambda-1}. \] Hence \[ dx = \frac{du}{\lambda\alpha\beta\,x^{\alpha-1} v(x)^{\beta-1} w(x)^{\lambda-1}}. \]

Substituting into the integral of \(f\), \[ \begin{align} \int_0^1 f(x;\boldsymbol{\theta})\,dx &= \frac{\lambda\alpha\beta}{B(\gamma,\delta+1)} \int_0^1 x^{\alpha-1}v(x)^{\beta-1}w(x)^{\gamma\lambda-1}z(x)^\delta\,dx\\ &= \frac{\lambda\alpha\beta}{B(\gamma,\delta+1)} \int_0^1 x^{\alpha-1}v(x)^{\beta-1}w(x)^{\gamma\lambda-\lambda}w(x)^{\lambda-1}z(x)^\delta\,dx\\ &= \frac{\lambda\alpha\beta}{B(\gamma,\delta+1)} \int_0^1 w(x)^{\lambda(\gamma-1)} z(x)^\delta\, \underbrace{x^{\alpha-1}v(x)^{\beta-1}w(x)^{\lambda-1} dx}_{du / (\lambda\alpha\beta)}\\ &= \frac{1}{B(\gamma,\delta+1)} \int_0^1 u^{\gamma-1}(1-u)^\delta\,du\\ &= \frac{B(\gamma,\delta+1)}{B(\gamma,\delta+1)}=1, \end{align} \tag{2.3} \] because \(w(x)^{\lambda(\gamma-1)} = (w(x)^\lambda)^{\gamma-1} = u^{\gamma-1}\) and \(z(x)=1-w(x)^\lambda = 1-u\). \(\square\)

The same change of variable yields the CDF.

Theorem 2.2 (Cumulative distribution function).
If \(X\sim \mathrm{GKw}(\boldsymbol{\theta})\), then for \(x\in(0,1)\), \[ \boxed{ F(x; \boldsymbol{\theta}) = I_{w(x)^\lambda}(\gamma,\delta+1) = I_{[1-(1-x^\alpha)^\beta]^\lambda}(\gamma,\delta+1). } \tag{2.4} \]

Proof. From the same substitution as in (2.3), \[ \begin{align} F(x) &= \int_0^x f(t;\boldsymbol{\theta})\,dt\\ &= \frac{1}{B(\gamma,\delta+1)} \int_0^{w(x)^\lambda} u^{\gamma-1}(1-u)^\delta\,du\\ &= I_{w(x)^\lambda}(\gamma,\delta+1). \end{align} \] At the endpoints, \(w(0^+)^\lambda = 0\) and \(w(1^-)^\lambda = 1\), so \(F(0)=0\) and \(F(1)=1\). \(\square\)

Corollary 2.1 (Quantile function and random generation).
Let \(0<p<1\) and let \[ y_p = I^{-1}_p(\gamma,\delta+1) \] denote the inverse, with respect to its argument, of the regularized incomplete beta function. Then \[ \boxed{ Q_{\mathrm{GKw}}(p) = \left\{1-\left[1-y_p^{1/\lambda}\right]^{1/\beta}\right\}^{1/\alpha}. } \tag{2.4a} \] Consequently, if \(Y\sim\mathrm{Beta}(\gamma,\delta+1)\), then \[ X=\left\{1-\left[1-Y^{1/\lambda}\right]^{1/\beta}\right\}^{1/\alpha} \tag{2.4b} \] has the GKw distribution.

Proof. Set \(F(x)=p\) in (2.4). Then \(w(x)^\lambda=y_p\), hence \(w(x)=y_p^{1/\lambda}\), \(v(x)^\beta=1-y_p^{1/\lambda}\), and inversion of \(v(x)=1-x^\alpha\) yields (2.4a). Equation (2.4b) follows by the probability integral transform. \(\square\)

2.2 The Hierarchical Structure

The GKw family exhibits a rich nested structure. Several well-known bounded distributions, or package-restricted portions of their usual parameter spaces, arise as particular choices (and mild reparameterizations) of \(\boldsymbol{\theta}\).

For a data point \(x\in(0,1)\) we will often write \[ v = 1-x^\alpha,\quad w = 1-(1-x^\alpha)^\beta,\quad z = 1-w^\lambda. \]

3.1 The Log-Likelihood Function

Let \(\mathbf{X} = (X_1,\dots,X_n)^\top\) be an i.i.d. sample from \(\mathrm{GKw}(\boldsymbol{\theta})\), with observed values \(\mathbf{x}=(x_1,\dots,x_n)^\top\). For each \(i\), define \[ v_i = v(x_i),\quad w_i = w(x_i),\quad z_i = z(x_i). \]

Definition 3.1 (Log-likelihood function).
The log-likelihood is \[ \ell(\boldsymbol{\theta};\mathbf{x}) = \sum_{i=1}^n \ln f(x_i;\boldsymbol{\theta}). \tag{3.1} \]

Theorem 3.1 (Decomposition of the log-likelihood).
The log-likelihood can be written as \[ \boxed{ \ell(\boldsymbol{\theta}) = n\ln(\lambda\alpha\beta) - n\ln B(\gamma,\delta+1) + \sum_{i=1}^n S_i(\boldsymbol{\theta}), } \tag{3.2} \] where \[ S_i(\boldsymbol{\theta}) = (\alpha-1)\ln x_i + (\beta-1)\ln v_i + (\gamma\lambda-1)\ln w_i + \delta\ln z_i. \tag{3.3} \]

Equivalently, \[ \ell(\boldsymbol{\theta}) = L_1 + L_2 + L_3 + L_4 + L_5 + L_6 + L_7 + L_8, \tag{3.4} \] where \[ \begin{align} L_1 &= n\ln\lambda, \tag{3.5}\\ L_2 &= n\ln\alpha, \tag{3.6}\\ L_3 &= n\ln\beta, \tag{3.7}\\ L_4 &= -n\ln B(\gamma,\delta+1), \tag{3.8}\\ L_5 &= (\alpha-1)\sum_{i=1}^n \ln x_i, \tag{3.9}\\ L_6 &= (\beta-1)\sum_{i=1}^n \ln v_i, \tag{3.10}\\ L_7 &= (\gamma\lambda-1)\sum_{i=1}^n \ln w_i, \tag{3.11}\\ L_8 &= \delta\sum_{i=1}^n \ln z_i. \tag{3.12} \end{align} \]

Proof. Take logarithms of (2.1) and sum over \(i\). \(\square\)

3.2 Maximum Likelihood Estimation

Definition 3.2 (Maximum likelihood estimator).
Whenever the likelihood maximum is attained, an MLE is any maximizer \[ \hat{\boldsymbol{\theta}}_n \in \operatorname*{arg\,max}_{\boldsymbol{\theta}\in\overline{\Theta}} \ell(\boldsymbol{\theta};\mathbf{x}). \tag{3.13} \] The set of maximizers need not be a singleton on non-identifiable manifolds. In a regular identifiable parameterization, the local maximizer considered below is assumed to be well defined with probability tending to one.

We will refer to the true parameter value as \(\boldsymbol{\theta}_0\in\overline{\Theta}\).

Theorem 3.2 (Consistency and asymptotic normality on a regular identifiable region).
Let \(\boldsymbol\theta_0\) belong to an open identifiable parameter region \(\Theta_R\subset\Theta^\circ\). Assume:

1. \(\boldsymbol\theta_0\) is the unique maximizer of \(M(\boldsymbol\theta)=\mathbb E_{\boldsymbol\theta_0}[\log f(X;\boldsymbol\theta)]\) on the parameter region considered, and the normalized log-likelihood obeys a uniform law of large numbers sufficient for M-estimator consistency;
2. \(\log f(x;\boldsymbol\theta)\) is twice continuously differentiable in a neighborhood of \(\boldsymbol\theta_0\), with domination conditions allowing differentiation under the integral sign and a law of large numbers for the Hessian;
3. the per-observation score \(U_1(\boldsymbol\theta)=\nabla_{\boldsymbol\theta}\log f(X;\boldsymbol\theta)\) has finite second moment and \[ \mathcal I(\boldsymbol\theta_0) =\mathbb E_{\boldsymbol\theta_0} \left[U_1(\boldsymbol\theta_0)U_1(\boldsymbol\theta_0)^\top\right] =-\mathbb E_{\boldsymbol\theta_0} \left[\nabla^2\log f(X;\boldsymbol\theta_0)\right] \] is positive definite.

Then any sequence of MLEs in this regular region satisfying the usual approximate-maximization condition is consistent, \[ \hat{\boldsymbol\theta}_n\xrightarrow{p}\boldsymbol\theta_0, \tag{3.14} \] and \[ \sqrt n(\hat{\boldsymbol\theta}_n-\boldsymbol\theta_0) \xrightarrow{d} \mathcal N_5\left(\mathbf 0,\mathcal I(\boldsymbol\theta_0)^{-1}\right). \tag{3.15} \]

Proof. The first assertion follows from the standard consistency theorem for M-estimators because \(n^{-1}\ell_n(\boldsymbol\theta)\) converges uniformly to \(M(\boldsymbol\theta)\), whose maximizer is unique at \(\boldsymbol\theta_0\). For asymptotic normality, expand the score around \(\boldsymbol\theta_0\): \[ \mathbf 0 =U_n(\hat{\boldsymbol\theta}_n) =U_n(\boldsymbol\theta_0) +H_n(\tilde{\boldsymbol\theta}_n) (\hat{\boldsymbol\theta}_n-\boldsymbol\theta_0), \] where \(\tilde{\boldsymbol\theta}_n\) lies between \(\hat{\boldsymbol\theta}_n\) and \(\boldsymbol\theta_0\). By the multivariate central limit theorem, \(n^{-1/2}U_n(\boldsymbol\theta_0)\Rightarrow\mathcal N_5(0,\mathcal I(\boldsymbol\theta_0))\), while consistency and the Hessian law of large numbers give \(-n^{-1}H_n(\tilde{\boldsymbol\theta}_n)\xrightarrow{p}\mathcal I(\boldsymbol\theta_0)\). Slutsky’s theorem yields (3.15). \(\square\)

Remark 3.1 (Why the qualification is essential).
Theorem 3.2 does not apply to the full five-parameter representation on the non-identifiable manifolds in Proposition 2.1. In particular, \(\delta=0\) makes \((\gamma,\lambda)\) confounded through \(\gamma\lambda\); \(\beta=1\) makes \((\alpha,\lambda)\) confounded through \(\alpha\lambda\); and \(\gamma=\lambda=1\) makes \((\beta,\delta)\) confounded through \(\beta(\delta+1)\). At such points the information matrix for the full embedding is singular and the usual five-dimensional Gaussian limit cannot hold. A reduced model can nevertheless possess regular asymptotics when it is written in an identifiable minimal parameterization.

3.3 Likelihood Ratio Tests

For nested models, likelihood ratio tests follow Wilks’ theorem.

Theorem 3.3 (Wilks’ theorem on a regular interior null).
Consider nested hypotheses \(H_0:\ \boldsymbol\theta\in\Theta_0\) versus \(H_1:\ \boldsymbol\theta\in\Theta_R\), where locally \(\Theta_0\) is a smooth identifiable submanifold of codimension \(r\), the true parameter is an interior point relative to the regular parameterization, and the Fisher information is nonsingular. Let \[ \Lambda_n =2\{\ell(\hat{\boldsymbol\theta}_n)-\ell(\hat{\boldsymbol\theta}_{0,n})\}. \tag{3.16} \] Then, under \(H_0\), \[ \Lambda_n\xrightarrow{d}\chi_r^2. \tag{3.17} \]

Proof. This is the standard Wilks expansion of the log-likelihood under local asymptotic normality and a regular interior constraint. The assumptions exclude boundary points and loss of identifiability. \(\square\)

Corollary 3.1 (LR tests within the GKw hierarchy).
The following regular reductions have the usual chi-squared limit provided the null point lies in an identifiable region with nonsingular information: \[ \begin{align} \text{GKw vs. BKw:}\quad &H_0:\ \lambda=1 &&\Rightarrow\ \Lambda_n\xrightarrow{d}\chi_1^2, \tag{3.18}\\[3pt] \text{GKw vs. KKw:}\quad &H_0:\ \gamma=1 &&\Rightarrow\ \Lambda_n\xrightarrow{d}\chi_1^2, \tag{3.19}\\[3pt] \text{BKw vs. Beta:}\quad &H_0:\ \alpha=\beta=1 &&\Rightarrow\ \Lambda_n\xrightarrow{d}\chi_2^2. \tag{3.20} \end{align} \] For the first two comparisons, this excludes, in particular, null points at which the full GKw representation loses identifiability.

The comparison KKw versus Kw is different. In distribution space, the Kw null inside KKw is the singular slice \[ H_0:\ \lambda=1, \tag{3.21} \] because, when \(\gamma=1\) is already fixed by KKw and \(\lambda=1\), Proposition 2.1(3) gives \[ f_{\mathrm{KKw}}(x)=\alpha\xi x^{\alpha-1}(1-x^\alpha)^{\xi-1}, \qquad \xi=\beta(\delta+1), \] which is exactly a Kw\((\alpha,\xi)\) density for every \(\delta\ge0\). The commonly used representative restriction \(\delta=0,\lambda=1\) selects one parameterization of the same Kw null but is not its full preimage in the redundant KKw coordinates. Consequently, under \(H_0\), the split of \(\xi\) into \(\beta\) and \(\delta\) is unidentified. The null is therefore a singular loss-of-identifiability problem (and the representative \(\delta=0\) also lies on the package boundary). Wilks’ ordinary chi-squared theorem does not apply, and one should not assume a standard chi-bar-square mixture without a model-specific derivation of the local experiment. In practice, a parametric bootstrap under the fitted two-parameter Kw null is a defensible calibration strategy, but its finite-sample performance should itself be checked by simulation for the intended sample size and parameter region. Self and Liang (1987) provide a general boundary-likelihood framework; the unidentified nuisance direction present here requires additional care.

---

4.1 The Score Vector

Definition 4.1 (Score function).
The score vector is \[ U(\boldsymbol{\theta}) = \nabla_{\boldsymbol{\theta}}\ell(\boldsymbol{\theta}) = \left( \frac{\partial\ell}{\partial\alpha}, \frac{\partial\ell}{\partial\beta}, \frac{\partial\ell}{\partial\gamma}, \frac{\partial\ell}{\partial\delta}, \frac{\partial\ell}{\partial\lambda} \right)^\top. \tag{4.1} \]

For the derivatives of \(v,w,z\) with respect to the parameters we will use \[ \begin{align} \frac{\partial v_i}{\partial\alpha} &= -x_i^\alpha\ln x_i,\\[2pt] \frac{\partial w_i}{\partial\alpha} &= \beta v_i^{\beta-1} x_i^\alpha\ln x_i,\\[2pt] \frac{\partial w_i}{\partial\beta} &= -v_i^\beta \ln v_i,\\[2pt] \frac{\partial z_i}{\partial\alpha} &= -\lambda w_i^{\lambda-1}\frac{\partial w_i}{\partial\alpha},\\[2pt] \frac{\partial z_i}{\partial\beta} &= -\lambda w_i^{\lambda-1}\frac{\partial w_i}{\partial\beta},\\[2pt] \frac{\partial z_i}{\partial\lambda} &= -w_i^{\lambda}\ln w_i. \end{align} \]

Theorem 4.1 (Score components).
The components of \(U(\boldsymbol{\theta})\) are

\[ \begin{align} \frac{\partial\ell}{\partial\alpha} &= \frac{n}{\alpha} + \sum_{i=1}^n \ln x_i - \sum_{i=1}^n 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], \tag{4.2}\\[6pt] \frac{\partial\ell}{\partial\beta} &= \frac{n}{\beta} + \sum_{i=1}^n \ln v_i - \sum_{i=1}^n v_i^\beta \ln v_i \left[ \frac{\gamma\lambda-1}{w_i} - \frac{\delta\lambda w_i^{\lambda-1}}{z_i} \right], \tag{4.3}\\[6pt] \frac{\partial\ell}{\partial\gamma} &= -n\bigl[\psi(\gamma) - \psi(\gamma+\delta+1)\bigr] + \lambda\sum_{i=1}^n \ln w_i, \tag{4.4}\\[6pt] \frac{\partial\ell}{\partial\delta} &= -n\bigl[\psi(\delta+1) - \psi(\gamma+\delta+1)\bigr] + \sum_{i=1}^n \ln z_i, \tag{4.5}\\[6pt] \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}. \tag{4.6} \end{align} \]

Proof. We differentiate the decomposition (3.4)–(3.12) term by term.

(i) Derivative with respect to \(\alpha\).

From (3.6) and (3.9), \[ \frac{\partial L_2}{\partial\alpha} = \frac{n}{\alpha}, \quad \frac{\partial L_5}{\partial\alpha} = \sum_{i=1}^n \ln x_i. \] Using \(\partial v_i/\partial\alpha = -x_i^\alpha\ln x_i\), \[ \frac{\partial L_6}{\partial\alpha} = (\beta-1)\sum_{i=1}^n \frac{1}{v_i}\frac{\partial v_i}{\partial\alpha} = -(\beta-1)\sum_{i=1}^n \frac{x_i^\alpha\ln x_i}{v_i}. \] Next, \(\partial w_i/\partial\alpha = \beta v_i^{\beta-1}x_i^\alpha\ln x_i\), so \[ \frac{\partial L_7}{\partial\alpha} = (\gamma\lambda-1)\sum_{i=1}^n \frac{1}{w_i}\frac{\partial w_i}{\partial\alpha} = (\gamma\lambda-1)\beta\sum_{i=1}^n \frac{v_i^{\beta-1}x_i^\alpha\ln x_i}{w_i}. \] Similarly, \[ \frac{\partial z_i}{\partial\alpha} = -\lambda w_i^{\lambda-1}\frac{\partial w_i}{\partial\alpha} = -\lambda\beta v_i^{\beta-1}w_i^{\lambda-1}x_i^\alpha\ln x_i, \] so \[ \frac{\partial L_8}{\partial\alpha} = \delta\sum_{i=1}^n \frac{1}{z_i}\frac{\partial z_i}{\partial\alpha} = -\delta\lambda\beta\sum_{i=1}^n \frac{v_i^{\beta-1}w_i^{\lambda-1}x_i^\alpha\ln x_i}{z_i}. \] Collecting terms gives (4.2).

(ii) Derivative with respect to \(\beta\).

From (3.7) and (3.10), \[ \frac{\partial L_3}{\partial\beta} = \frac{n}{\beta}, \quad \frac{\partial L_6}{\partial\beta} = \sum_{i=1}^n \ln v_i, \] since \(v_i\) does not depend on \(\beta\). Using \(\partial w_i/\partial\beta = -v_i^\beta\ln v_i\), \[ \frac{\partial L_7}{\partial\beta} = (\gamma\lambda-1)\sum_{i=1}^n \frac{1}{w_i}\frac{\partial w_i}{\partial\beta} = -(\gamma\lambda-1)\sum_{i=1}^n \frac{v_i^\beta\ln v_i}{w_i}. \] Furthermore, \[ \frac{\partial z_i}{\partial\beta} = -\lambda w_i^{\lambda-1}\frac{\partial w_i}{\partial\beta} = \lambda w_i^{\lambda-1}v_i^\beta\ln v_i, \] so \[ \frac{\partial L_8}{\partial\beta} = \delta\sum_{i=1}^n \frac{1}{z_i}\frac{\partial z_i}{\partial\beta} = \delta\lambda\sum_{i=1}^n \frac{w_i^{\lambda-1}v_i^\beta\ln v_i}{z_i}. \] Combining terms yields (4.3).

(iii) Derivative with respect to \(\gamma\).

Only \(L_4\) and \(L_7\) depend on \(\gamma\). From Lemma 1.1, \[ \frac{\partial L_4}{\partial\gamma} = -n\bigl[\psi(\gamma) - \psi(\gamma+\delta+1)\bigr], \qquad \frac{\partial L_7}{\partial\gamma} = \lambda\sum_{i=1}^n \ln w_i, \] giving (4.4).

(iv) Derivative with respect to \(\delta\).

Similarly, \[ \frac{\partial L_4}{\partial\delta} = -n\bigl[\psi(\delta+1) - \psi(\gamma+\delta+1)\bigr], \qquad \frac{\partial L_8}{\partial\delta} = \sum_{i=1}^n \ln z_i, \] giving (4.5).

(v) Derivative with respect to \(\lambda\).

We have \[ \frac{\partial L_1}{\partial\lambda} = \frac{n}{\lambda}, \qquad \frac{\partial L_7}{\partial\lambda} = \gamma\sum_{i=1}^n \ln w_i, \] and \[ \frac{\partial z_i}{\partial\lambda} = -w_i^\lambda\ln w_i \quad\Rightarrow\quad \frac{\partial L_8}{\partial\lambda} = \delta\sum_{i=1}^n \frac{1}{z_i}\frac{\partial z_i}{\partial\lambda} = -\delta\sum_{i=1}^n \frac{w_i^\lambda\ln w_i}{z_i}. \] Together these yield (4.6). \(\square\)

Implementation sign convention.
Equations (4.2)–(4.6) are derivatives of the log-likelihood \(\ell\). The optimization helpers in gkwdist are written for the negative log-likelihood; therefore their gradient and Hessian should satisfy \[ \nabla(-\ell)=-U(\boldsymbol\theta), \qquad \nabla^2(-\ell)=-H(\boldsymbol\theta)=\mathcal J(\boldsymbol\theta). \] This distinction is essential when comparing the formulas in this vignette with gr* and hs* implementation routines.

4.2 The Hessian and Observed Information Matrix

We now consider second-order derivatives. Let \[ H(\boldsymbol{\theta}) = \nabla^2\ell(\boldsymbol{\theta}) = \biggl[ \frac{\partial^2\ell}{\partial\theta_j\partial\theta_k} \biggr]_{j,k=1}^5 \] denote the Hessian matrix of the log-likelihood, where \(\boldsymbol{\theta} = (\alpha,\beta,\gamma,\delta,\lambda)^\top\).

Definition 4.2 (Observed information).
The observed information matrix is defined as \[ \mathcal{J}(\boldsymbol{\theta}) = -H(\boldsymbol{\theta}) = -\nabla^2\ell(\boldsymbol{\theta}). \tag{4.7} \]

To keep the formulas compact, for each observation \(i\) and each transformation \(u_i\in\{v_i,w_i,z_i\}\) we define, for parameters \(\theta_j,\theta_k\), \[ D^{u}_{jk}(i) = \frac{ \dfrac{\partial^2 u_i}{\partial\theta_j\partial\theta_k}u_i - \dfrac{\partial u_i}{\partial\theta_j} \dfrac{\partial u_i}{\partial\theta_k} }{u_i^2} = \frac{\partial^2}{\partial\theta_j\partial\theta_k}\ln u_i. \]

In particular, \[ \frac{\partial^2}{\partial\theta_j^2}\ln u_i = D^u_{jj}(i). \]

For implementation, all raw second derivatives entering \(D^u_{jk}(i)\) can be written explicitly. Let \[ A_i=x_i^\alpha,\qquad L_i=\ln x_i. \] Then \[ \begin{align} v_{i,\alpha} &=-A_iL_i, &v_{i,\alpha\alpha} &=-A_iL_i^2,\tag{4.7a}\\ w_{i,\alpha} &=\beta v_i^{\beta-1}A_iL_i, &w_{i,\beta} &=-v_i^\beta\ln v_i,\tag{4.7b}\\ w_{i,\alpha\alpha} &=\beta A_iL_i^2v_i^{\beta-2}\{v_i-(\beta-1)A_i\}, &w_{i,\alpha\beta} &=A_iL_iv_i^{\beta-1}(1+\beta\ln v_i),\tag{4.7c}\\ w_{i,\beta\beta} &=-v_i^\beta(\ln v_i)^2.\tag{4.7d} \end{align} \] For \(j,k\in\{\alpha,\beta\}\), \[ \begin{align} z_{i,j} &=-\lambda w_i^{\lambda-1}w_{i,j},\tag{4.7e}\\ z_{i,jk} &=-\lambda\left[(\lambda-1)w_i^{\lambda-2}w_{i,j}w_{i,k} +w_i^{\lambda-1}w_{i,jk}\right],\tag{4.7f}\\ z_{i,\lambda} &=-w_i^\lambda\ln w_i,\tag{4.7g}\\ z_{i,j\lambda} &=-w_i^{\lambda-1}w_{i,j}(1+\lambda\ln w_i),\tag{4.7h}\\ z_{i,\lambda\lambda} &=-w_i^\lambda(\ln w_i)^2.\tag{4.7i} \end{align} \] These formulas, together with the definition of \(D^u_{jk}(i)\), give every Hessian entry below without numerical differentiation.

4.3 Asymptotic Variance–Covariance Matrix

On a regular identifiable parameter region satisfying Theorem 3.2, the asymptotic variance–covariance matrix of the MLE is governed by the Fisher information.

Let \(\mathcal{I}(\boldsymbol{\theta}_0)\) denote the per-observation Fisher information, \[ \mathcal{I}(\boldsymbol{\theta}_0) = \mathbb{E}_{\boldsymbol{\theta}_0} \left[ -\nabla^2\ln f(X;\boldsymbol{\theta}) \right]. \]

For the full sample of size \(n\), the expected information is \(n\mathcal{I}(\boldsymbol{\theta}_0)\), while the observed information is \(\mathcal{J}(\hat{\boldsymbol{\theta}}_n)\).

Theorem 4.4 (Variance–covariance matrix of the MLE).
Under the regularity assumptions of Theorem 3.2, \[ \operatorname{Var}(\hat{\boldsymbol{\theta}}_n) \approx \frac{1}{n}\mathcal{I}(\boldsymbol{\theta}_0)^{-1} \approx \mathcal{J}(\hat{\boldsymbol{\theta}}_n)^{-1}, \tag{4.27} \] and the asymptotic standard error of \(\hat\theta_j\) is approximated by \[ \widehat{\operatorname{SE}}(\hat\theta_j) = \sqrt{\bigl[\mathcal{J}(\hat{\boldsymbol{\theta}}_n)^{-1}\bigr]_{jj}}. \tag{4.28} \]

Proof. The convergence \(\mathcal{J}(\hat{\boldsymbol{\theta}}_n)/n \xrightarrow{p} \mathcal{I}(\boldsymbol{\theta}_0)\) follows from a law of large numbers for the Hessian (Cox and Hinkley, 1974). Combining this with the asymptotic normality (3.15) yields (4.27)–(4.28). \(\square\)

---

5.1 Numerical Stability

Direct evaluation of \(v_i=1-x_i^\alpha\), \(w_i=1-v_i^\beta\), and \(z_i=1-w_i^\lambda\) can suffer catastrophic cancellation when the subtracted quantity is close to one. The implementation should therefore propagate the cascade on the log scale.

Algorithm 5.1 (Stable computation of \(\log(1-e^a)\) for \(a<0\)).

For \(a<0\), define \[ \operatorname{log1mexp}(a)= \begin{cases} \log1p(-e^a), & a\le -\log 2,\\[4pt] \log(-\operatorname{expm1}(a)), & -\log 2<a<0. \end{cases} \tag{5.1} \] The first branch avoids loss of relative accuracy when \(e^a\) is small, while the second avoids cancellation when \(a\) is close to zero. Thus \[ \begin{align} \log v_i &=\operatorname{log1mexp}(\alpha\log x_i),\tag{5.2a}\\ \log w_i &=\operatorname{log1mexp}(\beta\log v_i),\tag{5.2b}\\ \log z_i &=\operatorname{log1mexp}(\lambda\log w_i).\tag{5.2c} \end{align} \] This is the stable branch structure recommended by Mächler (2012).

The same principle should be used for ratios appearing in the score and Hessian. For example, \[ \frac{x_i^\alpha}{v_i} =\exp\{\alpha\log x_i-\log v_i\},\qquad \frac{v_i^\beta}{w_i} =\exp\{\beta\log v_i-\log w_i\},\qquad \frac{w_i^\lambda}{z_i} =\exp\{\lambda\log w_i-\log z_i\}. \tag{5.2d} \] Products such as \(v_i^{\beta-1}w_i^{\lambda-1}/z_i\) should likewise be assembled from sums and differences of logged factors whenever direct powers or divisions would overflow, underflow, or magnify cancellation.

5.2 Optimization

For the regular five-parameter model it is convenient to optimize an unconstrained negative log-likelihood. Define \[ \eta=(\eta_1,\ldots,\eta_5)^\top, \qquad (\alpha,\beta,\gamma,\delta,\lambda) =(e^{\eta_1},e^{\eta_2},e^{\eta_3},e^{\eta_4},e^{\eta_5}), \tag{5.3} \] so that \(\delta>0\). Boundary submodels such as EKw or Kw should be fitted as their own reduced parameterizations rather than approximated by driving \(\delta\) numerically to zero in the non-identifiable five-parameter representation. Let \[ m(\eta)=-\ell(\theta(\eta)). \]

Algorithm 5.2 (Maximum likelihood via BFGS applied to \(-\ell\)).

1. Initialization. Obtain a finite starting vector \(\eta^{(0)}\), preferably from a simpler identifiable submodel.
2. BFGS iteration. At iteration \(k\), evaluate \(m(\eta^{(k)})\) and its exact gradient. With \(H_k\) a positive-definite approximation to the inverse Hessian of \(m\), use the descent direction \[ p_k=-H_k\nabla m(\eta^{(k)}), \] choose a step length \(\rho_k>0\) satisfying standard line-search conditions (e.g. Wolfe conditions), and update \[ \eta^{(k+1)}=\eta^{(k)}+\rho_kp_k. \] Apply the standard BFGS inverse-Hessian update to \(H_k\).
3. Convergence checks. Require a small gradient norm together with stable parameter and objective changes; inspect the Hessian conditioning and profile likelihood when parameters are weakly identified.
4. Observed information. Transform the Hessian back to the original parameterization if standard errors are reported on the \(\theta\)-scale, and compute the observed information only at a regular local maximum.

This formulation avoids the sign inconsistency that arises if a positive-definite BFGS matrix is described as an approximation to the Hessian of the maximized log-likelihood, whose Hessian is negative definite at a strict local maximum.

Theorem 5.1 (Local superlinear convergence of BFGS).
Suppose \(m\) is twice continuously differentiable in a neighborhood of a local minimizer \(\hat\eta\), \(\nabla^2m(\hat\eta)\) is positive definite, the Hessian is locally Lipschitz continuous, the BFGS iterates converge to \(\hat\eta\), exact gradients are used, and the line search satisfies the standard conditions required by BFGS convergence theory. Then the BFGS iterates are q-superlinearly convergent: \[ \frac{\|\eta^{(k+1)}-\hat\eta\|}{\|\eta^{(k)}-\hat\eta\|}\longrightarrow0. \tag{5.4} \] This is a local convergence statement; exact gradients alone do not guarantee global convergence for an arbitrary non-concave likelihood.

Justification. Under the stated smoothness, curvature, convergence and line-search hypotheses, the standard BFGS theory yields the Dennis–Moré condition for the inverse-Hessian approximations, which implies q-superlinear convergence of the iterates. This is a theorem about the local minimization problem for \(m=-\ell\), not a claim that an arbitrary starting value reaches the global MLE; see Nocedal and Wright (2006), Chapter 6. \(\square\)

5.3 Gradient Accuracy

Numerical differentiation can be used to validate the analytic derivatives but is less efficient for routine computation.

Lemma 5.1 (Finite-difference error).
Consider the central finite-difference approximation to \(\partial\ell/\partial\theta_j\) with step size \(h>0\): \[ D_h = \frac{\ell(\boldsymbol{\theta}+h\mathbf{e}_j) - \ell(\boldsymbol{\theta}-h\mathbf{e}_j)}{2h}. \] Then \[ \left| D_h - \frac{\partial\ell}{\partial\theta_j} \right| = O(h^2) + O\!\left(\frac{\epsilon}{h}\right), \tag{5.5} \] where \(\epsilon\) is machine precision (approximately \(2.22\times10^{-16}\) in double precision). Balancing truncation and roundoff errors gives the familiar cube-root scaling \(h=O(\epsilon^{1/3})\) after the objective and parameter have been appropriately scaled. There is no universal numerical value of \(h\), because the constants depend on the local third derivative and on the scale of the likelihood evaluation.

Proof. Taylor expansion about \(\boldsymbol\theta\) gives \[ \ell(\boldsymbol\theta\pm h\mathbf e_j) =\ell(\boldsymbol\theta) \pm h\ell_j +\frac{h^2}{2}\ell_{jj} \pm\frac{h^3}{6}\ell_{jjj} +O(h^4), \] so subtraction and division by \(2h\) leave a truncation error \(O(h^2)\). If each floating-point objective evaluation has absolute perturbation of order \(O(\epsilon)\) on the working scale, subtracting two such values and dividing by \(2h\) contributes \(O(\epsilon/h)\). Balancing these two leading orders gives \(h^2\asymp\epsilon/h\), hence \(h=O(\epsilon^{1/3})\) after suitable scaling. See Nocedal and Wright (2006), Chapter 8. \(\square\)

In contrast, the analytical gradients of Theorem 4.1 can be evaluated with accuracy limited essentially only by floating-point roundoff and can be accumulated in one \(O(n)\) pass through the data, whereas a central-difference approximation to the full \(p\)-dimensional gradient requires \(2p\) likelihood evaluations (here \(p=5\)).

5.4 Practical Recommendations

Guideline 5.1 (Model selection within the GKw hierarchy).

1. Begin with simpler identifiable models such as the embedded Beta subfamily or the Kumaraswamy model when they are scientifically plausible.
2. Consider EKw or McDonald-type extensions when one additional shape mechanism is needed.
3. Use BKw or KKw when the additional flexibility is supported by the data and by stable likelihood geometry.
4. Fit the full five-parameter GKw model only when the likelihood is well identified: inspect profile likelihoods, the observed-information spectrum, sensitivity to starting values, and parametric-bootstrap behavior. There is no distribution-free sample-size cutoff such as \(n\ge500\) that guarantees identifiability or numerical stability.
5. Compare candidate models with likelihood criteria such as \[ \mathrm{AIC}=-2\ell(\hat{\boldsymbol\theta})+2p, \qquad \mathrm{BIC}=-2\ell(\hat{\boldsymbol\theta})+p\log n, \] while remembering that ordinary likelihood-ratio reference distributions can fail for boundary or non-identifiable reductions.

Guideline 5.2 (Diagnostics).

1. Q-Q plot. Use the quantile representation (2.4a), with the inverse beta function evaluated numerically.
2. Probability integral transform. The fitted values \(F(x_i;\hat{\boldsymbol\theta})\) should resemble a Uniform\((0,1)\) sample, but they are not exactly i.i.d. Uniform in finite samples because the same data were used to estimate \(\hat{\boldsymbol\theta}\). Formal calibration should therefore use a parametric bootstrap when needed.
3. Information conditioning. Examine the eigenvalues or a suitably scaled condition number of \(\mathcal J(\hat{\boldsymbol\theta})\). The raw condition number is parameterization-dependent, so fixed thresholds should be treated only as numerical warnings rather than universal inferential rules.
4. Positive definiteness. At a strict regular local maximum, the Hessian of \(\ell\) is negative definite and the observed information is positive definite. Failure of this property can indicate a saddle point, a flat/non-identifiable direction, a boundary optimum, or numerical error.
5. Profile likelihood. For a five-parameter model, profile plots are especially useful for detecting ridges involving \((\alpha,\lambda)\), \((\gamma,\lambda)\), or \((\beta,\delta)\) near the non-identifiable manifolds in Proposition 2.1.

5.5 Discussion

We have developed a rigorous mathematical framework for the Generalized Kumaraswamy (GKw) family, including:

1. Hierarchical embedding.
   The GKw family contains the McDonald, Kumaraswamy, exponentiated Kumaraswamy, Beta–Kumaraswamy and Kumaraswamy–Kumaraswamy distributions as submodels, together with the shifted Beta subfamily described in Theorem 2.8, with explicit parameter mappings.

2. Likelihood theory.
   We derived explicit expressions for the log-likelihood, the score vector and the full observed information matrix in terms of the cascade transformations \(v,w,z\), in a form suitable for stable numerical implementation.

3. Likelihood asymptotics and identifiability.
   Standard MLE asymptotics hold only on regular identifiable regions with nonsingular Fisher information. The full five-parameter representation is non-identifiable on important manifolds including \(\delta=0\), \(\beta=1\), and \(\gamma=\lambda=1\), so boundary and reduced models must be handled separately.

4. Computational considerations.
   Log-scale evaluations and carefully structured derivatives provide numerical stability and efficiency. In our C++ implementation via RcppArmadillo, analytical gradients and Hessians yield substantial speedups over finite-difference approximations, together with better numerical accuracy.

Natural extensions and computational developments include:

• Compact and numerically efficient representations for the full GKw moments. Carrasco et al. (2010) provide series representations; simple elementary closed forms are generally unavailable.
• Higher-order likelihood corrections and profile-likelihood inference in weakly identified regions.
• Boundary-aware likelihood-ratio calibration and parametric-bootstrap procedures for comparisons such as KKw versus Kw.
• Multivariate generalizations using copulas constructed from GKw marginals.
• Fully Bayesian treatments with priors and parameterizations designed to avoid or regularize non-identifiable directions.

The full GKw quantile is already available computationally through the inverse regularized beta function followed by the explicit cascade inversion in (2.4a); no elementary inverse-beta formula is required for simulation or Q-Q diagnostics.

The gkwdist R package provides the distributional and likelihood-computation routines underlying this vignette. The identifiability, singular-null and numerical-stability qualifications established here should be treated as part of the mathematical contract for interpreting those routines and for future package-level inference helpers.

---