The Lambert-F Distributions Class: An Alternative Family for Positive Data Analysis

In this article, we introduce a new probability distribution generator called the Lambert-F generator. For any continuous baseline distribution F, with positive support, the corresponding Lambert-F version is generated by using the new generator. The result is a new class of distributions with one extra parameter that generalizes the baseline distribution and whose quantile function can be expressed in closed form in terms of the Lambert W function. The hazard rate function of a Lambert-F distribution corresponds to a modification of the baseline hazard rate function, greatly increasing or decreasing the baseline hazard rate for earlier times. Herein, we study the main structural properties of the new class of distributions. Special attention is given to two particular cases that can be understood as two-parameter extensions of the well-known exponential and Rayleigh distributions. We discuss parameter estimation for the proposed models considering the moments and maximum likelihood methods. Finally, two applications were developed to illustrate the usefulness of the proposed distributions in the analysis of data from different real settings.


Introduction
In recent years, towards more flexibility, many studies have been developed in which new methods are proposed to add one or more parameters to a baseline probability distribution. These methods have given way to the generation of models with more complex parametric structures and more flexibility in aspects such as the shapes of the density and hazard rate functions and the asymmetry or kurtosis of the distribution.
One of the most popular methods is to propose a new cumulative distribution function (cdf) considering a suitable transformation of the cdf of a certain random variable of interest. More specifically, if X is a random variable with cdf F(x), then where r(·) is a nondecreasing function such that r(u) → 0 as u → 0 + , and r(u) → 1 as u → 1 − , for u ∈ (0, 1), is a cdf. Different analytical expressions for r(u) can be found in the literature. If r(u) = u α , α > 0, then the cdf in Equation (1) corresponds to the cdf of the exponentiated distributions class; see Gupta et al. [1], Nadarajah and Kotz [2], Al-hussaini [3], Castillo et al. [4] and Gómez-Déniz et al. [5], among others. If r(u) = I u (α, β), where α, β > 0 and I u (·, ·) denotes the incomplete beta function ratio,

Lambert-F Distribution Generator
Definition 1. A random variable X follows a Lambert-F distribution, denoted as X ∼ LF(η, α), if its cdf is given by where α ∈ (0, e) is an extra parameter and F(x; η) is the cdf of a continuous and positive baseline random variable with parameter vector η.
Note that the function F X (x; η, α) inherits the support of the baseline distribution F(x; η) and that F X (x; η, α) = F(x; η) when α = 1. We refer to the cdf presented in Equation (2) as the Lambert-F distribution generator. From now on, we will denote F(x; η) = F(x) and F X (x; η, α) = F X (x) to simplify the notation. Proposition 1. Let X ∼ LF(η, α). Then, the probability density function (pdf), the survival function (sf) and the hazard rate function (hrf) of X are given, respectively, by S X (x) = S(x)α F(x) and (4) where f (x), F(x), S(x) = 1 − F(x) and h(x) = f (x)/S(x) are the pdf, the cdf, the sf and the hrf of the baseline distribution, respectively.
From the above results, we note that the extra parameter α allows the pdf and the hrf of the Lambert-F distributions to present a wider range of shapes than those of the baseline distributions. These shapes are discussed in detail in the following section.

Properties
In this section, we study the pdf and hrf shapes of the Lambert-F distributions class, discuss the stochastic ordering of Lambert-F random variables and derive analytical expressions for the quantile function and the raw moments of the Lambert-F distributions class.

Shapes
The shapes of the pdf and the hrf presented in Corollary 1 can be described analytically. First of all, we see that: Secondly, under the assumption that f (x) and f (x) exist, the critical points of the pdf of X are the roots of the equation A root x = x 0 of Equation (6) will be a local maximum, a local minimum or an inflection point in the cases ϕ(x 0 ) < 0, ϕ(x 0 ) > 0 or ϕ(x 0 ) = 0, respectively, where The critical points of the hrf of X are the roots of the equation A root x = x 0 of Equation (7) will be a local maximum, a local minimum or an inflection point in the cases τ( While giving even greater attention to the hrf of the Lambert-F distributions, we point out the following: 1. The hrf of the Lambert-F distribution approximates the hrf of the baseline distribution F when x is large enough, that is, 2. The hrf of the Lambert-F distribution is greater than the hrf of the baseline distribution F if and only if α ∈ (1, e). 3. If α ∈ (1, e) and the hrf of the baseline distribution F is a nondecreasing monotonic function, then the hrf of the Lambert-F distribution is a nondecreasing monotonic function.
In view of the above, the application of the proposed Lambert transformation to a baseline distribution has a simple justification in terms of the hrf of the resulting model. The hrf induced by the Lambert-F generator is more distant from the baseline hrf for lower values of X (earlier times), as can be seen in Figures 1 and 2.
It should be taken into account that Equations (6) and (7) allow a detailed description of the shapes of the pdf and the hrf of a Lambert-F distribution once the definition functions of the baseline distribution are specified. The non-existence or existence of one or more solutions for these equations will depend jointly on the performance of the parameter α and on the analytical structure of the specified baseline distribution.

Stochastic Order
The ordering of continuous positive random variables is an important tool for judging comparative behavior. It is well known that a random variable X is smaller than a random variable Y in the stochastic are well known Shaked and Shanthikumar [20].
Proof. First, notice that is a nondecreasing function if and only if µ (x) ≥ 0, where Some calculations show that The remaining affirmations were derived from the implications in Equation (8).
A direct consequence of Proposition 2 is that the hrf of the Lambert-F distribution is less than the hrf of the baseline F distribution when α ∈ (1, e), which is consistent with the observation Equation (2) in Section 3.1.

Quantile Function
Proposition 3. Let X ∼ LF(η, α). Then, the quantile function of X is given by where F −1 (·; η) is the quantile function of the baseline distribution and W 0 (·) is the principal branch of the Lambert W function.
Proof. From Equation (2) it can be observed that the baseline cdf can be written as Then, solving this equation with respect to x, we obtain the expression in Equation (9) for α = 1. In the special case α = 1, the result is obtained directly by calculating the inverse function of the baseline cdf.
Proof. From the Lambert-F pdf in Corollary 1, we have that and by applying the change of variable u = F(x), the result is obtained.

Two Special Cases
In what follows, two new two-parameter models generated from the results in Corollary 1 are introduced. The well-known exponential and Rayleigh distributions are taken as baseline distributions in the generation of these new models.
Lambert-exponential model. The random variable X follows the Lambert-exponential distribution with scale parameter σ > 0, denoted as X ∼ LE(σ, α), if its pdf and hrf for x > 0 are given, respectively, by Lambert-Rayleigh model. The random variable X follows the Lambert-Rayleigh distribution with scale parameter σ > 0, denoted as X ∼ LR(σ, α), if its pdf and hrf for x > 0 are given, respectively, by The new two-parameter distributions described above are members of the well-known and popular shape and scale distributions family. The scale σ in both models is inherited from the respective baseline distribution, while the shape parameter α arises from the application of the Lambert transformation to the baseline distribution. Figures 1 and 2 display some plots of the pdfs and the hrfs of the above models for σ = 1 and different values of shape parameter.

Characterizing the LE and LR Distributions
This section describes the main structural properties of the LE and LR distributions.

Description and Comparison of Shapes
From the results in Section 3.1, it is possible to analytically describe the shapes of the pdf and the hrf of the above models. In the sequel, the pdf and the hrf of the LE and LR models can be analytically described.
, the pdf of X is a monotonically decreasing function when α ∈ (0, a] or α = 1, and it is a unimodal function for α ∈ (a, 1) or α ∈ (1, e). The mode of X is given by The hrf of X is a monotonically decreasing function when α ∈ (0, 1), a monotonically increasing function when α ∈ (1, e) and a constant function for α = 1.
The pdf of X is a unimodal function for α ∈ (0, e) (Mode without explicit analytic expression).
(c) The hrf of X is a monotonically increasing function when α ∈ [0.1, e) and presents the increasing-decreasing-increasing shape for α ∈ (0, 0.1). The local maximum and the local minimum are given by where W 0 (·) and W −1 (·) denote the principal and non-principal branches of the Lambert W function, respectively.
The shapes of the pdf and the hrf for the LE model are similar to those presented by other two-parameter models, such as Weibull (W) and gamma (G). It is important to note that the LE, W and G models are two-parameter extensions of the exponential model, so the comparison between these models is quite natural. A similar observation can be made for the LR and W models because both are two-parameter extensions of the Rayleigh model.
In Tables 1 and 2, we present a comparison of the shapes of density and hazard rate functions of the LE and LR models with those of the W and G models. In Table 1, it is seen that the pdfs of the LE distribution presents similar shapes to those of the W and G models. However, unlike the W and G models, the pdf of the LE distribution tends to σ −1 [1 − log(α)] as x → 0 + when it is unimodal; that is, it tends to a positive finite value. This led us to establish that the LE distribution can properly fit datasets whose frequency distributions are unimodal while having observations lumped around 0. On the other hand, the pdf of the LR model presents only the unimodal shape but (as will be seen later) with variations of asymmetry and kurtosis. In Table 2, it is seen that the shapes presented by the hrf of the LE distribution are similar to those of the W and G distributions. However, from the results presented in Table 2, and similarly to the behavior of the LE pdf, an important difference can be observed in the behavior of the LE hrf for lower values of x (times close to 0). On the other hand, the LR distribution is the only distribution (among the distributions considered) that has a hrf that can present the increasing-decreasing-increasing shape. Table 1. Comparison of the Lambert-exponential (LE), Lambert-Rayleigh (LR), Weibull (W) and gamma (G) models in terms of the pdf.

Model Shape Parameter (α) Interval Shape of the pdf, f X
If α = 1, the hrf of the LR model reduces to h(x) = x/σ 2 (the Rayleigh hrf) and the hrfs of the LE, W, and G models to h(x) = 1/σ (the exponential hrf).
Corollary 4. Let X 1 ∼ LE(σ, α) and X 2 ∼ LR(σ, α). Then, the asymmetry (β 1 (X i ) ) and kurtosis (β 2 (X i ) ) coefficients for X i , with i = 1, 2, are given by where m r,α (u i ), for r = 1, 2, 3, 4 and i = 1, 2, are as in Corollary 2. Figure 3 shows plots of the asymmetry and kurtosis coefficients for both the LE and LR distributions, and for the baseline distributions, exponential and Rayleigh, respectively. Note that the asymmetry and kurtosis values of the baseline distributions were extended to a range of values in the respective Lambert versions, showing greater flexibility of these latter distributions. The asymmetry and kurtosis ranges for the LE distribution were (1.456, 4.461) and (6.416, 48.814), respectively, and for the LR distribution (0.342, 1.274) and (3.027, 6.005), respectively. These ranges were calculated by minimizing and maximizing the asymmetry and kurtosis coefficients in Corollary 4. We used the integrate function of the R programming language R Core Team [21] to compute the m r,α (u i ) functions. The optimize function was used to minimize and maximize the coefficients.  (panels (a,b)) and for the Lambert-Rayleigh (solid line) and Rayleigh (circle) distributions (panels (c,d)).

Parameter Estimation
In this section, we discuss parameter estimation for the Lambert-F distributions using the moments and maximum likelihood methods.

Moment Estimators
For a random sample X 1 , . . . , X n of the random variable X ∼ LF(η 1 , η 2 , . . . , η k , α), the moment estimators are obtained by solving the equations generated by matching the first k + 1 raw moments with the first k + 1 sampling moments.
For a random sample X s,1 , . . . , X s,n of the random variable X s , for s = 1, 2, such that X 1 ∼ LE(σ, α) and X 2 ∼ LR(σ, α), the moment estimator α M for α is given by the root of the equation and the moment estimator σ M for σ is given by , for the LE distribution, , for the LR distribution, where m r,α (u s ), for r = 1, 2 and s = 1, 2, are as in Corollary 2 and X s and X 2 s are the first and second sampling moments, respectively, associated with the sample of the random variable X s . (11), it is seen that the moment estimators for α cannot be obtained in closed form. To obtain the moment estimate of this parameter, Equation (11) should be solved numerically. We solved Equation (11) using the uniroot function from the R programming language. Once obtained the estimate of α, we used this value to obtain the moment estimate for the σ parameter from Equation (12).

Maximum Likelihood Estimators
For a random sample X 1 , . . . , X n of a random variable X with LE or LR distribution (given in Section 4), the log-likelihood function is given by where f (x i ) = f (x i ; σ), F(x i ) = F(x i ; σ) and S(x i ) = 1 − F(x i ) are the baseline functions given in Table 3. Table 3. Baseline f (x), F(x), f 1 (x), f 11 (x), F 1 (x) and F 11 (x) functions (x > 0 and σ > 0) for the Lambert-exponential (LE) and Lambert-Rayleigh (LR) models.

LE LR
Since expressions of the ML estimators are not available in closed form, ML estimates are computed using the optim function in R via the L-BFGS-B method with the moments estimates, σ M and α M , as the initial values for the iterative process.
Under regularity conditions, the asymptotic distribution of the ML estimator θ = ( σ, α) of θ = (σ, α) is N 2 (θ, I(θ) −1 ), where I(θ) is the expected information matrix. Due to the analytical structure of Equation (13), it is not easy to obtain the analytical expression of this matrix. However, it can be approximated by minus the Hessian matrix evaluated at the ML estimate of the parameters. For the LE and LR models, the elements of the Hessian matrix are given by with f 11 (·) and F 11 (·) for the LE and LR models as in Table 3.
In the case of survival data, under the consideration of non-informative right-censoring and assuming independence between failure times (Y i ) and censoring times (C i ), i = 1, . . . , n, the observed time for the i-th individual is given by Thus, for a observed sample (x 1 , δ 1 ), . . . , (x n , δ n ), the log-likelihood function for the LE and LR models is given by where f (·), F(·) and S(·) are as in the Table 3. For δ i = 1, i = 1, . . . , n, Equation (14) is reduced to Equation (13). Inference based on Equation (14) can be performed in a similar manner, as was done in the uncensored case, as described above. Finally, note that the above procedures can be extended to the case where the baseline distribution has k parameters.

Simulation Study
In this section, we present a simulation study done to assess the performances of the moments and maximum likelihood estimates for the parameters of the models in Section 4. We generated N = 1000 random samples of sizes n = 50, 100, 200 from the LE and LR distributions, respectively, for different values of its parameters. The random numbers were generated by through the following steps: We used the W function of the LambertW package in R Goerg [22] to compute the principal branch of Lambert W function. Tables 4 and 5 show averages, empirical standard deviations (SD), averages of asymptotic standard errors (SE) and roots of the simulated mean square errors (RMSE) of the estimates of σ and α for the LE and LR distributions. Looking at Table 4, it can be seen that both the moments method and the ML method provide acceptable estimates of the parameters of the LE distribution. However, the ML method provides estimates with lower biases, and the SD and RMSE were smaller than those provided by the moments method. In addition, SD, SE and RMSE were closer for the ML method. The same held for SD and RMSE from the moments method. The estimators of the parameters of the LR distribution in Table 5 had similar behavior, with the exception that SD and RMSE were smaller for the moments method when the sample size was n = 50.

Application
In this section, we present two applications to real data using the Lambert-exponential (LE) and Lambert-Rayleigh (LR) distributions. With the first application, we provide evidence that the LR distribution may present a better fit than the Weibull (W), gamma (G), generalized Rayleigh (GR) Surles and Padgett [23] and Rayleigh (R) models. Similarly, with the second application we show that the LE distribution presents a better fit than the W, G, generalized exponential (GE) Gupta et al. [1] and exponential (E) models.

Nicotine Measurement Data
We considered a dataset of 346 nicotine measurements (milligrams per cigarette) collected in 1998 by the Federal Trade Commission (FTC), Washington, DC. For a recent analysis of this data, see Handique et al. [24]. Some descriptive statistics for these data are the following; average 0.852 milligrams per cigarette, standard deviation 0.334 milligrams per cigarette, sample asymmetry coefficient 0.172 and sample kurtosis coefficient 3.315.
Using results from Section 6.1, moment estimates were computed, leading to the following values: σ M = 0.519 and α M = 2.285. Using the moment estimates as initial values, ML estimates were computed, and they are presented in Table 6 with the standard errors in parenthesis. For each fitted model, the maximum value of the log-likelihood function is also reported in Table 6. Note that the LR model has a maximum value of the log-likelihood function larger than the other models. In order to compare the distributions, we computed the usual Akaike information criterion (AIC) Akaike [25] and Bayesian information criterion (BIC) Schwarz [26]. Table 6 reports AIC and BIC values for each fitted model. We can see that AIC and BIC show a better fit for the LR model. In Figure 4 (left panel), the empirical cdf for nicotine measurements is compared with the cdf of the fitted LR model, and we can see that the two curves are close. In the same figure (right panel), the histograms for the nicotine measurement data and the fitted density functions are presented. Here, we see that the LR distribution fits nicotine measurement data better than the other distributions.

Monoclonal Gammopathy Data
This dataset comprises survival times (days) from diagnosis to the last follow-up of 241 subjects diagnosed with apparently benign monoclonal gammopathy at Mayo Clinic (US). Of the 241 subjects, 16 survived until the end of the follow-up and three had monoclonal gammopathy of undetermined significance (MGUS) detected on the day of death. This dataset was previously analyzed in Kyle [27] and is currently available under the name mgus in the survival package in R Therneau [28].
Maximum likelihood estimates, with standard errors in parentheses; the maximum values of the log-likelihood functions; and AIC and BIC values for the E, GE, G, W and LE models are reported in Table 7. It can be noted that AIC and BIC show better fit of the LE model. In Table 7, for the LE model the estimate of α is α = 1.830 > 1 so that the hazard rate function monotonically increases to reach a plateau at 1/ σ = 1/4173.763 = 0.00024. In Figure 5 (left panel), the survival curves estimated by the fitted LE distribution and by the Kaplan-Meier estimator are close. In the right panel, the hazard rate function for each distribution fitted to the monoclonal gampathy data is presented. Here, we observe that the hazard rate of the LE distribution is lower than the hazard rates of the other distributions in the first 3600 days (approximately) after the diagnosis of MGUS. Opposite behavior was observed for a time greater than 3600 days.

Final Comments
In this paper, we proposed a new probability distribution generator called the Lambert-F generator. For any baseline distribution F, continuous and with positive support, the Lambert-F version is derived by applying the generator. The result is a new distributions class with one extra parameter and that generalizes the baseline distributions. The quantile function of the new class of distributions can be expressed in closed form in terms of the Lambert W function.
We proved that the hrf of a Lambert-F distribution corresponds to a perturbation of the baseline hrf, increasing or decreasing the baseline hrf for lower values of X (earlier times). We detailed two special cases corresponding to two-parameter extensions of the well-known exponential and Rayleigh distributions. We discussed moments and maximum likelihood estimators for the parameters of the proposed models. For both methods, we provided guidance on numerical procedures that might be used. Additionally, we carried out a simulation study to assess the behavior of the estimates. We found good performances for both estimators, but especially for maximum likelihood estimators, which yielded estimates with less bias. Finally, we developed two applications to real datasets, thereby providing evidence that the LE and LR distributions may present a better fit than other two-parameter distributions such as Weibull, gamma, generalized exponential and generalized Rayleigh.