A Logistic L-Moment-Based Analog for the Tukey gh , g , h , and hh System of Distributions

This paper introduces a standard logistic L-moment-based system of distributions. The proposed system is an analog to the standard normal conventional moment-based Tukey g-h, g, h, and h-h system of distributions. The system also consists of four classes of distributions and is referred to as i asymmetric γ-κ, ii log-logistic γ , iii symmetric κ, and iv asymmetric κL-κR. The system can be used in a variety of settings such as simulation or modeling events—most notably when heavytailed distributions are of interest. A procedure is also described for simulating γ-κ, γ , κ, and κL-κR distributions with specified L-moments and L-correlations. The Monte Carlo results presented in this study indicate that estimates of L-skew, L-kurtosis, and L-correlation associated with the γ-κ, γ , κ, and κL-κR distributions are substantially superior to their corresponding conventional productmoment estimators in terms of relative bias and relative standard error.


Introduction
Simulating or modeling phenomena that involve heavy-tailed distributions have increasingly become of interest in many areas. Some examples include investigations associated with modeling stream flow and flood frequency 1-3 , regional frequency analysis 4 , aircraft landing processes 5 , disaster analysis 6 , finance 7 , signal and image processing 8 , and latent traits in the social and behavioral sciences 9 . Further, it is also common practice for methodologists to investigate the type I error and power properties associated with inferential statistics using heavy-tailed distributions 10-14 . In many cases, these investigations may only require an elementary transformation to produce heavy-tailed distributions with specified values of conventional skew and kurtosis. A system of transformations that is 2 ISRN Probability and Statistics particularly well suited for this task is the Tukey g-h, g, h, and h-h classes, which can be used for simulating or modeling univariate and multivariate heavy-tailed distributions 14-20 . The quantile function associated with g-h distributions can be succinctly described as where Z ∼ N 0, 1 , the parameters g / 0 and h ≥ 0 control the skew and kurtosis of q Z , and the quantile functions for the one-parameter log-normal g or symmetric h distributions can be obtained by taking the lim h → 0 q Z or lim g → 0 q Z , respectively. The class of asymmetric h-h distributions can be obtained by a doubling technique of h distributions that is described in 18, 21 . Hence, an attractive feature of the Tukey system g-h, g, h, h-h classes is that it is computationally efficient because each class of distributions only requires the knowledge of one or two parameters and an algorithm that generates standard normal random deviates. One problem that arises in the context of heavy-tailed distributions is that conventional moment-based estimators of skew and kurtosis have unfavorable attributes insofar as they can be substantially biased, have high variance, or can be influenced by outliers 4, 22-25 . To obviate this problem, L-moment-based estimators such as L-skew and L-kurtosis were introduced to address the limitations associated with conventional estimators 22 . Specifically, some of the advantages that L-moments have over conventional moments are that they i exist whenever the mean of the distribution exists, ii are nearly unbiased for all sample sizes and distributions, and iii are more robust in the presence of outliers 4, 22-25 . Thus, it would be advantageous to have a L-moment-based Tukey system of transformations for the purpose of simulating or modeling univariate and multivariate heavy-tailed distributions. However, because of the complexities associated with 1.1 and the fact that the standard normal distribution function is not available in closed form, the derivation of a L-moment-based Tukey system would be problematic. We would note that the Tukey h and h-h classes of distributions have been characterized in the context of L-moments and the L-correlation see 21 . In view of the above, the present aim is to introduce an L-moment Tukey system analog that is based on the standard logistic distribution. More specifically, the system consists of four quantile functions that produce continuous symmetric and asymmetric distributions with specified values of L-skew, L-kurtosis, and L-correlation. The four classes of distributions are referred to as i asymmetric γ-κ, ii log-logistic γ, iii symmetric κ, and iv asymmetric κ L -κ R . Some of the advantages that the new system has is that the class of log-logistic γ distributions has demonstrated to be more efficient than its log-normal g counterpart in terms of computing its hazard function or when censored data are encountered in the context of survival analysis 26 . Note that the κ and κ L -κ R classes of distributions were recently introduced in an article by Headrick and Pant 27 .
The rest of the paper is outlined as follows. In Section 2, the cumulative distribution function and probability density function as well as other properties associated with the four classes of distributions are derived. A summary of univariate L-moment theory is also provided, and the derivations of the systems of equations for specifying values of L-skew and L-kurtosis for the four classes of distributions are subsequently provided. In Section 3, the coefficient of L-correlation is introduced, and the equations are developed for determining intermediate correlations for specified L-correlations associated with the four classes of distributions. In Section 4, the steps for implementing a simulation procedure are described. Numerical examples and the results of a simulation are also provided to confirm the derivations and compare the new methodology with its conventional moment-based counterpart. In Section 5, the results of the simulation are discussed.

Definitions and Properties for the System of
γ-κ, γ, κ, and κ L -κ R Distributions The derivation of the cumulative distribution function cdf , probability density function pdf , and other properties associated with the system of γ-κ, γ, κ, and κ L -κ R distributions begins with the following definitions.
Definition 2.1. Let X be a random variable that has a standard logistic distribution with cdf and pdf expressed as Definition 2.2. Let the quantile function q X associated with the system of γ-κ, γ, κ, and κ Lκ R distributions be defined as where q X in 2.3 -2.6 is a strictly monotone increasing function with real-valued parameters γ, κ, κ L , κ R . Equations 2.3 -2.6 produce distributions defined as i asymmetric γ-κ γ / 0, κ ≥ 0 , ii log-logistic γ γ / 0 , iii symmetric κ κ ≥ 0 , and iv asymmetric κ Lκ R κ L ≥ 0, κ R ≥ 0, κ L / κ R . The parameter ±γ in 2.3 and 2.4 controls the degree of skew associated with a distribution. Taking the negative of γ will change the direction of the skew but not its magnitude that is, q −γ,κ X −q γ,κ −X . The parameter κ in 2.3 and 2.5 controls the tail weight of a distribution where the function e κ|X| i preserves symmetry, ii is increasing for X ≥ 0 and κ ≥ 0, and iii produces increased tail weight as the value of κ becomes larger. Analogous to κ, the real-valued parameter κ L κ R in 2.6 controls the left right tail weight of a κ L -κ R distribution.

ISRN Probability and Statistics
Theorem 2.3. The cdf and pdf associated with the γ-κ, γ, κ, and κ L -κ R classes of distributions in 2.3 -2.6 are where q x > 0 and q x 0 1 for 2.9 .
Proof. The requirement that q X in 2.3 -2.6 is a strictly monotone increasing function implies that an inverse function q −1 exists and thus y F q X x Ψ x . Differentiating both sides with respect to q x yields dy/dq x f q x x . Hence, dy/dq x dy/dx / dq x / dx ψ x /q x . Whence, the pdf in 2.8 integrates to one because ψ x is the logistic pdf in 2.2 and will be nonnegative on its support q x for x ∈ −∞, ∞ given from Definition 2.2 that γ / 0, κ ≥ 0, κ L ≥ 0, κ R ≥ 0, and where lim x → ±∞ ψ x /q x 0 as lim x → ±∞ ψ x 0 and lim x → ±∞ 1/q x 0.
Definition 2.4. If the monotonicity assumption q x > 0 holds for the pdf in 2.8 for all x ∈ −∞, ∞ , then 2.8 is defined to be a global pdf.
Remark 2.5. Inspection of 2.2 and 2.3 -2.9 indicates that the height of a global pdf in 2.8 for any γ-κ, γ, κ, or κ L -κ R distribution at x 0 will be ψ x 1/4.
Remark 2.6. The mode associated with a global pdf in 2.8 is located at q x where x x is the critical number that solves d ψ x /q x /dq x 0 and globally maximizes ψ x /q x at the mode q x . It is noted that a global pdf in 2.8 will have a global maximum because the standard logistic pdf in 2.2 has a global maximum, and the transformations in 2.3 -2.6 are assumed to be strictly monotone increasing functions.
Remark 2.7. The median associated with a global pdf in 2.8 is located at q x 0 0. This can be shown by letting Ψ x 0.50 denote the 50th percentile. In general, we must have x 0 such that Ψ 0 Pr{X ≤ 0} 0.50 holds in 2.7 for the standard logistic distribution. As such, the limit of the quantile function q x locates the median at lim x → 0 q x 0 for the system of γ-κ, γ, κ, and κ L -κ R distributions. Remark 2.8. In terms of the κ L -κ R or κ class of distributions, the monotonicity assumption q x > 0 holds locally in 2.8 for cases where κ L,R < 0 and 1/κ L < x < 0 or 0 ≤ x < −1/κ R . This leads to Definition 2.9. Definition 2.9. A κ L -κ R or κ distribution is defined to have a local pdf if the monotonicity assumption holds in accordance with Remark 2.8 and 1 − Pr{q x ≤ |1/κ L,R | } ≤ ε, where ε is a user specified threshold probability based on the cdf in 2.7 e.g., ε ≤ 0.001 .
Examples of γ-κ, γ, κ, and κ L -κ R distributions based on the cdf and pdf in 2.7 and 2.8 are presented in Figure 1 through Figure 4, respectively. The heavy-tailed distributions in Figures 1 a , 2 a , 3 b , and 4 a are used in the simulation portion of this study in Section 4. In the next section, univariate L-moments are introduced, and the L-moments for the system of γ-κ, γ, κ, and κ L -κ R distributions are subsequently derived, and other properties are discussed.

Preliminaries on Univariate L-Moments
Let Y 1 , . . . , Y j , . . . , Y n be independent and identically distributed random variables each with continuous pdf f Y y , cdf F Y y , order statistics denoted as Y 1:n ≤ · · · ≤ Y j:n ≤ · · · ≤ Y n:n , and L-moments defined in terms of either linear combinations of i expectations of order statistics or ii probability weighted moments β i . For the purposes considered herein, the first four L-moments associated with Y j:n are expressed as 4, pages 20-22 where the β i are determined from where i 0, . . . , 3. The coefficients associated with β i in 2.14 are obtained from shifted orthogonal Legendre polynomials and are computed as shown in 4, page 20 or in 22 .
The L-moments λ 1 and λ 2 in 2.10 and 2.11 are measures of location and scale and are the arithmetic mean and one-half the coefficient of mean difference, respectively. Higher order L-moments are transformed to dimensionless quantities referred to as L-moment ratios defined as τ r λ r /λ 2 for r ≥ 3, and where τ 3 and τ 4 are the analogs to the conventional measures of skew and kurtosis. In general, L-moment ratios are bounded in the interval −1 < τ r < 1 as is the index of L-skew τ 3 where a symmetric distribution implies that all L-moment ratios with odd subscripts are zero. Other smaller boundaries can be found for more specific

Pdf characteristics
Parameters Percentiles
Taking the limit of 2.17 -2.20 as κ → 0 gives the class associated with the loglogistic γ distributions in 2.4 as Analogously, taking the limit of 2.17 -2.20 as γ → 0 gives the system associated with the symmetric κ family of distributions in 2.5 as in 27 with polygamma functions P · of p 1 P 0, 1/2 − κ/2 , p 2 P 0, 1 − κ/2 , p 3 P 1, 1/2 − κ/2 , and p 4 P 1, 1 − κ/2 . The L-moments for the class of asymmetric κ L -κ R distributions in 2.6 can be determined by separately evaluating and summing two integrals of the form in 2.16 as For an asymmetric κ L -κ R distribution with a global pdf to have finite values of λ 1 , λ 2 , τ 3 , τ 4 based on 2.29 will require κ L < 1 and κ R < 1. As such, integrating 2.29 to obtain β 0 , . . . , β 3 and subsequently substituting these terms into 2.10 -2.13 yield as in 27 As with the γ-κ class of distributions, given specified values of τ 3 and τ 4 , 2.32 and 2.33 can be numerically solved to obtain the corresponding values of κ L and κ R . Inspection of 2.32 and 2.33 reveals that interchanging values for the parameters of κ L and κ R reverses the direction of τ 3 and has no effect on τ 4 . As such, a graph of the region for feasible combinations of |τ 3 | and τ 4 for 2.32 and 2.33 is provided in Figure 5. Feasible combinations of L-skew and L-kurtosis will lie in the region above the curve graphed in the |τ 3 |, τ 4 plane of Figure 5. Note that the curve in Figure 5 was graphed by setting κ L 0 with κ R ∈ 0, 1 in 2.32 and 2.33 .
The conventional moment-based systems for the γ-κ and κ L -κ R classes of distributions are given in Appendices A and B, respectively. These systems were used to determine the values of skew and kurtosis associated with the distributions given in Figures 1-4. It is worthy to point out that the conventional moment-based systems have a disadvantage in terms of moment existence. That is, the integral in A.1 of Appendix A reveals that for the rth moment to exist we must have in general i γ κ < 1/r, ii κ < 1/r and iii 1/r γ < κ. For example, if the mean, variance, skew, and kurtosis exist, then we must have 1/r 0.25. And, analogously for the κ L -κ R class of distributions in Appendix B, the parameters are bounded in the range of 0 ≤ κ L < 0.25 and 0 ≤ κ R < 0.25. The advantage that the L-moment system has in this context is attributed to Hosking's Theorem 1 22 which states that if the mean λ 1 exists, then all other L-moments will have finite expectations. We note that the conventional moment-based systems for the log-logistic γ and symmetric κ classes of distributions can be obtained by simplifying the system as described in Appendix A. In the next section we first introduce the topic of the L-correlation and subsequently develop the methodology for simulating γ-κ, γ, κ, and κ L -κ R distributions with specified L-correlations.

L-Correlations for the System of γ-κ, γ, κ, and κ L -κ R Distributions
The coefficient of L-correlation 29 is introduced by considering two random variables Y j and Y k with distribution functions F Y j and F Y k , respectively. The second L-moments of Y j and Y k can alternatively be expressed as The second L-comoments of Y j toward Y k and Y k toward Y j are 3.4 As such, the L-correlations of Y j toward Y k and Y k toward Y j are expressed as The L-correlation in 3.5 or 3.6 is bounded such that −1 ≤ η jk ≤ 1 where a value of η jk 1 η jk −1 indicates a strictly increasing decreasing monotone relationship between the two variables. In general, we would also note that η jk / η kj . In the context of the L-moment-based γ-κ, γ, κ, and κ L -κ R classes of distributions, suppose that it is desired to simulate a T -variate distribution based on quantile functions of the forms in 2.3 -2.6 with a specified L-correlation matrix and where each distribution has its own specified values of τ 3 and τ 4 . Let Z 1 , . . . , Z T denote standard normal variables where the distribution functions and bivariate density function associated with Z j and Z k are expressed as Using 3.7 , it follows that the jth distribution associated with 2.3 -2.6 can be expressed as q j g Φ Z j , where g Φ Z j ln Φ Z j / 1 − Φ Z j is standard logistic because Φ Z j ∼ U 0, 1 . As such, using 3.5 , the L-correlation of q j g Φ Z j toward q k g Φ Z k can be evaluated using the solved value s of the parameter s i.e., γ j -κ j ; γ j ; κ j ; κ L j -κ R j for q j g Φ Z j , a specified intermediate correlation IC ρ jk in 3.9 , and the following integral generally expressed as We would point out that the purpose of the IC ρ jk in 3.9 and 3.10 is to adjust for the effect of the transformation q j g Φ Z j , which is induced by the parameters, such that q j g Φ Z j has its specified L-correlation η jk toward q k g Φ Z k . Further, to simplify the computation, the quantile function in 3.10 is standardized by a linear transformation such that it has a mean of zero and one-half the coefficient of mean difference equal to that of the unit-normal distribution as where λ 1 is a mean from 2.

The Procedure for Simulation and Monte Carlo Study
To implement the method for simulating γ-κ, γ, κ, and κ L -κ R distributions with specified L-moments and L-correlations we suggest the following six steps.
1 Specify the L-moments for T transformations of the forms in 2.3 -2.6 that is, and obtain the solutions for the parameters of γκ, γ, κ, or κ L -κ R distributions by solving 2.19 , 2.20 ; 2.23 , 2.24 ; 2.27 , 2.28 ; 2.32 , 2.33 using the specified values of L-skew τ 3 and L-kurtosis τ 4 for each distribution. Specify a T × T matrix of L-correlations η jk for q j g Φ Z j toward q k g Φ Z k , where j < k ∈ {1, 2, . . . , T}. 3 Assemble the ICs into a T × T matrix and decompose this matrix using a Cholesky factorization. Note that this step requires the IC matrix to be positive definite. * Intermediate Correlation for Distributions 1 and 4. See Table 2. * ρ jk 0.835142; q Rk X k * Exp κ Rk * Abs X k ; * Standardizing constants λ 1j , λ 1k 2.17 , 2.30 and α 2j , α 2k A.2 , B.2 from Appendices A and B * Sq γj ,κj q γj ,κj − λ 1j /α 2j ; Sq Lk q Lk − λ 1k /α 2k ; Sq Rk q Rk − λ 1k /α 2k ; * Compute the specified Pearson correlation in Table 1 . . . Z T a 1T V 1 a 2T V 2 · · · a iT V i · · · a jT V j · · · a TT V T , 4.1 Table 1: Specified correlation matrix for the distributions in Figures 1 a , 2 a , 3 b , and 4 a .  where V 1 , . . . , V T are independent standard normal random variables and where a ij represents the element in the ith row and the jth column of the matrix associated with the Cholesky factorization performed in Step 3 .
5 Substitute Z 1 , . . . , Z T from Step 4 into the following Taylor series-based expansion for the standard normal cdf 30 : where φ Z j denotes the standard normal pdf and where the absolute error associated with 4.2 is less than 8 × 10 −16 .
6 Substitute the zero-one uniform deviates, Φ Z j , generated from Step 5 into the T equations of the form of q j g Φ Z j , where g Φ Z j ln Φ Z j / 1 − Φ Z j is standard logistic to generate the γ-κ, γ, κ, and κ L -κ R distributions with the specified L-moments and L-correlations.
To demonstrate the steps above and evaluate the proposed method, a comparison between the proposed L-moment and conventional product-moment-based procedures is subsequently described. Specifically, the heavy-tailed distributions in Figures 1 a , 2 a , 3 b , and 4 a are used as a basis for a comparison using the specified correlation matrix in Table 1.  Tables 2 and 3 give the solved IC matrices for the conventional moment and L-moment-based methods, respectively. See Algorithm 2 for an example of computing ICs for the conventional method. Tables 4 and 5 give the results of the Cholesky decompositions on the IC matrices, which are then used to create Z 1 , . . . , Z 4 with the specified ICs by making use of the formulae given in 4.1 of Step 4 with T 4. The values of Z 1 , . . . , Z 4 are subsequently transformed to Φ Z 1 , . . . , Φ Z 4 using 4.2 and then substituted into equations of the forms in 2.3 -2.6 to produce q 1 g Φ Z 1 , . . . , q 4 g Φ Z 4 for both methods. In terms of the simulation, a Fortran algorithm was written for both methods to generate 25,000 independent sample estimates for the specified parameters of i conventional skew α 3 , kurtosis α 4 , and Pearson correlation ρ * jk and ii L-skew τ 3 , L-kurtosis τ 4 , and L-correlation η jk . All estimates were based on sample sizes of n 25 and n 1000.   The formulae used for computing estimates of α 3,4 were based on Fisher's k-statistics that is, the formulae currently used by most commercial software packages such as SAS, SPSS, and Minitab, for computing indices of skew and kurtosis where α 3,4 0 for the standard normal distribution . The formulae used for computing estimates of τ 3,4 were Headrick's Equations 2.4 and 2.6 25 . The estimate for ρ * jk was based on the usual formula for the Pearson product moment of correlation statistic, and the estimate for η jk was computed based on 3.5 using the empirical forms of the cdfs in 3.1 and 3.3 . The estimates for ρ * jk and η jk were both transformed using Fisher's z transformation. Bias-corrected accelerated bootstrapped average mean estimates, confidence intervals CIs , and standard errors were subsequently obtained for the estimates associated with the parameters α 3,4 , τ 3,4 , z ρ * jk , z η jk using 10,000 resamples via the commercial software package Spotfire S 31 . The bootstrap results for the estimates of the means and CIs associated with z ρ * jk and z η jk were transformed back to their original metrics i.e., estimates for ρ * jk and η jk . Further, if a parameter P was outside its associated bootstrap CI, then an index of relative bias RB was computed for the estimate E as RB E−P /P × 100. Note that the small amount of bias associated with any bootstrap CI containing a parameter was considered negligible and thus not reported. The results of the simulation are reported in Tables 6-13 and are discussed in the next section.

Discussion and Conclusion
One of the primary advantages that L-moments have over conventional moment-based estimators is that they can be far less biased when sampling is from distributions with more severe departures from normality e.g., 21, 25 . Inspection of the simulation results in Tables 6 and 7 of this study clearly indicates that this is also the case for the system of γ-κ, γ, κ, and κ L -κ R distributions. Specifically, the superiority that estimates of L-moment ratios τ 3 , τ 4 have over their corresponding conventional moment-based counterparts α 3 , α 4 is obvious. For example, with samples of size n 25 the estimates of skew and kurtosis for Distribution 1 Tables 6 and 8 were, on average, only 26.80% and 2.38% of their associated population parameters whereas the estimates of L-skew and L-kurtosis were 87.80% and 94.81% of their respective parameters. It is also evident from Tables 6 and 8 that L-skew and L-kurtosis are more efficient estimators as their relative standard errors RSE standard error/estimate × 100 are substantially smaller than the conventional estimators of skew and kurtosis. For example, in terms of Distribution 1, inspection of Tables 7 and 9 indicates RSE measures of  The results associated with the conventional Pearson and L-correlations are presented in Tables 10-13. Overall inspection of these tables indicates that the L-correlation is substantially superior to the Pearson correlation in terms of relative bias. For example, in terms of a moderate correlation Table 10, n 25, ρ * 12 0.40 the relative bias for Distributions 1 and 2 was 9.73% for the Pearson correlation compared to only 1.85% for the L-correlation Table 12, n 25, η 12 0.40 . Further, for large sample sizes Table 13, n 1000 , the L-correlation bootstrap CIs contained all of the population parameters whereas the Pearson correlation CIs contained none of the parameters Table 11 .
In summary, the proposed L-moment-based system of γ-κ, γ, κ, and κ L -κ R distributions is an attractive alternative to the traditional conventional-moment-based system. In particular, the L-moment-based system has distinct advantages when heavy-tailed distributions are of concern. Finally, we would note that Mathematica Version 8.0 30 source code is available from the authors for implementing the L-moment-based method.

A. System of Conventional Moment-Based Equations for γ-κ Distributions
The moments μ r 1,...,4 associated with the class of γ-κ distributions in 2.3 can be determined from      In general, to obtain defined values of μ r based on A.1 we must have i γ κ < 1/r, ii κ < 1/r, and iii 1/r γ < κ. The mean, variance, skew, and kurtosis are defined in general as in 17 A.2 The moments associated with the location and scale parameters in A.2 are