# A LogisticL-Moment-Based Analog for the Tukey g-h, g, h, and h-h System of Distributions

Authors: Todd C. Headrick; Mohan D. Pant
Journal: ISRN Probability and Statistics
 (2012)

---

## Abstract

This paper introduces a standard logisticL-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 κℒ-κℛ. The system can be used in a variety of settings such as simulation or modeling events—most notably when heavy-tailed distributions are of interest. A procedure is also described for simulating γ-κ, γ, κ, and κℒ-κℛ 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 κℒ-κℛ distributions are substantially superior to their corresponding conventional product-moment estimators in terms of relative bias and relative standard error.


---

## Body

## 1. 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 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 withg-h distributions can be succinctly described as
(1.1)q(Z)=g-1(egZ-1)ehZ2,
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 limh→0q(Z) or limg→0q(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 aL-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 anL-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 κℒ-κℛ. 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 κℒ-κℛ classes of distributions were recently introduced in an article by Headrick and Pant [27].The rest of the paper is outlined as follows. In Section2, 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.


## 2. Methodology

### 2.1. Definitions and Properties for the System ofγ-κ, γ, κ, and κℒ-κℛ Distributions

The derivation of the cumulative distribution function (cdf), probability density function (pdf), and other properties associated with the system ofγ-κ, γ, κ, and κℒ-κℛ distributions begins with the following definitions.Definition 2.1.
LetX be a random variable that has a standard logistic distribution with cdf and pdf expressed as
(2.1)FX(x)=Ψ(x)=Pr(X≤x)=1(1+e-x),(2.2)fX(x)=ψ(x)=e-x(1+e-x)2,-∞<x<+∞.Definition 2.2.
Let the quantile functionq(X) associated with the system of γ-κ, γ, κ, and κℒ-κℛ distributions be defined as
(2.3)q(X)=qγ,κ(X)=γ-1(eγX-1)eκ|X|,(2.4)q(X)=qγ,0(X)=limκ→0qγ,κ(X)=γ-1(eγX-1),(2.5)q(X)=q0,κ(X)=limγ→0qγ,κ(X)=Xeκ|X|,(2.6)q(X)=qκℒ,κℛ(X)={Xeκℒ|X|,forX≤0,Xeκℛ|X|,forX≥0,
where q(X) in (2.3)–(2.6) is a strictly monotone increasing function with real-valued parameters γ, κ, κℒ, κℛ. Equations (2.3)–(2.6) produce distributions defined as (i) asymmetric γ-κ (γ≠0, κ≥0), (ii) log-logistic γ (γ≠0), (iii) symmetric κ (κ≥0), and (iv) asymmetric κℒ-κℛ (κℒ≥0, κℛ≥0, κℒ≠κℛ). 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 κℒ (κℛ) in (2.6) controls the left (right) tail weight of a κℒ-κℛ distribution.Theorem 2.3.
The cdf and pdf associated with theγ-κ, γ, κ, and κℒ-κℛ classes of distributions in (2.3)–(2.6) are
(2.7)Fq(X)(x)=Ψ(x),(2.8)fq(X)(x)=ψ(x)q′(x)
with derivatives q′(x) expressed as
(2.9)q′(x)=qγ,κ′(x)=eγx+κ|x|+γ-1(eκ|x|(eγx-1)κd|x|dx),q′(x)=qγ,0′(x)=eγx,q′(x)=q0,κ′(x)=eκ|x|+xκeκ|x|d|x|dx,q′(x)=qκℒ,κℛ′(x)={eκℒ|x|-xκℒeκℒ|x|,forx≤0,eκℛ|x|+xκℛeκℛ|x|,forx≥0,
where q′(x)>0 and q′(x=0)=1 for (2.9).Proof.
The requirement thatq(X) in (2.3)–(2.6) is a strictly monotone increasing function implies that an inverse function (q-1) exists and thus y=Fq(X)(x)=Ψ(x). Differentiating both sides with respect to q(x) yields dy/dq(x)=fq(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, κℒ≥0, κℛ≥0, and where limx→±∞ψ(x)/q′(x)=0 as limx→±∞ψ(x)=0 and limx→±∞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 κℒ-κℛ 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 limx→0q(x)=0 for the system of γ-κ, γ, κ, and κℒ-κℛ distributions.Remark 2.8.
In terms of theκℒ-κℛ (or κ) class of distributions, the monotonicity assumption (q′(x)>0) holds locally in (2.8) for cases where κℒ,ℛ<0 and (1/κℒ<x<0 or 0≤x<-1/κℛ). This leads to Definition 2.9.Definition 2.9.
Aκℒ-κℛ (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/κℒ,ℛ|)}≤ε, where ε is a user specified threshold probability based on the cdf in (2.7) (e.g., ε≤0.001).
Examples ofγ-κ, γ, κ, and κℒ-κℛ 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 κℒ-κℛ distributions are subsequently derived, and other properties are discussed.Two asymmetricγ-κ pdfs and their cdfs. The parameters of skew (α3) and kurtosis (α4) are based on (A.2) in Appendix A. The parameters of L-skew (τ3) and L-kurtosis (τ4) are based on (2.19) and (2.20).
(a)
(b)Two log-logisticγ pdfs and their cdfs. The parameters of skew (α3) and kurtosis (α4) are based on (A.2) in Appendix A. The parameters of L-skew (τ3) and L-kurtosis (τ4) are based on (2.23) and (2.24).
(a)
(b)Two symmetricκ pdfs and their cdfs. The parameters of skew (α3) and kurtosis (α4) are based on (A.2) in Appendix A. The parameters of L-skew (τ3) and L-kurtosis (τ4) are based on (2.27) and (2.28).
(a)
(b)Two asymmetricκℒ-κℛ distributions pdfs and their cdfs. The parameters of skew (α3) and kurtosis (α4) are based on (A.2) in Appendix A. The parameters of L-skew (τ3) and L-kurtosis (τ4) are based on (2.32) and (2.33). The distribution in (b) is a local pdf in the range of q(x)∈(-8.085,+∞) where (2.7) yields 1-Pr{q(x)≤|-8.085|}=0.0003.
(a)
(b)


### 2.2. Preliminaries on UnivariateL-Moments

LetY1,…,Yj,…,Yn be independent and identically distributed random variables each with continuous pdf fY(y), cdf FY(y), order statistics denoted as Y1:n≤⋯≤Yj:n≤⋯≤Yn: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 Yj:n are expressed as [4, pages 20–22]
(2.10)λ1=E[Y1:1]=β0,(2.11)λ2=12E[Y2:2-Y1:2]=2β1-β0,(2.12)λ3=13E[Y3:3-2Y2:3+Y1:3]=6β2-6β1+β0,(2.13)λ4=14E[Y4:4-3Y3:4+3Y2:4-Y1:4]=20β3-30β2+12β1-β0,
where the βi are determined from
(2.14)βi=∫y{FY(y)}ifY(y)dy,
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].TheL-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 cases. For example, the index of L-kurtosis (τ4) has the boundary condition for continuous distributions of [28]
(2.15)5τ32-14<τ4<1.


### 2.3.L-Moments for the System of γ-κ, γ, κ, and κℒ-κℛ Distributions

The derivation of the first fourL-moments associated with the γ-κ class of distributions begins by defining the probability weighted moments based on (2.14) in terms of (2.1)–(2.3) as
(2.16)βi=∫-∞+∞q(x,γ,κ){Ψ(x)}iψ(x)dx.In general, to obtain finite values ofλ1, λ2, τ3, τ4 based on (2.16), we must have (i) γ+κ<1, (ii) κ<1, and (iii) 1+γ>κ. As such, integrating (2.16) for i=0,1,2,3 and using (2.10)–(2.13) gives λ1, λ2, τ3, and τ4 as
(2.17)λ1=((-γ-κ)h1+(γ-κ)h2+(-γ+κ)h3+2κh4+(γ+κ)h5-2κh6)(2γ),(2.18)λ2=(2γ-(γ+κ)2h1+(γ-κ)2(h2-h3)+(γ+κ)2h5)(2γ),(2.19)τ3=(2γ2-(γ+κ)3h1+(γ-κ)3h2-(γ-κ)3h3+2κ3h4+(γ+κ)3h5-2κ3h6)2γλ2,(2.20)τ4=(2γ(1+5γ2+5κ(1-3κ))-(γ+κ)2(1+5(γ+κ)2)h1+(1+5(γ-κ)2)(γ-κ)2(h2-h3)+(γ+κ)2(1+5(γ+κ)2)h5/(12γλ2))
with harmonic number functions H[·] of h1=H[(1/2)(-1-γ-κ)], h2=H[(1/2)(-1+γ-κ)], h3=H[(1/2)(γ-κ)], h4=H[(1/2)(-1-κ)], h5=H[(-1/2)(γ+κ)], and h6=H[-(1/2)κ]. As such, given user specied values of τ3 and τ4, (2.19) and (2.20) can be numerically solved to obtain the parameters for γ and κ.Taking the limit of (2.17)–(2.20) as κ→0 gives the class associated with the log-logistic γ distributions in (2.4) as
(2.21)λ1=πcsc(γπ)-1γ,(2.22)λ2=γπcsc(γπ),(2.23)τ3=γ,(2.24)τ4=16(1+5γ2).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]
(2.25)λ1=0,(2.26)λ2=12(2+κ(-4p1+4p2+κ(p3-p4))),(2.27)τ3=0,(2.28)τ4=16(1+5κ2+10κ(1+2κ+2κ2(h6-h4))2λ2)
with polygamma functions P[·] of p1=P[0,1/2-κ/2],  p2=P[0,1-κ/2],  p3=P[1,1/2-κ/2], and p4=P[1,1-κ/2].TheL-moments for the class of asymmetric κℒ-κℛ distributions in (2.6) can be determined by separately evaluating and summing two integrals of the form in (2.16) as
(2.29)βi=∫-∞0q(x,κℒ){Ψ(x)}iψ(x)dx+∫0+∞q(x,κℛ){Ψ(x)}iψ(x)dx.For an asymmetricκℒ-κℛ distribution with a global pdf to have finite values of λ1, λ2, τ3, τ4 based on (2.29) will require κℒ<1 and κℛ<1. As such, integrating (2.29) to obtain β0,…,β3 and subsequently substituting these terms into (2.10)–(2.13) yield as in [27]
(2.30)λ1=14(2p5-2p6-2p7+2p8-κℒp9+κℒp10+κℛp11-κℛp12),(2.31)λ2=14(4+κℒ(-4p5+4p6+κℒ(p9-p10))+4+κℛ(-4p7+4p8+κℛ(p11-p12))),(2.32)τ3={4κℛ-4κℒ+κℒ2(6p5-6p6+κℒ(p10-p9))+κℛ2(-6p7+6p8+κℛ(p11-p12))}(4λ2),(2.33)τ4={2(2+5κℒ(1+3κℒ)+5κℛ(1+3κℛ))-4(κℒ+10κℒ3)p5-4(κℛ+10κℛ3)p7+5κℒ((4+40κℒ2)p6+κℒ(1+5κℒ2)(p9-p10)+κℛ((4+40κℛ2)p8+κℛ(1+5κℛ2)(p11-p12)))}/(24λ2)
with polygamma functions P[·] of p5=P[0,1/2-κℒ/2], p6=P[0,1-κℒ/2], p7=P[0,1/2-κℛ/2], p8=P[0,1-κℛ/2], p9=P[1,1/2-κℒ/2], p10=P[1,1-κℒ/2], p11=P[1,1/2-κℛ/2], p12=P[1,1-κℛ/2].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 κℒ and κℛ. Inspection of (2.32) and (2.33) reveals that interchanging values for the parameters of κℒ and κℛ 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 κℒ=0 with κℛ∈[0,1] in (2.32) and (2.33).Figure 5
Graph of the region for feasible combinations of (absolute value)L-skew |τ3| and L-kurtosis τ4 for global κℒ-κℛ pdfs. Asymmetric κℒ-κℛ distributions will lie in the area above the curve graphed in the |τ3| and τ4 plane.The conventional moment-based systems for theγ-κ and κℒ-κℛ 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 κℒ-κℛ class of distributions in Appendix B, the parameters are bounded in the range of 0≤κℒ<0.25 and 0≤κℛ<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 κℒ-κℛ distributions with specified L-correlations.


## 2.1. Definitions and Properties for the System ofγ-κ, γ, κ, and κℒ-κℛ Distributions

The derivation of the cumulative distribution function (cdf), probability density function (pdf), and other properties associated with the system ofγ-κ, γ, κ, and κℒ-κℛ distributions begins with the following definitions.Definition 2.1.
LetX be a random variable that has a standard logistic distribution with cdf and pdf expressed as
(2.1)FX(x)=Ψ(x)=Pr(X≤x)=1(1+e-x),(2.2)fX(x)=ψ(x)=e-x(1+e-x)2,-∞<x<+∞.Definition 2.2.
Let the quantile functionq(X) associated with the system of γ-κ, γ, κ, and κℒ-κℛ distributions be defined as
(2.3)q(X)=qγ,κ(X)=γ-1(eγX-1)eκ|X|,(2.4)q(X)=qγ,0(X)=limκ→0qγ,κ(X)=γ-1(eγX-1),(2.5)q(X)=q0,κ(X)=limγ→0qγ,κ(X)=Xeκ|X|,(2.6)q(X)=qκℒ,κℛ(X)={Xeκℒ|X|,forX≤0,Xeκℛ|X|,forX≥0,
where q(X) in (2.3)–(2.6) is a strictly monotone increasing function with real-valued parameters γ, κ, κℒ, κℛ. Equations (2.3)–(2.6) produce distributions defined as (i) asymmetric γ-κ (γ≠0, κ≥0), (ii) log-logistic γ (γ≠0), (iii) symmetric κ (κ≥0), and (iv) asymmetric κℒ-κℛ (κℒ≥0, κℛ≥0, κℒ≠κℛ). 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 κℒ (κℛ) in (2.6) controls the left (right) tail weight of a κℒ-κℛ distribution.Theorem 2.3.
The cdf and pdf associated with theγ-κ, γ, κ, and κℒ-κℛ classes of distributions in (2.3)–(2.6) are
(2.7)Fq(X)(x)=Ψ(x),(2.8)fq(X)(x)=ψ(x)q′(x)
with derivatives q′(x) expressed as
(2.9)q′(x)=qγ,κ′(x)=eγx+κ|x|+γ-1(eκ|x|(eγx-1)κd|x|dx),q′(x)=qγ,0′(x)=eγx,q′(x)=q0,κ′(x)=eκ|x|+xκeκ|x|d|x|dx,q′(x)=qκℒ,κℛ′(x)={eκℒ|x|-xκℒeκℒ|x|,forx≤0,eκℛ|x|+xκℛeκℛ|x|,forx≥0,
where q′(x)>0 and q′(x=0)=1 for (2.9).Proof.
The requirement thatq(X) in (2.3)–(2.6) is a strictly monotone increasing function implies that an inverse function (q-1) exists and thus y=Fq(X)(x)=Ψ(x). Differentiating both sides with respect to q(x) yields dy/dq(x)=fq(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, κℒ≥0, κℛ≥0, and where limx→±∞ψ(x)/q′(x)=0 as limx→±∞ψ(x)=0 and limx→±∞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 κℒ-κℛ 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 limx→0q(x)=0 for the system of γ-κ, γ, κ, and κℒ-κℛ distributions.Remark 2.8.
In terms of theκℒ-κℛ (or κ) class of distributions, the monotonicity assumption (q′(x)>0) holds locally in (2.8) for cases where κℒ,ℛ<0 and (1/κℒ<x<0 or 0≤x<-1/κℛ). This leads to Definition 2.9.Definition 2.9.
Aκℒ-κℛ (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/κℒ,ℛ|)}≤ε, where ε is a user specified threshold probability based on the cdf in (2.7) (e.g., ε≤0.001).
Examples ofγ-κ, γ, κ, and κℒ-κℛ 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 κℒ-κℛ distributions are subsequently derived, and other properties are discussed.Two asymmetricγ-κ pdfs and their cdfs. The parameters of skew (α3) and kurtosis (α4) are based on (A.2) in Appendix A. The parameters of L-skew (τ3) and L-kurtosis (τ4) are based on (2.19) and (2.20).
(a)
(b)Two log-logisticγ pdfs and their cdfs. The parameters of skew (α3) and kurtosis (α4) are based on (A.2) in Appendix A. The parameters of L-skew (τ3) and L-kurtosis (τ4) are based on (2.23) and (2.24).
(a)
(b)Two symmetricκ pdfs and their cdfs. The parameters of skew (α3) and kurtosis (α4) are based on (A.2) in Appendix A. The parameters of L-skew (τ3) and L-kurtosis (τ4) are based on (2.27) and (2.28).
(a)
(b)Two asymmetricκℒ-κℛ distributions pdfs and their cdfs. The parameters of skew (α3) and kurtosis (α4) are based on (A.2) in Appendix A. The parameters of L-skew (τ3) and L-kurtosis (τ4) are based on (2.32) and (2.33). The distribution in (b) is a local pdf in the range of q(x)∈(-8.085,+∞) where (2.7) yields 1-Pr{q(x)≤|-8.085|}=0.0003.
(a)
(b)


## 2.2. Preliminaries on UnivariateL-Moments

LetY1,…,Yj,…,Yn be independent and identically distributed random variables each with continuous pdf fY(y), cdf FY(y), order statistics denoted as Y1:n≤⋯≤Yj:n≤⋯≤Yn: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 Yj:n are expressed as [4, pages 20–22]
(2.10)λ1=E[Y1:1]=β0,(2.11)λ2=12E[Y2:2-Y1:2]=2β1-β0,(2.12)λ3=13E[Y3:3-2Y2:3+Y1:3]=6β2-6β1+β0,(2.13)λ4=14E[Y4:4-3Y3:4+3Y2:4-Y1:4]=20β3-30β2+12β1-β0,
where the βi are determined from
(2.14)βi=∫y{FY(y)}ifY(y)dy,
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].TheL-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 cases. For example, the index of L-kurtosis (τ4) has the boundary condition for continuous distributions of [28]
(2.15)5τ32-14<τ4<1.


## 2.3.L-Moments for the System of γ-κ, γ, κ, and κℒ-κℛ Distributions

The derivation of the first fourL-moments associated with the γ-κ class of distributions begins by defining the probability weighted moments based on (2.14) in terms of (2.1)–(2.3) as
(2.16)βi=∫-∞+∞q(x,γ,κ){Ψ(x)}iψ(x)dx.In general, to obtain finite values ofλ1, λ2, τ3, τ4 based on (2.16), we must have (i) γ+κ<1, (ii) κ<1, and (iii) 1+γ>κ. As such, integrating (2.16) for i=0,1,2,3 and using (2.10)–(2.13) gives λ1, λ2, τ3, and τ4 as
(2.17)λ1=((-γ-κ)h1+(γ-κ)h2+(-γ+κ)h3+2κh4+(γ+κ)h5-2κh6)(2γ),(2.18)λ2=(2γ-(γ+κ)2h1+(γ-κ)2(h2-h3)+(γ+κ)2h5)(2γ),(2.19)τ3=(2γ2-(γ+κ)3h1+(γ-κ)3h2-(γ-κ)3h3+2κ3h4+(γ+κ)3h5-2κ3h6)2γλ2,(2.20)τ4=(2γ(1+5γ2+5κ(1-3κ))-(γ+κ)2(1+5(γ+κ)2)h1+(1+5(γ-κ)2)(γ-κ)2(h2-h3)+(γ+κ)2(1+5(γ+κ)2)h5/(12γλ2))
with harmonic number functions H[·] of h1=H[(1/2)(-1-γ-κ)], h2=H[(1/2)(-1+γ-κ)], h3=H[(1/2)(γ-κ)], h4=H[(1/2)(-1-κ)], h5=H[(-1/2)(γ+κ)], and h6=H[-(1/2)κ]. As such, given user specied values of τ3 and τ4, (2.19) and (2.20) can be numerically solved to obtain the parameters for γ and κ.Taking the limit of (2.17)–(2.20) as κ→0 gives the class associated with the log-logistic γ distributions in (2.4) as
(2.21)λ1=πcsc(γπ)-1γ,(2.22)λ2=γπcsc(γπ),(2.23)τ3=γ,(2.24)τ4=16(1+5γ2).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]
(2.25)λ1=0,(2.26)λ2=12(2+κ(-4p1+4p2+κ(p3-p4))),(2.27)τ3=0,(2.28)τ4=16(1+5κ2+10κ(1+2κ+2κ2(h6-h4))2λ2)
with polygamma functions P[·] of p1=P[0,1/2-κ/2],  p2=P[0,1-κ/2],  p3=P[1,1/2-κ/2], and p4=P[1,1-κ/2].TheL-moments for the class of asymmetric κℒ-κℛ distributions in (2.6) can be determined by separately evaluating and summing two integrals of the form in (2.16) as
(2.29)βi=∫-∞0q(x,κℒ){Ψ(x)}iψ(x)dx+∫0+∞q(x,κℛ){Ψ(x)}iψ(x)dx.For an asymmetricκℒ-κℛ distribution with a global pdf to have finite values of λ1, λ2, τ3, τ4 based on (2.29) will require κℒ<1 and κℛ<1. As such, integrating (2.29) to obtain β0,…,β3 and subsequently substituting these terms into (2.10)–(2.13) yield as in [27]
(2.30)λ1=14(2p5-2p6-2p7+2p8-κℒp9+κℒp10+κℛp11-κℛp12),(2.31)λ2=14(4+κℒ(-4p5+4p6+κℒ(p9-p10))+4+κℛ(-4p7+4p8+κℛ(p11-p12))),(2.32)τ3={4κℛ-4κℒ+κℒ2(6p5-6p6+κℒ(p10-p9))+κℛ2(-6p7+6p8+κℛ(p11-p12))}(4λ2),(2.33)τ4={2(2+5κℒ(1+3κℒ)+5κℛ(1+3κℛ))-4(κℒ+10κℒ3)p5-4(κℛ+10κℛ3)p7+5κℒ((4+40κℒ2)p6+κℒ(1+5κℒ2)(p9-p10)+κℛ((4+40κℛ2)p8+κℛ(1+5κℛ2)(p11-p12)))}/(24λ2)
with polygamma functions P[·] of p5=P[0,1/2-κℒ/2], p6=P[0,1-κℒ/2], p7=P[0,1/2-κℛ/2], p8=P[0,1-κℛ/2], p9=P[1,1/2-κℒ/2], p10=P[1,1-κℒ/2], p11=P[1,1/2-κℛ/2], p12=P[1,1-κℛ/2].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 κℒ and κℛ. Inspection of (2.32) and (2.33) reveals that interchanging values for the parameters of κℒ and κℛ 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 κℒ=0 with κℛ∈[0,1] in (2.32) and (2.33).Figure 5
Graph of the region for feasible combinations of (absolute value)L-skew |τ3| and L-kurtosis τ4 for global κℒ-κℛ pdfs. Asymmetric κℒ-κℛ distributions will lie in the area above the curve graphed in the |τ3| and τ4 plane.The conventional moment-based systems for theγ-κ and κℒ-κℛ 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 κℒ-κℛ class of distributions in Appendix B, the parameters are bounded in the range of 0≤κℒ<0.25 and 0≤κℛ<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 κℒ-κℛ distributions with specified L-correlations.


## 3.L-Correlations for the System of γ-κ, γ, κ, and κℒ-κℛ Distributions

The coefficient ofL-correlation [29] is introduced by considering two random variables Yj and Yk with distribution functions F(Yj) and F(Yk), respectively. The second L-moments of Yj and Yk can alternatively be expressed as
(3.1)λ2(Yj)=2Cov(Yj,F(Yj)),(3.2)λ2(Yk)=2Cov(Yk,F(Yk)).The secondL-comoments of Yj toward Yk and Yk toward Yj are
(3.3)λ2(Yj,Yk)=2Cov(Yj,F(Yk)),(3.4)λ2(Yk,Yj)=2Cov(Yk,F(Yj)).As such, theL-correlations of Yj toward Yk and Yk toward Yj are expressed as
(3.5)ηjk=λ2(Yj,Yk)λ2(Yj),(3.6)ηkj=λ2(Yk,Yj)λ2(Yk).TheL-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 theL-moment-based γ-κ, γ, κ, and κℒ-κℛ 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 Z1,…,ZT denote standard normal variables where the distribution functions and bivariate density function associated with Zj and Zk are expressed as
(3.7)Φ(zj)=Pr{Zj≤zj}=∫-∞zj(2π)-1/2exp{-wj22}dwj,(3.8)Φ(zk)=Pr{Zk≤zk}=∫-∞zk(2π)-1/2exp{-wk22}dwk,(3.9)fjk=(2π(1-ρjk2)1/2)-1exp{-(2(1-ρjk2))-1(zj2+zk2-2ρjkzjzk)}.Using (3.7), it follows that the jth distribution associated with (2.3)–(2.6) can be expressed as qj(g(Φ(Zj))), where g(Φ(Zj))=ln(Φ(Zj)/(1-Φ(Zj))) is standard logistic because Φ(Zj)~U(0,1). As such, using (3.5), the L-correlation of qj(g(Φ(Zj))) toward qk(g(Φ(Zk))) can be evaluated using the solved value(s) of the parameter(s) (i.e., γj-κj; γj; κj; κℒj-κℛj) for qj(g(Φ(Zj))), a specified intermediate correlation (IC) ρjk in (3.9), and the following integral generally expressed as
(3.10)ηjk=2π∬-∞+∞sj(qj(g(Φ(zj)),γj,κj,κℒj,κℛj))Φ(zk)fjkdzjdzk.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 qj(g(Φ(Zj))), which is induced by the parameters, such that qj(g(Φ(Zj))) has its specified L-correlation (ηjk) toward qk(g(Φ(Zk))). 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
(3.11)sj(qj(g(Φ(zj)),γj,κj,κℒj,κℛj))=δ×(qj(g(Φ(zj)),γj,κj,κℒj,κℛj)-λ1),
where λ1 is a mean from (2.17), (2.21), (2.25), or (2.30) and δ is a constant that scales λ2 in (2.18), (2.22), (2.26), or (2.31) and in the denominator of (3.5) to 1/π  as
(3.12)δ=δγj,κj=2γj(π(2γj-h1(γj+κj)2+(γj-κj)2(h2-h3)+h5(γj+κj)2)),δ=δγj=sin(γjπ)γjπ3/2,δ=δκj=2(π(2+κj(-4p1+4p2+κj(p3-p4)))),δ=δκℒj,κℛj=4/(π(4-4p5κℒj+4p6κℒj+p9κℒj2-p10κℒj2-4p7κℛj+4p8κℛj+p11κℛj2-p12κℛj2)).Analogously, theL-correlation of qk(g(Φ(Zk))) toward qj(g(Φ(Zj))) is expressed as
(3.13)ηkj=2π∬-∞+∞sk(qk(g(Φ(zk)),γk,κk,κℒk,κℛk))Φ(zj)fjkdzkdzj.Note also for the special case that ifqj(g(Φ(Zj))) in (3.10) and qk(g(Φ(Zk))) in (3.13) have the same parameters, that is, γj=γk; κj=κk; κℒj=κℒk; κℛj=κℛk, then ηjk=ηkj. Source code written in Mathematica [30] that implements the computation of an IC (ρjk) based on (3.10) is provided in Algorithm 1. The details for simulating γ-κ, γ, κ, and κℒ-κℛ distributions with specified values of L-skew, L-kurtosis, and L-correlations are described in the next section.Algorithm 1:Mathematica source code for computing intermediate correlations for specified L-correlations. The example is for distribution j=1 towards distribution k=2 (η12) in Figures 1(a) and 2(a). See Tables 1 and 3.
(*Intermediate Correlation for Distributions 1 and 2. See Table 3.*)
ρ
j
k
= 0.380206;
Needs[“MultivariateStatistics’”]
f
j
k = PDF[MultinormalDistribution[{0,0},{{1,ρjk},{ρjk,1}}],{Zj,Zk}];
F
j = CDF[NormalDistribution[0, 1],Zj];
F
k = CDF[NormalDistribution[0, 1],Zk];
(*Parameters for Figure 1(a)*)
γ
j
=
0.15;
κ
j
=
0.075;
X
j = Log[Fj/(1 -Fj)];
(*Quantile function from (2.3)*)
q
γ
j
,
κ
j
=
γ
j
-
1
*
(
Exp
[
γ
j
X
j
]
-
1
)
*
Exp
[
κ
j
*
Abs
[
X
j
]
];
(*Standardizing constants λ1 from (2.17) and δγj,κj from (3.12)*)
S
q
γ
j
,
κ
j
=
δ
γ
j
,
κ
j
*
(
q
γ
j
,
κ
j
-
λ
1
);
(*Compute the specified L-correlation in Table 1 for Distributions 1 and 2*)
η
j
k
=
2
π
*
NIntegrate
[
S
q
γ
j
,
κ
j
*
F
k
*
f
j
k
,
{
Z
j
,
-
10
,
10
}
,
{
Z
k
,
-
10
,
10
}
Method
→
MultiDimensional]
0.40


## 4. The Procedure for Simulation and Monte Carlo Study

To implement the method for simulatingγ-κ, γ, κ, and κℒ-κℛ distributions with specified L-moments and L-correlations we suggest the following six steps.(1)
Specify theL-moments for T transformations of the forms in (2.3)–(2.6) that is, q1(g(Φ(Z1))),…,qT(g(Φ(ZT))) and obtain the solutions for the parameters of γ-κ, γ, κ, or κℒ-κℛ 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 qj(g(Φ(Zj))) toward qk(g(Φ(Zk))), where j<k∈{1,2,…,T}.(2)
Compute the (Pearson) intermediate correlations (ICs)ρjk by substituting the solutions of the parameters from Step (1) into (3.10) and then numerically integrate to solve for ρjk (see Algorithm 1 for an example). Repeat this step separately for all T(T-1)/2 pairwise combinations of correlations.(3)
Assemble the ICs into aT×T matrix and decompose this matrix using a Cholesky factorization. Note that this step requires the IC matrix to be positive definite.(4)
Use the results of the Cholesky factorization from Step (3) to generateT standard normal variables (Z1,…,ZT) correlated at the intermediate levels as follows:
(4.1)Z1=a11V1,Z2=a12V1+a22V2,⋮Zj=a1jV1+a2jV2+⋯+aijVi+⋯+ajjVj,⋮ZT=a1TV1+a2TV2+⋯+aiTVi+⋯+ajTVj+⋯+aTTVT,
where V1,…,VT are independent standard normal random variables and where aij represents the element in the ith row and the jth column of the matrix associated with the Cholesky factorization performed in Step (3).(5)
SubstituteZ1,…,ZT from Step (4) into the following Taylor series-based expansion for the standard normal cdf [30]:
(4.2)Φ(Zj)=(12)+ϕ(Zj){Zj+Zj33+Zj5(3·5)+Zj7(3·5·7)+⋯},
where ϕ(Zj) 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,Φ(Zj), generated from Step (5) into the T equations of the form of qj(g(Φ(Zj))), where g(Φ(Zj))=ln(Φ(Zj)/(1-Φ(Zj))) is standard logistic to generate the γ-κ, γ, κ, and κℒ-κℛ distributions with the specified L-moments and L-correlations.To demonstrate the steps above and evaluate the proposed method, a comparison between the proposedL-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 Z1,…,Z4 with the specified ICs by making use of the formulae given in (4.1) of Step (4) with T=4. The values of Z1,…,Z4 are subsequently transformed to Φ(Z1),…,Φ(Z4) using (4.2) and then substituted into equations of the forms in (2.3)–(2.6) to produce q1(g(Φ(Z1))),…,q4(g(Φ(Z4))) for both methods.Table 1
Specified correlation matrix for the distributions in Figures1(a), 2(a), 3(b), and 4(a).
1
2
3
4
1
1
2
0.40
1
3
0.60
0.50
1
4
0.80
0.70
0.60
1Table 2
Intermediate correlations for the conventional moment procedure.
1
2
3
4
1
1
2
0.446226
1
3
0.671890
0.566802
1
4
0.835142
0.740484
0.653826
1Table 3
Intermediate correlations for theL-moment procedure.
1
2
3
4
1
1
2
0.380206
1
3
0.576862
0.486201
1
4
0.782049
0.686737
0.568733
1Table 4
Cholesky decompositions for the conventional moment procedure.
a
11
=
1
a
12
=
0.446226
a
13
=
0.671890
a
14
=
0.835142
0
a
22
=
0.894920
a
23
=
0.298336
a
24
=
0.411011
0
0
a
33
=
0.677908
a
34
=
-
0.044131
0
0
0
a
44
=
0.362849Table 5
Cholesky decompositions for theL-moment procedure.
a
11
=
1
a
12
=
0.380206
a
13
=
0.576862
a
14
=
0.782049
0
a
22
=
0.924902
a
23
=
0.288544
a
24
=
0.421015
0
0
a
33
=
0.764181
a
34
=
-
0.005081
0
0
0
a
44
=
0.459478Algorithm 2:Mathematica source code for computing intermediate correlations for specified conventional Pearson correlations. The example is for distributions j=1 and k=4 (ρ14*) in Figures 1(a) and 4(a). See Tables 1 and 2.
(
*Intermediate Correlation for Distributions 1 and 4. See Table 2.*)
ρ
j
k
= 0.835142;
Needs[“MultivariateStatistics’”]
f
j
k
=
PDF[MultinormalDistribution[{0,0},{{1,ρjk},{ρjk,1}}],{Zj,Zk}];
F
j
=
CDF[NormalDistribution[0,1],Zj];
F
k
=
CDF[NormalDistribution[0,1],Zk];
(*Parameters for Figure 1(a) and Figure 4(a) in Figure 4*)
γ
j
=
0.15;
κ
j
=
0.075;
κ
ℒ
k
=
0.0475;
κ
ℛ
k
=
0.1358;
X
j
=
Log[Fj/(1 -Fj)];
X
k
=
Log[Fk/(1-Fk)];
(*Quantile functions from (2.3) and (2.6)*)
q
γ
j
,
κ
j
=
γ
j
-
1
*
(
Exp
[
γ
j
X
j
]
-
1
)
*
Exp
[
κ
j
*
Abs
[
X
j
]
];
q
ℒ
k
=
X
k
*
Exp
[
κ
ℒ
k
*
Abs
[
X
k
]
];
q
ℛ
k
=
X
k
*
Exp
[
κ
ℛ
k
*
Abs
[
X
k
]
];
(*Standardizing constants λ1j, λ1k ((2.17), (2.30)) and α2j, α2k ((A.2), (B.2)) from
AppendicesA and B*)
S
q
γ
j
,
κ
j
=
(
q
γ
j
,
κ
j
-
λ
1
j
)
/
α
2
j;
S
q
ℒ
k
=
(
q
ℒ
k
-
λ
1
k
)
/
α
2
k;
S
q
ℛ
k
=
(
q
ℛ
k
-
λ
1
k
)
/
α
2
k;
(*Compute the specified Pearson correlation in Table 1 for Distributions 1 and 4.*)
ρ
j
k
* = NIntegrate[Sqγj,κj*Piecewise[{{Sqℒk, Xk≤0}, {Sqℛk, Xk>0}}]*fjk, {Zj, −10, 10}, {Zk,
−10, 10}, Method→MultiDimensional]
0.80In 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.Table 6
Skew (α3) and Kurtosis (α4) results for the conventional moment procedure. Sample size of n=25.
Dist.
Parameter
Estimate
95% bootstrap C.I.
Standard error
Relative bias %
1
α
3
=
3.488
0.9348
0.9229,
0.9495
0.00676
−73.2
α
4
=
120.0
2.861
2.8107,
2.9087
0.02499
−97.6
2
α
3
=
2.485
1.074
1.0641, 1.0853
0.00541
−56.8
α
4
=
26.57
2.240
2.1963,  2.2856
0.02272
−91.6
3
α
3
=
0.0
0.0125
−0.0052, 0.0294
0.00875
—
α
4
=
39.83
3.365
3.3157,
3.4122
0.02464
−91.6
4
α
3
=
1.83
0.5777
0.5633,
0.5907
0.00693
−68.4
α
4
=
23.48
2.5360
2.4915,
2.5778
0.02190
−89.2Table 7
Skew (α3) and Kurtosis (α4) results for the conventional moment procedure. Sample size of n=1000.
Dist.
Parameter
Estimate
95% bootstrap CI
Standard error
Relative bias %
1
α
3
=
3.488
2.562
2.5383, 2.5823
0.01117
−26.5
α
4
=
120.0
22.15
21.6873, 22.6698
0.24850
−81.5
2
α
3
=
2.485
2.180
2.1668,  2.1944
0.00697
−12.3
α
4
=
26.57
13.36
13.0936, 13.6467
0.14100
−49.7
3
α
3
=
0.0
−0.0051
−0.0265, 0.0163
0.01100
—
α
4
=
39.83
18.57
18.2203, 18.9412
0.18330
−53.4
4
α
3
=
1.83
1.54
1.5246, 1.5539
0.00743
−15.8
α
4
=
23.48
12.91
12.6537, 13.1903
0.13610
−45.0Table 8
L-skew (τ3) and L-kurtosis (τ4) results. Sample size of  n=25.
Dist.
Parameter
Estimate
95% bootstrap CI
Standard error
Relative bias %
1
τ
3
=
0.1664
0.1461
0.1443, 0.1480
0.00095
−12.2
τ
4
=
0.2506
0.2376
0.2362, 0.2389
0.00069
−5.19
2
τ
3
=
0.20
0.1833
0.1819, 0.1849
0.00079
−8.35
τ
4
=
0.20
0.1911
0.1898, 0.1923
0.00064
−3.35
3
τ
3
=
0.0
0.0014
−0.0008, 0.0037
0.00115
—
τ
4
=
0.2958
0.2809
0.2795, 0.2823
0.00070
−5.04
4
τ
3
=
0.10
0.0879
0.0859, 0.0897
0.00096
−12.1
τ
4
=
0.25
0.2400
0.2388, 0.2414
0.00065
−4.00Table 9
L-skew (τ3) and L-kurtosis (τ4) results. Sample size of n=1000.
Dist.
Parameter
Estimate
95% bootstrap CI
Standard error
Relative bias %
1
τ
3
=
0.1664
0.1657
0.1653, 0.1660
0.00018
−0.421
τ
4
=
0.2506
0.2502
0.2499, 0.2504
0.00012
−0.160
2
τ
3
=
0.20
0.1995
0.1992, 0.1997
0.00013
−0.250
τ
4
=
0.20
0.1997
0.1995, 0.1999
0.00011
−0.150
3
τ
3
=
0.0
−0.0001
−0.0005, 0.0003
0.00022
—
τ
4
=
0.2958
0.2954
0.2951, 0.2956
0.00012
−0.135
4
τ
3
=
0.10
0.0995
0.0992, 0.0999
0.00017
−0.500
τ
4
=
0.25
0.2497
0.2495, 0.2499
0.00011
−0.120Table 10
Correlation results for the conventional moment procedure. Sample size ofn=25.
Parameter
Estimate
95% bootstrap CI
Standard error
Relative bias %
ρ
12
*
=
0.40
0.4389
0.4365, 0.4413
0.00153
9.73
ρ
13
*
=
0.60
0.6522
0.6504, 0.6540
0.00153
8.70
ρ
14
*
=
0.80
0.8296
0.8287, 0.8308
0.00165
3.70
ρ
23
*
=
0.50
0.5441
0.5420, 0.5460
0.00142
8.82
ρ
24
*
=
0.70
0.7306
0.7293, 0.7321
0.00157
4.37
ρ
34
*
=
0.60
0.6417
0.6400, 0.6435
0.00154
6.95Table 11
Correlation results for the conventional moment procedure. Sample size ofn=1000.
Parameter
Estimate
95% bootstrap CI
Standard error
Relative bias %
ρ
12
*
=
0.40
0.4040
0.4035, 0.4044
0.00029
1.00
ρ
13
*
=
0.60
0.6067
0.6063, 0.6072
0.00036
1.12
ρ
14
*
=
0.80
0.8041
0.8038, 0.8044
0.00042
0.51
ρ
23
*
=
0.50
0.5037
0.5032, 0.5040
0.00029
0.74
ρ
24
*
=
0.70
0.7027
0.7024, 0.7030
0.00033
0.39
ρ
34
*
=
0.60
0.6041
0.6036, 0.6044
0.00032
0.68Table 12
Correlation results for theL-moment procedure. Sample size of n=25.
Parameter
Estimate
95% bootstrap CI
Standard error
Relative bias %
η
12
=
0.40
0.4074
0.4049, 0.4099
0.00154
1.85
η
13
=
0.60
0.6068
0.6049, 0.6087
0.00153
1.13
η
14
=
0.80
0.8058
0.8046, 0.8067
0.00152
0.73
η
23
=
0.50
0.5078
0.5056, 0.5100
0.00151
1.56
η
24
=
0.70
0.7076
0.7061, 0.7091
0.00154
1.09
η
34
=
0.60
0.6063
0.6045, 0.6083
0.00155
1.05Table 13
Correlation results for theL-moment procedure. Sample size of n=1000.
Parameter
Estimate
95% bootstrap CI
Standard error
Relative bias %
η
12
=
0.40
0.4000
0.3996, 0.4004
0.00023
—
η
13
=
0.60
0.6002
0.5999, 0.6005
0.00023
—
η
14
=
0.80
0.8001
0.7999, 0.8002
0.00023
—
η
23
=
0.50
0.4997
0.4994, 0.5001
0.00023
—
η
24
=
0.70
0.7002
0.6999, 0.7004
0.00023
—
η
34
=
0.60
0.5999
0.5996, 0.6002
0.00023
—


## 5. Discussion and Conclusion

One of the primary advantages thatL-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 κℒ-κℛ 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 RSE(α^3)=0.436 and RSE(α^4)=1.12 compared with RSE(τ^3)=0.109 and RSE(τ^4)=0.048. This demonstrates that L-skew and L-kurtosis have more precision because they have less variance around their estimates.The results associated with the conventional Pearson andL-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 proposedL-moment-based system of γ-κ, γ, κ, and κℒ-κℛ 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.



---

*Source: 245986-2012-10-22.xml*