The Beta Generalized Half-Normal Distribution: New Properties

Abstract

We study some mathematical properties of the beta generalized half-normal distribution recently proposed by Pescim et al. (2010). This model is quite flexible for analyzing positive real data since it contains as special models the half-normal, exponentiated half-normal, and generalized half-normal distributions. We provide a useful power series for the quantile function. Some new explicit expressions are derived for the mean deviations, Bonferroni and Lorenz curves, reliability, and entropy. We demonstrate that the density function of the beta generalized half-normal order statistics can be expressed as a mixture of generalized half-normal densities. We obtain two closed-form expressions for their moments and other statistical measures. The method of maximum likelihood is used to estimate the model parameters censored data. The beta generalized half-normal model is modified to cope with long-term survivors may be present in the data. The usefulness of this distribution is illustrated in the analysis of four real data sets.

1. Introduction

Cooray and Ananda [1] pioneered the generalized half-normal (GHN) distribution with shape parameter 𝛼>0 and scale parameter 𝜃>0 defined by the cumulative distribution function (cdf) 𝐺_𝛼,𝜃(𝑥)𝑥=2Φ𝜃^𝛼()−1=erf𝑥/𝜃^𝛼√2, (1) where the standard normal cdf Φ(𝑥) and the error function erf(𝑥) are given by Φ(𝑥)=12𝑥1+erf√2(𝑥)=2,erf√𝜋^𝑥₀𝑒^−𝑡2𝑑𝑡. (2)

Following an idea due to Eugene et al. [2], Pescim et al. [3] proposed the beta generalized half-normal (BGHN) distribution, which seems to be superior over the GHN model for some applications. The justification for the practicability of this model is based on the fatigue crack growth under variable stress or cyclic load. In this paper, we study several mathematical properties of the BGHN model with the hope that it will attract wider applications in reliability, engineering and in other areas of research.

The four-parameter BGHN cdf is defined from (1) by (for 𝑥>0), 𝐹(𝑥)=𝐼_{2Φ[(𝑥/𝜃)}^𝛼_]−1()=1𝑎,𝑏𝐵(𝑎,𝑏)^{2Φ[(𝑥/𝜃)𝛼}₀^]−1𝑤^𝑎−1(1−𝑤)^𝑏−1𝑑𝑤, (3) where 𝐵(𝑎,𝑏)=[Γ(𝑎)Γ(𝑏)]/Γ(𝑎+𝑏) is the beta function, 𝐼_𝑦(𝑎,𝑏)=𝐵(𝑎,𝑏)⁻¹∫^𝑦₀𝑤^𝑎−1(1−𝑤)^𝑏−1𝑑𝑤 is the incomplete beta function ratio, and 𝑎>0 and 𝑏>0 are two additional shape parameters.

The probability density function (pdf) and the hazard rate function (hrf) corresponding to (3) are 𝑓(𝑥)=2^𝑏−12/𝜋(𝛼/𝑥)(𝑥/𝜃)^𝛼𝑒^{−(1/2)(𝑥/𝜃)2𝛼}𝐵()×𝑥𝑎,𝑏2Φ𝜃^𝛼−1^𝑎−1𝑥1−Φ𝜃^𝛼𝑏−1,ℎ(𝑥)=2^𝑏−12𝜋𝛼𝑥𝑥𝜃^𝛼𝑒^{−(1/2)(𝑥/𝜃)2𝛼}×𝑥2Φ𝜃^𝛼−1^𝑎−1𝑥1−Φ𝜃^𝛼𝑏−1×𝐵(𝑎,𝑏)1−𝐼_{2Φ[(𝑥/𝜃)}^𝛼_]−1(𝑎,𝑏)⁻¹, (4)(5) respectively. Hereafter, a random variable with pdf (4) is denoted by 𝑋∼BGHN(𝛼,𝜃,𝑎,𝑏).

Pescim et al. [3] demonstrated that the cdf and pdf of 𝑋 can be expressed as infinite power series of the GHN cumulative distribution. Here, all expansions in power series are around the point zero. If 𝑏>0 is a real noninteger, we can expand the binomial term in (3) to obtain 𝐹(𝑥)=^∞_𝑗=0𝑤_𝑗𝑥2Φ𝜃^𝛼−1^𝑎+𝑗, (6) where 𝑤_𝑗=𝑤_𝑗(𝑎,𝑏)=(−1)^𝑗𝐵(𝑎,𝑏)⁻¹(𝑎+𝑗)⁻¹_𝑗^𝑏−1.

The pdf corresponding to (6) can be expressed as 𝑓(𝑥)=2𝜋𝛼𝑥𝑥𝜃^𝛼𝑒^{−(1/2)(𝑥/𝜃)2𝛼}×^∞_𝑗=0𝑤_𝑗𝑥2Φ𝜃^𝛼−1^{𝑎+𝑗−1}. (7) If 𝑎>0 is an integer, (7) provides the BGHN density function as an infinite power series of the GHN cumulative distribution. If 𝑏 is an integer, the index 𝑗 in the previous sum stops at 𝑏−1. Otherwise, if 𝑎 is a real noninteger, we can expand {2Φ[(𝑥/𝜃)^𝛼]−1}^{𝑎+𝑗−1} as follows: 𝑥2Φ𝜃^𝛼−1^{𝑎+𝑗−1}=^∞_𝑟=0𝑠_𝑟𝑥𝑎+𝑗−12Φ𝜃^𝛼−1^𝑟, (8) where 𝑠_𝑟𝑎+𝑗−1=^∞_𝑘=𝑟(−1)^𝑟+𝑘𝑘𝑘𝑟𝑎+𝑗𝑎+𝑗−1. (9) Hence, 𝑓(𝑥)=2𝜋𝛼𝑥𝑥𝜃^𝛼𝑒^{−(1/2)(𝑥/𝜃)2𝛼∞}_𝑟=0𝑏_𝑟()erf𝑥/𝜃^𝛼√2^𝑟, (10) whose coefficients are 𝑏_𝑟=𝑏_𝑟()=1𝑎,𝑏𝐵(𝑎,𝑏)^∞_𝑗=0(−1)^𝑗𝑗𝑠𝑏−1_𝑟𝑎+𝑗−1. (11)

The BGHN density function (4) allows for greater flexibility of its tails and can be widely applied in many areas of engineering and biology. We study the mathematical properties of this distribution because it extends some important distributions previously considered in the literature. In fact, the GHN model, with parameters 𝛼 and 𝜃, is clearly an important special case for 𝑎=𝑏=1, with a continuous crossover towards models with different shapes, for example, a particular combination of skewness and kurtosis. The BGHN distribution also contains the exponentiated generalized half-normal (EGHN) and half-normal (HN) distributions as submodels for 𝑏=1 and 𝑎=𝑏=𝛼=1, respectively. Moreover, while the transformation (3) is not analytically tractable in the general case, the formulas related to the BGHN distribution turn out manageable, and with the use of modern computer resources with analytic and numerical capabilities, may turn into adequate tools comprising the arsenal of applied statisticians.

The paper is outlined as follows. We derive an expansion for the quantile function in Section 2. Some statistical measures for the BGHN distribution such as moments, generating function, mean deviations, Rényi entropy and reliability are studied in Section 3. The computational issues relating to the infinite series for structural properties of the BGHN distribution are discussed in Section 4. In Section 5, we derive algebraic expressions for the moments, moment generating function (mgf), mean deviations and Rényi entropy for the order statistics. The estimation by maximum likelihood—including the case of censoring—is presented in Section 6. In Section 7, we propose a BGHN mixture model for survival data with long-term survivors. Section 8 illustrates the importance of the BGHN distribution applied to four real data sets. Finally, concluding remarks are given in Section 9.

2. Power Series for the Quantile Function

Power series methods are at the heart of many aspects of applied mathematics and statistics. Quantile functions are in widespread use in probability distributions and general statistics and often find representations in terms of power series. The quantile function for a probability distribution has many uses in both the theory and application of probability. It may be used to generate values of a random variable having 𝐹(𝑥) as its distributions function. This fact serves as the basis of a method for simulating a sample from an arbitrary distribution with the aid of a uniform random number generator.

The quantile function of 𝑋, say 𝑥=𝑄_BGHN(𝑢)=𝐹⁻¹(𝑢), can be obtained by inverting the cumulative function (3). Now, we provide a power series expansion for 𝑄_BGHN(𝑢) that can be useful to determine some mathematical measures of the BGHN distribution. First, an expansion for the inverse of the incomplete beta function 𝐼_𝑥(𝑎,𝑏)=𝑢 can be found in wolfram website (http://functions.wolfram.com/06.23.06.0004.01) as 𝑧=𝑄_𝐵(𝑢)=^∞_𝑖=0𝑑_𝑖𝑢^𝑖/𝑎, (12) where 𝑄_𝐵(𝑢) is the beta quantile function, 𝑑_𝑖=𝑒_𝑖[𝑎𝐵(𝑎,𝑏)]^1/𝑎 for 𝑎>0 and 𝑒₀=0, 𝑒₁=1, 𝑒₂=(𝑏−1)/(𝑎+1),…. The coefficients 𝑒^′_𝑖s for 𝑖≥2 can be obtained from a cubic recursion of the form 𝑒_𝑖=1𝑖²+()(×𝛼−2𝑖+1−𝛼)1−𝛿_𝑖,2^𝑖−1_𝑟=2𝑒_𝑟𝑒_{𝑖+1−𝑟}×[𝑟()(+1−𝑎)(𝑖−𝑟−𝑟𝑟−1)]_{𝑖−1𝑟=1𝑖−𝑟𝑠=1}𝑒_𝑟𝑒_𝑠𝑒_{𝑖+1−𝑟−𝑠}×[𝑟(𝑟−𝑎)+𝑠(𝑎+𝑏−2)(𝑖+1−𝑟−𝑠)], (13) where 𝛿_𝑖,2=1 if 𝑖=2 and 𝛿_𝑖,2=0 if 𝑖≠2. In the last equation, we note that the quadratic term only contributes for 𝑖≥3.

Following Steinbrecher [4], the quantile function of the standard normal distribution, say 𝑄_𝑁(𝑢)=Φ⁻¹(𝑢), can be expanded as 𝑄_𝑁(𝑢)=Φ⁻¹(𝑢)=^∞_𝑘=0𝑏_𝑘𝜁^2𝑘+1, (14) where 𝜁=√2𝜋(𝑢−1/2) and the 𝑏^′_𝑘s can be calculated from 𝑏₀=1 and 𝑏_𝑘+1=12(2𝑘+3)^𝑘_𝑟=0(2𝑟+1)(2𝑘−2𝑟+1)𝑏_𝑟𝑏_𝑘−𝑟().𝑟+1)(2𝑟+1 (15) Here, 𝑏₁=1/6, 𝑏₂=7/120, 𝑏₃=127/7560,….

The function 𝑄_𝑁(𝑢) can be expressed as a power series given by 𝑄_𝑁(𝑢)=^∞_𝑗=0𝑓_𝑗𝑢^𝑗, (16) where 𝑓_𝑗=∑^∞_𝑘=𝑗(−1/2)^{𝑘−𝑗𝑘}_𝑗𝑔_𝑘 and the quantities 𝑔_𝑘 are defined from the coefficients in (14) by 𝑔_𝑘=0 for 𝑘=0,2,4,… and 𝑔_𝑘=(2𝜋)^𝑘/2𝑏_(𝑘−1)/2 for 𝑘=1,3,5,….

By inverting the normal cumulative function in 𝑧=2Φ[(𝑥/𝜃)^𝛼]−1, and using (16), we can express the quantile function of 𝑋 in terms of 𝑧 as 𝑥=𝜃𝑄_𝑁𝑧+12^1/𝛼=𝜃^∞_𝑗=0𝑓_𝑗𝑧+12^𝑗1/𝛼=𝜃^∞𝑗_𝑗=0𝑝=0𝑓_𝑗2^𝑗𝑗𝑝𝑧^𝑝1/𝛼. (17) Replacing ∑^∞_𝑗=0∑^𝑗_𝑝=0 by ∑^∞_𝑝=0∑^∞_𝑗=𝑝 in last equation, we obtain 𝑥=𝜃^∞_𝑝=0ℎ_𝑝𝑧^𝑝1/𝛼, (18) where ℎ_𝑝=∑^∞_𝑗=𝑝2^−𝑗𝑓_𝑗^𝑗_𝑝. Using the same steps of (8), we can write ^∞_𝑝=0ℎ_𝑝𝑧^𝑝1/𝛼=^∞_𝑟=0𝑠_𝑟𝛼^−1∞_𝑝=0ℎ_𝑝𝑧^𝑝𝑟, (19) where 𝑠_𝑟(𝛼⁻¹)=∑^∞_𝑗=𝑟(−1)^{𝑟+𝑗𝛼−1}_𝑗^𝑗_𝑟.

We use throughout an equation of Gradshteyn and Ryzhik [5] for a power series raised to a positive integer 𝑗^∞_𝑖=0𝑎_𝑖𝑥^𝑖𝑗=^∞_𝑖=0𝑐_𝑗,𝑖𝑥^𝑖, (20) where the coefficients 𝑐_𝑗,𝑖 (for 𝑖=1,2,…) are easily obtained from the recurrence equation 𝑐_𝑗,𝑖=𝑖𝑎₀^𝑖−1_𝑚=1𝑚𝑗+1−𝑖𝑎_𝑚𝑐_{𝑗,𝑖−𝑚} (21) and 𝑐_𝑗,0=𝑎^𝑗₀. From (19) and (20), we have 𝑥=^∞_𝑝=0ℎ^⋆_𝑝𝑧^𝑝, (22) where ℎ^⋆_𝑝=𝜃∑^∞_𝑟=0𝑠_𝑟(𝛼⁻¹)ℎ_𝑟,𝑝 for 𝑝≥0, ℎ_𝑟,𝑝=(𝑝ℎ₀)⁻¹∑^𝑝_𝑚=1[𝑚(𝑟+1)−𝑝]ℎ_𝑚ℎ_{𝑟,𝑝−𝑚} for 𝑝≥1 and ℎ_𝑟,0=ℎ^𝑟₀. By inserting (12) in (22), we obtain 𝑥=𝑄_BGHN(𝑢)=^∞_𝑝=0ℎ^⋆_𝑝^∞_𝑖=0𝑑_𝑖𝑢^{𝑖/𝑎𝑝}. (23) From (20), it follows (∑^∞_𝑖=0𝑑_𝑖𝑢^𝑖/𝑎)^𝑝=∑^∞_𝑖=0𝑑_𝑝,𝑖𝑢^𝑖/𝑎, where the quantities 𝑑_𝑝,𝑖 are determined recursively from 𝑑_0,𝑖=𝑑^𝑖₀ and 𝑑_𝑝,𝑖=(𝑖𝑑₀)⁻¹∑^𝑖_𝑚=1[𝑚(𝑝+1)−𝑖]𝑑_𝑚𝑑_{𝑝,𝑖−𝑚} (for 𝑖=1,2,…). Finally, we obtain 𝑥=𝑄_BGHN(𝑢)=^∞_𝑖=0𝑣_𝑖𝑢^𝑖/𝑎, (24) where 𝑣_𝑖=∑^∞_𝑝=0ℎ^⋆_𝑝𝑑_𝑝,𝑖. Equation (24) gives the BGHN quantile function as a power series and represents the main result of this section.

3. BGHN Properties

3.1. Moments and Generating Function

Here, we provide new expressions for the moments and mgf of 𝑋 based upon the power series for its quantile function as alternative results form those obtained by Pescim et al. [3].

We can write from (24) 𝜇^′_𝑠(𝛼,𝜃,𝑎,𝑏)=𝐸𝑋^𝑠=^∞₀𝑥^𝑠𝑓(𝑥)=𝑑𝑥_10∞𝑖=0𝑣_𝑖𝑢^{𝑖/𝑎𝑠}𝑑𝑢, (25) and then using (20), we obtain 𝜇^′_𝑠(𝛼,𝜃,𝑎,𝑏)=𝑎^∞_𝑖=0𝑣_𝑠,𝑖(),𝑖+𝑎 (26) where the quantities 𝑣_𝑠,𝑖 (for 𝑖=1,2,…) are easily determined from the recurrence equation 𝑣_𝑠,𝑖=𝑖𝑣₀^𝑖−1_𝑚=1[𝑚(𝑠+1)−𝑖]𝑣_𝑚𝑣_{𝑠,𝑖−𝑚}, (27) with 𝑣_𝑠,0=𝑣^𝑠₀.

We now provide a new alternative representation for the mgf of 𝑋, say 𝑀_{𝛼,𝜃,𝑎,𝑏}(𝑡)=𝐸(𝑒^𝑡𝑋), based on the quantile power series (24). We can write 𝑀_{𝛼,𝜃,𝑎,𝑏}(𝑡)=^∞₀𝑒^𝑡𝑥𝑓(𝑥)=𝑑𝑥₁₀exp𝑡^∞_𝑖=0𝑣_𝑖𝑢^𝑖/𝑎𝑑𝑢. (28) We expand the exponential function and use the same algebra that leads to (26) 𝑀_{𝛼,𝜃,𝑎,𝑏}(𝑡)=^∞_𝑟=0𝑡^𝑟𝑟!¹₀^∞_𝑖=0𝑣_𝑖𝑢^{𝑖/𝑎𝑟}=𝑑𝑢_{∞𝑟,𝑖=0}𝑡^𝑟𝑣_𝑟,𝑖𝑟!¹₀𝑢^𝑖/𝑎𝑑𝑢 (29) and then 𝑀_{𝛼,𝜃,𝑎,𝑏}(𝑡)=𝑎^∞_𝑟,𝑖=0𝑡^𝑟𝑣_𝑟,𝑖().𝑎+𝑖𝑟! (30) Equations (26) and (30) are the main results of this section.

3.2. Mean and Median Deviations

The amount of scatter in a population is evidently measured to some extent by the mean deviations in relation to the mean and the median defined by 𝛿₁(𝑋)=^∞₀𝑥−𝜇^′₁𝑓(𝑥)𝛿𝑑𝑥,₂(𝑋)=^∞₀|𝑥−𝑀|𝑓(𝑥)𝑑𝑥, (31) respectively, where 𝜇^′₁=𝐸(𝑋) and 𝑀=Median(𝑋) denotes the median. Here, 𝑀 is calculated as the solution of the nonlinear equation 𝐼_{2Φ[(𝑀/𝜃)}^𝛼_]−1(𝑎,𝑏)=1/2. We define 𝑇(𝑞)=∫^∞_𝑞𝑥𝑓(𝑥)𝑑𝑥, which is determined below. The measures 𝛿₁(𝑋) and 𝛿₂(𝑋) can be written in terms of 𝜇^′₁ and 𝑇(𝑞) as 𝛿₁(𝑋)=2𝜇^′₁𝐹𝜇^′₁−2𝜇^′₁+2𝑇𝜇^′₁,𝛿₂(𝑋)=2𝑇(𝑀)−𝜇^′₁. (32) For more details, see Paranaíba et al. [6]. Clearly, 𝐹(𝑀) and 𝐹(𝜇^′₁) are determined from (3). From (10), we have 𝑇𝑞=𝛼2𝜋^∞_𝑟=0𝑏_𝑟^∞_𝑞𝑥𝜃^𝛼𝑒^{−(1/2)(𝑥/𝜃)2𝛼}()erf𝑥/𝜃^𝛼√2^𝑟𝑑𝑥. (33) Setting 𝑢=(𝑥/𝜃)^𝛼 in the last equation gives 𝑇𝑞=𝜃2𝜋^∞_𝑟=0𝑏_𝑟^∞_{𝑞/𝜃𝛼}𝑢^1/𝛼𝑒^−𝑢2/2𝑢erf√2^𝑟𝑑𝑢. (34) Using the power series for the error function erf(𝑥)=(2/√𝜋)∑^∞_𝑚=0((−1)^𝑚𝑥^2𝑚+1/(2𝑚+1)𝑚!) (see, e.g., [7]), we obtain after some algebra 𝑇𝑞=𝜃2𝜋^∞_𝑟=0𝑏_𝑟2√𝜋^𝑟×^∞_𝑚1=0⋯^∞_𝑚𝑟=0(−1)^{𝑚1+⋯+𝑚𝑟}2𝑚₁+1⋯2𝑚_𝑟+1𝑚₁!⋯𝑚_𝑟!×Γ𝑚₁+⋯+𝑚_𝑟+𝑟2+1+12𝛼2,12𝑞𝜃^𝛼, (35) where Γ(𝑝,𝑥)=∫^∞_𝑥𝑣^𝑝−1𝑒^−𝑣𝑑𝑣 denotes the complementary incomplete gamma function for 𝑝>0. The measures 𝛿₁(𝑋) and 𝛿₂(𝑋) are immediately calculated from (35).

Bonferroni and Lorenz curves have applications not only in economics to study income and poverty, but also in other fields such as reliability, demography, insurance and medicine. They are defined by 𝐵𝜌=1𝜌𝜇^′₁^𝑞0(𝑥)𝐿𝜌=1𝑥𝑓𝑑𝑥,𝜇^′₁^𝑞0𝑥𝑓(𝑥)𝑑𝑥, (36) respectively, where 𝑞=𝑄_BGHN(𝜌)=𝑄_𝐵(𝜈) and 𝜈=2Φ[(𝜌/𝜃)^𝛼]−1 (Section 2) for a given probability 𝜌. From ∫^𝑞₀𝑥𝑓(𝑥)𝑑𝑥=𝜇^′₁−𝑇(𝜌), we obtain 𝐵(𝜌)=𝜌⁻¹[1−𝑇(𝜌)/𝜇^′₁] and 𝐿(𝜌)=1−𝑇(𝜌)/𝜇^′₁.

3.3. Rényi Entropy

The Rényi information of order 𝜉 for a continuous random variable with density function 𝑓(𝑥) is defined as 𝚥_𝑅(𝜉)=1𝐼(𝜉,1−𝜉log) (37) where 𝐼(𝜉)=∫𝑓^𝜉(𝑥)𝑑𝑥, 𝜉>0 and 𝜉≠1. Applications of the Rényi entropy can be found in several areas such as physics, information theory and engineering to describe many nonlinear dynamical or chaotic systems [8], and in statistics as certain appropriately scaled test statistics (relative Rényi information) for testing hypotheses in parametric models [9]. Rényi [10] generalized the concept of information theory which allows for different averaging of probabilities via 𝜉.

For the BGHN distribution (4), the Rényi entropy is defined by 𝚥_𝑅(𝜉)=1𝐼(𝜉,1−𝜉log) (38) where 𝐼(𝜉)=∫𝑓^𝜉(𝑥)𝑑𝑥, 𝜉>0 and 𝜉≠1. From (4), we have 𝐼(𝜉)=2^{𝜉(𝑏−1)}𝛼2/𝜋^𝜉[𝐵(𝑎,𝑏)]^𝜉×^∞₀𝑥^−𝜉𝑥𝜃^𝛼𝜉𝑒^{−(𝜉/2)(𝑥/𝜃)2𝛼}𝑥2Φ𝜃^𝛼−1^{𝜉(𝑎−1)}×𝑥1−Φ𝜃^{𝛼𝜉(𝑏−1)}𝑑𝑥. (39) For |𝑧|<1 and 𝑏 is a real noninteger, the power series holds (1−𝑧)^𝑏−1=^∞_𝑗=0(−1)^𝑗𝑗𝑧𝑏−1_𝑗, (40) where the binomial coefficient is defined for any real. Using (40) in (39) twice, 𝐼(𝜉) can be expressed as 𝐼(𝜉)=2^{𝜉(𝑏−1)}𝛼2/𝜋^𝜉[𝐵(𝑎,𝑏)]^𝜉×^∞_𝑗,𝑘=0(−1)^𝑗+𝑘2^𝑗𝜉()𝑗×𝜉()𝑘𝑎−1+1𝑏−1+𝑗+1×^∞₀𝑥^−𝜉𝑥𝜃^𝛼𝜉𝑒^{−(𝜉/2)(𝑥/𝜃)2𝛼}Φ𝑥𝜃^𝛼𝑘𝑑𝑥. (41) Substituting Φ(𝑥) by the error function and setting 𝑢=(𝑥/𝜃)^𝛼, 𝐼(𝜉) reduces to 𝐼(𝜉)=𝛼^𝜉−1𝜃^1−𝜉2^{𝜉(𝑏−1)}×2𝜋^𝜉∞𝑘_{𝑗,𝑘=0𝑙=0}𝑤_{𝑗,𝑘,𝑙}()×𝑎,𝑏_∞0𝑢^{((𝜉(𝛼−1)+1)/𝛼)−1}𝑒^{−(𝜉/2)𝑢2}𝑢erf√2^𝑙𝑑𝑢. (42) Following similar algebra that lead to (35), we obtain 𝐼(𝜉)=𝛼^𝜉−1𝜃^1−𝜉2^{[𝜉(𝑏−1)+(𝜉(𝛼−1)+1)/2𝛼−1]}×2𝜋^𝜉∞𝑘_{𝑗,𝑘=0𝑙=0}𝑤_{𝑗,𝑘,𝑙}()×𝑎,𝑏_∞𝑚1=0⋯^∞_𝑚𝑙=0(−1)^{𝑚1+⋯+𝑚𝑙}2𝑚₁+1⋯2𝑚_𝑙+1𝑚₁!⋯𝑚_𝑙!×𝜉^{−[𝑚1+⋯+𝑚𝑙+𝑙/2+(𝜉(𝛼−1)+1)/2𝛼]}×Γ𝑚₁+⋯+𝑚_𝑙+𝑙2+𝜉(𝛼−1)+1.2𝛼 (43) Finally, the Rényi entropy reduces to 𝚥_𝑅(𝜉)=(1−𝜉)⁻¹×)(𝛼)+()(𝜃)+𝜉()+𝜉()(𝜉−1log1−𝜉log𝑏−1𝛼−1+1(2)2𝛼−1log+𝜉log2𝜋+log^∞𝑘_{𝑗,𝑘=0𝑙=0}𝑤_{𝑗,𝑘,𝑙}(𝑎,𝑏)+log^∞_𝑚1=0⋯^∞_𝑚𝑙=0(−1)^{𝑚1+⋯+𝑚𝑙}×2𝑚₁+1⋯2𝑚_𝑙+1×𝑚₁!⋯𝑚_𝑙!⁻¹−𝑚₁+⋯+𝑚_𝑙+𝑙2+𝜉(𝛼−1)+1(𝜉)Γ𝑚2𝛼log+log₁+⋯+𝑚_𝑙+𝑙2+𝜉(𝛼−1)+1.2𝛼 (44)

3.4. Reliability

In the context of reliability, the stress-strength model describes the life of a component which has a random strength 𝑋₁ that is subjected to a random stress 𝑋₂. The component fails at the instant that the stress applied to it exceeds the strength, and the component will function satisfactorily whenever 𝑋₁>𝑋₂. Hence, 𝑅=Pr(𝑋₂<𝑋₁) is a measure of component reliability. Here, we derive 𝑅 when 𝑋₁ and 𝑋₂ have independent BGHN(𝛼,𝜃,𝑎₁,𝑏₁) and BGHN(𝛼,𝜃,𝑎₂,𝑏₂) distributions with the same shape parameters 𝛼 and 𝜃. The reliability 𝑅 becomes 𝑅=^∞₀𝑓₁(𝑥)𝐹₂(𝑥)𝑑𝑥, (45) where the cdf of 𝑋₂ and the density of 𝑋₁ are obtained from (6) and (10) as 𝐹₂(𝑥)=^∞_𝑗=0𝑤_𝑗𝑎₂,𝑏₂𝑥2Φ𝜃^𝛼−1^𝑎2+𝑗,𝑓₁(𝑥)=2𝜋𝛼𝑥𝑥𝜃^𝛼𝑒^{−(1/2)(𝑥/𝜃)2𝛼}×^∞_𝑟=0𝑏_𝑟𝑎₁,𝑏₁𝑥2Φ𝜃^𝛼−1^𝑟, (46) respectively, where 𝑤_𝑗𝑎₂,𝑏₂=(−1)^𝑗𝐵𝑎₂,𝑏₂𝑏₂𝑗,𝑏−1_𝑟𝑎₁,𝑏₁=^∞_𝑗=0(−1)^𝑗𝐵𝑎₁,𝑏₁𝑏₁𝑗𝑠−1_𝑟𝑎₁+𝑗−1 (47) refer to 𝑋₂ and 𝑋₁, respectively. Hence, 𝑅=𝛼2𝜋^∞_𝑗,𝑟=0𝑤_𝑗𝑎₂,𝑏₂𝑏_𝑟𝑎₁,𝑏₁×^∞₀𝑥⁻¹𝑥𝜃^𝛼𝑒^{−(1/2)(𝑥/𝜃)2𝛼}𝑥2Φ𝜃^𝛼−1^{𝑎2+𝑗+𝑟}𝑑𝑥. (48)

Setting 𝑢=2Φ[(𝑥/𝜃)^𝛼]−1 in the last integral, the reliability of 𝑋 reduces to 𝑅=^∞_𝑗,𝑟=0𝑤_𝑗𝑎₂,𝑏₂𝑏_𝑟𝑎₁,𝑏₁𝑎₂.+𝑗+𝑟+1 (49)

4. Computational Issues

Here, we show the practical use of (24), (26), (30), (35), and (49). If we truncate these infinite power series, we have a simple way to compute the quantile function, moments, mgf, mean deviations and reliability of the BGHN distribution. The question is how large the truncation limit should be?

We now provide evidence that each infinite summation can be truncated at twenty to yield sufficient accuracy. Let 𝐷₁ denote the absolute difference between the integrated version, 1𝐵(𝑎,𝑏)^{2Φ[(𝑥/𝜃)𝛼}₀^]−1𝑡^𝑎−1(1−𝑡)^𝑏−1𝑑𝑡=𝑈, (50) where 𝑈∼𝑈(0,1) is a uniform variate on the unit interval, and the truncated version of (24) averaged over (all with step 0.01) 𝑥=0.01,…,5, 𝑎=0.01,…,10, 𝑏=0.01,…,10, 𝛼=0.01,…,10 and 𝜃=0.01,…,10. Let 𝐷₂ denote the absolute difference between integrated version 𝜇^′_𝑠=2^𝑏−1𝐵(𝑎,𝑏)2𝜋^∞₀𝑥^𝑠𝛼𝑥𝑥𝜃^𝛼𝑒^{−(1/2)(𝑥/𝜃)2𝛼}×𝑥2Φ𝜃^𝛼−1^𝑎−1×𝑥1−Φ𝜃^𝛼𝑏−1𝑑𝑥 (51) and the truncated version of (26) averaged over 𝑥=0.01,…,5, 𝑎=0.01,…,10, 𝑏=0.01,…,10, 𝛼=0.01,…,10 and 𝜃=0.01,…,10. Let 𝐷₃ denote the absolute difference between the truncated version of (30) and the integrated version, 𝑀_{𝛼,𝜃,𝑎,𝑏}(𝑡)=2^𝑏−1𝐵(𝑎,𝑏)2𝜋^∞₀𝛼𝑥𝑥𝜃^𝛼1exp𝑡𝑥−2𝑥𝜃^2𝛼×𝑥2Φ𝜃^𝛼−1^𝑎−1×𝑥1−Φ𝜃^𝛼𝑏−1𝑑𝑥, (52) averaged over 𝑥=0.01,…,5, 𝑎=0.01,…,10, 𝑏=0.01,…,10, 𝛼=0.01,…,10 and 𝜃=0.01,…,10. Let 𝐷₄ denote the absolute difference between the truncated version of (35) and the integrated version, 𝑇𝑞=𝛼2^𝑏−1𝐵(𝑎,𝑏)2𝜋^𝑞₀𝑥𝜃^𝛼−1×exp2𝑥𝜃^2𝛼𝑥2Φ𝜃^𝛼−1^𝑎−1×𝑥1−Φ𝜃^𝛼𝑏−1𝑑𝑥, (53) averaged over 𝑥=0.01,…,5, 𝑎=0.01,…,10, 𝑏=0.01,…,10, 𝛼=0.01,…,10 and 𝜃=0.01,…,10. Let 𝐷₅ denote the absolute difference between the truncated version of (49) and the integrated version (45), averaged over 𝑥=0.01,…,5, 𝑎=0.01,…,10, 𝑏=0.01,…,10, 𝛼=0.01,…,10 and 𝜃=0.01,…,10.

We obtain the following estimates after extensive computations: 𝐷₁=2.31×10⁻²⁰, 𝐷₂=8.47×10⁻¹⁸, 𝐷₃=1.22×10⁻²¹, 𝐷₄=1.51×10⁻²² and 𝐷₅=9.41×10⁻¹⁹. These estimates are small enough to suggest that the truncated versions of (24), (26), (30), (35) and (49) are reasonable practical use.

It would be important to verify that each (untruncated) infinite series (such as (24), (26), (30), (35) and (49)) is convergent and provide valid values for its arguments. However, this will be a difficult mathematical problem that could be investigated in future work.

5. Properties of the BGHN Order Statistics

Order statistics have been used in a wide range of problems, including robust statistical estimation and detection of outliers, characterization of probability distributions and goodness-of-fit tests, entropy estimation, analysis of censored samples, reliability analysis, quality control, and strength of materials.

5.1. Mixture Form

Suppose 𝑋₁,…,𝑋_𝑛 is a random sample of size 𝑛 from a continuous distribution and let 𝑋_1∶𝑛<𝑋_2∶𝑛<⋯<𝑋_𝑛∶𝑛 denote the corresponding order statistics. There has been a large amount of work relating to moments of order statistics 𝑋_𝑖∶𝑛. See Arnold et al. [11], David and Nagaraja [12], and Ahsanullah and Nevzorov [13] for excellent accounts. It is well known that the density function of 𝑋_𝑖∶𝑛 is given by 𝑓_𝑖∶𝑛(𝑥)=𝑓(𝑥)𝐵()𝐹(𝑥)𝑖,𝑛−𝑖+1^𝑖−1[1−𝐹(𝑥)]^𝑛−𝑖. (54)

For the BGHN distribution, Pescim et al. [3] obtained 𝑓_𝑖∶𝑛(𝑥)=^{𝑛−𝑖∞}_{𝑘=0𝑚1=0}⋯^∞_{𝑚𝑖+𝑘−1=0}𝜂_{𝑘,𝑀𝑘}𝑓_{𝑘,𝑀𝑘}(𝑥), (55) where 𝑀_𝑘 denotes a sequence (𝑚₁,…,𝑚_{𝑖+𝑘−1}) of 𝑖+𝑘−1 nonnegative integers, 𝑓_{𝑘,𝑀𝑘}(𝑥) denotes the BGHN(𝛼,𝜃,(𝑖+𝑘)𝑎+∑^{𝑖+𝑘−1}_𝑗=1𝑚_𝑗,𝑏) density function defined under 𝑀_𝑘 and the constants 𝜂_{𝑘,𝑀𝑘} are given by 𝜂_{𝑘,𝑀𝑘}=(−1)^{𝑘+∑𝑖+𝑘−1𝑗=1𝑚𝑗}_𝑘^𝑛−𝑖𝐵𝑎(𝑖+𝑘)+∑^{𝑖+𝑘−1}_𝑗=1𝑚_𝑗,𝑏Γ(𝑏)^{𝑖+𝑘−1}𝐵(𝑎,𝑏)^𝑖+𝑘𝐵(𝑖,𝑛−𝑖+1)∏^{𝑖+𝑘−1}_𝑗=1Γ𝑏−𝑚_𝑗𝑚_𝑗!𝑎+𝑚_𝑗. (56) The quantities 𝜂_{𝑘,𝑀𝑘} are easily obtained given 𝑘 and a sequence 𝑀_𝑘 of indices 𝑚₁,…,𝑚_{𝑖+𝑘−1}. The sums in (55) extend over all (𝑖+𝑘)-tuples (𝑘,𝑚₁,…,𝑚_{𝑖+𝑘−1}) and can be implementable on a computer. If 𝑏>0 is an integer, (55) holds but the indices 𝑚₁,…,𝑚_{𝑖+𝑘−1} vary from zero to 𝑏−1. Equation (55) reveals that the density function of the BGHN order statistics can be expressed as an infinite mixture of BGHN densities.

5.2. Moments of Order Statistics

The moments of the BGHN order statistics can be written directly in terms of the moments of BGHN distributions from the mixture form (55). We have 𝐸𝑋^𝑠_𝑖∶𝑛=^{𝑛−𝑖∞}_{𝑘=0𝑚1=0}⋯^∞_{𝑚𝑖+𝑘−1=0}𝜂_{𝑘,𝑀𝑘}𝐸𝑋^𝑠_𝑖,𝑘, (57) where 𝑋_𝑖,𝑘∼BGHN(𝛼,𝜃,(𝑖+𝑘)𝑎+∑^{𝑖+𝑘−1}_𝑗=1𝑚_𝑗,𝑏) and 𝐸(𝑋^𝑠_𝑖,𝑘) can be determined from (26). Inserting (26) in (57) and changing indices, we can write 𝐸𝑋^𝑠_𝑖∶𝑛=^{𝑛−𝑖∞}_{𝑘=0𝑚1=0}⋯^∞_{𝑚𝑖+𝑘−1=0}𝜂_{𝑘,𝑀𝑘}^∞_𝑞=0𝑎^⋆𝑣_𝑠,𝑞𝑞+𝑎^⋆, (58) where 𝑎^⋆=(𝑖+𝑘)𝑎+^{𝑖+𝑘−1}_𝑗=1𝑚_𝑗. (59)

The moments 𝐸(𝑋^𝑠_𝑖∶𝑛) can be determined based on the explicit sum (58) using the software Mathematica. They are also computed from (55) by numerical integration using the statistical software R [14]. The values from both techniques are usually close when ∞ is replaced by a moderate number such as 20 in (58). For some values 𝑎=2, 𝑏=3 and 𝑛=5, Table 1 lists the numerical values for the moments 𝐸(𝑋^𝑠_𝑖∶𝑛), (𝑠=1,…,4; 𝑖=1,…,5) and for the variance, skewness and kurtosis.

Table 1 

Moments of the BGHN order statistics for 𝜃=8.5, 𝛼=3, 𝑎=2, and 𝑏=3.

5.3. Generating Function

The mgf of the BGHN order statistics can be written directly in terms of the BGHN mgf’s from the mixture form (55) and (30). We obtain 𝑀_𝑖∶𝑛(𝑡)=^{𝑛−𝑖∞}_{𝑘=0𝑚1=0}⋯^∞_{𝑚𝑖+𝑘−1=0}𝜂_{𝑘,𝑀𝑘}𝑀_{𝛼,𝜃,𝑎⋆,𝑏}(𝑡), (60) where 𝑀_{𝛼,𝜃,𝑎⋆,𝑏}(𝑡) is the mgf of the BGHN(𝛼,𝜃,𝑎^⋆,𝑏) distribution obtained from (30).

5.4. Mean Deviations

The mean deviations of the order statistics in relation to the mean and the median are defined by 𝛿₁𝑋_𝑖∶𝑛=^∞₀𝑥−𝜇_𝑖𝑓_𝑖∶𝑛(𝑥)𝛿𝑑𝑥,₂𝑋_𝑖∶𝑛=^∞₀𝑥−𝑀_𝑖𝑓_𝑖∶𝑛(𝑥)𝑑𝑥, (61) respectively, where 𝜇_𝑖=𝐸(𝑋_𝑖∶𝑛) and 𝑀_𝑖=Median(𝑋_𝑖∶𝑛) denotes the median. Here, 𝑀_𝑖 is obtained as the solution of the nonlinear equation ^𝑛_𝑟=𝑖𝑛𝑟𝐼_{2Φ[(𝑀𝑖/𝜃)}^𝛼_]−1(𝑎,𝑏)^𝑟×1−𝐼_{2Φ[(𝑀𝑖/𝜃)}^𝛼_]−1(𝑎,𝑏)^𝑛−𝑟=12. (62) The measures 𝛿₁(𝑋_𝑖∶𝑛) and 𝛿₂(𝑋_𝑖∶𝑛) follow from 𝛿₁𝑋_𝑖∶𝑛=2𝜇_𝑖𝐹_𝑖∶𝑛𝜇_𝑖−2𝜇_𝑖+2𝐽_𝑖𝜇_𝑖,𝛿₂𝑋_𝑖∶𝑛=2𝐽_𝑖𝑀_𝑖−𝜇_𝑖, (63) where 𝐽_𝑖(𝑞)=∫^∞_𝑞𝑥𝑓_𝑖∶𝑛(𝑥)𝑑𝑥. Using (55), we have 𝐽_𝑖𝑞=^{𝑛−𝑖∞}_{𝑘=0𝑚1=0}⋯^∞_{𝑚𝑖+𝑘−1=0}𝜂_{𝑖,𝑀𝑘}𝑇𝑞, (64) where 𝑇(𝑞) is obtained from (35) using 𝑎^⋆ given by (59) and 𝑏^⋆_𝑟=𝑏_𝑟𝑎^⋆=1,𝑏𝐵(𝑎^⋆,𝑏)^∞_𝑗=0(−1)^𝑗𝑗𝑠𝑏−1_𝑟𝑎^⋆+𝑗−1, (65) where 𝑠_𝑟𝑎^⋆+𝑗−1=^∞_𝑘=𝑟(−1)^𝑟+𝑘𝑎^⋆𝑎+𝑗^⋆𝑘𝑘𝑟+𝑗−1. (66)

Bonferroni and Lorenz curves of the order statistics are given by 𝐵_𝑖∶𝑛𝜌=1𝜌𝜇_𝑖^𝑞₀𝑥𝑓_𝑖∶𝑛(𝑥)𝐿𝑑𝑥,_𝑖∶𝑛𝜌=1𝜇_𝑖^𝑞₀𝑥𝑓_𝑖∶𝑛(𝑥)𝑑𝑥, (67) respectively, where 𝑞=𝐹⁻¹_𝑖∶𝑛(𝜌) for a given probability 𝜌. From ∫^𝑞₀𝑥𝑓_𝑖∶𝑛(𝑥)𝑑𝑥=𝜇_𝑖−𝐽_𝑖(𝜌), we obtain 𝐵_𝑖∶𝑛𝜌=1𝜌−𝐽_𝑖𝜌𝜌𝜇_𝑖,𝐿_𝑖∶𝑛𝜌𝐽=1−^𝑖𝜌𝜇_𝑖. (68)

5.5. Rényi Entropy

The Rényi entropy of the order statistics is defined by 𝚥_𝑅(𝜆)=1[𝐻(𝜆,1−𝜆log)] (69) where 𝐻(𝜆)=∫𝑓^𝜆_𝑖∶𝑛(𝑥)𝑑𝑥, 𝜆>0 and 𝜆≠1. From (54), it follows that 𝐻(𝜆)=1[𝐵(𝑖,𝑛−𝑖+1)]^𝜆∞0𝑓^𝜆(𝑥)×[𝐹(𝑥)]^{𝜆(𝑖−1)}[1−𝐹(𝑥)]^{𝜆(𝑛−𝑖)}𝑑𝑥. (70) Using (40) in (70), we obtain 𝐻(𝜆)=1[𝐵(𝑖,𝑛−𝑖+1)]^𝜆∞_𝑘1=0(−1)^𝑘1𝜆()𝑘𝑖−1₁×^∞₀𝑓^𝜆(𝑥)[𝐹(𝑥)]^{𝜆(𝑖−1)+𝑘1}𝑑𝑥. (71) For 𝜆>0 and 𝜆≠1, we can expand [𝐹(𝑥)]^{𝜆(𝑖−1)+𝑘1} as 𝐹(𝑥)^{𝜆(𝑖−1)+𝑘1}={1−[1−𝐹(𝑥)]}^{𝜆(𝑖−1)+𝑘1}=^∞_𝑝1=0(−1)^𝑝1𝜆(𝑖−1)+𝑘₁𝑝+1₁[1−𝐹(𝑥)]^𝑝1, (72) and then 𝐹(𝑥)^{𝜆(𝑖−1)+𝑘1}=^∞_𝑝1^𝑝₌₀¹_𝑙1=0(−1)^𝑝1+𝑙1×𝜆(𝑖−1)+𝑘₁𝑝+1₁𝑝₁𝑙₁𝐹(𝑥)^𝑙1. (73) Hence, from (70), we can write 𝐻(𝜆)=1[𝐵(𝑖,𝑛−𝑖+1)]^𝜆×^∞_𝑘1,𝑝1^𝑝₌₀¹_𝑙1=0(−1)^{𝑘1+𝑝1+𝑙1}×𝜆()𝑘𝑖−1₁𝜆(𝑖−1)+𝑘₁𝑝+1₁𝑝₁𝑙₁×^∞₀𝑓^𝜆(𝑥)𝐹(𝑥)^𝑙1𝑑𝑥. (74) By replacing 𝑓(𝑥) and 𝐹(𝑥) by (10) and (16) given by Pescim et al. [3], we obtain 𝐻(𝜆)=2^{𝜆(𝑏−1)}𝛼2/𝜋^𝜆[𝐵(𝑖,𝑛−𝑖+1)]^𝜆×^∞_𝑘1,𝑝1^𝑝₌₀¹_𝑙1=0[Γ(𝑏)]^𝑙1(−1)^{𝑘1+𝑝1+𝑙1}[𝐵(𝑎,𝑏)]^𝑙1×𝜆()𝑘𝑖−1₁𝜆(𝑖−1)+𝑘₁𝑝+1₁𝑝₁𝑙₁×^∞₀𝑥^−𝜆𝑥𝜃^𝛼𝜆𝑒^{−(𝜆/2)(𝑥/𝜃)2𝛼}×𝑥2Φ𝜃^𝛼−1^{𝜆(𝑎−1)}×𝑥1−Φ𝜃^{𝛼𝜆(𝑏−1)}×^∞_𝑗=0(−1)^𝑗2Φ(𝑥/𝜃)^𝛼−1^𝑗Γ𝑏−𝑗𝑗!𝑎+𝑗^𝑙1𝑑𝑥. (75)

Using the identity (∑^∞_𝑖=0𝑎_𝑖)^𝑘=∑^∞_{{𝑚1,…,𝑚𝑘}=0}𝑎_𝑚1⋯𝑎_𝑚𝑘 (for 𝑘 positive integer) in (4), we have 𝐻(𝜆)=2^{𝜆(𝑏−1)}𝛼2/𝜋^𝜆[𝐵(𝑖,𝑛−𝑖+1)]^𝜆×^∞_𝑘1,𝑝1^𝑝₌₀¹_𝑙1=0𝑐_{𝑘1,𝑝1,𝑙1}()×𝑎,𝑏_∞𝑚1=0⋯^∞_𝑚𝑙1=0(−1)^{𝑚1+⋯+𝑚𝑙1}∏^𝑙1_𝑗=1Γ𝑏−𝑚_𝑗𝑚_𝑗!𝑎+𝑚_𝑗×^∞₀𝑥^−𝜆𝑥𝜃^𝛼𝜆𝑒^{−(𝜆/2)(𝑥/𝜃)2𝛼}×𝑥2Φ𝜃^𝛼−1^{𝑎(𝜆+𝑙1)−𝜆+∑𝑙1𝑗=1𝑚𝑗}×𝑥1−Φ𝜃^{𝛼𝜆(𝑏−1)}𝑑𝑥, (76) where 𝑐_{𝑘1,𝑝1,𝑙1}()=()𝑎,𝑏−1^{𝑘1+𝑝1+𝑙1}[Γ(𝑏)]^𝑙1[𝐵(𝑎,𝑏)]^𝑙1×𝜆()𝑘𝑖−1₁𝜆(𝑖−1)+𝑘₁𝑝+1₁𝑝₁𝑙₁. (77) Using the power series (40) in (76) twice and replacing Φ(𝑥) by the error function erf(𝑥) (defined in Section 3.2) we obtain after some algebra 𝐻(𝜆)=𝛼^𝜆−1𝜃^1−𝜆2^{[𝜆(𝑏−1)+(𝜆(𝛼−1)+1)/2𝛼−1]}2/𝜋^𝜆[𝐵(𝑖,𝑛−𝑖+1)]^𝜆×^∞_𝑘1,𝑝1^𝑝₌₀¹_𝑙1=0𝑐_{𝑘1,𝑝1,𝑙1}()×𝑎,𝑏_∞𝑚1=0⋯^∞_𝑚𝑙1=0(−1)^{𝑚1+⋯+𝑚𝑙1}∏^𝑙1_𝑗=1Γ𝑏−𝑚_𝑗𝑚_𝑗!𝑎+𝑚_𝑗×^∞_𝑟1,𝑠1^𝑠₌₀¹_𝑣1=0𝑑_{𝑟1,𝑠1,𝑣1}×^∞_𝑚1=0⋯^∞_𝑚𝑠1=0(−1)^{𝑚1+⋯+𝑚𝑠1}2𝑚₁+1⋯2𝑚_𝑠1+1𝑚₁!⋯𝑚_𝑠1!×𝜆^{−[𝑚1+⋯+𝑚𝑠1+𝑠1/2+(𝜆(𝛼−1)+1)/2𝛼]}×Γ𝑚₁+⋯+𝑚_𝑠1+𝑠₁2+𝜆(𝛼−1)+1,2𝛼 (78) where the quantity 𝑑_{𝑟,𝑠,𝑡} is well defined by Pescim et al. [3]. (More details, see equation (30) in [3]).

Finally, the entropy of the order statistics can be expressed as 𝚥_𝑅(𝜆)=(1−𝜆)⁻¹×()(𝛼)+()(𝜃)+𝜆()+𝜆()𝜆−1log1−𝜆log𝑏−1𝛼−1+1(2)2𝛼−1log+𝜆log2𝜋−𝜆log[𝐵(𝑖,𝑛−𝑖+1)]+log^∞_𝑘1,𝑝1^𝑝₌₀¹_𝑙1=0𝑐_{𝑘1,𝑝1,𝑙1}(𝑎,𝑏)+log^∞_𝑚1=0⋯^∞_𝑚𝑙1=0(−1)^{𝑚1+⋯+𝑚𝑙1}∏^𝑙1_𝑗=1Γ𝑏−𝑚_𝑗𝑚_𝑗!𝑎+𝑚_𝑗+log^∞_𝑟1,𝑠1^𝑠₌₀¹_𝑣1=0𝑑_{𝑟1,𝑠1,𝑣1}+log^∞_𝑚1=0⋯^∞_𝑚𝑠1=0(−1)^{𝑚1+⋯+𝑚𝑠1}×2𝑚₁+1⋯2𝑚_𝑠1+1×𝑚₁!⋯𝑚_𝑠1!⁻¹−𝑚₁+⋯+𝑚_𝑠1+𝑠₁2+𝜆(𝛼−1)+1(𝜆)Γ𝑚2𝛼log+log₁+⋯+𝑚_𝑠1+𝑠₁2+𝜆(𝛼−1)+1.2𝛼 (79) Equation (79) is the main result of this section.

An alternative expansion for the density function of the order statistics derived from the identity (20) was proposed by Pescim et al. [3]. A second density function representation and alternative expressions for some measures of the order statistics are given in Appendix A.

6. Lifetime Analysis

Let 𝑋_𝑖 be a random variable having density function (4), where 𝝎=(𝛼,𝜃,𝑎,𝑏)^𝑇 is the vector of parameters. The data encountered in survival analysis and reliability studies are often censored. A very simple random censoring mechanism that is often realistic is one in which each individual 𝑖 is assumed to have a lifetime 𝑋_𝑖 and a censoring time 𝐶_𝑖, where 𝑋_𝑖 and 𝐶_𝑖 are independent random variables. Suppose that the data consist of 𝑛 independent observations 𝑥_𝑖=min(𝑋_𝑖,𝐶_𝑖) for 𝑖=1,…,𝑛. The distribution of 𝐶_𝑖 does not depend on any of the unknown parameters of the distribution of 𝑋_𝑖. Parametric inference for such data are usually based on likelihood methods and their asymptotic theory. The censored log-likelihood 𝑙(𝝎) for the model parameters is 𝑙(𝝎)=𝑟log2𝜋+_𝑖∈𝐹𝛼log𝑥_𝑖+𝛼_𝑖∈𝐹𝑥log^𝑖𝜃−12_𝑖∈𝐹𝑥_𝑖𝜃^2𝛼()(2)[𝐵(+()+𝑟𝑏−1log−𝑟log𝑎,𝑏)]𝑎−1_𝑖∈𝐹𝑥log2Φ^𝑖𝜃^𝛼+()−1𝑏−1_𝑖∈𝐹𝑥log1−Φ^𝑖𝜃^𝛼+_𝑖∈𝐶log1−𝐼_{2Φ[(𝑥𝑖/𝜃)}^𝛼_]−1(𝑎,𝑏), (80) where 𝑟 is the number of failures and 𝐹 and 𝐶 denote the uncensored and censored sets of observations, respectively.

The score functions for the parameters 𝛼, 𝜃, 𝑎, and 𝑏 are given by 𝑈_𝛼(𝝎)=𝑟2+_𝑖∈𝐹𝑥log^𝑖𝜃𝑥1−^𝑖𝜃^2𝛼+2(𝑎−1)√2𝜋_𝑖∈𝐹𝑣_𝑖log𝑥_𝑖/𝜃𝑃𝑥_𝑖+(1−𝑏)√𝜋_𝑖∈𝐹𝑣_𝑖log𝑥_𝑖/𝜃𝐷𝑥_𝑖−2^𝑏−1_𝑖∈𝐶𝑣_𝑖log𝑥_𝑖/𝜃𝑃𝑥_𝑖^𝑎−1𝐷𝑥_𝑖^𝑏−11−𝑄𝑥_𝑖,𝑈_𝜃(𝝎)𝛼=−𝑟𝜃+𝛼𝜃_𝑖∈𝐹𝑥_𝑖𝜃^2𝛼+2(𝑎−1)√2𝜋×_𝑖∈𝐹𝑣_𝑖(𝛼/𝜃)𝑃𝑥_𝑖+(1−𝑏)√𝜋_𝑖∈𝐹𝑣_𝑖(𝛼/𝜃)𝐷𝑥_𝑖−2^𝑏−1_𝑖∈𝐶𝑣_𝑖(𝛼/𝜃)𝑃𝑥_𝑖^𝑎−1𝐷𝑥_𝑖^𝑏−11−𝑄𝑥_𝑖,𝑈_𝑎(𝝎)=𝑟𝜓(𝑎+𝑏)−𝜓(𝑎)+_𝑖∈𝐹log𝑃𝑥_𝑖−_𝑖∈𝐶𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑎1−𝑄𝑥_𝑖,𝑈_𝑏(𝝎)=𝑟log(2)+𝑟𝜓(𝑎+𝑏)−𝜓(𝑏)+_𝑖∈𝐹log𝐷𝑥_𝑖−_𝑖∈𝐶𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑏1−𝑄𝑥_𝑖, (81) where 𝑣_𝑖−1=exp2𝑥_𝑖𝜃^2𝛼𝑥_𝑖𝜃^𝛼,𝑃𝑥_𝑖𝑥=2Φ^𝑖𝜃^𝛼𝐷𝑥−1,_𝑖𝑥=1−Φ^𝑖𝜃^𝛼,𝑄𝑥_𝑖=𝐼_{2Φ[(𝑥𝑖/𝜃)}^𝛼_]−1(),𝐼̇𝑎,𝑏_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑎=𝜕1𝜕𝑎𝐵(𝑎,𝑏)^{𝑃(𝑥𝑖)}₀𝑤^𝑎−1(1−𝑤)^𝑏−1,𝐼̇𝑑𝑤_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑏=𝜕1𝜕𝑏𝐵(𝑎,𝑏)^{𝑃(𝑥𝑖)}₀𝑤^𝑎−1(1−𝑤)^𝑏−1𝑑𝑤 (82) and 𝜓(⋅) is the digamma function.

The maximum likelihood estimate (MLE) 𝝎 of 𝝎 is obtained by solving numerically the nonlinear equations 𝑈_𝛼(𝝎)=0, 𝑈_𝜃(𝝎)=0, 𝑈_𝑎(𝝎)=0, and 𝑈_𝑏(𝝎)=0. We can use iterative techniques such as the Newton-Raphson type algorithm to obtain the estimate 𝝎. We employ the subroutine NLMixed in SAS [15].

For interval estimation of (𝛼,𝜃,𝑎,𝑏) and tests of hypotheses on these parameters we obtain the observed information matrix since the expected information matrix is very complicated and will require numerical integration. The 4×4 observed information matrix 𝐉(𝝎) is 𝐉(𝝎)𝐋=−^𝛼𝛼𝐋_𝛼𝜃𝐋_𝛼𝑎𝐋_𝛼𝑏⋅𝐋_𝜃𝜃𝐋_𝜃𝑎𝐋_𝜃𝑏⋅⋅𝐋_𝑎𝑎𝐋_𝑎𝑏⋅⋅⋅𝐋_𝑏𝑏, (83) whose elements are given in Appendix B.

Under conditions that are fulfilled for parameters in the interior of the parameter space but not on the boundary, the asymptotic distribution of √𝑛(𝝎−𝝎) is 𝑁₄(0,𝐊(𝝎)⁻¹), where 𝐊(𝝎) is the expected information matrix. The normal approximation is valid if 𝐊(𝝎) is replaced by 𝐉(𝝎), that is, the observed information matrix evaluated at 𝝎. The multivariate normal 𝑁₄(0,𝐉(𝝎)⁻¹) distribution can be used to construct approximate confidence intervals for the model parameters. The likelihood ratio (LR) statistic can be used for testing the goodness of fit of the BGHN distribution and for comparing this distribution with some of its special submodels.

We can compute the maximum values of the unrestricted and restricted log-likelihoods to construct LR statistics for testing some submodels of the BGHN distribution. For example, we may use the LR statistic to check if the fit using the BGHN distribution is statistically “superior” to a fit using the modified Weibull or exponentiated Weibull distribution for a given data set. In any case, hypothesis testing of the type 𝐻₀:𝝎=𝝎₀ versus 𝐻:𝝎≠𝝎₀ can be performed using LR statistics. For example, the test of 𝐻₀:𝑏=1 versus 𝐻:𝑏≠1 is equivalent to compare the EGHN and BGHN distributions. In this case, the LR statistic becomes 𝑤=2{𝑙(𝛼,𝜃,𝑎,𝑏̂)−𝑙(𝛼,𝜃,𝑎,1)}, where 𝛼, 𝜃, 𝑎, and 𝑏̂ are the MLEs under 𝐻 and 𝛼, 𝜃 and 𝑎 are the estimates under 𝐻₀. We emphasize that first-order asymptotic results for LR statistics can be misleading for small sample sizes. Future research can be conducted to derive analytical or bootstrap Bartlett corrections.

7. A BGHN Model for Survival Data with Long-Term Survivors

In population based cancer studies, cure is said to occur when the mortality in the group of cancer patients returns to the same level as that expected in the general population. The cure fraction is of interest to patients and also a useful measure when analyzing trends in cancer patient survival. Models for survival analysis typically assume that every subject in the study population is susceptible to the event under study and will eventually experience such event if the follow-up is sufficiently long. However, there are situations when a fraction of individuals are not expected to experience the event of interest, that is, those individuals are cured or not susceptible. For example, researchers may be interested in analyzing the recurrence of a disease. Many individuals may never experience a recurrence; therefore, a cured fraction of the population exists. Cure rate models have been used to estimate the cured fraction. These models are survival models which allow for a cured fraction of individuals. These models extend the understanding of time-to-event data by allowing for the formulation of more accurate and informative conclusions. These conclusions are otherwise unobtainable from an analysis which fails to account for a cured or insusceptible fraction of the population. If a cured component is not present, the analysis reduces to standard approaches of survival analysis. Cure rate models have been used for modeling time-to-event data for various types of cancers, including breast cancer, non-Hodgkins lymphoma, leukemia, prostate cancer, and melanoma. Perhaps the most popular type of cure rate models is the mixture model introduced by Berkson and Gage [16] and Maller and Zhou [17]. In this model, the population is divided into two sub-populations so that an individual either is cured with probability 𝑝, or has a proper survival function 𝑆(𝑥) with probability 1−𝑝. This gives an improper population survivor function 𝑆^∗(𝑥) in the form of the mixture, namely, 𝑆^∗(𝑥)=𝑝+1−𝑝𝑆(𝑥),𝑆(∞)=0,𝑆^∗(∞)=𝑝. (84) Common choices for 𝑆(𝑥) in (84) are the exponential and Weibull distributions. Here, we adopt the BGHN distribution. Mixture models involving these distributions have been studied by several authors, including Farewell [18], Sy and Taylor [19], and Ortega et al. [20]. The book by Maller and Zhou [17] provides a wide range of applications of the long-term survivor mixture model. The use of survival models with a cure fraction has become more and more frequent because traditional survival analysis do not allow for modeling data in which nonhomogeneous parts of the population do not represent the event of interest even after a long follow-up. Now, we propose an application of the BGHN distribution to compose a mixture model for cure rate estimation. Consider a sample 𝑥₁,…,𝑥_𝑛, where 𝑥_𝑖 is either the observed lifetime or censoring time for the 𝑖th individual. Let a binary random variable 𝑧_𝑖, for 𝑖=1,…,𝑛, indicating that the 𝑖th individual in a population is at risk or not with respect to a certain type of failure, that is, 𝑧_𝑖=1 indicates that the 𝑖th individual will eventually experience a failure event (uncured) and 𝑧_𝑖=0 indicates that the individual will never experience such event (cured).

For an individual, the proportion of uncured 1−𝑝 can be specified such that the conditional distribution of 𝑧 is given by Pr(𝑧_𝑖=1)=1−𝑝. Suppose that the 𝑋_𝑖’s are independent and identically distributed random variables having the BGHN distribution with density function (4).

The maximum likelihood method is used to estimate the parameters. Thus, the contribution of an individual that failed at 𝑥_𝑖 to the likelihood function is given by 1−𝑝2^𝑏−12/𝜋𝛼/𝑥_𝑖𝑥_𝑖/𝜃^𝛼exp−(1/2)𝑥_𝑖/𝜃^2𝛼𝐵()×𝑥𝑎,𝑏2Φ^𝑖𝜃^𝛼−1^𝑎−1𝑥1−Φ^𝑖𝜃^𝛼𝑏−1 (85) and the contribution of an individual that is at risk at time 𝑥_𝑖 is 𝑝+1−𝑝1−𝐼_{2Φ[(𝑥𝑖/𝜃)}^𝛼_]−1(𝑎,𝑏). (86) The new model defined by (85) and (86) is called the BGHN mixture model with long-term survivors. For 𝑎=𝑏=1, we obtain the GHN mixture model (a new model) with long-term survivors.

Thus, the log-likelihood function for the parameter vector 𝜽=(𝑝,𝛼,𝜃,𝑎,𝑏)^𝑇 follows from (85) and (86) as 𝑙(𝜽)2=𝑟log1−𝑝^𝑏−12/𝜋𝐵(𝑎,𝑏)+_𝑖∈𝐹𝛼log𝑥_𝑖+𝛼_𝑖∈𝐹𝑥_𝑖𝜃−12_𝑖∈𝐹𝑥log^𝑖𝜃^2𝛼+(𝑎−1)_𝑖∈𝐹𝑥log2Φ^𝑖𝜃^𝛼+()−1𝑏−1_𝑖∈𝐹𝑥log1−Φ^𝑖𝜃^𝛼+_𝑖∈𝐶log𝑝+1−𝑝1−𝐼_{2Φ[(𝑥𝑖/𝜃)}^𝛼_]−1(𝑎,𝑏), (87) where 𝐹 and 𝐶 denote the sets of individuals corresponding to lifetime observations and censoring times, respectively, and 𝑟 is the number of uncensored observations (failures).

8. Applications

In this section, we give four applications using well-known data sets (three with censoring and one uncensored data set) to demonstrate the flexibility and applicability of the proposed model. The reason for choosing these data is that they allow us to show how in different fields it is necessary to have positively skewed distributions with nonnegative support. These data sets present different degrees of skewness and kurtosis.

8.1. Censored Data: Transistors

The first data set consists of the lifetimes of 𝑛=34 transistors in an accelerated life test. Three of the lifetimes are censored, so the censoring fraction is 3/34 (or 8.8%). Wilk et al. [21] stated that “there is reason, from past experience, to expect that the gamma distribution might reasonably approximate the failure time distribution.” Wilk et al. [21] and also Lawless [22] fitted a gamma distribution. Alternatively, we fit the BGHN model to these data. We estimate the parameters 𝛼, 𝜃, 𝑎, and 𝑏 by maximum likelihood. The MLEs of 𝛼 and 𝜃 for the GHN distribution are taken as starting values for the iterative procedure. The computations were performed using the NLMixed procedure in SAS [15]. Table 2 lists the MLEs (and the corresponding standard errors in parentheses) of the model parameters and the values of the following statistics for some models: Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), and Consistent Akaike Information Criterion (CAIC). The results indicate that the BGHN model has the lowest AIC, BIC, and CAIC values, and, therefore, it could be chosen as the best model. Further, we calculate the maximum unrestricted and restricted log-likelihoods and the LR statistics for testing some submodels. For example, the LR statistic for testing the hypotheses 𝐻₀: 𝑎=𝑏=1 versus 𝐻₁: 𝐻₀  is not true, that is, to compare the BGHN and GHN models, is 𝑤=2{−120.1−(−126.0)}=11.8 (𝑃 value =0.0027), which gives favorable indication towards the BGHN model. In special, the Weibull density function is given by 𝑓(𝑥)=𝛾𝜆^𝛾𝑥^𝛾−1−𝑥exp𝜆^𝛾,𝑥,𝜆,𝛾>0. (88) Figure 1(a) displays plots of the empirical survival function and the estimated survival functions of the BGHN, GHN and Weibull distributions. We obtain a good fit of the BGHN model to these data.

Table 2 

MLEs of the model parameters for the transistors data, the corresponding SEs (given in parentheses), and the AIC, CAIC, and BIC statistics.

(a)

(b)

(c)

Figure 1 

Estimated survival function for the BGHN, GHN, and Weibull distributions and the empirical survival; (a) for transistors data, (b) for radiotherapy data, and (c) USS Halfbeak diesel engine data.

8.2. Censored Data: Radiotherapy

The data refer to the survival time (days) for cancer patients (𝑛=51) undergoing radiotherapy [23]. The percentage of censored observations was 17.65%. Fitting the BGHN distribution to these data, we obtain the MLEs of the model parameters listed in Table 3. The values of the AIC, CAIC, and BIC statistics are smaller for the BGHN distribution as compared to those values of the GHN and Weibull distributions. The LR statistic for testing the hypotheses 𝐻₀: 𝑎=𝑏=1 versus 𝐻₁: 𝐻₀  is not true, that is, to compare the BGHN and GHN models, is 𝑤=2{−293.3−(−299.4)}=12.2 (𝑃 value =0.0022), which yields favorable indication towards the BGHN distribution. Thus, this distribution seems to be a very competitive model for lifetime data analysis. Figure 1(b) displays the plots of the empirical survival function and the estimated survival functions for the BGHN, GHN, and Weibull distributions. We note a good fit of the BGHN distribution.

Table 3 

MLEs of the model parameters for the radiotherapy data, the corresponding SEs (given in parentheses), and the AIC, CAIC, and BIC statistics.

8.3. Uncensored Data: USS Halfbeak Diesel Engine

Here, we compare the fits of the BGHN, GHN, and Weibull distributions to the data set presented by Ascher and Feingold [24] from a USS Halfbeak (submarine) diesel engine. The data denote 73 failure times (in hours) of unscheduled maintenance actions for the USS Halfbeak number 4 main propulsion diesel engine over 25.518 operating hours. Table 4 gives the MLEs of the model parameters. The values of the AIC, CAIC, and BIC statistics are smaller for the BGHN distribution when compared to those values of the GHN and Weibull distributions. The LR statistic for testing the hypotheses 𝐻₀: 𝑎=𝑏=1 versus 𝐻₁: 𝐻₀  is not true, that is, to compare the BGHN and GHN models, is 𝑤=2{−196.1−(−217.2)}=42.1 (𝑃 value <0.0001), which indicates that the first distribution is superior to the second in terms of model fitting. Figure 1(c) displays plots of the empirical survival function and the estimated survival function of the BGHN, GHN, and Weibull distributions. In fact, the BGHN distribution provides a better fit to these data.

Table 4 

MLEs of the model parameters for the USS Halfbeak diesel engine data, the corresponding SEs (given in parentheses), and the AIC, CAIC, and BIC statistics.

8.4. Melanoma Data with Long-Term Survivors

In this section, we provide an application of the BGHN model to cancer recurrence. The data are part of a study on cutaneous melanoma (a type of malignant cancer) for the evaluation of postoperative treatment performance with a high dose of a certain drug (interferon alfa-2b) in order to prevent recurrence. Patients were included in the study from 1991 to 1995, and follow-up was conducted until 1998. The data were collected by Ibrahim et al. [25]. The survival time 𝑇 is defined as the time until the patient’s death. The original sample size was 𝑛=427 patients, 10 of whom did not present a value for explanatory variable tumor thickness. When such cases were removed, a sample of size 𝑛=417 patients was retained. The percentage of censored observations was 56%. Table 5 lists the MLEs of the model parameters. The values of the AIC, CAIC and BIC statistics are smaller for the BGHN mixture model when compared to those values of the GHN mixture model. The LR statistic for testing the hypotheses 𝐻₀: 𝑎=𝑏=1 versus 𝐻₁: 𝐻₀  is not true, that is, to compare the BGHN and GHN mixture models, becomes 𝑤=2{−524.75−(−534.05)}=18.6 (𝑃 value <0.0001), which indicates that the BGHN mixture model is superior to the GHN mixture model in terms of model fitting. Figure 2 provides plots of the empirical survival function and the estimated survival functions of the BGHN, and GHN mixture models. Note that the cured proportion estimated by the BGHN mixture model (𝑝_BGHN=0.4871) is more appropriate than that one estimated by the GHN mixture model (𝑝_GHN=0.5150). Further, the BGHN mixture model provides a better fit to these data.

Table 5 

MLEs of the model parameters for the melanoma data, the corresponding SEs (given in parentheses), and the AIC, CAIC, and BIC statistics.

Figure 2 

Estimated survival functions for the BGHN and GHN mixture models and the empirical survival for the melanoma data.

9. Concluding Remarks

In this paper, we discuss some mathematical properties of the beta generalized half normal (BGHN) distribution introduced recently by Pescim et al. [3]. The incorporation of additional parameters provides a very flexible model. We derive a power series expansion for the BGHN quantile function. We provide some new structural properties such as moments, generating function, mean deviations, Rényi entropy, and reliability. We investigate properties of the order statistics. The method of maximum likelihood is used for estimating the model parameters for uncensored and censored data. We also propose a BGHN model for survival data with long-term survivors. Applications to four real data sets indicate that the BGHN model provides a flexible alternative to (and a better fit than) some classical lifetimes models.

Appendices

A. Appendix A

Here, we derive some properties of the BGHN order statistics. Using the identity in Gradshteyn and Ryzhik [5] and after some algebraic manipulations, we obtain an alternative form for the density function of the 𝑖th order statistic, namely, 𝑓_𝑖∶𝑛(𝑥)=^{𝑛−𝑖∞}_𝑘=0𝑗=0(−1)^𝑘_𝑘^𝑛−𝑖Γ(𝑏)^{𝑖+𝑘−1}𝐵𝑎(𝑖+𝑘)+𝑗,𝑏𝑑_{𝑖,𝑗,𝑘}𝐵(𝑎,𝑏)^𝑖+𝑘𝐵(𝑖,𝑛−1+𝑖)×𝑓_{𝑖,𝑗,𝑘}(𝑥), (A.1) where 𝑓_{𝑖,𝑗,𝑘}(𝑥) denotes the density of the BGHN (𝛼,𝜃,𝑎(𝑖+𝑘)+𝑗,𝑏) distribution and the quantity 𝑑_{𝑖,𝑗,𝑘} is defined by Pescim et al. [3].

From (A.1), we obtain alternative expressions for the moments and mgf of the BGHN order statistics given by 𝐸𝑋^𝑠_𝑖∶𝑛=^{𝑛−𝑖∞}_𝑘=0𝑗=0(−1)^𝑘_𝑘^𝑛−𝑖Γ(𝑏)^{𝑖+𝑘−1}𝐵𝑎(𝑖+𝑘)+𝑗,𝑏𝑑_{𝑖,𝑗,𝑘}𝐵(𝑎,𝑏)^𝑖+𝑘𝐵()𝑋𝑖,𝑛−𝑖+1×𝐸^𝑠_{𝑖,𝑗,𝑘},𝑀_𝑖∶𝑛(𝑡)=^{𝑛−𝑖∞}_𝑘=0𝑗=0(−1)^𝑘_𝑘^𝑛−𝑖Γ(𝑏)^{𝑖+𝑘−1}𝐵𝑎(𝑖+𝑘)+𝑗,𝑏𝑑_{𝑖,𝑗,𝑘}𝐵(𝑎,𝑏)^𝑖+𝑘𝐵(𝑖,𝑛−𝑖+1)×𝑀_{𝛼,𝜃,𝑎,𝑏}(𝑡). (A.2)

The mean deviations of the order statistics in relation to the mean and to the median can be expressed as 𝛿₁𝑋_𝑖∶𝑛=2𝜇_𝑖𝐹_𝑖∶𝑛𝜇_𝑖−2𝜇_𝑖+2𝐽_𝑖𝜇_𝑖,𝛿₂𝑋_𝑖∶𝑛=2𝐽_𝑖𝑀_𝑖−𝜇_𝑖, (A.3) respectively, where 𝜇_𝑖=𝐸(𝑋_𝑖∶𝑛), 𝜇_𝑖=Median(𝑋_𝑖∶𝑛) and 𝐽_𝑖(𝑞) are defined in Section 5.4. We use (A.1) to obtain 𝐽_𝑖(𝑞)𝐽_𝑖𝑞=^{𝑛−𝑖∞}_𝑘=0𝑗=0(−1)^𝑘_𝑘^𝑛−𝑖Γ(𝑏)^{𝑖+𝑘−1}𝐵𝑎(𝑖+𝑘)+𝑗,𝑏𝑑_{𝑖,𝑗,𝑘}𝐵(𝑎,𝑏)^𝑖+𝑘𝐵()𝑞,𝑖,𝑛−𝑖+1×𝑇 (A.4) where 𝑇(𝑞) comes from (35).

Bonferroni and Lorenz curves of the order statistic are defined, in Section 5.4, by 𝐵_𝑖∶𝑛𝜌=1𝜌𝜇_𝑖^𝑞₀𝑥𝑓_𝑖∶𝑛(𝑥)𝐿𝑑𝑥,_𝑖∶𝑛𝜌=1𝜇_𝑖^𝑞₀𝑥𝑓_𝑖∶𝑛(𝑥)𝑑𝑥, (A.5) respectively, where 𝑞=𝐹⁻¹_𝑖∶𝑛(𝜌). From ∫^𝑞₀𝑥𝑓_𝑖∶𝑛(𝑥)𝑑𝑥=𝜇_𝑖−𝐽_𝑖(𝜌), these curves also can be expressed as 𝐵_𝑖∶𝑛𝜌=1𝜌−𝐽_𝑖𝜌𝜌𝜇_𝑖,𝐿_𝑖∶𝑛𝜌𝐽=1−^𝑖𝜌𝜇_𝑖. (A.6) Using (A.4), we obtain the Bonferroni and Lorenz curves for the BGHN order statistics.

B. Appendix B

Here, we give the necessary formulas to obtain the second-order partial derivatives of the log-likelihood function. After some algebraic manipulations, we obtain 𝐋_𝛼𝛼=2_𝑖∈𝐹𝑥_𝑖𝜃^2𝛼log²𝑥_𝑖𝜃+2(𝑎−1)√×2𝜋_𝑖∈𝐹𝑣_𝑖log²𝑥_𝑖/𝜃1−𝑥_𝑖/𝜃^2𝛼𝑃𝑥_𝑖−2√2𝜋_𝑖∈𝐹𝑣_𝑖log𝑥_𝑖/𝜃𝑃𝑥_𝑖²+(1−𝑏)√𝜋_𝑖∈𝐹𝑣_𝑖log²𝑥_𝑖/𝜃1−𝑥_𝑖/𝜃^2𝛼𝐷𝑥_𝑖+1√𝜋_𝑖∈𝐹𝑣_𝑖log𝑥_𝑖/𝜃𝐷𝑥_𝑖²−2^𝑏−1_𝑖∈𝐶𝑣_𝑖log²𝑥_𝑖𝜃𝑥1−^𝑖𝜃^2𝛼×𝑃𝑥_𝑖^𝑎−1𝐷𝑥_𝑖^𝑏−1×1−𝑄𝑥_𝑖⁻¹+2(𝑎−1)√2𝜋_𝑖∈𝐶𝑣²_𝑖log²𝑥_𝑖𝜃𝑥1−^𝑖𝜃^2𝛼×𝑃𝑥_𝑖^𝑎−2𝐷𝑥_𝑖^𝑏−1×[1−𝑄𝑥_𝑖]⁻¹+(1−𝑏)√𝜋_𝑖∈𝐶𝑣²_𝑖log²𝑥_𝑖𝜃𝑥1−^𝑖𝜃^2𝛼×𝑃(𝑥_𝑖)^𝑎−1𝐷𝑥_𝑖^𝑏−2×[1−𝑄𝑥_𝑖]⁻¹+2^𝑏−1_𝑖∈𝐶𝑣_𝑖𝑥log^𝑖𝜃𝑃𝑥_𝑖^(𝑎−1)×𝐷𝑥_𝑖^(𝑏−1)×1−𝑄𝑥_𝑖⁻¹²,𝐋_𝛼𝜃𝑟=−𝜃𝛼+2𝜃_𝑖∈𝐹𝑥_𝑖𝜃^2𝛼𝑥log^𝑖𝜃+1𝜃_𝑖∈𝐹𝑥_𝑖𝜃^2𝛼+2(𝑎−1)√2𝜋_𝑖∈𝐹𝑣_𝑖log²𝑥_𝑖/𝜃1−𝑥_𝑖/𝜃^2𝛼+𝑣_𝑖/𝜃𝑃𝑥_𝑖−2√2𝜋_𝑖∈𝐹𝑣_𝑖(𝛼/𝜃)𝑃𝑥_𝑖²+(1−𝑏)√𝜋_𝑖∈𝐹𝑣_𝑖log𝑥_𝑖/𝜃1−𝑥_𝑖/𝜃^2𝛼+𝑣_𝑖/𝜃𝐷𝑥_𝑖+1√𝜋_𝑖∈𝐹𝑣_𝑖(𝛼/𝜃)𝐷𝑥_𝑖²−2^𝑏−1_𝑖∈𝐶𝑣_𝑖log²𝑥_𝑖𝜃𝑥1−^𝑖𝜃^2𝛼+𝑣_𝑖𝜃×𝑃𝑥_𝑖^𝑎−1𝐷𝑥_𝑖^𝑏−1×1−𝑄𝑥_𝑖⁻¹+2(𝑎−1)√2𝜋_𝑖∈𝐶𝑣²_𝑖𝛼𝜃𝑥log^𝑖𝜃𝑃𝑥_𝑖^𝑎−2×𝐷𝑥_𝑖^𝑏−1×1−𝑄𝑥_𝑖⁻¹+(1−𝑏)√𝜋_𝑖∈𝐶𝑣²_𝑖𝛼𝜃𝑥log^𝑖𝜃𝑃𝑥_𝑖^𝑎−1×𝐷𝑥_𝑖^𝑏−2×1−𝑄𝑥_𝑖⁻¹+2^𝑏−1_𝑖∈𝐶𝑣²_𝑖𝛼𝜃𝑥log^𝑖𝜃𝑃𝑥_𝑖^2(𝑎−1)×𝐷𝑥_𝑖^2(𝑏−1)×1−𝑄𝑥_𝑖²⁻¹,𝐋_𝛼𝑎=2(𝑎−1)√2𝜋_𝑖∈𝐹𝑣_𝑖log𝑥_𝑖/𝜃𝑃𝑥_𝑖−_𝑖∈𝐶𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑎𝛼1−𝑄𝑥_𝑖−2^𝑏−1_𝑖∈𝐶𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑎𝑣_𝑖𝑥log^𝑖𝜃𝑃𝑥_𝑖^𝑎−1×𝐷𝑥_𝑖^𝑏−1×1−𝑄𝑥_𝑖²⁻¹,𝐋_𝛼𝑏=(1−𝑏)√𝜋_𝑖∈𝐹𝑣_𝑖log𝑥_𝑖/𝜃𝐷𝑥_𝑖−_𝑖∈𝐶𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑏𝛼1−𝑄𝑥_𝑖−2^𝑏−1_𝑖∈𝐶𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑏𝑣_𝑖𝑥log^𝑖𝜃𝑃𝑥_𝑖^𝑎−1×𝐷𝑥_𝑖^𝑏−1×1−𝑄𝑥_𝑖²⁻¹,𝐋_𝜃𝜃𝛼=𝑟𝜃²−𝛼𝜃²_𝑖∈𝐹𝑥_𝑖𝜃^2𝛼𝛼−2𝜃²_𝑖∈𝐹𝑥_𝑖𝜃^2𝛼+2(𝑎−1)√2𝜋_𝑖∈𝐹𝑣_𝑖(𝛼/𝜃)𝑥_𝑖/𝜃^2𝛼−1/𝜃−1𝑃𝑥_𝑖−2√2𝜋_𝑖∈𝐹𝑣_𝑖(𝛼/𝜃)𝑃𝑥_𝑖²+(1−𝑏)√𝜋_𝑖∈𝐹𝑣_𝑖(𝛼/𝜃)𝑥_𝑖/𝜃^2𝛼−1/𝜃−1𝐷𝑥_𝑖+1√𝜋_𝑖∈𝐹𝑣_𝑖(𝛼/𝜃)𝐷𝑥_𝑖²−2^𝑏−1_𝑖∈𝐶𝑣_𝑖𝛼𝜃𝑥_𝑖𝜃^2𝛼−1𝜃−1𝑃𝑥_𝑖^𝑎−1×𝐷𝑥_𝑖^𝑏−1×[1−𝑄𝑥_𝑖]⁻¹+2(𝑎−1)√2𝜋_𝑖∈𝐶𝑣²_𝑖(𝛼/𝜃)²𝑃𝑥_𝑖^𝑎−2𝐷𝑥_𝑖^𝑏−11−𝑄𝑥_𝑖+(1−𝑏)√𝜋_𝑖∈𝐶𝑣²_𝑖(𝛼/𝜃)²𝑃𝑥_𝑖^𝑎−1𝐷𝑥_𝑖^𝑏−21−𝑄𝑥_𝑖+2^𝑏−1_𝑖∈𝐶𝑣²_𝑖(𝛼/𝜃)²𝑃𝑥_𝑖^2(𝑎−1)𝐷𝑥_𝑖^2(𝑏−1)1−𝑄𝑥_𝑖²,𝐋_𝜃𝑎=2(𝑎−1)√2𝜋_𝑖∈𝐹𝑣_𝑖(𝛼/𝜃)𝑃𝑥_𝑖−_𝑖∈𝐶𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑎𝜃1−𝑄𝑥_𝑖2^𝑏−1×_𝑖∈𝐶𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑎𝑣_𝑖(𝛼/𝜃)𝑃𝑥_𝑖^(𝑎−1)𝐷𝑥_𝑖^(𝑏−1)1−𝑄𝑥_𝑖²,𝐋_𝜃𝑏=(1−𝑏)√𝜋_𝑖∈𝐹𝑣_𝑖(𝛼/𝜃)𝐷𝑥_𝑖−_𝑖∈𝐶𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑏𝜃1−𝑄𝑥_𝑖2^𝑏−1×_𝑖∈𝐶𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑏𝑣_𝑖(𝛼/𝜃)𝑃𝑥_𝑖^(𝑎−1)𝐷𝑥_𝑖^(𝑏−1)1−𝑄𝑥_𝑖²,𝐋_𝑎𝑎=𝑟𝜓^′(𝑎+𝑏)−𝜓^′(𝑎)−_𝑖∈𝐶𝐼̈_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑎1−𝑄𝑥_𝑖−_𝑖∈𝐶𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑎1−𝑄𝑥_𝑖²,𝐋_𝑎𝑏=𝑟𝜓^′(𝑎+𝑏)−_𝑖∈𝐶𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑎𝑏1−𝑄𝑥_𝑖−_𝑖∈𝐶𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑎𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑏1−𝑄𝑥_𝑖²,𝐋_𝑏𝑏=𝑟𝜓^′(𝑎+𝑏)−𝜓^′(𝑏)−_𝑖∈𝐶𝐼̈_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑏1−𝑄𝑥_𝑖−_𝑖∈𝐶𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑏1−𝑄𝑥_𝑖², (B.1) where 𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑎𝛼=𝜕²𝑄𝑥𝜕𝑎𝜕𝛼^𝑖,𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑏𝛼=𝜕²𝑄𝑥𝜕𝑏𝜕𝛼^𝑖,𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑎𝜃=𝜕²𝑄𝑥𝜕𝑎𝜕𝜃^𝑖,𝐼̇_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑏𝜃=𝜕²𝑄𝑥𝜕𝑏𝜕𝜃^𝑖,𝐼̈_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑎=𝜕²𝜕𝑎²𝑄𝑥_𝑖,𝐼̈_{𝑃(𝑥𝑖)}(𝑎,𝑏)|_𝑏=𝜕²𝜕𝑏²𝑄𝑥_𝑖. (B.2)

Acknowledgments

The authors would like to thank the editor, the associate editor, and the referees for carefully reading the paper and for their comments, which greatly improved the paper.