A simple and efficientmethod for finding the closest generalized lambda distribution to a specific model

The four-parameter Generalized Lambda distribution (GLD) can be used to approximate many probability distributions. We present a simple and efficient two-stage process for finding optimal GLD parameters to approximate a specified distribution. The probability density quantile function is first used to find the best GLD shape parameters. Given those shape parameters, it is then straightforward to find the best location and scale parameters. We highlight the excellent performance of our approach with comparisons to two existing and popular methods for a wide choice of distributions. Finally, we show that this is method can be used with other distributions by providing applications also to the Generalized Beta distribution. Subjects: Science; Mathematics & Statistics; Statistics & Probability; Statistics; Mathematical Statistics; Statistical Computing; Statistical Theory & Methods


Introduction
introduced the Generalized Lambda distribution (GLD), a flexible fourparameter family of distributions that can be used to approximate a very large number of probability distributions. It is defined by its quantile function, which depends on a location parameter λ 1 , inverse scale parameter λ 2 , and shape parameters λ 3 and λ 4 . There are several re-parameterizations of the GLD, but the FKML parameterization (Freimer, Kollia, Mudholkar, & Lin, 1988) shown in (1) is simplest to work with, being defined for all real parameter values subject only to λ 2 > 0. This GLD is used extensively throughout the text by Karian and Dudewicz (2016) regarding fitting statistical distributions to data and the authors of a book (Pfaff, 2016)

PUBLIC INTEREST STATEMENT
The Generalized Lambda distribution (GLD) is a commonly used distribution in many fields, including finance and economics. Due to its flexibility, the GLD is widely used to approximate other distributions. Motivated by this fact, this article presents a new method which can be simply executed to obtain the optimal parameters for the GLD to approximate another distribution. Further, the results are compared with two other methods that are generally used to estimate GLD parameters to show that this new method has favourable qualities.
Given the correct parameters, the GLD distribution equates to several well-known distributions (e.g. the uniform, exponential, logistic) and for many others, a close approximation is possible (e.g. normal, Cauchy, log-normal). While the use of two or more shape parameters increases the flexibility of generalized distributions such as the GLD, an exhaustive search for optimal parameters becomes difficult due to the dimension over which the search need be carried out. We propose to use the probability density quantile (Staudte, 2017) which depends only on the shape parameters and which is defined on the interval ½0; 1 making the task of finding optimal shape parameters a simple one. It is then straightforward to find the location and inverse scale parameters. Our investigations to be detailed later reveal excellent performance of this new approach when compared to commonly used methods. Additionally, the method we propose is applicable to other distributions such as the generalized beta distribution which we explore in Section 5. Consequently, researchers exploring new generalized distributions will also find our method useful.
An outline of this paper is as follows: In Section 2 we introduce definitions used throughout. We also work through some simple examples that help to motivate our work. In Section 3 we briefly discuss two existing methods that are commonly used to identify GLD parameters before introducing our procedure. Examples are considered in Section 4, where our method is shown to perform relatively well in identifying optimal GLD parameters. Further, a web application is also provided. We then show in Section 5 that the method is suitable also for fitting generalized beta distributions to a given model. A discussion and further work are proposed in Section 6.

Definitions and some motivating examples
In this section we introduce definitions and notation to be used in the sequel.

Quantile, quantile density and density quantile functions
For u 2 ð0; 1Þ, let QðuÞ ¼ F À1 ðuÞ denote the quantile function associated with distribution function F. Let f ¼ F 0 denote the probability density function and assume f ðxÞ > 0 for all x in its domain. The quantile density function (Parzen, 1979) is defined as qðuÞ ¼ Q 0 ðuÞ ¼ 1=f QðuÞ ½ . It was earlier called the sparsity index by Tukey (1965). The quantile density function is often useful in nonparametric modeling and inference. For example, for d QðuÞ ¼ F À1 n ðuÞ denoting the estimated pth quantile arising from the empirical distribution function F n for a simple random sample of n observations, an asymptotic variance for d QðuÞ is uð1 À uÞq 2 ðuÞ=n, e.g. page 52 of Staudte and Sheather (1990). Parzen (1979) called the reciprocal of the quantile density function the density quantile function, and we will denote it here by f Q ðuÞ ¼ f ½QðuÞ:

Probability density quantile functions
The probability density quantile, introduced recently by Staudte (2017), is shown to be useful in the study of shape in a location-scale free environment; see also Staudte and Xia (2018).
The pdQ is defined for all lattice distributions and continuous distributions having square-integrable densities; it is free of location and scale parameters and defined on the finite domain ½0; 1. Further, it retains information regarding shape such as asymmetry and tail behaviour (Staudte, 2017). It, therefore, can provide a simple means to identify suitable shape parameters. Once we introduce the GLD distribution next, we will provide some simple examples using the pdQ to find shape parameters.

The generalized lambda distribution
The Generalized Lambda Distribution (GLD) of Freimer et al. (1988) is defined in terms of its quantile function: where λ 1 is a location parameter, 1=λ 2 > 0 a scale parameter and λ 3 ; λ 4 are shape parameters.
It is easy to see that the quantile density function for the GLD is qðuÞ h i so that the density quantile function is In general, no closed-form solution for the integral κ ¼ ò 1 0 f Q ðuÞdu exists, but it can be evaluated computationally and quite efficiently since integration occurs only between the finite bounds zero and one.

Some motivating examples
The following examples are used to motivate our approach that follows. While for these examples exact solutions are possible, our method provides a way to determine parameter choices for the GLD to closely approximate the true density when exact solutions are not possible.

The uniform distribution
For a simple motivating example, we consider the Uniform ða; bÞ distribution with quantile function Q 1 ðuÞ ¼ a þ uðb À aÞ and quantile density function q 1 ðuÞ ¼ Q 0 1 ðuÞ ¼ ðb À aÞ. Of course a simple inspection of the quantile functions for the uniform and GLD distributions reveals that λ 1 ¼ ða þ bÞ=2, λ 2 ¼ 2=ðb À aÞ, λ 3 ¼ 1 and λ 4 ¼ 1 results in the Uniformð0; 1Þ distribution. Below we show how the same results can be arrived at via the GLD.
The density quantile function for the uniform is f Q1 ðuÞ ¼ 1=ðb À aÞ where ò 1 0 f Q1 ðuÞdu ¼ 1=ðb À aÞ leading to the uniform pdQ f Ã Q 1 ðuÞ ¼ 1: Now by choosing λ 3 ¼ 1 and λ 4 ¼ 1 we can assure that the GLD pdQ is also free of u, as it is for the uniform distribution. For the GLD distribution f Q ðuÞ ¼ λ and f Ã Q ðuÞ ¼ 1 so that for these choices of λ 3 and λ 4 we have equality between the uniform and GLD pdQs.

Existing methods
We start by briefly describing two commonly used methods for obtaining GLD parameters. These methods are available in the bda package (Wang, 2015) available for the R statistical software (R Core Team, 2017) and will be used later to assess our new approach which is also detailed here.
3.1.1. The moments method Karian and Dudewicz (2003) and Lakhany and Mausser (2000) describe the moments matching method for GLD on the basis of the first four central moments. To estimate the GLD parameters, they match mean, variance, skewness and kurtosis of the GLD with the moments resulting in a system of non-linear equations that can be solved using optimization methods to obtain parameters.

The percentile matching method
Rather than moments, the percentile matching method equates a selected number of percentiles with the GLD theoretical percentiles. By optimizing this set of non-linear equations, parameters can be obtained for the given Generalized Lambda Distribution as discussed by Karian and Dudewicz (1999) and Tarsitano (2005). This method is preferred over the moments method as it applies less weight to the extremes and also could be used in situations where moments do not exist.

Choosing parameters based on the pdQ
Throughout let Q denote the GLD quantile function given in (1) and let f Ã Q denote the corresponding GLD pdQ where we write, for a given u, f Ã Q ðu; λ 3 ; λ 4 Þ to emphasize that the GLD is a function of the two shape parameters. Further, let Q 1 and f Ã Q 1 be the quantile function and pdQ for another distribution to be approximated.

Step 1: choosing shape parameters
In the first step, we seek the GLD shape parameters λ 3 and λ 4 that minimizes the expected squared distance between the pdQs f Ã Q and f Ã Q 1 given as arg min While this need not be an onerous task, since the integral needs to be evaluated over the finite bounds 0 to 1, it is simple to approximate this step using a discrete set of us given as for a positive integer J. Further, it may be appropriate to apply less weighting to the extreme tails to ensure that the GLD distribution is close to the distribution to be approximated in regions of high density mass. Therefore, we suggest choosing the optimal shape parameters using where w is a weight function to be chosen. We have chosen simple choices of w to be wðuÞ ¼ 1 (no weighting), wðuÞ ¼ u (weighting down the left tail), wðuÞ ¼ 1 À u (weighting down the right tail) and wðuÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uð1 À uÞ p (decreasing weights towards both the left and right tails). In what follows when considering weighting, we refer to these as left, right and both. Even for J ¼ 1000 standard optimization functionality in packages such as R can find the optimal λ 3 and λ 4 quickly.
There may be some difficulty in identifying shape parameters near zero due to the limiting behavior of ðu λ 3 À 1Þ=λ 3 and ðð1 À uÞ λ 4 À 1Þ=λ 4 which tend to logðuÞ and logð1 À uÞ respectively as the shape parameters tend to zero. An example of this is given in Section 2.4.2 for the exponential distribution. We, therefore, suggest a small tolerance value for the shape parameters such that when the parameters are close to zero such that ðu λ 3 À 1Þ=λ 3 and ðð1 À uÞ λ4 À 1Þ=λ 4 are replaced with their respective log forms. We choose a tolerance of 0.02 which is also used in the Box-Cox transformation, which has a similar form, in the R MASS package (Venables & Ripley, 2002).

Step 2: choosing location and scale parameters
Following the selection of optimal λ 3 and λ 4 in Step 1, we write the GLD quantile function as Qðu; λ 1 ; λ 2 ; λ Ã 3 ; λ Ã 4 Þ which is now a function of u and the unknown λ 1 and λ 2 . We consider two methods to identify suitable location and scale parameters.
Distributional least squares. The optimal λ 1 and λ 2 can be chosen as which is the method of distributional least squares for the linear regression of the GLD quantiles on the quantiles from the distribution to be approximated. Consequently the closed-form solution for this minimization exists in the form of usual least squares where the intercept is the optimal choice for λ 1 and the inverse of the slope parameter is the optimal choice for the scale parameter λ 2 . In this step we have found that only a few u k 's are needed to obtain suitable parameters and we set In what follows we choose K ¼ 5 although other choices may also be considered. For more on distributional least squares see, e.g., page 203 of Gilchrist (2000).

Examples
In this section, we assess the performance of our pdQ method with different weight functions and compare them with the existing moments and percentile matching methods. We use the Hellinger distance (Hellinger, 1909) to measure the distance between the true and GLD-approximated density.
As can be seen in Table 1, similar results are obtained for symmetric distributions regardless of the weight function used. For the skewed distributions, choosing the appropriate weight function improves the GLD approximation. This is due to the fact that the application of less weight to the extreme tails would enhance the GLD approximation by applying more focus in regions with higher density. An example of this could be found for cases like right skewed exponential and Pareto distributions whereby weighting down the right tail provides the close-to-exact approximation. In many cases, the quartile matching approach to identify the best parameter choices is the superior performer and will be our preferred choice moving forward. For the least squares approach, we tried increasing K but found that this typically returned poorer results.
In Table 2 we compare our pdQ method with the moments and percentile matching methods. The pdQ method typically outperforms the existing methods, returning excellent approximations for all of the distributions considered, especially for those that are highly skewed. This is apparent in Figure 1 where the approximated density curves from the pdQ method mimics the true density extremely well compared to the other methods for the lognormal, Gamma and exponential distributions. All methods do very well at approximating the normal density. The non-existence of the first four moments for the selected parameter choice for the Cauchy, Frechet and Dagum distributions makes the moment method inapplicable for those instances. However, both the percentile and the pdQ methods could be used for these distributions.
We further investigated the pdQ method and the percentile method by considering the scale-and location-free pdQ plots given in Figure 2. The pdQ curves produced by the pdQ method shows a better approximation to the true pdQ curve compared to those generated by the percentile matching method.
In Table 3 we provide the chosen parameters of the GLD distribution by the pdQ method. A large value of λ 3 for the exponential is due to the fact that the theoretically optimal value is 1 (see Section 2.4.2). The pdQ method also chooses the correct optimal choices for the uniform distribution that we considered earlier in Section 2.4.1.
We have also made a Shiny (Chang, Cheng, Allaire, Xie, & McPherson, 2017) application that can be used to find parameters of the GLD for several distributions. The app can be found at https:// lukeprendergast.shinyapps.io/GLDlambdas/ We are happy to include other distributions on request and improvements to this application will continue to be made.    Figure 1. Plots of the true and approximating GLD probability densities for several distributions using the moment matching (green), percentile matching (blue) and pdQ (red) methods. The true probability density is the black curve.

Extensions
As an another possibility, we consider the Generalized Beta Distribution (GB) with the parameterization below as expressed in the R bda package (Wang, 2015). This parameterization form includes a location parameter (a), a scale parameter (b) and two shape parameters α and β. The quantile function is Figure 2. Plots of the pdQ for several distributions using the percentile matching (blue) and pdQ (red) methods. The true pdQ is the black curve.
where Sðu; α; βÞ is the quantile function of the usual beta distribution, also commonly known as the inverse regularized beta function. A closed-form expression for Sðu; α; βÞ exists only if α ¼ 1, β ¼ 1 or α ¼ β ¼ 1. Further the density quantile function for the GB can be obtained as follows, where qðuÞ denotes the quantile density function of the beta distribution. The pdQ for the GB can be obtained similarly using numerical integration as a closed-form solution is not available Table 4 summarizes the performance of the pdQ approach in choosing optimal parameters for the GB distributions in approximating several distributions. Similar results can be found by Karian and Dudewicz (2016) for parameter choices of the GB using the moments method. Excellent approximations can be obtained for most of the distributions such as the Uniform, Normal, exponential, χ 2 , gamma, beta and triangular distributions using the pdQ approach. The parameterization in 5 indicates that the location and scale are the lower and upper bounds of support for the GB. Therefore, approximations for those parameters get larger for distributions which have infinite support in an effort to capture enough of the density being approximated. Further, shape parameters also can be very large due to the theoretically optimal value being infinite. As an example, both the pdQ method here and the moments matching method in Karian and Dudewicz (2016) provide similar approximations for the normal distribution. The moments matching method approximates the standard normal with parameter choices −1414. 214, 1414.214, 999998.848, 999998.848 while the pdQ approach provides similar choices in Table 4. Moreover, Karian and Dudewicz (2016) report that it is possible to further improve this approximation by letting the two shape parameters get even larger. For the standard beta distribution, the theoretically optimal choices for the four GB distribution are the first and second shape parameter values and a location and scale of zero and one respectively. As can be seen in the table, the pdQ approach correctly identifies these parameters.
We now further look at two specific cases where an exact solution can be obtained for the GB using the pdQ method.

The uniform distribution
We revisit Section 2.4.1 for the uniform distribution. As described below, we can use the same process to identify the theoretically optimal parameters for the GB. The pdQ for the uniform can be derived as f Ã Q1 ðuÞ ¼ 1. When α ¼ 1 and β ¼ 1, the closed-form expression for the quantile function of the beta distribution is found to be Sðu; α ¼ 1; β ¼ 1Þ ¼ u (Sharma & Chakrabarty, 2017). Hence, for these parameter choices, we have the GB pdQ as f Ã Q ðuÞ ¼ 1 and equality between the uniform and GB pdQs. As with the beta distribution, the pdQ approach has identified the optimal parameters.

Discussion
By concentrating on first identifying suitable shape parameters using the probability density quantile function, we have introduced a method for choosing optimal GLD parameters to approximate other distributions. Our method typically outperforms the approaches based on moment and percentile matching, and it is simple and efficient to implement. We also showed that the approach is also suitable for other distributions such as the generalized beta distribution. We are currently investigating estimation methods using some of the ideas developed here in the population setting.