Comparing Estimation Methods for the FPLD

Abstract

We use the quantile function to define statistical models. In particular, we present a five-parameter version of the generalized lambda distribution (FPLD). Three alternative methods for estimating its parameters are proposed and their properties are investigated and compared by making use of real and simulated datasets. It will be shown that the proposed model realistically approximates a number of families of probability distributions, has feasible methods for its parameter estimation, and offers an easier way to generate random numbers.

1. Introduction

Statistical distributions can be used to summarize, in a small number of parameters, the patterns observed in empirical work and to uncover existing features of data, which are not immediately apparent. They may be used to characterize variability or uncertainty in a dataset in a compact way. Furthermore, statistical distributions may both facilitate the mathematical analysis of the basic structure of empirical observations and interrelate information from two or more sources. The aim of this paper is to make a contribution to the use of quantile statistical methods in fitting a probability distribution to data, following the line of thought indicated by Parzen [1] and Gilchrist [2, Section 1.10]. In particular, we adopt a five-parameter version of the generalized lambda distribution𝑋(𝑝,𝝀)=𝜆₁+𝜆₂2{(1−𝜆₃𝑝)[^𝜆4−1𝜆₄]−(1+𝜆₃𝑞)[^𝜆5−1𝜆₅]}, (1.1) where 0≤𝑝≤1,𝑞=1−𝑝. 𝑋(𝑝,𝝀) is the quantile function and 𝑝∈[0,1]. If 𝜆₂≥0 and 𝜆₃∈[−1,1], then (1.1) is a continuous and increasing function of 𝑝. Here, 𝜆₁ controls, albeit not exclusively, the location of an FPLD; the parameter 𝜆₂ acts for a multiplier to the translated quantile function [𝑋(𝑝,𝝀)−𝜆₁] and is, thus, a scale parameter; 𝜆₃,𝜆₄,𝜆₅ influence the shape of [𝑋(𝑝,𝝀)]. The limiting cases 𝜆₄→0 and/or 𝜆₅→0 must be interpreted according to L'Hopital's rule. The references [2–7] are illustrative regarding the history of various parametrizations of the distribution (see also [8]). The (1.1) version was proposed by Gilchrist [2, page 163] as an extension of the generalized lambda distribution (GLD) in the parametrization given by Freimer et al. [5]𝑋(𝑝,𝝀)=𝜆₁+1𝜆₂𝑝{[^𝜆4−1𝜆₄𝑞]−[^𝜆5−1𝜆₅]}, (1.2) where the scale parameter 𝜆₂ appears as a reciprocal and 𝜆₃=0.

The FPLD can be useful in problem solving when there are difficulties with specifying the precise form of an error distribution, because it can assume a wide variety of curve shapes and, more importantly, it uses only one general formula over the entire range of data. Its versatility, however, is obtained at the cost of an unusually large number of parameters if compared with the array of studies referred to in the literature on fitting distributions to data. Gilchrist [9] notes the rather strange fact that, throughout the history of statistics, researchers seem to have been quite happy to use many parameters in, for example, regression, yet they have kept to two or three parameters in most distributional models. Clearly, a high number of parameters are a less challenging task in our age of readily available computer power. Moreover, the number of individual data in a statistical sample has considerably increased from what was customary a few decades ago, so there should be no special reason for not using five- (or more) parameter distributional models.

Regrettably, the richness of shape in (1.1) is not accompanied by a comparable development in parameter estimation techniques. For instance, one of the classical methods for computing 𝝀=(𝜆₁,…,𝜆₅) can be based on matching the first five moments of the sample observations. In this case, though, the problem of modelling the data would be restricted to distributions possessing fairly light tails because the fifth moment must be finite. In addition, the percentiles estimates (e.g., Karian and Dudewicz [10]) and methods using location and scale-free shape functionals (e.g., King and MacGillivray [11]) depend markedly on the particular choice of percentage points. On the other hand, maximum likelihood estimation is computationally demanding because of the complex and unwieldy objective function (see Tarsitano [12]). Asquith [13] as well as Karvanen and Nuutinen [14] equated sample L-moments to those of the fitted distribution. However, higher L-moments are inelastic to changes in the parameters, leading to difficulties in estimation. These questions have been the focus of intense academic research in the past few years, for example, Ramberg et al. [4], King and MacGilllivray [15], Karian and Dudewicz [7], Lakhany and Mausser [16], Rayner and MacGillivray [17], Fournier et al. [18], and Su [19].

In the present paper, we start from the fact that the vector of parameters 𝝀 can be divided

into linear and nonlinear parameters, so that we can apply the nonlinear least squares (NLS) method, proposed by Öztürk and Dale [20], which is driven by the minimization of the sum of squared differences between the observed order statistics and their expected values. Furthermore, we will develop a criterion driven by the minimization of the sum of the least absolute deviations between the observed order statistics and their medians. The two estimation procedures are based on a two-stage inner/outer optimization process. The first stage is the inner step in which the linear parameters are estimated for fixed values of the nonlinear parameters. The second stage is the outer step in which a controlled random search (CRS) technique is applied in order to determine the best companion value for the nonlinear parameters. The two stages are carried out consecutively, and the process is stopped when the difference between the values of the minima of two successive steps becomes smaller than a prefixed tolerance threshold.

The rest of the paper is organized as follows. The Section 2 describes the density of the FPLD random variables and underlines the advantage of its use in practical work. The Section 3 derives the nonlinear least squares for the estimation of 𝝀. In Section 4, we outline a procedure for applying the criterion of least absolute deviations to the FPLD. The CRS, an optimization algorithm based on a purely heuristic, derivative-free technique is presented in the fifth section. In Section 6, the effectiveness of the proposed methods is demonstrated and compared using real and simulated data. Finally, we conclude and outline future research directions in Section 7.

2. Overview of the Five-Parameter Lambda Distribution (FPLD)

The FPLD is a highly flexible and adaptable tool for modelling empirical and theoretical distributions and, hence, might have many potential applications. In this section, we summarize its general features. The density-quantile function of a FPLD random variable is𝜆𝑓[𝑋(𝑝,𝝀)]=²⁻¹𝛾₁𝑝^𝜆4−1+𝛾₂𝑞^𝜆5−1,𝛾₁=(1−𝜆₃)2,𝛾₂=(1+𝜆₃)2. (2.1)

The range of a FPLD is given by [𝑋(0,𝝀),𝑋(1,𝝀)] and is different for different values of 𝝀; in particular, the extremes are [𝜆₁−𝜆₂𝛾₁,𝜆₁+𝜆₂𝛾₂] if 𝜆₄,𝜆₅>0; [𝜆₁−𝜆₂𝛾₁,∞) if 𝜆₄>0, 𝜆₅≤0; and (−∞,𝜆₁+𝜆₂𝛾₂] if 𝜆₄≤0, 𝜆₅>0. Hence, the FPLD can model distributions whose support extends to infinity in either or both directions as well as those with finite support.

The rth moment of the linear transformation 𝑍=2(𝑋−𝜆₁)/𝜆₂+𝑘₁−𝑘₂ where 𝑘₁=(1−𝜆₃)/𝜆₄, 𝑘₂=(1+𝜆₃)/𝜆₅ is𝐸(𝑍^𝑟)=^𝑟_𝑗=0(−1)^𝑗𝑟𝑗𝑘₁^{(𝑟−𝑗)}𝑘^𝑗₂𝐵[1+(𝑟−𝑗)𝜆₄,1+𝑗𝜆₅]. (2.2) The symbol 𝐵(𝛼,𝛽) with 𝛼,𝛽>0 denotes the complete beta function. As a consequence, the 𝑟th moment of 𝑍 exists and is finite if 𝜆₄ and 𝜆₅ are greater than 𝑟⁻¹.

The central moments of 𝑋 which can be derived from (2.2) are𝜎²=𝛼₁−2𝑘₁𝑘₂𝐵(𝜆₄+1,𝜆₅𝛾+1),₁=𝛼₂𝜎³−3(𝑘₁𝑘₂)²𝜎³[𝐵(2𝜆₄+1,𝜆₅+1)𝑘₂−𝐵(𝜆₄+1,2𝜆₅+1)𝑘₁𝛾],₂=𝛼₃𝜎⁴−4(𝑘₁𝑘₂)³𝜎⁴[𝐵(3𝜆₄+1,𝜆₅+1)𝑘²₂+𝐵(𝜆₄+1,3𝜆₅+1)𝑘²₁−+,3𝐵(2𝜆₄+1,2𝜆₅+1)2𝑘₁𝑘₂𝛼],_𝑖=𝑘₁^𝑖+1[(𝑖+1)𝜆₄++1](−1)^𝑖+1𝑘₂^𝑖+1[(𝑖+1)𝜆₅+1],𝑖=1,2,3. (2.3)

The shape of the FPLD is sensitive to its parameter values; in fact, the moments of 𝑋(𝑝,𝝀) are determined by a combination of all of the vector 𝝀 elements. If 𝜆₄=𝜆₅ and 𝜆₃=0, then the FPLD is symmetric about 𝜆₁ because, in this case, (1.1) satisfies the condition 𝑋(𝑝,𝝀)=−𝑋(𝑞,𝝀). If 𝜆₄=𝜆₅, then (2.1) is skewed to the left (right) if 𝜆₃<0 (𝜆₃>0), which suggests a natural interpretation of 𝜆₃ as a parameter which is prevalently related to the asymmetry of 𝑋(𝑝,𝝀). For a FPLD, interchanging 𝜆₄ and 𝜆₅ and simultaneously changing the sign of 𝜆₃ reverses the type of skewness. If 𝜆₃=−1 (𝜆₃=1), then 𝜆₄ (𝜆₅) indicates the left (right) tail flatness in the sense that the smaller 𝜆₄ (𝜆₅) the flatter the left (right) tail of the FPLD density. In practice, 𝜆₄ and 𝜆₅ capture the tail order on the two sides of the support of 𝑋(𝑝,𝝀).

The densities (2.1) are zeromodal if {max(𝜆₄,𝜆₅)>1}∧{min(𝜆₄,𝜆₅)<1}, unimodal with continuous tails if {max(𝜆₄,𝜆₅)<1}, unimodal with truncated tails if {min(𝜆₄,𝜆₅)>2}, U-shaped if 1<𝜆₄,𝜆₅<2, and S-shaped if {max(𝜆₄,𝜆₅)>2}∧{min(𝜆₄,𝜆₅)>1}. Curves corresponding to large positive values of 𝜆₄ and 𝜆₅ have extreme peakedness and short high tails. For 𝜆₂=0 the FPLD degenerates 𝑋(𝑝,𝝀)=𝜆₁ to a one-point distribution.

The FPLD family is a good model for use in Monte Carlo simulation and robustness studies because it contains densities which cover a large spectrum of arbitrarily shaped densities. For example, if 𝜆₄→0,𝜆₅→0, then 𝑋(𝑝,𝝀) converges on a potentially asymmetric form of the logistic distribution; if 𝜆₄→∞ and 𝜆₅→0, then 𝑋(𝑝,𝝀) is an exponential distribution whereas for 𝜆₄→0,𝜆₅→∞, 𝑋(𝑝,𝝀) becomes a reflected exponential distribution. Furthermore, FPLD fits data containing extreme values well; in fact, for 𝜆₄→∞,|𝜆₅|<∞, the FPLD corresponds to the generalized Pareto distribution. For 𝜆₅→∞,|𝜆₄|<∞, the FPLD generates the power-function distribution. The rectangular random variable is present in four versions: (𝜆₃=−1,𝜆₄=1), (𝜆₃=1,𝜆₅=1), (𝜆₄=1,𝜆₅=1), and (𝜆₃=0,𝜆₄=2,𝜆₅=2).

The FPLD could have an important role as a general (not necessarily Gaussian) distribution if it is able to produce good approximations to many of the densities commonly encountered in applications. There are several measures of closeness and overlapping of two statistical distributions based on the distance between density functions or between distributions functions (see, e.g., [21] and [7, pages 196–201]). On the other hand, the FPLD is expressed in the quantile domain, so, for a reliable assessment of the agreement between the exact distribution and that provided by the FPLD, we should use a measure based on quantiles. Statistical literature is not very helpful in providing criteria to quantify the proximity between distributions naturally arising from a quantile framework. In the present paper we have used the Tchebycheff metric, that is, the closeness between the theoretical quantile function 𝑋(𝑝) and the FPLD is quantified by the maximum absolute difference (Maxd) between observed and fitted quantiles:max_{1≤𝑖≤500}|𝑋(𝑝_𝑖)−𝑋(𝑝_𝑖,𝝀)|,where𝑝_𝑖=𝑖501,𝑖=1,2,…,500. (2.4) The Tchebycheff metric is attractive because it leads to the minimization of the worst case and does not imply a heavy computational burden.

The five parameters of 𝑋(𝑝,𝝀) have been estimated using the controlled random search method described in Section 5. Table 1 shows the accuracy of the FPLD fit to some standard symmetric distributions.

Table 1 

Comparison of FPLD approximations for symmetric distributions.

Our findings suggest that the FPLD fit is reasonably good for all of the models included in the table, with the only possible exception being the Cauchy random variable.

One of the major motives for generalizing the four-parameter GLD to the five-parameter FPLD is to have a more flexible, and, hence, more suitable, distribution for fitting purposes. How much better, though, is the FPLD when compared to the GLD? Are the fitting improvements that the FPLD provides significant enough to offset the complexity that is introduced through a fifth parameter? Table 2 shows the results of the fitting procedure for some asymmetric distributions. The first row presents the parameter estimates and the Tchebycheff metric for the FPLD approximation; the second row presents the same quantities for the GLD approximation using the FMKL parametrization (see Freimer et al. [5]). The values in the table indicate that the upper bound on the discrepancy between the true and the approximated distribution is systematically lower for the FPLD than it is for the GLD. The amounts appear to be large enough to justify the inclusion of an additional linear parameter in the GLD so as to form the FPLD.

Table 2 

Comparison of FPLD and GLD approximations for asymmetric distributions.

The most important point to note in Tables 1 and 2 is that the FPLD is a flexible class of very rich distributional models which cover the Gaussian and other common distributions. In this sense, the FPLD is a valid candidate for being fitted to data when the experimenter does not want to be committed to the use of a particular distribution.

3. Nonlinear Least Squares Estimation

Öztürk and Dale [20] proposed a nonlinear least squares (NLS) framework to estimate vector 𝝀 for a four-parameter generalized distribution. The method can easily be generalized to handle a FPLD.

For a given random sample 𝑆={𝑥₁,𝑥₂,…,𝑥_𝑛} of 𝑛-independent observations from 𝑋(𝑝,𝝀), the 𝑖th order statistic may be expressed as𝑥_𝑖∶𝑛=𝐸(𝑥_𝑖∶𝑛)+𝑒_𝑖,𝑖=1,2,…,𝑛, (3.1) in which the ordered observations are compared with their expected values under the hypothesized distribution 𝑋(𝑝,𝝀) and 𝑒_𝑖 is a measure of the discrepancy between the observed and modelled 𝑖th value. The error terms {𝑒_𝑖,𝑖=1,2,…,𝑛} are such that 𝐸(𝑒_𝑖)=0 and 𝜎²(𝑒_𝑖)=𝜎²_𝑖. The heteroskedasticity in (3.1) is not the only violation of standard assumptions of the classical regression model. In fact, the error terms are not independent and do not come from a symmetrical distribution (except for the median). However, since the purpose was to obtain an approximate solution, Öztürk and Dale [20] ignored these inadequacies.

The expected value of the 𝑖th order statistic from a FPLD is available in closed form𝐸(𝑥_𝑖∶𝑛)=𝜆₁+(𝛽₁/𝜆₄)∫¹₀𝑝^{𝜆4+𝑖−1}𝑞^𝑛−𝑖𝑑𝑝+(𝛽₂/𝜆₅)[1−∫¹₀𝑝^𝑖−1𝑞^{𝜆5+𝑛−𝑖}𝑑𝑝],𝐵(𝑖,𝑛+1−𝑖) (3.2) where 𝛽₁=0.5(1−𝜆₃)𝜆₂,𝛽₂=0.5(1+𝜆₃)𝜆₂. It follows that the deterministic component of (3.1) can be written as𝑥_𝑖∶𝑛(𝝀)=𝛽₀𝑈_0,𝑖+𝛽₁𝑈_1,𝑖+𝛽₂𝑈_2,𝑖,𝑖=1,2,…,𝑛, (3.3) with 𝛽₀=𝜆₁ and𝑈_0,𝑖=1,𝑈_1,𝑖=𝜆₄⁻¹[Γ(𝑛+1)Γ(𝑖+𝜆₄)Γ(𝑖)Γ(𝑛+1+𝜆₄)𝑈−1],_2,𝑖=𝜆₅⁻¹[1−Γ(𝑛+1)Γ(𝑛+1−𝑖+𝜆₅)Γ(𝑛+1−𝑖)Γ(𝑛+1+𝜆₅)]. (3.4)

The least squares approach calls for 𝜷 in (3.3) to be chosen so as to minimize𝑆₂(𝝀)=^𝑛_𝑖=1[𝑥_𝑖∶𝑛−𝑥_𝑖∶𝑛(𝝀)]². (3.5) The 𝑆₂(𝝀) is linear in 𝛽₁,𝛽₂,𝛽₃ but not in (𝜆₄,𝜆₅). According to Lawton and Sylvestre [22] the solution to (3.5) can be obtained by fixing the value of (𝜆₄,𝜆₅) and, then, applying the ordinary least squares to solve the linear problem𝜷(𝜆₄,𝜆₅)=(𝐔^𝑡_𝝀𝐔_𝝀)⁻¹𝐔^𝑡_𝝀𝐱 (3.6) for 𝜷, where 𝐔_𝝀 denotes the (𝑛×3) matrix with elements given by (3.4) and 𝐱 is the (𝑛×1) vector of the ordered observations. The conversion of 𝛽₀,𝛽₁,𝛽₂ into 𝜆₁,𝜆₂,𝜆₃ is straightforward:𝜆₁=𝛽₀,𝜆₂=𝛽₁0.5(1−𝜆₃)=𝛽₂0.5(1+𝜆₃),𝜆₃=1−(𝛽₁/𝛽₂)1+(𝛽₁/𝛽₂). (3.7)Öztürk and Dale [20] used an approximation of 𝐸(𝑥_𝑖∶𝑛) to overcome the problems created by repeated use of the gamma function in (3.4). More specifically, they used 𝐸(𝑥_𝑖∶𝑛)=𝑋(𝑝⁺_𝑖,𝝀), with the plotting positions 𝑝⁺_𝑖=𝑖(𝑛+1)⁻¹,𝑖=1,2,…,𝑛. In this case, the pseudoregressors are𝑉_0,𝑖=1,𝑉_1,𝑖=𝜆₄⁻¹[(𝑝⁺_𝑖)^𝜆4−1],𝑉_2,𝑖=𝜆₅⁻¹[(𝑝⁺_{𝑛+1−𝑖})^𝜆5−1]. (3.8) Parameter estimates are then obtained as in (3.5) through (3.7) with 𝑉_𝑗,𝑖 in place of 𝑈_𝑗,𝑖 for 𝑗=1,2,3. We can now define the reduced form of the quantile function (3.5) in which 𝜷 is substituted into (3.3). Clearly, the reduced quantile function only depends on the pair (𝜆₄,𝜆₅). We might try to find a better estimate of (𝜆₄,𝜆₅) by using one of the function minimization techniques to solve (3.5) and then go through exactly the same procedure as the one described above. The two-stage process is repeated until the correction for 𝑆₂(𝝀) becomes sufficiently small.

4. Least Absolute Deviations Estimation

Filliben [23] noted that order statistics means have three undesirable properties: (1) no uniform technique exists for generating the 𝐸(𝑥_𝑖∶𝑛) for all distributions; (2) for almost all distributions 𝐸(𝑥_𝑖∶𝑛) is difficult or time consuming to compute and so must be stored or approximated; (3) for other distributions (e.g., Cauchy), the expected value of the order statistics 𝐸(𝑥_𝑖∶𝑛) may not always be defined. All three of these drawbacks are avoided in general by choosing to measure the location of the 𝑖th order statistic by its median rather than by its mean.

The use of the median in preference to the mean leads naturally to the absolute error being the loss function instead of the squared error. Gilchrist [2] called distributional residuals the deviations of the observed data from a given quantile function 𝑒_𝑖=𝑥_𝑖∶𝑛−𝑋(𝑝_𝑖,𝝀),𝑖=1,2,…,𝑛. In this context, the method of nonlinear least squares discussed in the previous section is defined as distributional least squares since it involves using the ordered data as a dependent variable, whereas the deterministic component of the model is defined by a quantile function.

The least squares criterion is reasonable when data come from a normal distribution or when linear estimates of the parameters are required. While data may contain outliers or highly influential values, the method of least squares has a clear disadvantage in that it may be pulled by extremely large errors. The least absolute deviation estimation of parameters is more suitable in such cases.

In this section we consider the problem of finding a vector of constants 𝐲∈𝐑^𝑛 which minimizes the objective functionmin_{𝑦∈𝑅𝑛}^𝑛_𝑖=1|𝑥_𝑖∶𝑛−𝑦_𝑖|, (4.1) where 𝑥_𝑖∶𝑛 is the 𝑖th order statistic of the sample 𝑆={𝑥₁,𝑥₂,…,𝑥_𝑛} from the quantile function 𝑋(𝑝,𝝀). Rao [24] observed that the solution to this problem is the vector of marginal medians, that is, the 𝑖th element of 𝐲 is the median of the order statistics 𝑥_𝑖∶𝑛.

The quantile function (see, e.g., Gilchrist [2, page 86]) of the 𝑖th order statistic for a sample of size 𝑛 from an FPLD is 𝑋[𝑝^∗_𝑖(𝑝),𝝀] with𝑝^∗_𝑖(𝑝)=𝐵⁻¹(𝑝,𝑖,𝑛+1−𝑖), (4.2) where 𝐵⁻¹(⋅) is the inverse of the incomplete beta function. Consequently, the median of the 𝑖th order statistic can be computed using𝑝^∗_𝑖=𝑋[𝐵⁻¹(0.5,𝑖,𝑛+1−𝑖),𝝀],𝑖=1,2,…,𝑛. (4.3) The quantities in (4.3) are a solution of (4.1) in the sense that, for a fixed value of 𝝀, the following expression^𝑛_𝑖=1|𝑥_𝑖∶𝑛−𝑋[𝑝^∗_𝑖,𝝀]| (4.4) attains its minimum value. It follows that the estimate of 𝝀 that minimizes (4.4) can be obtained by the use of distributional absolutes (DLAs).

There is nothing in the previous section that cannot be extended to the criterion (4.4). In the light of the decomposition of 𝝀 into two groups, linear parameters (𝜆₁,𝜆₂,𝜆₃), and nonlinear parameters (𝜆₄,𝜆₅), another two-stage procedure can be devised. In the inner stage, we determine the vector 𝜷 which minimizes𝑆₁(𝝀)=^𝑛_𝑖=1|𝑥_𝑖∶𝑛−𝑥_𝑖∶𝑛(𝝀)|,with𝑥_𝑖∶𝑛(𝝀)=²_𝑗=0𝛽_𝑗𝑊_𝑖,𝑗, (4.5) with respect of the linear parameters 𝜷 and for a given value of the nonlinear parameters (𝜆₄,𝜆₅). In this case, the pseudo regressors are𝑊_0,𝑖=1,𝑊_1,𝑖=𝜆₄⁻¹{[𝑝^∗_𝑖]^𝜆4−1},𝑊_2,𝑖=𝜆₅⁻¹{[1−𝑝^∗_𝑖]^𝜆5−1}. (4.6)

The computation of the LAD estimates for 𝜷 in (4.5) can be formulated as a linear programming optimization and the standard simplex method can be employed to solve this problem. For a detailed treatment of the method of least absolutes see, for example, Koenker and D'Orey [25]. The conversion of 𝛽₀,𝛽₁,𝛽₂ into 𝜆₁,𝜆₂,𝜆₃ takes place as shown in (3.7). Once (𝜆₁,𝜆₂,𝜆₃) have been estimated, we can obtain a new value for (𝜆₄,𝜆₅) with the same outer step as that used in the previous section. The inner and outer steps can be alternated until a satisfactory result is found.

Least absolute deviations have, undoubtedly, greater sophistication than least squares, although two important points should be made regarding this issue. First, the annoying computation of 𝑝^∗_𝑖 has to be executed just once. In addition, the fact that 𝐵(𝑝,𝑖,𝑛+1−𝑖)=1−𝐵(1−𝑝,𝑛+1−𝑖,𝑖) implies that 𝑝^∗_{𝑛+1−𝑖}=1−𝑝^∗_𝑖 and, consequently, the number of medians to be computed is halved.

5. Controlled Random Search (CRS)

The inner/outer scheme outlined in previous sections has the merit of reducing the dimension of the nonlinear problem. To apply it, however, we must minimize two relatively complex objective functions 𝑆₁(𝝀) and 𝑆₂(𝝀) of two unknowns, 𝜆₄,𝜆₅, which generally have more than one minimum. The solutions to these equations are difficult to obtain through traditional optimization approaches. In order to perform this special task, we resorted to a direct search optimization technique (CRS) proposed by Price [26], which is suitable for searching the global minima for continuous functions that are not differentiable, or whose derivatives are difficult to compute or to approximate (see [27, 28]).

Some well-known direct search procedures such as the downhill simplex method of Nelder and Mead [29] and the pattern-search method of Hooke and Jeeves [30] are familiar to statisticians. However, these techniques are really local optimization methods; in practice, they are designed so as to converge on a single optimum point and, therefore, they, unavoidably, discard information related to all other possible optima. In contrast, a wise application of the CRS overcomes this drawback since this method is especially suited to multimodal functions. In effect, CRS procedures have more randomness in their decision process, thus increasing the chances of finding a global optimum. Of course, global optimality cannot be guaranteed in finite time.

A CRS procedure employs a storage 𝐴 of points, the number of which is determined by the user for the particular problem to be solved. Repeated evaluation of the function to be optimized, 𝑓(𝐱), is performed at points randomly chosen from the storage of points and the search region is progressively contracted by substituting the worst point with a better one. The search continues until an iteration limit is reached, or until a desired tolerance between minimum and maximum values in the stored 𝑓(𝐱) values is attained. Apparently, these methods are less known and less frequently used in statistics (two works that have appeared in literature are those by Shelton et al. [31] as well as Křivý and Tvrdík [32]).

For practical purposes, it is necessary to confine the search within a prescribed bounded domain. Let 𝐶⊂𝐑^𝑘 be the set of admissible values of the 𝑘=2 nonlinear parameters of the FPLD model. The most promising region that we have found in our applications is the rectangle delimited by𝐶={(𝜆₄,𝜆₅)∣−0.999≤𝜆₄≤3;−0.999≤𝜆₅≤3}. (5.1) Furthermore, 𝐶 ensures that the mean of the order statistics from a FPLD always exists and the ranges of 𝜆₄,𝜆₅ cover most of the distributional shapes commonly encountered.

The theory and practice of global optimization have progressed rapidly over the last few years, and a wide variety of modifications of the basic CRS algorithm are now available. In particular, we refer to the scheme discussed in Brachetti et al. [33].

Any CRS implies many subjective choices of constants and weights. The version used in our paper is as follows. (1)Generate a set 𝐴 of 𝑚₀>(𝑘+1) random admissible points in 𝐂. Find the minimum and the maximum in 𝐴 and their function values: (𝐱_𝐿,𝑓_𝐿),(𝐱_𝐻,𝑓_𝐻). (2)If |𝑓_𝐻−𝑓_𝐿|/(10⁻⁷+|𝑓_𝐻|+|𝑓_𝐿|)<10⁻⁸, then stop. Also stop if the limit to the number of function evaluations has been reached. The point with the lowest objective function value found is returned. If required, return all the points in 𝐴. (3)Choose at random, without duplication, a set 𝐵 of 𝑘 points in (𝐴−𝐱_𝐿), include 𝐱_𝐿 in 𝐵, and compute the centroid 𝐱=(𝑘+1)⁻¹∑_𝐱∈𝐵𝐱. Define 𝐴₁=𝐴−𝐵. (4)Select randomly a point 𝐱_𝑟∈𝐴₁. (5)Generate uniform random numbers 𝑑₁,𝑑₂,…,𝑑_𝑘 in [0,2] and compute the trial point 𝑥_𝑖,𝑇=(1+𝑑_𝑖)𝐱_𝑖−𝑑_𝑖𝐱_𝑟,𝑖,for𝑖=1,2,…,𝑘; thus, the exploration of the space around the centroid is based on a randomized reflection in the simplex originally used by Price [26]. (6)If 𝐱_𝑇∈𝐌, then compute 𝑓_𝑇=𝑓(𝐱_𝑇) and go to step (7). If 𝐱_𝑇∉𝐌, then repeat step (5) a fixed number 𝑚₁ of times. If 𝐱_𝑇 is still unfeasible, then exclude 𝐱_𝑟 from 𝐴₁, redefine the set 𝐴₁ and, if it is not empty, repeat step (4) for at most 𝑚₂ times; otherwise go to step (9). (7)If 𝑓_𝐿<𝑓_𝑇<𝑓_𝐻, then replace 𝑓_𝐻 with 𝑓_𝑇 and 𝐱_𝐻 with 𝐱_𝑇 and then redetermine the maximum 𝐱_𝐻 and 𝑓_𝐻. (8)If 𝑓_𝑇<𝑓_𝐿, then replace 𝑓_𝐻 with 𝑓_𝐿, 𝐱_𝐻 with 𝐱_𝐿, 𝑓_𝐿 with 𝑓_𝑇, 𝐱_𝐿 with 𝐱_𝑇, and 𝐱_𝐻 and 𝑓_𝐻. (9)Randomly choose two distinct points 𝐱₂,𝐱₃ in (𝐴−𝐱_𝐿) and compute 𝑑_𝑄,𝑖=(𝑥_2,𝑖−𝑥_3,𝑖)𝑓_𝐿+(𝑥_3,𝑖−𝑥_𝐿,𝑖)𝑓₂+(𝑥_𝐿,𝑖−𝑥_2,𝑖)𝑓₃, (5.2) for 𝑖=1,2,…,𝑘.(10)If |𝑑_𝑄,𝑖|>10⁻⁷ for each 𝑖, then build the following quadratic interpolation of the objective function: 𝑥_𝑄,𝑖=[(𝑥_2,𝑖−𝑥_3,𝑖)²𝑓_𝐿+(𝑥_3,𝑖−𝑥_𝐿,𝑖)²𝑓₂+(𝑥_𝐿,𝑖−𝑥_2,𝑖)²𝑓₃2𝑑_𝑄,𝑖] (5.3) for 𝑖=1,2,…,𝑘. If some of the |𝑑_𝑄,𝑖|≤10⁻⁷, then generate 𝑘 uniform random numbers 𝑑_𝑖,𝑖=1,2,…,𝑘, in [0,1] and compute 𝑥_𝑖,𝑄=𝑑_𝑖𝑥_𝑖,𝐿+(1−𝑑_𝑖)𝑥_𝑖,𝐻,for𝑖=1,2,…,𝑘. (5.4)(11)If 𝐱_𝑄∉𝐌, then repeat steps (9)-(10) no more than 𝑚₃ times. If 𝐱_𝑄 is still unfeasible, then return to step (2).(12)Execute steps (7)-(9) with 𝐱_𝑄 in place of 𝐱_𝑇 and 𝑓_𝑄 in place of 𝑓_𝑇 and then return to step (2).

We experimented by using the CRS algorithm, with 𝐱=𝝀, 𝑓(𝐱)=𝑆₂(𝝀), 𝑚₀=30𝑘, 𝑚₁=𝑚₂=𝑚₃=4(𝑘+1), and 𝑘=2. For the random component of our implementation we used quasirandom Faure sequences (Faure [34]) rather than pseudorandom numbers. This choice was made in order to cover the search region 𝐂 more evenly (see, e.g., Morokoff and Caflisch [35]).

The CRS procedure has several desirable features. Boundaries and objective function can be arbitrarily complex because the procedure does not require the calculation of derivatives. Shelton et al. [31] observed that, if several points within the optimization region are equally good, the procedure will make this clear to the user.

The general outline presented here should give some appreciation of the philosophy underlying the design of CRS. It is important not to forget that, although this method proved very robust in many applications, typical nonlinear problems which have a complicated bound-constrained, possibly multiextremal objective function are challenging tasks and the success of an automated algorithm is still dependent upon the researcher's guidance. In particular, intervention may be needed to contract the box constraints on the parameters in order to perform a finer search in a final stage.

6. Comparison of the Procedures

In this section we present results of some numerical experiments carried out in order to study the behavior of the FPLD model. Moreover, we perform a comparison between estimated and theoretical values to test the accuracy of the estimators proposed in the previous sections and investigate their properties through a Monte Carlo sampling process.

6.1. Fitting to Data

To illustrate the efficacy of the FPLD as a modeling device we applied it to four sets of real data. To judge the overall agreement between estimated and observed quantiles we used the correlation between the empirical values 𝑥_𝑖∶𝑛 centered on the arithmetic mean 𝑥 and the hypothesized values∑𝑟(𝝀)=^𝑛𝑖=1[𝑥_𝑖∶𝑛−𝑥][𝑋(𝑝_𝑖,𝝀)−𝜇(𝝀)]{∑^𝑛_𝑖=1[𝑥_𝑖∶𝑛−𝑥]²∑ⁿ_𝑖=1[𝑋(𝑝_𝑖,𝝀)−𝜇(𝝀)]²}^0.5, (6.1) where 𝑝_𝑖=𝑖/(𝑛+1), and𝜇(𝝀)=¹₀𝑋(𝑝,𝝀)𝑑𝑝=𝜆₁+𝜆₂2[1+𝜆₃1+𝜆₅−1−𝜆₃1+𝜆₄]. (6.2) From another point of view, (6.1) can be considered a measurement of the linearity of the Q-Q plot for empirical versus theoretical standardized quantiles. The coefficient (6.1) will always be positive since 𝑥_𝑖∶𝑛 and 𝑋(𝑝_𝑖,𝝀) are both arranged in ascending order of magnitude. More specifically, 𝑟(𝝀) takes values in the [0,1] interval. The case 𝑟(𝝀)=0 corresponds to a noninformative model 𝑋(𝑝_𝑖,𝝀)=𝜆₁,𝑖=1,2,…,𝑛; a value of the coefficient 𝑟(𝝀) close to one will suggest a good fit of the FPLD. The model that causes data to be most like a straight line on its Q-Q plot is the FPLD that most closely resembles the underlying distribution of the data. Correlations that are too small indicate a lack of fit. Finally, 𝑟(𝝀)=1 means that 𝑥_𝑖∶𝑛=𝑋(𝑝_𝑖,𝝀),𝑖=1,2,…,𝑛, that is, the observed data correspond exactly to their expected values.

Table 3 reports our findings on the following datasets:(i)recorded annual drought of Iroquois River recorded near Chebanse (IL), 𝑛=32, Ghosh [36],(ii)times in seconds to serve light vans of the Severn Bridge River crossing in Britain, 𝑛=47, Gupta and Parzen [37],(iii)peak concentrations of accidental releases of toxic gases, 𝑛=100, Hankin and Lee [38],(iv)observations on overall length of bolts, 𝑛=200, Pal [39].

Table 3 

Comparison of FPLD interpolations.

The symbol 𝜈 denotes the number of function evaluations required to find the global optimum of the criterion. The acronym OD refers to the plotting positions 𝑝⁺_𝑖 proposed by Öztürk and Dale [20].

The best overall performance is achieved by NLS; comparably good results are also obtained for DLA and they are only slightly inferior for OD. Based on these few values, the computational difficulties encountered in NLS may be worth the effort. The estimates provided by DLA can be very different from the others because they are based on a different distance metric. The complexity of DLA is higher than that of NLS or OD, due to its more intricate search. The limited impact on the modelling adequacy does not compensate for the extra energy expended for the DLA procedure.

6.2. Simulations

In this section we carry out a simulation study to evaluate the performance of the proposed methods mainly with respect to their biases and squared errors for different sample sizes. In particular, we consider the parameter combination that allows the FPLD to approximate to the standard Gaussian distribution: 𝝀₀=(0.00000,1.35921,0.00000,0.13312,0.13312). (6.3) The vector 𝝀₀ will act as our reference point. The reason we use the Gaussian model is that the normality assumption is one of the most extensively used in statistics and a new method should at least work in this default case.

Two simple coefficients of performance have been considered for comparingMRB(𝜆_𝑗)=𝑁^𝑁−1_𝑖=1|𝜆_𝑖,𝑗−𝜆_0,𝑗𝛼_𝑗|,where𝛼_𝑗=𝜆_0,𝑗if𝜆_0,𝑗≠01if𝜆_0,𝑗,=0RMSE(𝜆_𝑗)=[𝑁^𝑁−1_𝑖=1(𝜆_𝑖,𝑗−𝜆_0,𝑗)²]^0.5. (6.4) The mean relative bias (MRB) quantifies the average relative magnitude of the accuracy for each parameter and the root mean squared error (RMSE) represents the mean Euclidean distance between simulation and measurement. Methods that yield estimates which are closer to the true parameter value have lower bias, higher precision, and lower RMSEs.

The statistics in (6.4) were calculated by generating 𝑁=10,000 different random samples of size 𝑛∈(30,60,120,250,500) from 𝑋(𝑝,𝝀₀). The results are summarized in Table 4

Table 4 

Simulation results.

For all methods, MRB and RMSE decrease as the sample size increases, an indication that all of the estimators are consistent. The sign of the bias for the exponential parameters is almost always negative showing a systematic underestimation of (𝜆₄,𝜆₅) probably due to the fact that the Gaussian has an infinite range (which would imply a FPLD with negative values of the exponential parameters), but it is well approximated by a FPLD with a limited support.

As might have been expected, the bias and the standard error of the exponential parameters (𝜆₄,𝜆₅) yielded little or no discernible differences, particularly for the largest samples. The bias of the NLS technique is generally smaller than those of the other methods. This effect is particularly noticeable in the exponential parameters. The standard error of OD and DLA are generally larger than those for the NLS so that this method can be considered as the best of those looked at in this paper. However, for moderate sample sizes (𝑛>120), OD is nearly as good as NLS, but takes less time to compute, so the choice of method depends on the particular dataset. The DLA solution does not seem to be a real improvement over methods based on the least squares criterion.

7. Conclusions

Statistics professionals are often faced with the problem of fitting a set of data by means of a stochastic approximation. For this purpose, it is essential to select a single parametric family of distributions that offers simple but still realistic descriptions of most of the potential models of the underlying phenomenon. A very useful and tractable parametric model is the five-parameter generalized lambda distribution. By means of the FPLD model we have obtained some satisfactory fitting of empirical data and at least one of the many subordinate models can provide a good approximation to several common distributions.

More specifically, we have analyzed three methods for estimating the parameters of the FPLD (two of which are completely new). The common denominator of all of the methods is the replacement of linear parameters by their least squares or least absolutes estimates, given the value of the nonlinear parameters that, in turn, are determined by using a controlled random search technique. The proposed procedures are simple to apply as they only require a derivative-free optimization over a bounded region. Thus, many obstacles to the use of the FPLD are resolved. The efficiency of the new algorithms has been verified by applications to real and simulated data.

Our findings indicate that NLS estimators outperform least absolute deviation estimators with respect to bias and mean square error. For large sample, however, the difference between NLS and the approximated, but more manageable, scheme OD proposed by Öztürk and Dale [20] becomes imperceptible.

As for future research, we plan to compare different diagnostic tools used to analyze the fit of a quantile model. It would also be interesting to incorporate in the methods described in this paper a weighting scheme (e.g., weights proportional to the local density) so that the tails and/or the middle portion of the samples become more detectable.