Statistical Inferences and Applications of the Half Exponential Power Distribution

Abstract

We investigate the statistical inferences and applications of the half exponential power distribution for the first time. The proposed model defined on the nonnegative reals extends the half normal distribution and is more flexible. The characterizations and properties involving moments and some measures based on moments of this distribution are derived. The inference aspects using methods of moment and maximum likelihood are presented. We also study the performance of the estimators using the Monte Carlo simulation. Finally, we illustrate it with two real applications.

1. Introduction

The well-known exponential power (EP) distribution or the generalized normal distribution has the following density function: 𝑓(𝑥)=𝑝^1−1/𝑝𝑒2Γ1/𝑝^{−|𝑥|𝑝/𝑝},−∞<𝑥<∞, (1) where 𝑝>0 is the shape parameter. This family consists of a wide range of symmetric distributions and allows continuous variation from normality to nonnormality. It includes the normal distribution 𝑍∼𝑁(0,1) as the special case when 𝑝=2 and the Laplace distribution when 𝑝=1. Nadarajah [1] provided a comprehensive treatment of its mathematical properties.

Its tails can be more platykurtic (𝑝>2) or more leptokurtic (𝑝<2) than the normal distribution (𝑝=2). The distribution has been widely used in the Bayes analysis and robustness studies (see Box and Tiao [2], Genc [3], Goodman and Kotz [4], and Tiao and Lund [5].)

On the other hand, since the most popular models used to describe the lifetime process are defined on nonnegative measurements, which motivate us to take a positive truncation in the model (1) and develop a half exponential power (HEP) distribution. As far as we know, this model has not been previously studied although, we believe, it plays an important role in data analysis. The resulting nonnegative half exponential power distribution generalizes the half normal (HN) distribution, and it is more flexible. In our work, we aim to investigate the statistical features of the nonnegative model and apply them to fit the lifetime data.

The rest of this paper is organized as follows: in Section 2, we present the new distribution and study its properties. Section 3 discusses the inference, moments, and maximum likelihood estimation for the parameters. In Section 4, we discuss a useful technique, a half normal plot with a simulated envelope, to assess the model adequacy. Simulation studies are performed in Section 5. Section 6 gives two illustrative examples and reports the results. Section 7 concludes our work.

2. The Half Exponential Power Distribution

2.1. The Density and Hazard Function

Definition 1. A random variable 𝑋 has a half exponential power slash distribution if its density function with scale parameter 𝜎>0 takes 𝑓(𝑥)=𝑝^1−1/𝑝𝑒𝜎Γ1/𝑝^{−𝑥𝑝/𝑝𝜎𝑝},𝑥≥0, (2) where 𝜎>0 and 𝑝>0. We denote it as 𝑋∼HEP(𝜎,𝑝).
Figure 1(a) displays some plots of the density function of the half exponential power distribution with various parameters.
The cumulative distribution function of the half exponential power distribution 𝑋∼HEP(𝜎,𝑝) is given as follows. For 𝑥≥0, 𝐹(𝑥)=^𝑥₀𝑓_𝑋(𝑢)𝑑𝑢=^𝑥₀𝑝^1−1/𝑝𝑒𝜎Γ1/𝑝^{−𝑢𝑝/𝑝𝜎𝑝}=𝛾𝑑𝑢1/𝑝,𝑥^𝑝/𝑝𝜎^𝑝Γ,1/𝑝 (3) where 𝛾(,) is the lower incomplete gamma function, defined as 𝛾(𝑠,𝑥)=∫^𝑥₀𝑡^𝑠−1𝑒^−𝑡𝑑𝑡.
The hazard rate function (also known as the failure rate function) of the half exponential power distribution is given by, for 𝑥≥0, ℎ(𝑥)=𝑓(𝑥)(𝑥)=𝑝1−𝐹^1−1/𝑝𝑒^{−𝑥𝑝/𝑝𝜎𝑝}𝜎Γ1/𝑝−𝛾1/𝑝,𝑥^𝑝/𝑝𝜎^𝑝. (4)
Since Γ(𝑠)−𝛾(𝑠,𝑥)∼𝑥^𝑠−1𝑒^−𝑥, as 𝑥→∞, we obtain ℎ(𝑥)∼𝑥^𝑝−1/𝜎^𝑝. Therefore, the hazard rate function is increasing for 𝑝≥1 and decreasing for 0<𝑝<1. Figure 1(b) displays some plots of the hazard rate function of the half exponential power distribution with various parameters.

(a) Density function

(b) Hazard function

Figure 1 

The density and hazard rate functions of HEP(𝜎,𝑝) for 𝜎=1.

2.2. Moments and Measures Based on Moments

Proposition 2. Let 𝑋∼𝐻𝐸𝑃(𝜎,𝑝), for 𝑘=1,2,3,…; the 𝑘th noncentral moments are given by 𝜇_𝑘=𝔼𝑋^𝑘=𝑝^𝑘/𝑝𝜎^𝑘ΓΓ1/𝑝𝑘+1𝑝. (5)

The following results are immediate consequences of (5).

Corollary 3. Let 𝑋∼𝐻𝐸𝑃(𝜎,𝑝). The mean and variance of 𝑋 are given by 𝑝𝔼𝑋=^1/𝑝𝜎ΓΓ21/𝑝𝑝,(𝑋)=𝑝Var^2/𝑝𝜎²Γ1/𝑝Γ3/𝑝−Γ2/𝑝²Γ1/𝑝². (6)

Corollary 4. Let 𝑋∼𝐻𝐸𝑃(𝜎,𝑝). The skewness and kurtosis coefficients of 𝑋 are given by 𝛽₁=2Γ2/𝑝³−3Γ1/𝑝Γ2/𝑝Γ3/𝑝Γ1/𝑝Γ3/𝑝−Γ(2/𝑝)^23/2+Γ1/𝑝²Γ4/𝑝Γ1/𝑝Γ3/𝑝−Γ(2/𝑝)^23/2,𝛽₂=−3Γ2/𝑝⁴+6Γ1/𝑝Γ2/𝑝²Γ3/𝑝Γ1/𝑝Γ3/𝑝−Γ(2/𝑝)²²−4Γ1/𝑝²Γ2/𝑝Γ4/𝑝+Γ1/𝑝³Γ5/𝑝Γ1/𝑝Γ3/𝑝−Γ(2/𝑝)²². (7)

Figure 2 shows the skewness and kurtosis coefficients with various parameters for the HEP model.

(a) Skewness coefficient

(b) Skewness coefficient in log scale

(c) Kurtosis coefficient

(d) Kurtosis coefficient in log scale

Figure 2 

The plot for the skewness and kurtosis coefficients with various parameters.

3. Inference

3.1. Moment Estimation

Let 𝑋₁,𝑋₂,…,𝑋_𝑛 be a random sample from the distribution HEP(𝜎,𝑝). From (5), we have 𝔼𝑋=(𝑝^1/𝑝𝜎/Γ(1/𝑝))Γ(2/𝑝) and 𝔼𝑋²=(𝑝^2/𝑝𝜎²/Γ(1/𝑝))Γ(3/𝑝). Replacing 𝔼𝑋 and 𝔼𝑋² with the corresponding sample estimators, we obtain the moment equations 1𝑋=𝑛^𝑛_𝑖=1𝑋_𝑖=𝑝^1/𝑝𝜎ΓΓ21/𝑝𝑝,𝑋²=1𝑛^𝑛_𝑖=1𝑋²_𝑖=𝑝^2/𝑝𝜎²ΓΓ31/𝑝𝑝. (8)

The estimate 𝑝 is the solution to Γ1/𝑝Γ3/𝑝Γ2/𝑝²=𝑋²𝑋², (9) which can be solved numerically. And the estimate 𝜎 is given by 𝜎=𝑋Γ1/𝑝𝑝^1/𝑝Γ𝑝.2/ (10)

It is clear that, for the special case when 𝑝 is known, estimator 𝜎 is unbiased and its mean squared error (MSE) is given by (𝜎)=𝜎MSE²Γ1/𝑝Γ3/𝑝−Γ2/𝑝²𝑛Γ2/𝑝². (11)

In the following proposition, we present the asymtotic property of the moment estimators.

Proposition 5. Let 𝑋₁,𝑋₂,…,𝑋_𝑛 be a random sample of size 𝑛 from the distribution 𝐻𝐸𝑃(𝜎,𝑝), and let 𝜽=(𝜎,𝑝); then, if 𝜇₆=𝔼𝑋⁶<∞ and 𝜽 is the moment estimator of 𝜽, one has √𝑛𝜽−𝜽⟶^𝑑𝑁₂𝟎,𝐇⁻¹Σ𝐇^−1𝑇 (12) as 𝑛→∞, where Σ=({𝜇_𝑖+𝑗−𝜇_𝑖𝜇_𝑗}_𝑖𝑗) and 𝐇 is given by 𝐇=𝐇(𝜽)=𝜕𝜇₁𝜕𝜎𝜕𝜇₁𝜕𝑝𝜕𝜇₂𝜕𝜎𝜕𝜇₂𝜕𝑝 (13)
whose entries are given by 𝜕𝜇₁=𝑝𝜕𝜎^1/𝑝Γ2/𝑝Γ,1/𝑝𝜕𝜇₁𝑝𝜕𝑝=−^−2+1/𝑝𝜎Γ2/𝑝−1+log𝑝−𝜓1/𝑝+2𝜓2/𝑝Γ,1/𝑝𝜕𝜇₂=𝜕𝜎2𝑝^2/𝑝𝜎Γ3/𝑝Γ,1/𝑝𝜕𝜇₂𝑝𝜕𝑝=−^−2+2/𝑝𝑠²Γ3/𝑝−2+2log𝑝−𝜓1/𝑝+3𝜓3/𝑝Γ,1/𝑝 (14)

where 𝜓() is the digamma function defined as the logarithmic derivative of the gamma function, 𝜓(𝑥)=(𝑑/𝑑𝑥)logΓ(𝑥)=Γ^′(𝑥)/Γ(𝑥).

Remark 6. A consistent estimator for the asymptotic covariance matrix 𝐇⁻¹Σ[𝐇⁻¹]^𝑇 can be obtained by replacing parameters with their corresponding moment estimators.

3.2. Maximum Likelihood Estimation

In this section, we consider the maximum likelihood estimation about the parameter 𝜽=(𝜎,𝑝) of the HEP model defined in (2). The log likelihood for a random sample 𝑥₁,𝑥₂,…,𝑥_𝑛 is 𝑙(𝜽)=log^𝑛_𝑖=1𝑓𝑥_𝑖1=𝑛1−𝑝1log𝑝−𝑛log𝜎−𝑛logΓ𝑝−1𝑝𝜎^𝑝𝑛_𝑖=1𝑥^𝑝_𝑖. (15)

By taking the partial derivatives of the log-likelihood function with respect to 𝜎  and 𝑝, respectively, and equalizing the obtained expressions to zero, the following maximum likelihood estimating equations are obtained: 𝑙_𝜎𝑛=−𝜎+1𝜎^𝑛𝑝+1_𝑖=1𝑥^𝑝_𝑖𝑙=0,_𝑝=𝑛log𝑝+𝑝−1𝑝²+𝑛𝜓1/𝑝𝑝²+1+𝑝log𝜎𝜎^𝑝𝑝^2𝑛_𝑖=1𝑥^𝑝_𝑖−1𝑝𝜎^𝑝𝑛_𝑖=1𝑥^𝑝_𝑖log𝑥_𝑖=0. (16)

In general, there are no explicit solutions for the above maximum likelihood estimating equations. The estimates can be obtained by means of numerical procedures such as the Newton-Raphson method. The program 𝑅 provides the nonlinear optimization routine optim for solving such problems.

For asymptotic inference of 𝜽=(𝜎,𝑝), we need the Fisher information matrix 𝐈(𝜽). It is known that its inverse is the asymptotic variance matrix of the maximum likelihood estimators. For the case of a single observation (𝑛=1), we take the second-order derivatives of the log-likelihood function in (15).

Consider, 𝑙_𝜎𝜎=1𝜎²−𝑝+1𝜎^𝑝+2𝑥^𝑝,𝑙_𝜎𝑝=1𝜎^𝑝+1𝑥^𝑝,𝑙log𝑥−log𝜎_𝑝𝑝1=−𝑝⁴𝜎^𝑝−3𝑝𝜎^𝑝+𝑝²𝜎^𝑝+2𝑝𝑥^𝑝+2𝑝𝜎^𝑝log𝑝+2𝑝²𝑥^𝑝log𝜎+𝑝³𝑥^𝑝log𝜎²−2𝑝²𝑥^𝑝log𝑥−2𝑝³𝑥^𝑝log𝜎log𝑥+𝑝³𝑥^𝑝log𝑥²+2𝑝𝜎^𝑝𝜓1𝑝+𝜎^𝑝𝜓^′1𝑝. (17) Using the facts 𝔼𝑥^𝑝=𝜎^𝑝,𝔼𝑥^𝑝=𝜎log𝑥^𝑝𝑝log𝜎+log𝑝+𝜓1+1/𝑝𝑝,𝔼𝑥^𝑝log𝑥²=𝜎^𝑝𝑝log𝜎+log𝑝+𝜓1+1/𝑝²+𝜓^′1+1/𝑝𝑝², (18) we can obtain the elements of the Fisher information matrix: 𝐼₁₁=−𝔼𝑙_𝜎𝜎=𝑝𝜎²,𝐼₁₂=−𝔼𝑙_𝜎𝑝=log𝑝+𝜓1+1/𝑝,𝐼𝜎𝑝₂₁=−𝔼𝑙_𝑝𝜎=log𝑝+𝜓1+1/𝑝,𝐼𝜎𝑝₂₂=−𝔼𝑙_𝑝𝑝=−𝑝−𝑝²+𝑝log𝑝+𝜓1+1/𝑝²𝑝⁴+𝑝𝜓^′1+1/𝑝+𝜓^′1/𝑝𝑝⁴. (19)

Proposition 7. Let 𝑋₁,𝑋₂,…,𝑋_𝑛 be a random sample of size 𝑛 from the distribution 𝐻𝐸𝑃(𝜎,𝑝), let 𝜽=(𝜎,𝑝), and 𝜽 is the maximum likelihood estimator of 𝜽, one has √𝑛𝜽−𝜽⟶^𝑑𝑁₂𝟎,𝐈(𝜽)⁻¹. (20)

4. Assessment of Model Adequacy

In this section, we introduce a useful tool, a half normal plot with a simulated envelope which will be used to evaluate the HEP model in Section 6. The advantage of this technique is its ease of interpretation without knowing the distribution of the residuals.

Atkinson [6] proposed this diagnostic plot to detect potential outliers and influential observations in linear regression models. A simulated envelope is added to the plot to aid overall assessment, whereby the observed residuals are expected to lie within the boundary of the envelope if the presumed model has been correctly specified.

The method of simulated envelope and its corresponding transformations have been widely applied in many applications (see Flack and Flores [7], Ferrari and Cribari-Neto [8], da Silva Ferreira et al. [9], and so forth.) The simulated envelope technique compares the observed statistics with those of the data generated from the proposed model. Any sizeble departure of the observed residuals from the simulated quantities may be thought as evidence against the adequacy of the proposed model. Here is the procedure to produce the half normal plot with simulated envelopes. (1)Fit the model to the observed data (sample size = 𝑛). (2)Generate a sample of 𝑛 observations based on the fitted model. (3)Fit the model to the above generated sample and compute the ordered absolute values of the standard residuals. (4)Repeat the above steps 𝑘 times.(5)Consider the 𝑛 sets of the 𝑘-ordered statistics; calculate the average, minimum, and maximum values across each set. (6)Plot these values together with the ordered residuals from the original data against the half normal scores Φ⁻¹((𝑖+𝑛−1/8)/(2𝑛+1/2)).

The minimum and maximum values of the 𝑘-ordered statistics constitute a simulated envelope to guide assessment of the model adequacy. Atkinson [6] suggested using 𝑘=19 since there is a 5% chance to detect the largest residual being outside the boundary of the simulated envelope. Moreover, other types of residuals such as deviance or score residual may be used in the procedure. For example, da Silva Ferreira et al. [9] used the Mahalanobis distance to assess their models. The horizontal axis can also show other variables such as index.

5. Simulation Study

In this section, we conduct some simulations and study the properties of the estimators numerically.

We perform a simulation to illustrate the behaviors of the moment and MLE estimators for parameters 𝜽=(𝜎,𝑝), respectively. The simulation is conducted by the software 𝑅. We generate 1000 samples of size 𝑛=100, 𝑛=150, and 𝑛=200 from the HEP(𝜎,𝑝) distribution for fixed parameters 𝜎 and 𝑝.

The random numbers can be generated as follows. We first generate random numbers 𝑌 from an exponential power distribution with 𝜇=0, 𝜎, and 𝑝, the procedures can be found in Chiodi [10]; then we take the absolute value of the random numbers, 𝑋=|𝑌|. It follows that 𝑋∼HEP(𝜎,𝑝).

The estimators are computed using the results in Section 3. The empirical means and standard deviations of the estimators are presented in Tables 1 and 2, respectively. The simulation studies show that the parameters are well estimated, and the estimates are asymptotically unbiased. The empirical MSEs decrease as sample size increases as expected. Further, MLEs are more efficient than moment estimators.

Table 1 

Empirical means and SD for the moment estimators of 𝜎 and 𝑝.

Table 2 

Empirical means and SD for the MLE estimators of 𝜎 and 𝑝.

6. Real Data Illustration

In this section, we analyze two real datasets to fit with the proposed model. The applications demonstrate that the HEP model fits the data better than the HN model.

6.1. Application 1

The data are the plasma ferritin concentration measurements of 202 athletes collected at the Australian Institute of Sport. This dataset has been studied by several authors (see Azzalini and Dalla Valle [11], Cook and Weisberc [12], and Elal-Olivero et al. [13].)

The descriptive statistics for the dataset are shown in Table 3, where 𝑏₁ and 𝑏₂ are the sample skewness and kurtosis coefficients. Notice that the dataset presents nonnegative measurements.

Table 3 

Summary of the plasma ferritin concentration measurements.

We fit the dataset with the half normal and the half exponential power distribution, respectively, using maximum likelihood method. The MLE estimators are computed using 𝑅, and the results are reported in Table 4. The usual Akaike information criterion (AIC) and Bayesian information criterion (BIC) to measure of the goodness of fit are also computed: AIC=2𝑘−2log𝐿 and BIC=𝑘log𝑛−2log𝐿, where, 𝑘 is the number of parameters in the distribution and 𝐿 is the maximized value of the likelihood function. The results indicate that HEP model has the lower values for the AIC and BIC statistics, and thus it is a better model. Figures 3(a) and 3(b) display the fitted models using the MLE estimates.

Table 4 

Maximum likelihood parameter estimates (with (SD)) of the HN and HEP models for the plasma ferritin concentration data.

(a) Histogram and fitted curves

(b) Empirical and fitted CDF

Figure 3 

Models fitted for the plasma ferritin concentration dataset.

The diagnostic procedure introduced in Section 4 is implemented for both models. The simulated envelope plots are shown in Figures 4(a) and 4(b). Most of the observed residuals are either near or outside the boundary of the envelope, indicating inadequacy of the fitted HN model. On the other hand, the observed residuals corresponding to the HEP model in Figure 4(b) are well within the simulated envelope, indicating that the HEP model provides a better fit to the data.

(a) Half normal

(b) Half exponential power

Figure 4 

Simulated envelopes for on HN and HEP models.

6.2. Application 2

We consider the stress-rupture dataset and the life of fatigue fracture of Kevlar 49/epoxy that are subject to the pressure at the 90% level. The dataset has been previously studied by Andrews and Herzberg [14], Barlow et al. [15], and Olmos et al. [16].

Table 5 summarizes the dataset. This dataset also shows nonnegative asymmetry. Same as before, we fit the dataset with the half normal and the half exponential power distribution, respectively, using maximum likelihood method. The results are reported in Table 6. The AIC and BIC are presented as well, and the results show that HEP model fits better. Figures 5(a) and 5(b) display the fitted models using the MLE estimates.

Table 5 

Summaryofthe life of fatigue fracture.

Table 6 

Maximum likelihood parameter estimates (with (SD)) of the HN and HEP models for the life of fatigue fracture data.

(a) Histogram and fitted curves

(b) Empirical and fitted CDF

Figure 5 

Models fitted for the life of fatigue fracture dataset.

The diagnostic procedure introduced in Section 4 is implemented for both models. The simulated envelope plots are shown in Figures 6(a) and 6(b). The observed residuals corresponding to the HEP model in Figure 6(b) are well within the simulated envelope, indicating that the HEP model provides a better fit to the data.

(a) Half normal

(b) Half exponential power

Figure 6 

Simulated envelopes for on HN and HEP models.

7. Concluding Remarks

In this paper, we have studied the half exponential power distribution HEP(𝜎,𝑝) in detail. This nonnegative distribution contains the half normal distribution as its special case. Probabilistic and inferential properties are studied. A simulation is conducted and demonstrates the good performance of the moment and maximum likelihood estimators. We apply the model to two real datasets, illustrating that the proposed model is appropriate and flexible in real applications. There are a number of possible extensions of the current work. Mixture modeling using the proposed distributions is the most natural extension. Other extensions of the current work include a generalization of the distribution to multivariate settings.

Appendix

Proofs of Propositions

Proof of Proposition 2. Consider, 𝔼𝑋^𝑘=^∞₀𝑥^𝑘𝑝^1−1/𝑝𝑒𝜎Γ1/𝑝^{−𝑥𝑝/𝑝𝜎𝑝}=𝑝𝑑𝑥_1−1/𝑝𝜎Γ1/𝑝^∞₀𝑥^𝑘𝑒^{−𝑥𝑝/𝑝𝜎𝑝}=𝑝𝑑𝑥_1−1/𝑝𝜎𝜎Γ1/𝑝^𝑘+1𝑝^{(𝑘+1)/𝑝−1}Γ𝑘+1𝑝=𝑝^𝑘/𝑝𝜎^𝑘ΓΓ1/𝑝𝑘+1𝑝. (A.1)

Proof of Proposition 5. This result follows directly by using standard large sample theory for moment estimators, as discussed in Sen and Singer [17].

Proof of Proposition 7. It follows directly by using the large sample theory for maximum likelihood estimators and the Fisher information matrix given above.