# The Beta Generalized Half-Normal Distribution: New Properties

Authors: Gauss M. Cordeiro; Rodrigo R. Pescim; Edwin M. M. Ortega; Clarice G. B. Demétrio
Journal: Journal of Probability and Statistics
 (2013)

---

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


---

## Body

## 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)
(1)Gα,θ(x)=2Φ[(xθ)α]-1=erf[(x/θ)α2],
where the standard normal cdf Φ(x) and the error function erf(x) are given by
(2)Φ(x)=12[1+erf(x2)],erf(x)=2π∫0x‍e-t2dt.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 x>0),
(3)F(x)=I2Φ[(x/θ)α]-1(a,b)=1B(a,b)∫02Φ[(x/θ)α]-1‍wa-1(1-w)b-1dw,
where B(a,b)=[Γ(a)Γ(b)]/Γ(a+b) is the beta function, Iy(a,b)=B(a,b)-1∫0y‍wa-1(1-w)b-1dw is the incomplete beta function ratio, and a>0 and b>0 are two additional shape parameters.The probability density function (pdf) and the hazard rate function (hrf) corresponding to (3) are
(4)f(x)=2b-12/π(α/x)(x/θ)αe-(1/2)(x/θ)2αB(a,b)×{2Φ[(xθ)α]-1}a-1{1-Φ[(xθ)α]}b-1,(5)h(x)=(2b-12π(αx)(xθ)αe-(1/2)(x/θ)2α×{2Φ[(xθ)α]-1}a-1{1-Φ[(xθ)α]}b-1)×(B(a,b){1-I2Φ[(x/θ)α]-1(a,b)})-1,
respectively. Hereafter, a random variable with pdf (4) is denoted by X~BGHN(α,θ,a,b).Pescim et al. [3] demonstrated that the cdf and pdf of X can be expressed as infinite power series of the GHN cumulative distribution. Here, all expansions in power series are around the point zero. If b>0 is a real noninteger, we can expand the binomial term in (3) to obtain
(6)F(x)=∑j=0∞‍wj{2Φ[(xθ)α]-1}a+j,
where wj=wj(a,b)=(-1)jB(a,b)-1(a+j)-1(b-1j).The pdf corresponding to (6) can be expressed as
(7)f(x)=2π(αx)(xθ)αe-(1/2)(x/θ)2α×∑j=0∞‍wj{2Φ[(xθ)α]-1}a+j-1.
If a>0 is an integer, (7) provides the BGHN density function as an infinite power series of the GHN cumulative distribution. If b is an integer, the index j in the previous sum stops at b-1. Otherwise, if a is a real noninteger, we can expand {2Φ[(x/θ)α]-1}a+j-1 as follows:
(8){2Φ[(xθ)α]-1}a+j-1=∑r=0∞‍sr(a+j-1){2Φ[(xθ)α]-1}r,
where
(9)sr(a+j-1)=∑k=r∞‍(-1)r+ka+j(a+j-1k)(kr).
Hence,
(10)f(x)=2π(αx)(xθ)αe-(1/2)(x/θ)2α∑r=0∞‍br{erf[(x/θ)α2]}r,
whose coefficients are
(11)br=br(a,b)=1B(a,b)∑j=0∞‍(-1)j(b-1j)sr(a+j-1).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 a=b=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 b=1 and a=b=α=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 Section2. 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 havingF(x) 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 ofX, say x=QBGHN(u)=F-1(u), can be obtained by inverting the cumulative function (3). Now, we provide a power series expansion for QBGHN(u) that can be useful to determine some mathematical measures of the BGHN distribution. First, an expansion for the inverse of the incomplete beta function Ix(a,b)=u can be found in wolfram website (http://functions.wolfram.com/06.23.06.0004.01) as
(12)z=QB(u)=∑i=0∞‍diui/a,
where QB(u) is the beta quantile function, di=ei[aB(a,b)]1/a for a>0 and e0=0, e1=1, e2=(b-1)/(a+1),…. The coefficients ei′s for i≥2 can be obtained from a cubic recursion of the form
(13)ei=1[i2+(α-2)i+(1-α)]×{(1-δi,2)∑r=2i-1‍erei+1-r×[r(1-a)(i-r)-r(r-1)]+∑r=1i-1‍∑s=1i-r‍eresei+1-r-s×[r(r-a)+s(a+b-2)(i+1-r-s)]∑r=2i-1‍erei+1-r},
where δi,2=1 if i=2 and δi,2=0 if i≠2. In the last equation, we note that the quadratic term only contributes for i≥3.Following Steinbrecher [4], the quantile function of the standard normal distribution, say QN(u)=Φ-1(u), can be expanded as
(14)QN(u)=Φ-1(u)=∑k=0∞‍bkζ2k+1,
where ζ=2π(u-1/2) and the bk′s can be calculated from b0=1 and
(15)bk+1=12(2k+3)∑r=0k‍(2r+1)(2k-2r+1)brbk-r(r+1)(2r+1).
Here, b1=1/6, b2=7/120, b3=127/7560,….The functionQN(u) can be expressed as a power series given by
(16)QN(u)=∑j=0∞‍fjuj,
where fj=∑k=j∞‍(-1/2)k-j(kj)gk and the quantities gk are defined from the coefficients in (14) by gk=0 for k=0,2,4,… and gk=(2π)k/2b(k-1)/2 for k=1,3,5,….By inverting the normal cumulative function inz=2Φ[(x/θ)α]-1, and using (16), we can express the quantile function of X in terms of z as
(17)x=θ[QN(z+12)]1/α=θ[∑j=0∞‍fj(z+12)j]1/α=θ[∑j=0∞‍∑p=0j‍fj2j(jp)zp]1/α.
Replacing ∑j=0∞‍∑p=0j‍ by ∑p=0∞‍∑j=p∞‍ in last equation, we obtain
(18)x=θ(∑p=0∞‍hpzp)1/α,
where hp=∑j=p∞‍2-jfj(jp). Using the same steps of (8), we can write
(19)(∑p=0∞‍hpzp)1/α=∑r=0∞‍sr(α-1)(∑p=0∞‍hpzp)r,
where sr(α-1)=∑j=r∞‍(-1)r+j(α-1j)(jr).We use throughout an equation of Gradshteyn and Ryzhik [5] for a power series raised to a positive integer j(20)(∑i=0∞‍aixi)j=∑i=0∞‍cj,ixi,
where the coefficients cj,i (for i=1,2,…) are easily obtained from the recurrence equation
(21)cj,i=(ia0)-1∑m=1i‍[m(j+1)-i]amcj,i-m
and cj,0=a0j. From (19) and (20), we have
(22)x=∑p=0∞‍hp⋆zp,
where hp⋆=θ∑r=0∞‍sr(α-1)hr,p for p≥0, hr,p=(ph0)-1∑m=1p[m(r+1)-p]hmhr,p-m for p≥1 and hr,0=h0r. By inserting (12) in (22), we obtain
(23)x=QBGHN(u)=∑p=0∞‍hp⋆(∑i=0∞‍diui/a)p.
From (20), it follows (∑i=0∞‍diui/a)p=∑i=0∞‍dp,iui/a, where the quantities dp,i are determined recursively from d0,i=d0i and dp,i=(id0)-1∑m=1i‍[m(p+1)-i]dmdp,i-m (for i=1,2,…). Finally, we obtain
(24)x=QBGHN(u)=∑i=0∞‍viui/a,
where vi=∑p=0∞‍hp⋆dp,i. 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 ofX 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)
(25)μs′(α,θ,a,b)=E(Xs)=∫0∞‍xsf(x)dx=∫01‍(∑i=0∞‍viui/a)sdu,
and then using (20), we obtain
(26)μs′(α,θ,a,b)=a∑i=0∞‍vs,i(i+a),
where the quantities vs,i (for i=1,2,…) are easily determined from the recurrence equation
(27)vs,i=(iv0)-1∑m=1i‍[m(s+1)-i]vmvs,i-m,
with vs,0=v0s.We now provide a new alternative representation for the mgf ofX, say Mα,θ,a,b(t)=E(etX), based on the quantile power series (24). We can write
(28)Mα,θ,a,b(t)=∫0∞‍etxf(x)dx=∫01‍exp[t(∑i=0∞‍viui/a)]du.
We expand the exponential function and use the same algebra that leads to (26)
(29)Mα,θ,a,b(t)=∑r=0∞‍trr!∫01‍(∑i=0∞‍viui/a)rdu=∑r,i=0∞‍trvr,ir!∫01‍ui/adu
and then
(30)Mα,θ,a,b(t)=a∑r,i=0∞‍trvr,i(a+i)r!.
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)δ1(X)=∫0∞‍|x-μ1′|f(x)dx,δ2(X)=∫0∞‍|x-M|f(x)dx,
respectively, where μ1′=E(X) and M=Median(X) denotes the median. Here, M is calculated as the solution of the nonlinear equation I2Φ[(M/θ)α]-1(a,b)=1/2. We define T(q)=∫q∞‍xf(x)dx, which is determined below. The measures δ1(X) and δ2(X) can be written in terms of μ1′ and T(q) as
(32)δ1(X)=2μ1′F(μ1′)-2μ1′+2T(μ1′),δ2(X)=2T(M)-μ1′.
For more details, see Paranaíba et al. [6]. Clearly, F(M) and F(μ1′) are determined from (3). From (10), we have
(33)T(q)=α2π∑r=0∞‍br∫q∞‍(xθ)αe-(1/2)(x/θ)2α{erf[(x/θ)α2]}rdx.
Setting u=(x/θ)α in the last equation gives
(34)T(q)=θ2π∑r=0∞‍br∫(q/θ)α∞‍u1/αe-u2/2[erf(u2)]rdu.
Using the power series for the error function erf(x)=(2/π)∑m=0∞‍((-1)mx2m+1/(2m+1)m!) (see, e.g., [7]), we obtain after some algebra
(35)T(q)=θ2π∑r=0∞‍br(2π)r×∑m1=0∞‍⋯∑mr=0∞‍(-1)m1+⋯+mr(2m1+1)⋯(2mr+1)m1!⋯mr!×Γ[(m1+⋯+mr+r2+12α+12),12(qθ)α],
where Γ(p,x)=∫x∞‍vp-1e-vdv denotes the complementary incomplete gamma function for p>0. The measures δ1(X) and δ2(X) 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(36)B(ρ)=1ρμ1′∫0qxf(x)dx,L(ρ)=1μ1′∫0qxf(x)dx,
respectively, where q=QBGHN(ρ)=QB(ν) and ν=2Φ[(ρ/θ)α]-1 (Section 2) for a given probability ρ. From ∫0q‍xf(x)dx=μ1′-T(ρ), we obtain B(ρ)=ρ-1[1-T(ρ)/μ1′] and L(ρ)=1-T(ρ)/μ1′.


### 3.3. Rényi Entropy

The Rényi information of orderξ for a continuous random variable with density function f(x) is defined as
(37)ȷR(ξ)=11-ξlog[I(ξ)],
where I(ξ)=∫‍fξ(x)dx, ξ>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
(38)ȷR(ξ)=11-ξlog[I(ξ)],
where I(ξ)=∫fξ(x)dx, ξ>0 and ξ≠1. From (4), we have
(39)I(ξ)=2ξ(b-1)(α2/π)ξ[B(a,b)]ξ×∫0∞‍x-ξ(xθ)αξe-(ξ/2)(x/θ)2α{2Φ[(xθ)α]-1}ξ(a-1)×{1-Φ[(xθ)α]}ξ(b-1)dx.
For |z|<1 and b is a real noninteger, the power series holds
(40)(1-z)b-1=∑j=0∞‍(-1)j(b-1j)zj,
where the binomial coefficient is defined for any real. Using (40) in (39) twice, I(ξ) can be expressed as
(41)I(ξ)=2ξ(b-1)(α2/π)ξ[B(a,b)]ξ×∑j,k=0∞‍(-1)j+k2j(ξ(a-1)+1j)×(ξ(b-1)+j+1k)×∫0∞‍x-ξ(xθ)αξe-(ξ/2)(x/θ)2α{Φ[(xθ)α]}kdx.
Substituting Φ(x) by the error function and setting u=(x/θ)α, I(ξ) reduces to
(42)I(ξ)=αξ-1θ1-ξ2ξ(b-1)×(2π)ξ∑j,k=0∞‍∑l=0k‍wj,k,l(a,b)×∫0∞‍u((ξ(α-1)+1)/α)-1e-(ξ/2)u2[erf(u2)]ldu.
Following similar algebra that lead to (35), we obtain
(43)I(ξ)=αξ-1θ1-ξ2[ξ(b-1)+(ξ(α-1)+1)/2α-1]×(2π)ξ∑j,k=0∞‍∑l=0k‍wj,k,l(a,b)×∑m1=0∞‍⋯∑ml=0∞‍(-1)m1+⋯+ml(2m1+1)⋯(2ml+1)m1!⋯ml!×ξ-[m1+⋯+ml+l/2+(ξ(α-1)+1)/2α]×Γ(m1+⋯+ml+l2+ξ(α-1)+12α).
Finally, the Rényi entropy reduces to
(44)ȷR(ξ)=(1-ξ)-1×{(m1+⋯+ml+l2+ξ(α-1)+12α)(ξ-1)log(α)+(1-ξ)log(θ)+[ξ(b-1)+ξ(α-1)+12α-1]log(2)+ξlog(2π)+log[∑j,k=0∞‍∑l=0k‍wj,k,l(a,b)]+log[∑m1=0∞‍⋯∑ml=0∞‍(-1)m1+⋯+ml×((2m1+1)⋯(2ml+1)×m1!⋯ml!)-1∑ml=0∞‍]-[m1+⋯+ml+l2+ξ(α-1)+12α]log(ξ)+log[Γ(m1+⋯+ml+l2+ξ(α-1)+12α)]}.


### 3.4. Reliability

In the context of reliability, the stress-strength model describes the life of a component which has a random strengthX1 that is subjected to a random stress X2. The component fails at the instant that the stress applied to it exceeds the strength, and the component will function satisfactorily whenever X1>X2. Hence, R=Pr(X2<X1) is a measure of component reliability. Here, we derive R when X1 and X2 have independent BGHN(α,θ,a1,b1) and BGHN(α,θ,a2,b2) distributions with the same shape parameters α and θ. The reliability R becomes
(45)R=∫0∞‍f1(x)F2(x)dx,
where the cdf of X2 and the density of X1 are obtained from (6) and (10) as
(46)F2(x)=∑j=0∞‍wj(a2,b2){2Φ[(xθ)α]-1}a2+j,f1(x)=2π(αx)(xθ)αe-(1/2)(x/θ)2α×∑r=0∞‍br(a1,b1){2Φ[(xθ)α]-1}r,
respectively, where
(47)wj(a2,b2)=(-1)jB(a2,b2)(b2-1j),br(a1,b1)=∑j=0∞‍(-1)jB(a1,b1)(b1-1j)sr(a1+j-1)
refer to X2 and X1, respectively. Hence,
(48)R=α2π∑j,r=0∞‍wj(a2,b2)br(a1,b1)×∫0∞‍x-1(xθ)αe-(1/2)(x/θ)2α{2Φ[(xθ)α]-1}a2+j+rdx.Settingu=2Φ[(x/θ)α]-1 in the last integral, the reliability of X reduces to
(49)R=∑j,r=0∞‍wj(a2,b2)br(a1,b1)a2+j+r+1.


## 3.1. Moments and Generating Function

Here, we provide new expressions for the moments and mgf ofX 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)
(25)μs′(α,θ,a,b)=E(Xs)=∫0∞‍xsf(x)dx=∫01‍(∑i=0∞‍viui/a)sdu,
and then using (20), we obtain
(26)μs′(α,θ,a,b)=a∑i=0∞‍vs,i(i+a),
where the quantities vs,i (for i=1,2,…) are easily determined from the recurrence equation
(27)vs,i=(iv0)-1∑m=1i‍[m(s+1)-i]vmvs,i-m,
with vs,0=v0s.We now provide a new alternative representation for the mgf ofX, say Mα,θ,a,b(t)=E(etX), based on the quantile power series (24). We can write
(28)Mα,θ,a,b(t)=∫0∞‍etxf(x)dx=∫01‍exp[t(∑i=0∞‍viui/a)]du.
We expand the exponential function and use the same algebra that leads to (26)
(29)Mα,θ,a,b(t)=∑r=0∞‍trr!∫01‍(∑i=0∞‍viui/a)rdu=∑r,i=0∞‍trvr,ir!∫01‍ui/adu
and then
(30)Mα,θ,a,b(t)=a∑r,i=0∞‍trvr,i(a+i)r!.
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)δ1(X)=∫0∞‍|x-μ1′|f(x)dx,δ2(X)=∫0∞‍|x-M|f(x)dx,
respectively, where μ1′=E(X) and M=Median(X) denotes the median. Here, M is calculated as the solution of the nonlinear equation I2Φ[(M/θ)α]-1(a,b)=1/2. We define T(q)=∫q∞‍xf(x)dx, which is determined below. The measures δ1(X) and δ2(X) can be written in terms of μ1′ and T(q) as
(32)δ1(X)=2μ1′F(μ1′)-2μ1′+2T(μ1′),δ2(X)=2T(M)-μ1′.
For more details, see Paranaíba et al. [6]. Clearly, F(M) and F(μ1′) are determined from (3). From (10), we have
(33)T(q)=α2π∑r=0∞‍br∫q∞‍(xθ)αe-(1/2)(x/θ)2α{erf[(x/θ)α2]}rdx.
Setting u=(x/θ)α in the last equation gives
(34)T(q)=θ2π∑r=0∞‍br∫(q/θ)α∞‍u1/αe-u2/2[erf(u2)]rdu.
Using the power series for the error function erf(x)=(2/π)∑m=0∞‍((-1)mx2m+1/(2m+1)m!) (see, e.g., [7]), we obtain after some algebra
(35)T(q)=θ2π∑r=0∞‍br(2π)r×∑m1=0∞‍⋯∑mr=0∞‍(-1)m1+⋯+mr(2m1+1)⋯(2mr+1)m1!⋯mr!×Γ[(m1+⋯+mr+r2+12α+12),12(qθ)α],
where Γ(p,x)=∫x∞‍vp-1e-vdv denotes the complementary incomplete gamma function for p>0. The measures δ1(X) and δ2(X) 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(36)B(ρ)=1ρμ1′∫0qxf(x)dx,L(ρ)=1μ1′∫0qxf(x)dx,
respectively, where q=QBGHN(ρ)=QB(ν) and ν=2Φ[(ρ/θ)α]-1 (Section 2) for a given probability ρ. From ∫0q‍xf(x)dx=μ1′-T(ρ), we obtain B(ρ)=ρ-1[1-T(ρ)/μ1′] and L(ρ)=1-T(ρ)/μ1′.


## 3.3. Rényi Entropy

The Rényi information of orderξ for a continuous random variable with density function f(x) is defined as
(37)ȷR(ξ)=11-ξlog[I(ξ)],
where I(ξ)=∫‍fξ(x)dx, ξ>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
(38)ȷR(ξ)=11-ξlog[I(ξ)],
where I(ξ)=∫fξ(x)dx, ξ>0 and ξ≠1. From (4), we have
(39)I(ξ)=2ξ(b-1)(α2/π)ξ[B(a,b)]ξ×∫0∞‍x-ξ(xθ)αξe-(ξ/2)(x/θ)2α{2Φ[(xθ)α]-1}ξ(a-1)×{1-Φ[(xθ)α]}ξ(b-1)dx.
For |z|<1 and b is a real noninteger, the power series holds
(40)(1-z)b-1=∑j=0∞‍(-1)j(b-1j)zj,
where the binomial coefficient is defined for any real. Using (40) in (39) twice, I(ξ) can be expressed as
(41)I(ξ)=2ξ(b-1)(α2/π)ξ[B(a,b)]ξ×∑j,k=0∞‍(-1)j+k2j(ξ(a-1)+1j)×(ξ(b-1)+j+1k)×∫0∞‍x-ξ(xθ)αξe-(ξ/2)(x/θ)2α{Φ[(xθ)α]}kdx.
Substituting Φ(x) by the error function and setting u=(x/θ)α, I(ξ) reduces to
(42)I(ξ)=αξ-1θ1-ξ2ξ(b-1)×(2π)ξ∑j,k=0∞‍∑l=0k‍wj,k,l(a,b)×∫0∞‍u((ξ(α-1)+1)/α)-1e-(ξ/2)u2[erf(u2)]ldu.
Following similar algebra that lead to (35), we obtain
(43)I(ξ)=αξ-1θ1-ξ2[ξ(b-1)+(ξ(α-1)+1)/2α-1]×(2π)ξ∑j,k=0∞‍∑l=0k‍wj,k,l(a,b)×∑m1=0∞‍⋯∑ml=0∞‍(-1)m1+⋯+ml(2m1+1)⋯(2ml+1)m1!⋯ml!×ξ-[m1+⋯+ml+l/2+(ξ(α-1)+1)/2α]×Γ(m1+⋯+ml+l2+ξ(α-1)+12α).
Finally, the Rényi entropy reduces to
(44)ȷR(ξ)=(1-ξ)-1×{(m1+⋯+ml+l2+ξ(α-1)+12α)(ξ-1)log(α)+(1-ξ)log(θ)+[ξ(b-1)+ξ(α-1)+12α-1]log(2)+ξlog(2π)+log[∑j,k=0∞‍∑l=0k‍wj,k,l(a,b)]+log[∑m1=0∞‍⋯∑ml=0∞‍(-1)m1+⋯+ml×((2m1+1)⋯(2ml+1)×m1!⋯ml!)-1∑ml=0∞‍]-[m1+⋯+ml+l2+ξ(α-1)+12α]log(ξ)+log[Γ(m1+⋯+ml+l2+ξ(α-1)+12α)]}.


## 3.4. Reliability

In the context of reliability, the stress-strength model describes the life of a component which has a random strengthX1 that is subjected to a random stress X2. The component fails at the instant that the stress applied to it exceeds the strength, and the component will function satisfactorily whenever X1>X2. Hence, R=Pr(X2<X1) is a measure of component reliability. Here, we derive R when X1 and X2 have independent BGHN(α,θ,a1,b1) and BGHN(α,θ,a2,b2) distributions with the same shape parameters α and θ. The reliability R becomes
(45)R=∫0∞‍f1(x)F2(x)dx,
where the cdf of X2 and the density of X1 are obtained from (6) and (10) as
(46)F2(x)=∑j=0∞‍wj(a2,b2){2Φ[(xθ)α]-1}a2+j,f1(x)=2π(αx)(xθ)αe-(1/2)(x/θ)2α×∑r=0∞‍br(a1,b1){2Φ[(xθ)α]-1}r,
respectively, where
(47)wj(a2,b2)=(-1)jB(a2,b2)(b2-1j),br(a1,b1)=∑j=0∞‍(-1)jB(a1,b1)(b1-1j)sr(a1+j-1)
refer to X2 and X1, respectively. Hence,
(48)R=α2π∑j,r=0∞‍wj(a2,b2)br(a1,b1)×∫0∞‍x-1(xθ)αe-(1/2)(x/θ)2α{2Φ[(xθ)α]-1}a2+j+rdx.Settingu=2Φ[(x/θ)α]-1 in the last integral, the reliability of X reduces to
(49)R=∑j,r=0∞‍wj(a2,b2)br(a1,b1)a2+j+r+1.


## 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. LetD1 denote the absolute difference between the integrated version,
(50)1B(a,b)∫02Φ[(x/θ)α]-1‍ta-1(1-t)b-1dt=U,
where U~U(0,1) is a uniform variate on the unit interval, and the truncated version of (24) averaged over (all with step 0.01) x=0.01,…,5, a=0.01,…,10, b=0.01,…,10, α=0.01,…,10 and θ=0.01,…,10. Let D2 denote the absolute difference between integrated version
(51)μs′=2b-1B(a,b)2π∫0∞‍xs(αx)(xθ)αe-(1/2)(x/θ)2α×{2Φ[(xθ)α]-1}a-1×{1-Φ[(xθ)α]}b-1dx
and the truncated version of (26) averaged over x=0.01,…,5, a=0.01,…,10, b=0.01,…,10, α=0.01,…,10 and θ=0.01,…,10. Let D3 denote the absolute difference between the truncated version of (30) and the integrated version,
(52)Mα,θ,a,b(t)=2b-1B(a,b)2π∫0∞‍(αx)(xθ)αexp[tx-12(xθ)2α]×{2Φ[(xθ)α]-1}a-1×{1-Φ[(xθ)α]}b-1dx,
averaged over x=0.01,…,5, a=0.01,…,10, b=0.01,…,10, α=0.01,…,10 and θ=0.01,…,10. Let D4 denote the absolute difference between the truncated version of (35) and the integrated version,
(53)T(q)=α2b-1B(a,b)2π∫0q‍(xθ)α×exp[-12(xθ)2α]{2Φ[(xθ)α]-1}a-1×{1-Φ[(xθ)α]}b-1dx,
averaged over x=0.01,…,5, a=0.01,…,10, b=0.01,…,10, α=0.01,…,10 and θ=0.01,…,10. Let D5 denote the absolute difference between the truncated version of (49) and the integrated version (45), averaged over x=0.01,…,5, a=0.01,…,10, b=0.01,…,10, α=0.01,…,10 and θ=0.01,…,10.We obtain the following estimates after extensive computations:D1=2.31×10-20, D2=8.47×10-18, D3=1.22×10-21, D4=1.51×10-22 and D5=9.41×10-19. 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

SupposeX1,…,Xn is a random sample of size n from a continuous distribution and let X1:n<X2:n<⋯<Xn:n denote the corresponding order statistics. There has been a large amount of work relating to moments of order statistics Xi:n. 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 Xi:n is given by
(54)fi:n(x)=f(x)B(i,n-i+1)F(x)i-1[1-F(x)]n-i.For the BGHN distribution, Pescim et al. [3] obtained
(55)fi:n(x)=∑k=0n-i∑m1=0∞‍⋯∑mi+k-1=0∞ηk,Mkfk,Mk(x),
where Mk denotes a sequence (m1,…,mi+k-1) of i+k-1 nonnegative integers, fk,Mk(x) denotes the BGHN(α,θ,(i+k)a+∑j=1i+k-1‍mj,b) density function defined under Mk and the constants ηk,Mk are given by
(56)ηk,Mk=(-1)k+∑j=1i+k-1mj(n-ik)B(a(i+k)+∑j=1i+k-1mj,b)Γ(b)i+k-1B(a,b)i+kB(i,n-i+1)∏j=1i+k-1Γ(b-mj)mj!(a+mj).
The quantities ηk,Mk are easily obtained given k and a sequence Mk of indices m1,…,mi+k-1. The sums in (55) extend over all (i+k)-tuples (k,m1,…,mi+k-1) and can be implementable on a computer. If b>0 is an integer, (55) holds but the indices m1,…,mi+k-1 vary from zero to b-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
(57)E(Xi:ns)=∑k=0n-i∑m1=0∞‍⋯∑mi+k-1=0∞ηk,MkE(Xi,ks),
where Xi,k~BGHN(α,θ,(i+k)a+∑j=1i+k-1‍mj,b) and E(Xi,ks) can be determined from (26). Inserting (26) in (57) and changing indices, we can write
(58)E(Xi:ns)=∑k=0n-i∑m1=0∞‍⋯∑mi+k-1=0∞ηk,Mk∑q=0∞‍a⋆vs,q(q+a⋆),
where
(59)a⋆=(i+k)a+∑j=1i+k-1mj.The momentsE(Xi:ns) 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 a=2, b=3 and n=5, Table 1 lists the numerical values for the moments E(Xi:ns), (s=1,…,4; i=1,…,5) and for the variance, skewness and kurtosis.Table 1
Moments of the BGHN order statistics forθ=8.5, α=3, a=2, and b=3.
Order statistics→
X
1
:
5
X
2
:
5
X
3
:
5
X
4
:
5
X
5
:
5
s
=
1
2.67942
5.81946
4.73975
1.71571
0.23289
s
=
2
18.10007
39.31172
32.01807
11.59006
1.57328
s
=
3
127.8683
277.7184
226.1923
81.87821
11.11451
s
=
4
938.314
2037.933
1659.828
600.8327
81.55966
Variance
10.92078
5.44561
9.55281
8.64636
1.51904
Skewness
0.57767
−1.13596
−0.54601
1.27135
5.36291
Kurtosis
1.61751
4.06406
1.89737
2.91095
31.07642


### 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
(60)Mi:n(t)=∑k=0n-i‍∑m1=0∞‍⋯∑mi+k-1=0∞ηk,MkMα,θ,a⋆,b(t),
where Mα,θ,a⋆,b(t) is the mgf of the BGHN(α,θ,a⋆,b) 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)δ1(Xi:n)=∫0∞|x-μi|fi:n(x)dx,δ2(Xi:n)=∫0∞|x-Mi|fi:n(x)dx,
respectively, where μi=E(Xi:n) and Mi=Median(Xi:n) denotes the median. Here, Mi is obtained as the solution of the nonlinear equation
(62)∑r=in(nr){I2Φ[(Mi/θ)α]-1(a,b)}r×{1-I2Φ[(Mi/θ)α]-1(a,b)}n-r=12.
The measures δ1(Xi:n) and δ2(Xi:n) follow from
(63)δ1(Xi:n)=2μiFi:n(μi)-2μi+2Ji(μi),δ2(Xi:n)=2Ji(Mi)-μi,
where Ji(q)=∫q∞xfi:n(x)dx. Using (55), we have
(64)Ji(q)=∑k=0n-i‍∑m1=0∞‍⋯∑mi+k-1=0∞ηi,MkT(q),
where T(q) is obtained from (35) using a⋆ given by (59) and
(65)br⋆=br(a⋆,b)=1B(a⋆,b)∑j=0∞‍(-1)j(b-1j)sr(a⋆+j-1),
where
(66)sr(a⋆+j-1)=∑k=r∞‍(-1)r+ka⋆+j(a⋆+j-1k)(kr).Bonferroni and Lorenz curves of the order statistics are given by(67)Bi:n(ρ)=1ρμi∫0qxfi:n(x)dx,Li:n(ρ)=1μi∫0qxfi:n(x)dx,
respectively, where q=Fi:n-1(ρ) for a given probability ρ. From ∫0qxfi:n(x)dx=μi-Ji(ρ), we obtain
(68)Bi:n(ρ)=1ρ-Ji(ρ)ρμi,Li:n(ρ)=1-Ji(ρ)μi.


### 5.5. Rényi Entropy

The Rényi entropy of the order statistics is defined by(69)ȷR(λ)=11-λlog[H(λ)],
where H(λ)=∫‍fi:nλ(x)dx, λ>0 and λ≠1. From (54), it follows that
(70)H(λ)=1[B(i,n-i+1)]λ∫0∞‍fλ(x)×[F(x)]λ(i-1)[1-F(x)]λ(n-i)dx.
Using (40) in (70), we obtain
(71)H(λ)=1[B(i,n-i+1)]λ∑k1=0∞‍(-1)k1(λ(i-1)k1)×∫0∞‍fλ(x)[F(x)]λ(i-1)+k1dx.
For λ>0 and λ≠1, we can expand [F(x)]λ(i-1)+k1 as
(72)F(x)λ(i-1)+k1={1-[1-F(x)]}λ(i-1)+k1=∑p1=0∞(-1)p1(λ(i-1)+k1+1p1)[1-F(x)]p1,
and then
(73)F(x)λ(i-1)+k1=∑p1=0∞‍∑l1=0p1‍(-1)p1+l1×(λ(i-1)+k1+1p1)(p1l1)F(x)l1.
Hence, from (70), we can write
(74)H(λ)=1[B(i,n-i+1)]λ×∑k1,p1=0∞‍∑l1=0p1‍(-1)k1+p1+l1×(λ(i-1)k1)(λ(i-1)+k1+1p1)(p1l1)×∫0∞‍fλ(x)F(x)l1dx.
By replacing f(x) and F(x) by (10) and (16) given by Pescim et al. [3], we obtain
(75)H(λ)=2λ(b-1)[α(2/π)]λ[B(i,n-i+1)]λ×∑k1,p1=0∞‍∑l1=0p1‍[Γ(b)]l1(-1)k1+p1+l1[B(a,b)]l1×(λ(i-1)k1)(λ(i-1)+k1+1p1)(p1l1)×∫0∞‍x-λ(xθ)αλe-(λ/2)(x/θ)2α×{2Φ[(xθ)α]-1}λ(a-1)×{1-Φ[(xθ)α]}λ(b-1)×[∑j=0∞‍(-1)j{2Φ[(x/θ)α]-1}jΓ(b-j)j!(a+j)]l1dx.Using the identity(∑i=0∞‍ai)k=∑{m1,…,mk}=0∞‍am1⋯amk (for k positive integer) in (4), we have
(76)H(λ)=2λ(b-1)[α(2/π)]λ[B(i,n-i+1)]λ×∑k1,p1=0∞‍∑l1=0p1ck1,p1,l1(a,b)×∑m1=0∞‍⋯∑ml1=0∞‍(-1)m1+⋯+ml1∏j=1l1Γ(b-mj)mj!(a+mj)×∫0∞x-λ(xθ)αλe-(λ/2)(x/θ)2α×{2Φ[(xθ)α]-1}a(λ+l1)-λ+∑j=1l1mj×{1-Φ[(xθ)α]}λ(b-1)dx,
where
(77)ck1,p1,l1(a,b)=(-1)k1+p1+l1[Γ(b)]l1[B(a,b)]l1×(λ(i-1)k1)(λ(i-1)+k1+1p1)(p1l1).
Using the power series (40) in (76) twice and replacing Φ(x) by the error function erf(x) (defined in Section 3.2) we obtain after some algebra
(78)H(λ)=αλ-1θ1-λ2[λ(b-1)+(λ(α-1)+1)/2α-1](2/π)λ[B(i,n-i+1)]λ×∑k1,p1=0∞∑l1=0p1‍ck1,p1,l1(a,b)×∑m1=0∞‍⋯∑ml1=0∞‍(-1)m1+⋯+ml1∏j=1l1Γ(b-mj)mj!(a+mj)×∑r1,s1=0∞∑v1=0s1‍dr1,s1,v1×∑m1=0∞‍⋯∑ms1=0∞‍(-1)m1+⋯+ms1(2m1+1)⋯(2ms1+1)m1!⋯ms1!×λ-[m1+⋯+ms1+s1/2+(λ(α-1)+1)/2α]×Γ(m1+⋯+ms1+s12+λ(α-1)+12α),
where the quantity dr,s,t 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(79)ȷR(λ)=(1-λ)-1×{∑k1,p1=0∞‍(λ-1)log(α)+(1-λ)log(θ)+[λ(b-1)+λ(α-1)+12α-1]log(2)+λlog(2π)-λlog[B(i,n-i+1)]+log[∑k1,p1=0∞‍∑l1=0p1‍ck1,p1,l1(a,b)]+log[∑m1=0∞‍⋯∑ml1=0∞‍(-1)m1+⋯+ml1∏j=1l1Γ(b-mj)mj!(a+mj)]+log(∑r1,s1=0∞‍∑v1=0s1‍dr1,s1,v1)+log[∑m1=0∞‍⋯∑ms1=0∞‍((-1)m1+⋯+ms1)×(ms1!(2m1+1)⋯(2ms1+1)×m1!⋯ms1!)-1∑ms1=0∞‍]-[m1+⋯+ms1+s12+λ(α-1)+12α]log(λ)+log[Γ(m1+⋯+ms1+s12+λ(α-1)+12α)]∑k1,p1=0∞‍}.
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.


## 5.1. Mixture Form

SupposeX1,…,Xn is a random sample of size n from a continuous distribution and let X1:n<X2:n<⋯<Xn:n denote the corresponding order statistics. There has been a large amount of work relating to moments of order statistics Xi:n. 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 Xi:n is given by
(54)fi:n(x)=f(x)B(i,n-i+1)F(x)i-1[1-F(x)]n-i.For the BGHN distribution, Pescim et al. [3] obtained
(55)fi:n(x)=∑k=0n-i∑m1=0∞‍⋯∑mi+k-1=0∞ηk,Mkfk,Mk(x),
where Mk denotes a sequence (m1,…,mi+k-1) of i+k-1 nonnegative integers, fk,Mk(x) denotes the BGHN(α,θ,(i+k)a+∑j=1i+k-1‍mj,b) density function defined under Mk and the constants ηk,Mk are given by
(56)ηk,Mk=(-1)k+∑j=1i+k-1mj(n-ik)B(a(i+k)+∑j=1i+k-1mj,b)Γ(b)i+k-1B(a,b)i+kB(i,n-i+1)∏j=1i+k-1Γ(b-mj)mj!(a+mj).
The quantities ηk,Mk are easily obtained given k and a sequence Mk of indices m1,…,mi+k-1. The sums in (55) extend over all (i+k)-tuples (k,m1,…,mi+k-1) and can be implementable on a computer. If b>0 is an integer, (55) holds but the indices m1,…,mi+k-1 vary from zero to b-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
(57)E(Xi:ns)=∑k=0n-i∑m1=0∞‍⋯∑mi+k-1=0∞ηk,MkE(Xi,ks),
where Xi,k~BGHN(α,θ,(i+k)a+∑j=1i+k-1‍mj,b) and E(Xi,ks) can be determined from (26). Inserting (26) in (57) and changing indices, we can write
(58)E(Xi:ns)=∑k=0n-i∑m1=0∞‍⋯∑mi+k-1=0∞ηk,Mk∑q=0∞‍a⋆vs,q(q+a⋆),
where
(59)a⋆=(i+k)a+∑j=1i+k-1mj.The momentsE(Xi:ns) 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 a=2, b=3 and n=5, Table 1 lists the numerical values for the moments E(Xi:ns), (s=1,…,4; i=1,…,5) and for the variance, skewness and kurtosis.Table 1
Moments of the BGHN order statistics forθ=8.5, α=3, a=2, and b=3.
Order statistics→
X
1
:
5
X
2
:
5
X
3
:
5
X
4
:
5
X
5
:
5
s
=
1
2.67942
5.81946
4.73975
1.71571
0.23289
s
=
2
18.10007
39.31172
32.01807
11.59006
1.57328
s
=
3
127.8683
277.7184
226.1923
81.87821
11.11451
s
=
4
938.314
2037.933
1659.828
600.8327
81.55966
Variance
10.92078
5.44561
9.55281
8.64636
1.51904
Skewness
0.57767
−1.13596
−0.54601
1.27135
5.36291
Kurtosis
1.61751
4.06406
1.89737
2.91095
31.07642


## 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
(60)Mi:n(t)=∑k=0n-i‍∑m1=0∞‍⋯∑mi+k-1=0∞ηk,MkMα,θ,a⋆,b(t),
where Mα,θ,a⋆,b(t) is the mgf of the BGHN(α,θ,a⋆,b) 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)δ1(Xi:n)=∫0∞|x-μi|fi:n(x)dx,δ2(Xi:n)=∫0∞|x-Mi|fi:n(x)dx,
respectively, where μi=E(Xi:n) and Mi=Median(Xi:n) denotes the median. Here, Mi is obtained as the solution of the nonlinear equation
(62)∑r=in(nr){I2Φ[(Mi/θ)α]-1(a,b)}r×{1-I2Φ[(Mi/θ)α]-1(a,b)}n-r=12.
The measures δ1(Xi:n) and δ2(Xi:n) follow from
(63)δ1(Xi:n)=2μiFi:n(μi)-2μi+2Ji(μi),δ2(Xi:n)=2Ji(Mi)-μi,
where Ji(q)=∫q∞xfi:n(x)dx. Using (55), we have
(64)Ji(q)=∑k=0n-i‍∑m1=0∞‍⋯∑mi+k-1=0∞ηi,MkT(q),
where T(q) is obtained from (35) using a⋆ given by (59) and
(65)br⋆=br(a⋆,b)=1B(a⋆,b)∑j=0∞‍(-1)j(b-1j)sr(a⋆+j-1),
where
(66)sr(a⋆+j-1)=∑k=r∞‍(-1)r+ka⋆+j(a⋆+j-1k)(kr).Bonferroni and Lorenz curves of the order statistics are given by(67)Bi:n(ρ)=1ρμi∫0qxfi:n(x)dx,Li:n(ρ)=1μi∫0qxfi:n(x)dx,
respectively, where q=Fi:n-1(ρ) for a given probability ρ. From ∫0qxfi:n(x)dx=μi-Ji(ρ), we obtain
(68)Bi:n(ρ)=1ρ-Ji(ρ)ρμi,Li:n(ρ)=1-Ji(ρ)μi.


## 5.5. Rényi Entropy

The Rényi entropy of the order statistics is defined by(69)ȷR(λ)=11-λlog[H(λ)],
where H(λ)=∫‍fi:nλ(x)dx, λ>0 and λ≠1. From (54), it follows that
(70)H(λ)=1[B(i,n-i+1)]λ∫0∞‍fλ(x)×[F(x)]λ(i-1)[1-F(x)]λ(n-i)dx.
Using (40) in (70), we obtain
(71)H(λ)=1[B(i,n-i+1)]λ∑k1=0∞‍(-1)k1(λ(i-1)k1)×∫0∞‍fλ(x)[F(x)]λ(i-1)+k1dx.
For λ>0 and λ≠1, we can expand [F(x)]λ(i-1)+k1 as
(72)F(x)λ(i-1)+k1={1-[1-F(x)]}λ(i-1)+k1=∑p1=0∞(-1)p1(λ(i-1)+k1+1p1)[1-F(x)]p1,
and then
(73)F(x)λ(i-1)+k1=∑p1=0∞‍∑l1=0p1‍(-1)p1+l1×(λ(i-1)+k1+1p1)(p1l1)F(x)l1.
Hence, from (70), we can write
(74)H(λ)=1[B(i,n-i+1)]λ×∑k1,p1=0∞‍∑l1=0p1‍(-1)k1+p1+l1×(λ(i-1)k1)(λ(i-1)+k1+1p1)(p1l1)×∫0∞‍fλ(x)F(x)l1dx.
By replacing f(x) and F(x) by (10) and (16) given by Pescim et al. [3], we obtain
(75)H(λ)=2λ(b-1)[α(2/π)]λ[B(i,n-i+1)]λ×∑k1,p1=0∞‍∑l1=0p1‍[Γ(b)]l1(-1)k1+p1+l1[B(a,b)]l1×(λ(i-1)k1)(λ(i-1)+k1+1p1)(p1l1)×∫0∞‍x-λ(xθ)αλe-(λ/2)(x/θ)2α×{2Φ[(xθ)α]-1}λ(a-1)×{1-Φ[(xθ)α]}λ(b-1)×[∑j=0∞‍(-1)j{2Φ[(x/θ)α]-1}jΓ(b-j)j!(a+j)]l1dx.Using the identity(∑i=0∞‍ai)k=∑{m1,…,mk}=0∞‍am1⋯amk (for k positive integer) in (4), we have
(76)H(λ)=2λ(b-1)[α(2/π)]λ[B(i,n-i+1)]λ×∑k1,p1=0∞‍∑l1=0p1ck1,p1,l1(a,b)×∑m1=0∞‍⋯∑ml1=0∞‍(-1)m1+⋯+ml1∏j=1l1Γ(b-mj)mj!(a+mj)×∫0∞x-λ(xθ)αλe-(λ/2)(x/θ)2α×{2Φ[(xθ)α]-1}a(λ+l1)-λ+∑j=1l1mj×{1-Φ[(xθ)α]}λ(b-1)dx,
where
(77)ck1,p1,l1(a,b)=(-1)k1+p1+l1[Γ(b)]l1[B(a,b)]l1×(λ(i-1)k1)(λ(i-1)+k1+1p1)(p1l1).
Using the power series (40) in (76) twice and replacing Φ(x) by the error function erf(x) (defined in Section 3.2) we obtain after some algebra
(78)H(λ)=αλ-1θ1-λ2[λ(b-1)+(λ(α-1)+1)/2α-1](2/π)λ[B(i,n-i+1)]λ×∑k1,p1=0∞∑l1=0p1‍ck1,p1,l1(a,b)×∑m1=0∞‍⋯∑ml1=0∞‍(-1)m1+⋯+ml1∏j=1l1Γ(b-mj)mj!(a+mj)×∑r1,s1=0∞∑v1=0s1‍dr1,s1,v1×∑m1=0∞‍⋯∑ms1=0∞‍(-1)m1+⋯+ms1(2m1+1)⋯(2ms1+1)m1!⋯ms1!×λ-[m1+⋯+ms1+s1/2+(λ(α-1)+1)/2α]×Γ(m1+⋯+ms1+s12+λ(α-1)+12α),
where the quantity dr,s,t 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(79)ȷR(λ)=(1-λ)-1×{∑k1,p1=0∞‍(λ-1)log(α)+(1-λ)log(θ)+[λ(b-1)+λ(α-1)+12α-1]log(2)+λlog(2π)-λlog[B(i,n-i+1)]+log[∑k1,p1=0∞‍∑l1=0p1‍ck1,p1,l1(a,b)]+log[∑m1=0∞‍⋯∑ml1=0∞‍(-1)m1+⋯+ml1∏j=1l1Γ(b-mj)mj!(a+mj)]+log(∑r1,s1=0∞‍∑v1=0s1‍dr1,s1,v1)+log[∑m1=0∞‍⋯∑ms1=0∞‍((-1)m1+⋯+ms1)×(ms1!(2m1+1)⋯(2ms1+1)×m1!⋯ms1!)-1∑ms1=0∞‍]-[m1+⋯+ms1+s12+λ(α-1)+12α]log(λ)+log[Γ(m1+⋯+ms1+s12+λ(α-1)+12α)]∑k1,p1=0∞‍}.
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

LetXi be a random variable having density function (4), where ω=(α,θ,a,b)T 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 i is assumed to have a lifetime Xi and a censoring time Ci, where Xi and Ci are independent random variables. Suppose that the data consist of n independent observations xi=min(Xi,Ci) for i=1,…,n. The distribution of Ci does not depend on any of the unknown parameters of the distribution of Xi. Parametric inference for such data are usually based on likelihood methods and their asymptotic theory. The censored log-likelihood l(ω) for the model parameters is
(80)l(ω)=rlog(2π)+∑i∈F‍log(αxi)+α∑i∈F‍log(xiθ)-12∑i∈F‍(xiθ)2α+r(b-1)log(2)-rlog[B(a,b)]+(a-1)∑i∈F‍log{2Φ[(xiθ)α]-1}+(b-1)∑i∈F‍log{1-Φ[(xiθ)α]}+∑i∈C‍log{1-I2Φ[(xi/θ)α]-1(a,b)},
where r is the number of failures and F and C denote the uncensored and censored sets of observations, respectively.The score functions for the parametersα, θ, a, and b are given by
(81)Uα(ω)=r2+∑i∈F‍log(xiθ)[1-(xiθ)2α]+2(a-1)2π∑i∈F‍vilog(xi/θ)P(xi)+(1-b)π∑i∈F‍vilog(xi/θ)D(xi)-2b-1∑i∈C‍vilog(xi/θ)[P(xi)]a-1[D(xi)]b-1[1-Q(xi)],Uθ(ω)=-r(αθ)+(αθ)∑i∈F‍(xiθ)2α+2(a-1)2π×∑i∈F‍vi(α/θ)P(xi)+(1-b)π∑i∈F‍vi(α/θ)D(xi)-2b-1∑i∈C‍vi(α/θ)[P(xi)]a-1[D(xi)]b-1[1-Q(xi)],Ua(ω)=r[ψ(a+b)-ψ(a)]+∑i∈F‍log[P(xi)]-∑i∈C‍I˙P(xi)(a,b)|a[1-Q(xi)],Ub(ω)=rlog(2)+r[ψ(a+b)-ψ(b)]+∑i∈F‍log{D(xi)}-∑i∈C‍I˙P(xi)(a,b)|b[1-Q(xi)],
where
(82)vi=exp[-12(xiθ)2α](xiθ)α,P(xi)=2Φ[(xiθ)α]-1,D(xi)=1-Φ[(xiθ)α],Q(xi)=I2Φ[(xi/θ)α]-1(a,b),I˙P(xi)(a,b)|a=∂∂a[1B(a,b)∫0P(xi)‍wa-1(1-w)b-1dw],I˙P(xi)(a,b)|b=∂∂b[1B(a,b)∫0P(xi)‍wa-1(1-w)b-1dw]
and ψ(·) is the digamma function.The maximum likelihood estimate (MLE)ω^ of ω is obtained by solving numerically the nonlinear equations Uα(ω)=0, Uθ(ω)=0, Ua(ω)=0, and Ub(ω)=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(α,θ,a,b) 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 J(ω) is
(83)J(ω)=-(LααLαθLαaLαb·LθθLθaLθb··LaaLab···Lbb),
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 ofn(ω^-ω) is N4(0,K(ω)-1), where K(ω) is the expected information matrix. The normal approximation is valid if K(ω) is replaced by J(ω^), that is, the observed information matrix evaluated at ω^. The multivariate normal N4(0,J(ω^)-1) 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 typeH0:ω=ω0 versus H:ω≠ω0 can be performed using LR statistics. For example, the test of H0:b=1 versus H:b≠1 is equivalent to compare the EGHN and BGHN distributions. In this case, the LR statistic becomes w=2{l(α^,θ^,a^,b^)-l(α~,θ~,a~,1)}, where α^, θ^, a^, and b^ are the MLEs under H and α~, θ~ and a~ are the estimates under H0. 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 p, or has a proper survival function S(x) with probability 1-p. This gives an improper population survivor function S*(x) in the form of the mixture, namely,
(84)S*(x)=p+(1-p)S(x),S(∞)=0,S*(∞)=p.
Common choices for S(x) 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 x1,…,xn, where xi is either the observed lifetime or censoring time for the ith individual. Let a binary random variable zi, for i=1,…,n, indicating that the ith individual in a population is at risk or not with respect to a certain type of failure, that is, zi=1 indicates that the ith individual will eventually experience a failure event (uncured) and zi=0 indicates that the individual will never experience such event (cured).For an individual, the proportion of uncured1-p can be specified such that the conditional distribution of z is given by Pr(zi=1)=1-p. Suppose that the Xi’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 atxi to the likelihood function is given by
(85)(1-p)2b-12/π(α/xi)(xi/θ)αexp[-(1/2)(xi/θ)2α]B(a,b)×{2Φ[(xiθ)α]-1}a-1{1-Φ[(xiθ)α]}b-1
and the contribution of an individual that is at risk at time xi is
(86)p+(1-p){1-I2Φ[(xi/θ)α]-1(a,b)}.
The new model defined by (85) and (86) is called the BGHN mixture model with long-term survivors. For a=b=1, we obtain the GHN mixture model (a new model) with long-term survivors.Thus, the log-likelihood function for the parameter vectorθ=(p,α,θ,a,b)T follows from (85) and (86) as
(87)l(θ)=rlog[(1-p)2b-12/πB(a,b)]+∑i∈F‍log(αxi)+α∑i∈F‍(xiθ)-12∑i∈F‍log[(xiθ)2α]+(a-1)∑i∈F‍log{2Φ[(xiθ)α]-1}+(b-1)∑i∈F‍log{1-Φ[(xiθ)α]}+∑i∈C‍log{p+(1-p)[1-I2Φ[(xi/θ)α]-1(a,b)]},
where F and C denote the sets of individuals corresponding to lifetime observations and censoring times, respectively, and r 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 ofn=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 α, θ, a, and b 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 H0: a=b=1 versus H1: H0  is not true, that is, to compare the BGHN and GHN models, is w=2{-120.1-(-126.0)}=11.8 (P value =0.0027), which gives favorable indication towards the BGHN model. In special, the Weibull density function is given by
(88)f(x)=γλγxγ-1exp[-(xλ)γ],x,λ,γ>0.
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.
Model
α
θ
a
b
AIC
CAIC
BIC
BGHN
0.0479 (0.0062)
19.3753 (2.8892)
400.87 (0.4000)
192.13 (0.2480)
248.2
249.6
254.3
GHN
0.9223 (0.1361)
25.7408 (3.8627)
1 (—)
1 (—)
256.0
256.4
259.1
λ
γ
Weibull
20.7819 (2.8215)
1.3414 (0.1721)
263.4
267.8
270.4Estimated 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.
(a)
(b)
(c)


### 8.2. Censored Data: Radiotherapy

The data refer to the survival time (days) for cancer patients (n=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 H0: a=b=1 versus H1: H0  is not true, that is, to compare the BGHN and GHN models, is w=2{-293.3-(-299.4)}=12.2 (P 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.
Model
α
θ
a
b
AIC
CAIC
BIC
BGHN
0.0456 (0.0051)
540.40 (108.19)
223.66 (0.4313)
112.38 (0.1560)
594.6
595.5
602.3
GHN
0.7074 (0.0879)
533.45 (88.6917)
1 (—)
1 (—)
602.8
603.1
606.7
λ
γ
Weibull
361.76 (52.4678)
1.0239 (0.1062)
705.7
706.0
709.6


### 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 H0: a=b=1 versus H1: H0  is not true, that is, to compare the BGHN and GHN models, is w=2{-196.1-(-217.2)}=42.1 (P 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.
Model
α
θ
a
b
AIC
CAIC
BIC
BGHN
7.0649 (0.0027)
18.9581 (0.0027)
0.2368 (0.0416)
0.1015 (0.0134)
400.2
401.6
409.3
GHN
3.8208 (0.4150)
21.8838 (0.4969)
1 (—)
1 (—)
438.4
439.0
442.9
λ
γ
Weibull
21.1972 (0.6034)
4.2716 (0.4627)
455.5
456.1
460.1


### 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 T is defined as the time until the patient’s death. The original sample size was n=427 patients, 10 of whom did not present a value for explanatory variable tumor thickness. When such cases were removed, a sample of size n=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 H0: a=b=1 versus H1: H0  is not true, that is, to compare the BGHN and GHN mixture models, becomes w=2{-524.75-(-534.05)}=18.6 (P 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 (pBGHN=0.4871) is more appropriate than that one estimated by the GHN mixture model (pGHN=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.
Model
p
α
θ
a
b
AIC
CAIC
BIC
BGHN Mixture
0.4871 (0.0349)
0.2579 (0.0197)
54.8258 (17.2126)
19.4970 (1.0753)
40.9300 (2.2931)
1059.5
1059.6
1079.6
GHN Mixture
0.5150 (0.0279)
1.2553 (0.0799)
2.5090 (0.1475)
1 (—)
1 (—)
1074.1
1074.2
1086.2Figure 2
Estimated survival functions for the BGHN and GHN mixture models and the empirical survival for the melanoma data.


## 8.1. Censored Data: Transistors

The first data set consists of the lifetimes ofn=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 α, θ, a, and b 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 H0: a=b=1 versus H1: H0  is not true, that is, to compare the BGHN and GHN models, is w=2{-120.1-(-126.0)}=11.8 (P value =0.0027), which gives favorable indication towards the BGHN model. In special, the Weibull density function is given by
(88)f(x)=γλγxγ-1exp[-(xλ)γ],x,λ,γ>0.
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.
Model
α
θ
a
b
AIC
CAIC
BIC
BGHN
0.0479 (0.0062)
19.3753 (2.8892)
400.87 (0.4000)
192.13 (0.2480)
248.2
249.6
254.3
GHN
0.9223 (0.1361)
25.7408 (3.8627)
1 (—)
1 (—)
256.0
256.4
259.1
λ
γ
Weibull
20.7819 (2.8215)
1.3414 (0.1721)
263.4
267.8
270.4Estimated 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.
(a)
(b)
(c)


## 8.2. Censored Data: Radiotherapy

The data refer to the survival time (days) for cancer patients (n=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 H0: a=b=1 versus H1: H0  is not true, that is, to compare the BGHN and GHN models, is w=2{-293.3-(-299.4)}=12.2 (P 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.
Model
α
θ
a
b
AIC
CAIC
BIC
BGHN
0.0456 (0.0051)
540.40 (108.19)
223.66 (0.4313)
112.38 (0.1560)
594.6
595.5
602.3
GHN
0.7074 (0.0879)
533.45 (88.6917)
1 (—)
1 (—)
602.8
603.1
606.7
λ
γ
Weibull
361.76 (52.4678)
1.0239 (0.1062)
705.7
706.0
709.6


## 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 H0: a=b=1 versus H1: H0  is not true, that is, to compare the BGHN and GHN models, is w=2{-196.1-(-217.2)}=42.1 (P 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.
Model
α
θ
a
b
AIC
CAIC
BIC
BGHN
7.0649 (0.0027)
18.9581 (0.0027)
0.2368 (0.0416)
0.1015 (0.0134)
400.2
401.6
409.3
GHN
3.8208 (0.4150)
21.8838 (0.4969)
1 (—)
1 (—)
438.4
439.0
442.9
λ
γ
Weibull
21.1972 (0.6034)
4.2716 (0.4627)
455.5
456.1
460.1


## 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 T is defined as the time until the patient’s death. The original sample size was n=427 patients, 10 of whom did not present a value for explanatory variable tumor thickness. When such cases were removed, a sample of size n=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 H0: a=b=1 versus H1: H0  is not true, that is, to compare the BGHN and GHN mixture models, becomes w=2{-524.75-(-534.05)}=18.6 (P 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 (pBGHN=0.4871) is more appropriate than that one estimated by the GHN mixture model (pGHN=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.
Model
p
α
θ
a
b
AIC
CAIC
BIC
BGHN Mixture
0.4871 (0.0349)
0.2579 (0.0197)
54.8258 (17.2126)
19.4970 (1.0753)
40.9300 (2.2931)
1059.5
1059.6
1079.6
GHN Mixture
0.5150 (0.0279)
1.2553 (0.0799)
2.5090 (0.1475)
1 (—)
1 (—)
1074.1
1074.2
1086.2Figure 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 thebeta 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.



---

*Source: 491628-2013-12-31.xml*