# Comparing Estimation Methods for the FPLD

Authors: Agostino Tarsitano
Journal: Journal of Probability and Statistics
 (2010)

---

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


---

## Body

## 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(1.1)X(p,λ)=λ1+λ22{(1-λ3)[pλ4-1λ4]-(1+λ3)[qλ5-1λ5]},
where 0≤p≤1,q=1-p. X(p,λ) is the quantile function and p∈[0,1]. If λ2≥0 and λ3∈[-1,1], then (1.1) is a continuous and increasing function of p. Here, λ1 controls, albeit not exclusively, the location of an FPLD; the parameter λ2 acts for a multiplier to the translated quantile function [X(p,λ)-λ1] and is, thus, a scale parameter; λ3,λ4,λ5 influence the shape of [X(p,λ)]. The limiting cases λ4→0 and/or λ5→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.2)X(p,λ)=λ1+1λ2{[pλ4-1λ4]-[qλ5-1λ5]},
where the scale parameter λ2 appears as a reciprocal and λ3=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 λ=(λ1,…,λ5) 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 dividedinto 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 Section2 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(2.1)f[X(p,λ)]=λ2-1γ1pλ4-1+γ2qλ5-1,γ1=(1-λ3)2,γ2=(1+λ3)2.The range of a FPLD is given by[X(0,λ),X(1,λ)] and is different for different values of λ; in particular, the extremes are [λ1-λ2γ1,λ1+λ2γ2] if λ4,λ5>0; [λ1-λ2γ1,∞) if λ4>0, λ5≤0; and (-∞,λ1+λ2γ2] if λ4≤0, λ5>0. Hence, the FPLD can model distributions whose support extends to infinity in either or both directions as well as those with finite support.Therth moment of the linear transformation Z=2(X-λ1)/λ2+k1-k2 where k1=(1-λ3)/λ4, k2=(1+λ3)/λ5 is(2.2)E(Zr)=∑j=0r(-1)j(rj)k1(r-j)k2jB[1+(r-j)λ4,1+jλ5].
The symbol B(α,β) with α,β>0 denotes the complete beta function. As a consequence, the rth moment of Z exists and is finite if λ4 and λ5 are greater than r-1.The central moments ofX which can be derived from (2.2) are(2.3)σ2=α1-2k1k2B(λ4+1,λ5+1),γ1=α2σ3-3(k1k2)2σ3[B(2λ4+1,λ5+1)k2-B(λ4+1,2λ5+1)k1],γ2=α3σ4-4(k1k2)3σ4[B(3λ4+1,λ5+1)k22+B(λ4+1,3λ5+1)k12+,-3B(2λ4+1,2λ5+1)2k1k2],αi=k1i+1[(i+1)λ4+1]+(-1)i+1k2i+1[(i+1)λ5+1],i=1,2,3.The shape of the FPLD is sensitive to its parameter values; in fact, the moments ofX(p,λ) are determined by a combination of all of the vector λ elements. If λ4=λ5 and λ3=0, then the FPLD is symmetric about λ1 because, in this case, (1.1) satisfies the condition X(p,λ)=-X(q,λ). If λ4=λ5, then (2.1) is skewed to the left (right) if λ3<0 (λ3>0), which suggests a natural interpretation of λ3 as a parameter which is prevalently related to the asymmetry of X(p,λ). For a FPLD, interchanging λ4 and λ5 and simultaneously changing the sign of λ3 reverses the type of skewness. If λ3=-1 (λ3=1), then λ4 (λ5) indicates the left (right) tail flatness in the sense that the smaller λ4 (λ5) the flatter the left (right) tail of the FPLD density. In practice, λ4 and λ5 capture the tail order on the two sides of the support of X(p,λ).The densities (2.1) are zeromodal if {max(λ4,λ5)>1}∧{min(λ4,λ5)<1}, unimodal with continuous tails if {max(λ4,λ5)<1}, unimodal with truncated tails if {min(λ4,λ5)>2}, U-shaped if 1<λ4,λ5<2, and S-shaped if {max(λ4,λ5)>2}∧{min(λ4,λ5)>1}. Curves corresponding to large positive values of λ4 and λ5 have extreme peakedness and short high tails. For λ2=0 the FPLD degenerates X(p,λ)=λ1 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λ4→0,λ5→0, then X(p,λ) converges on a potentially asymmetric form of the logistic distribution; if λ4→∞ and λ5→0, then X(p,λ) is an exponential distribution whereas for λ4→0,λ5→∞, X(p,λ) becomes a reflected exponential distribution. Furthermore, FPLD fits data containing extreme values well; in fact, for λ4→∞,|λ5|<∞, the FPLD corresponds to the generalized Pareto distribution. For λ5→∞,|λ4|<∞, the FPLD generates the power-function distribution. The rectangular random variable is present in four versions: (λ3=-1,λ4=1), (λ3=1,λ5=1), (λ4=1,λ5=1), and  (λ3=0,λ4=2,λ5=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 X(p) and the FPLD is quantified by the maximum absolute difference (Maxd) between observed and fitted quantiles:(2.4)max1≤i≤500|X(pi)-X(pi,λ)|,wherepi=i501,i=1,2,…,500.
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 ofX(p,λ) 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.
Modelλ1λ2λ3λ4λ5MaxdGaussian(0,1)0.000001.359210.000000.133120.133120.0065Laplace:(0,1)0.000001.323970.00000-0.09556-0.095560.0957Student's t:ν=20.000000.910140.00000-0.45155-0.451550.0973Cauchy(0,1)0.000000.662670.00000-0.99250-0.992500.1673Beta(0.5,0.5)0.000001.632020.000001.454561.454560.0593Our 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? Table2 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.
Modelλ1λ2λ3λ4λ5MaxdChi-squareν=31.100102.843130.615050.485790.013380.00462.018221.075700.97596-0.033940.0827Lognormal(4,2)3.445601.073050.161470.210080.047760.00483.565913.640670.267150.015500.0193Weibull(1,3)0.925650.50134-0.153890.308750.125790.00130.874538.107800.263470.164780.0052Gumbel(0,1)-0.000501.502390.329840.19211-0.000460.00200.344062.447440.33889-0.048450.0549Inv. Gaussian(0.5,6)0.455210.186060.174940.257860.029610.00140.4772920.644100.328170.000390.0044The most important point to note in Tables1 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 sampleS={x1,x2,…,xn} of n-independent observations from X(p,λ), the ith order statistic may be expressed as(3.1)xi:n=E(xi:n)+ei,i=1,2,…,n,
in which the ordered observations are compared with their expected values under the hypothesized distribution X(p,λ) and ei is a measure of the discrepancy between the observed and modelled ith value. The error terms {ei,i=1,2,…,n} are such that E(ei)=0 and σ2(ei)=σi2. 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 theith order statistic from a FPLD is available in closed form(3.2)E(xi:n)=λ1+(β1/λ4)∫01pλ4+i-1qn-idp+(β2/λ5)[1-∫01pi-1qλ5+n-idp]B(i,n+1-i),
where β1=0.5(1-λ3)λ2,β2=0.5(1+λ3)λ2. It follows that the deterministic component of (3.1) can be written as(3.3)x¯i:n(λ)=β0U0,i+β1U1,i+β2U2,i,i=1,2,…,n,
with β0=λ1 and(3.4)U0,i=1,U1,i=λ4-1[Γ(n+1)Γ(i+λ4)Γ(i)Γ(n+1+λ4)-1],U2,i=λ5-1[1-Γ(n+1)Γ(n+1-i+λ5)Γ(n+1-i)Γ(n+1+λ5)].The least squares approach calls forβ in (3.3) to be chosen so as to minimize(3.5)S2(λ)=∑i=1n[xi:n-x¯i:n(λ)]2.
The S2(λ) is linear in β1,β2,β3 but not in (λ4,λ5). According to Lawton and Sylvestre [22] the solution to (3.5) can be obtained by fixing the value of (λ4,λ5) and, then, applying the ordinary least squares to solve the linear problem(3.6)β̂(λ4,λ5)=(UλtUλ)-1Uλtx
for β, where Uλ denotes the (n×3) matrix with elements given by (3.4) and x is the (n×1) vector of the ordered observations. The conversion of β0,β1,β2 into λ1,λ2,λ3 is straightforward:(3.7)λ̂1=β̂0,λ̂2=β̂10.5(1-λ̂3)=β̂20.5(1+λ̂3),λ̂3=1-(β̂1/β̂2)1+(β̂1/β̂2).Öztürk and Dale [20] used an approximation of E(xi:n) to overcome the problems created by repeated use of the gamma function in (3.4). More specifically, they used E(xi:n)=X(pi+,λ), with the plotting positions pi+=i(n+1)-1,i=1,2,…,n. In this case, the pseudoregressors are(3.8)V0,i=1,V1,i=λ4-1[(pi+)λ4-1],V2,i=λ5-1[(pn+1-i+)λ5-1].
Parameter estimates are then obtained as in (3.5) through (3.7) with Vj,i in place of Uj,i for j=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 (λ4,λ5). We might try to find a better estimate of (λ4,λ5) 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 S2(λ) 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 E(xi:n) for all distributions; (2) for almost all distributions E(xi:n) 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 E(xi:n) may not always be defined. All three of these drawbacks are avoided in general by choosing to measure the location of the ith 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 ei=xi:n-X(pi,λ),i=1,2,…,n. 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 constantsy∈Rn which minimizes the objective function(4.1)miny∈Rn∑i=1n|xi:n-yi|,
where xi:n is the ith order statistic of the sample S={x1,x2,…,xn} from the quantile function X(p,λ). Rao [24] observed that the solution to this problem is the vector of marginal medians, that is, the ith element of y is the median of the order statistics xi:n.The quantile function (see, e.g., Gilchrist [2, page 86]) of the ith order statistic for a sample of size n from an FPLD is X[pi*(p),λ] with(4.2)pi*(p)=B-1(p,i,n+1-i),
where B-1(·) is the inverse of the incomplete beta function. Consequently, the median of the ith order statistic can be computed using(4.3)pi*=X[B-1(0.5,i,n+1-i),λ],i=1,2,…,n.
The quantities in (4.3) are a solution of (4.1) in the sense that, for a fixed value of λ, the following expression(4.4)∑i=1n|xi:n-X[pi*,λ]|
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 (λ1,λ2,λ3), and nonlinear parameters (λ4,λ5), another two-stage procedure can be devised. In the inner stage, we determine the vector β̃ which minimizes(4.5)S1(λ)=∑i=1n|xi:n-x̃i:n(λ)|,withx̃i:n(λ)=∑j=02βjWi,j,
with respect of the linear parameters β and for a given value of the nonlinear parameters (λ4,λ5). In this case, the pseudo regressors are(4.6)W0,i=1,W1,i=λ4-1{[pi*]λ4-1},W2,i=λ5-1{[1-pi*]λ5-1}.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 β0,β1,β2 into λ1,λ2,λ3 takes place as shown in (3.7). Once (λ1,λ2,λ3) have been estimated, we can obtain a new value for (λ4,λ5) 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 ofpi* has to be executed just once. In addition, the fact that B(p,i,n+1-i)=1-B(1-p,n+1-i,i) implies that pn+1-i*=1-pi* 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 functionsS1(λ) and S2(λ) of two unknowns, λ4,λ5, 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 storageA 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, f(x), 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 f(x) 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. LetC⊂Rk be the set of admissible values of the k=2 nonlinear parameters of the FPLD model. The most promising region that we have found in our applications is the rectangle delimited by(5.1)C={(λ4,λ5)∣-0.999≤λ4≤3;-0.999≤λ5≤3}.
Furthermore, C ensures that the mean of the order statistics from a FPLD always exists and the ranges of λ4,λ5 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 setA of m0>(k+1) random admissible points in C. Find the minimum and the maximum in A and their function values: (xL,fL),(xH,fH).(2)
If|fH-fL|/(10-7+|fH|+|fL|)<10-8, 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 A.(3)
Choose at random, without duplication, a setB of k points in (A-xL), include xL in B, and compute the centroid x¯=(k+1)-1∑x∈Bx. Define A1=A-B.(4)
Select randomly a pointxr∈A1.(5)
Generate uniform random numbersd1,d2,…,dk in [0,2] and compute the trial point xi,T=(1+di)xi¯-dixr,i,fori=1,2,…,k; thus, the exploration of the space around the centroid is based on a randomized reflection in the simplex originally used by Price [26].(6)
IfxT∈M, then compute fT=f(xT) and go to step (7). If xT∉M, then repeat step (5) a fixed number m1 of times. If xT is still unfeasible, then exclude xr from A1, redefine the set A1 and, if it is not empty, repeat step (4) for at most m2 times; otherwise go to step (9).(7)
IffL<fT<fH, then replace fH with fT and xH with xT and then redetermine the maximum xH and fH.(8)
IffT<fL, then replace fH with fL, xH with xL, fL with fT, xL with xT, and xH and fH.(9)
Randomly choose two distinct pointsx2,x3 in (A-xL) and compute
(5.2)dQ,i=(x2,i-x3,i)fL+(x3,i-xL,i)f2+(xL,i-x2,i)f3,
for i=1,2,…,k.(10)
If|dQ,i|>10-7 for each i, then build the following quadratic interpolation of the objective function:
(5.3)xQ,i=[(x2,i-x3,i)2fL+(x3,i-xL,i)2f2+(xL,i-x2,i)2f32dQ,i]
for i=1,2,…,k. If some of the |dQ,i|≤10-7, then generate k uniform random numbers di,i=1,2,…,k, in [0,1] and compute
(5.4)xi,Q=dixi,L+(1-di)xi,H,fori=1,2,…,k.(11)
IfxQ∉M, then repeat steps (9)-(10) no more than m3 times. If xQ is still unfeasible, then return to step (2).(12)
Execute steps(7)–(9) with xQ in place of xT and fQ in place of fT and then return to step (2).We experimented by using the CRS algorithm, withx=λ, f(x)=S2(λ), m0=30k, m1=m2=m3=4(k+1), and k=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 C 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 valuesxi:n centered on the arithmetic mean x¯ and the hypothesized values(6.1)r(λ)=∑i=1n[xi:n-x¯][X(pi,λ)-μ(λ)]{∑i=1n[xi:n-x¯]2∑i=1n[X(pi,λ)-μ(λ)]2}0.5,
where pi=i/(n+1), and(6.2)μ(λ)=∫01X(p,λ)dp=λ1+λ22[1+λ31+λ5-1-λ31+λ4].
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 xi:n and X(pi,λ) are both arranged in ascending order of magnitude. More specifically, r(λ) takes values in the [0,1] interval. The case r(λ)=0 corresponds to a noninformative model X(pi,λ)=λ1,i=1,2,…,n; a value of the coefficient r(λ) 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, r(λ)=1 means that xi:n=X(pi,λ),i=1,2,…,n, that is,  the observed data correspond exactly to their expected values.Table3 reports our findings on the following datasets:(i)
recorded annual drought of Iroquois River recorded near Chebanse (IL),n=32, Ghosh [36],(ii)
times in seconds to serve light vans of the Severn Bridge River crossing in Britain,n=47, Gupta and Parzen [37],(iii)
peak concentrations of accidental releases of toxic gases,n=100, Hankin and Lee [38],(iv)
observations on overall length of bolts,n=200, Pal [39].Table 3
Comparison of FPLD interpolations.
Methodλ1λ2λ3λ4λ5νr(λ)Ghosh [36]NLS13.555348.53610.9978-0.54010.185112830.7812OD14.708646.50570.9917-0.63840.114011580.3821DLA15.221751.2328091600.95850.203692870.3771Gupta and Parzen [37]NLS5.43702.9152-0.15440.6544-0.109710970.8468OD6.32233.4372-0.60180.8189-0.403111390.9100DLA4.98642.65330.10580.5644-0.089526730.7579Hankin and Lee [38]NLS5.48824.62590.6731-0.01520.202510950.6982OD5.54274.44530.6759-0.08960.162010490.5605DLA5.93924.02340.5977-0.07640.093553950.6147Pal [39]NLS6.39810.19820.7108-0.05180.21849610.6771OD6.39710.19330.7330-0.12660.198310820.5337DLA6.41830.18030.5849-0.00870.132922720.6142The symbolν denotes the number of function evaluations required to find the global optimum of the criterion. The acronym OD refers to the plotting positions pi+ 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:(6.3)λ0=(0.00000,1.35921,0.00000,0.13312,0.13312).
The vector λ0 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 comparing(6.4)MRB(λj)=N-1∑i=1N|λi,j-λ0,jαj|,whereαj={λ0,jifλ0,j≠01ifλ0,j=0,RMSE(λj)=[N-1∑i=1N(λi,j-λ0,j)2]0.5.
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 N=10,000 different random samples of size n∈(30,60,120,250,500) from X(p,λ0). The results are summarized in Table 4Table 4
Simulation results.
MethodIndexnλ1λ2λ3λ4λ5NLSMRB300.00490.5616-0.00420.36270.338260-0.00610.19940.0037-0.1199-0.07561200.00040.0789-0.0015-0.0565-0.06872500.00140.0342-0.0017-0.0205-0.0265500-0.00210.01760.0027-0.01100.0035OD300.00330.5080-0.0055-0.7632-0.769760-0.01490.17820.0140-1.0219-0.89671200.00040.0636-0.0013-0.5128-0.5173250-0.00470.01900.0055-0.2752-0.23955000.00210.0035-0.0016-0.1402-0.1441DLA300.00830.5739-0.0047-0.0693-0.079160-0.00460.20450.0020-0.4196-0.44491200.00720.0702-0.0082-0.2173-0.26102500.00580.0280-0.0058-0.0898-0.11745000.00080.0111-0.0014-0.0495-0.0535NLSRMSE300.88142.14700.70790.75720.7578600.67290.98530.59840.40330.40261200.47280.36860.46020.21400.21282500.32250.15670.33400.13000.13105000.22080.10040.23690.08790.0875OD300.86102.12560.74460.86170.8664600.72570.93550.66280.51210.51041200.52900.32620.52260.28320.28232500.35670.16300.37650.16360.16305000.24100.10130.26300.10270.1029DLA300.87332.17300.74050.85940.8580600.70401.01390.63580.47750.47841200.50580.30630.49510.24120.24262500.33780.16060.35370.14470.14555000.23140.10110.24970.09440.0946For 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(λ4,λ5) 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(λ4,λ5) 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 (n>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.


## 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 valuesxi:n centered on the arithmetic mean x¯ and the hypothesized values(6.1)r(λ)=∑i=1n[xi:n-x¯][X(pi,λ)-μ(λ)]{∑i=1n[xi:n-x¯]2∑i=1n[X(pi,λ)-μ(λ)]2}0.5,
where pi=i/(n+1), and(6.2)μ(λ)=∫01X(p,λ)dp=λ1+λ22[1+λ31+λ5-1-λ31+λ4].
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 xi:n and X(pi,λ) are both arranged in ascending order of magnitude. More specifically, r(λ) takes values in the [0,1] interval. The case r(λ)=0 corresponds to a noninformative model X(pi,λ)=λ1,i=1,2,…,n; a value of the coefficient r(λ) 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, r(λ)=1 means that xi:n=X(pi,λ),i=1,2,…,n, that is,  the observed data correspond exactly to their expected values.Table3 reports our findings on the following datasets:(i)
recorded annual drought of Iroquois River recorded near Chebanse (IL),n=32, Ghosh [36],(ii)
times in seconds to serve light vans of the Severn Bridge River crossing in Britain,n=47, Gupta and Parzen [37],(iii)
peak concentrations of accidental releases of toxic gases,n=100, Hankin and Lee [38],(iv)
observations on overall length of bolts,n=200, Pal [39].Table 3
Comparison of FPLD interpolations.
Methodλ1λ2λ3λ4λ5νr(λ)Ghosh [36]NLS13.555348.53610.9978-0.54010.185112830.7812OD14.708646.50570.9917-0.63840.114011580.3821DLA15.221751.2328091600.95850.203692870.3771Gupta and Parzen [37]NLS5.43702.9152-0.15440.6544-0.109710970.8468OD6.32233.4372-0.60180.8189-0.403111390.9100DLA4.98642.65330.10580.5644-0.089526730.7579Hankin and Lee [38]NLS5.48824.62590.6731-0.01520.202510950.6982OD5.54274.44530.6759-0.08960.162010490.5605DLA5.93924.02340.5977-0.07640.093553950.6147Pal [39]NLS6.39810.19820.7108-0.05180.21849610.6771OD6.39710.19330.7330-0.12660.198310820.5337DLA6.41830.18030.5849-0.00870.132922720.6142The symbolν denotes the number of function evaluations required to find the global optimum of the criterion. The acronym OD refers to the plotting positions pi+ 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:(6.3)λ0=(0.00000,1.35921,0.00000,0.13312,0.13312).
The vector λ0 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 comparing(6.4)MRB(λj)=N-1∑i=1N|λi,j-λ0,jαj|,whereαj={λ0,jifλ0,j≠01ifλ0,j=0,RMSE(λj)=[N-1∑i=1N(λi,j-λ0,j)2]0.5.
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 N=10,000 different random samples of size n∈(30,60,120,250,500) from X(p,λ0). The results are summarized in Table 4Table 4
Simulation results.
MethodIndexnλ1λ2λ3λ4λ5NLSMRB300.00490.5616-0.00420.36270.338260-0.00610.19940.0037-0.1199-0.07561200.00040.0789-0.0015-0.0565-0.06872500.00140.0342-0.0017-0.0205-0.0265500-0.00210.01760.0027-0.01100.0035OD300.00330.5080-0.0055-0.7632-0.769760-0.01490.17820.0140-1.0219-0.89671200.00040.0636-0.0013-0.5128-0.5173250-0.00470.01900.0055-0.2752-0.23955000.00210.0035-0.0016-0.1402-0.1441DLA300.00830.5739-0.0047-0.0693-0.079160-0.00460.20450.0020-0.4196-0.44491200.00720.0702-0.0082-0.2173-0.26102500.00580.0280-0.0058-0.0898-0.11745000.00080.0111-0.0014-0.0495-0.0535NLSRMSE300.88142.14700.70790.75720.7578600.67290.98530.59840.40330.40261200.47280.36860.46020.21400.21282500.32250.15670.33400.13000.13105000.22080.10040.23690.08790.0875OD300.86102.12560.74460.86170.8664600.72570.93550.66280.51210.51041200.52900.32620.52260.28320.28232500.35670.16300.37650.16360.16305000.24100.10130.26300.10270.1029DLA300.87332.17300.74050.85940.8580600.70401.01390.63580.47750.47841200.50580.30630.49510.24120.24262500.33780.16060.35370.14470.14555000.23140.10110.24970.09440.0946For 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(λ4,λ5) 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(λ4,λ5) 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 (n>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.



---

*Source: 295042-2010-06-14.xml*