Modeling extreme value data with an upside down bathtub-shaped failure rate model

1 Introduction

The heavy right-tailed Pareto model which is characterized by the distribution function

and the probability distribution function (PDF)

occurs in a diverse range of physical phenomena. Generally, it is useful when there is an equilibrium in distribution of “small” to “large” values, e.g., the size of transmitted files on a computer network consisting many small files and few large ones, or the size of human settlements consisting of many small and few large cities villages/hamlets. Moreover, the sizes of solar flares, oil reserves in oil fields, earthquakes, corporations, and lunar craters have similar property which is referred to as “power law” property. Newman [1] reviewed some power law forms and theories explaining them. The Pareto model is recognized by its heavy right tail in the literature and shows a decreasing failure rate function. It is useful in biology, reliability engineering, survival analysis, quality control, economics, computer science, geophysics, and many other scientific fields. For detailed information about Pareto and related distributions and their features see Arnold [2], Zhang et al. [3], and Zhang et al. [4].

Bak and Sneppen [5], Sornette [6], and Carlson and Doyle [7] among many others used the Pareto as a power law model in their research. Also, Burroughs and Tebbens [8] fitted the Pareto model to earthquake and wildfire observations, and Schroeder et al. [9] described plate fault data by the Pareto model. Moreover, some researchers defined modified versions of the Pareto model and applied them in their studies. These modified Pareto models are more flexible for data generated in various phenomena. For example, Akinsete et al. [10] introduced a beta Pareto model; Nassar and Nada [11] and Mahmoudi [12] proposed a beta generalization of the Pareto distribution, Alzaatreh et al. [13] used the gamma distribution to propose a modified Pareto model; and Zea et al. [14], Elbatal [15] and Bourguignon et al. [16] defined extensions of the Pareto distribution. Papastathopoulos and Tawn [17] applied one extended Pareto model for tail estimation. Moreover, Mead [18], Elbatal and Aryal [19], Korkmaz et al. [20], Ghitany et al. [21], Tahir et al. [22], Ihtisham et al. [23], Haj Ahmad and Almetwally [24], Jayakumar et al. [25], and recently Jayakumar et al. [26] defined and studied a new model with heavier right tail than Pareto.

In this article, a new flexible generalized Pareto distribution with heavy right tail and upside down bathtub-shaped (UBT) FR function is introduced and studied. The novelty of the model is that it gathers the heavy right tail same as the Pareto model and UBT FR form in one model. Thus, the main advantage of the proposed model is that it is useful when the data show a fat right tail and UBT FR function. Such data can be observed in hydrology or other situations with extreme values. The remaining of the article is organized as follows. In Section 2, the new model is defined and its basic attributes like the moments and quantiles are discussed. Then, some important dynamic measures of it like FR and mean residual life (MRL) functions are studied. The aim of Section 3 is to estimate the parameters of the proposed model. In Section 4, the efficiency and consistency of the considered models are investigated by simulations. Then, the proposed model and some alternatives are fitted to consecutive flood discharges of the Floyd river located in James, Iowa, USA, during 1935 to 1973 to show its applicability.

2 The new modified Pareto distribution

The new generalized Pareto, GP ( a , b , c ) , model is defined by the distribution function,

and the PDF

In Figure 1, the PDF is drawn for some values of the parameters and show a unimodal form for it. It seems that the coefficient e − b x changes the form of the PDF from decreasing to increasing in an early period and the mode of the model increases with b . If b = 0 , the GP shows the baseline Pareto model and if a tends to zero, it converges to a model with the distribution function

which is a special case of the modified Weibull model defined by Kayid and Djemili [27].

Figure 1

Left: The PDF of GP for some parameter values. Right: The FR function of GP for some parameter values.

Like the baseline Pareto, the proposed GP model has a heavy right tail. For example, in comparison with the well-known Weibull distribution we can write

Moreover, in comparison with the baseline Pareto model

which indicates that GP has a heavy right tail like the baseline Pareto model.

The expectation of a random variable X following the Pareto distribution function (1) is finite and equals c ( 1 − a ) − 1 for a < 1 and is infinite for a ≥ 1 . Moreover, the k th moment of this random variable X is finite for a < k − 1 and infinite for a ≥ k − 1 . The following proposition proves similar result for the k th moment of the proposed GP model. It shows that when a < 1 , the expectation of GP ( a , b , c ) is finite and otherwise it is infinite.

One good tendency measure which may be applied in place of the moments is the quantile function which at point p equals q ( p ) = F − 1 ( p ) . This function is useful in generating simulations, estimating the parameters and describing the model characteristics like skewness and kurtosis. For example, Bowley [28] and MacGillivray [29] defined the skewness based on the quantile function. Moreover, Moors [30] presented a measure of kurtosis in terms of quantiles. For GP ( a , b , c ) , the quantile function at p could be obtained by solving the following equation in terms of x .

where A = c a ( ( 1 − p ) − a − 1 ) .

2.1 Dynamic measures

The FR function of the proposed model is

(7) λ ( x ) = 1 c e − b x 1 + b x 1 + a c x e − b / x − 1 , a > 0 , b > 0 , c > 0 , x ≥ 0 .

Proposition 2

The FR function of GP ( a , b , c ) is UBT. Moreover, the point maximizing the FR function is the root of the following equation:

(8) − 2 a b c x 2 − a c x 3 + b 2 e b x = 0 .

Proof

By differentiation of the FR function, we found that the sign of h ′ ( x ) is the same as sign of the following function:

(9) γ ( x ) = − 2 a b c x 2 − a c x 3 + b 2 e b x ,

which is a decreasing function and γ ( x ) → ∞ as x tends to zero, and γ ( x ) → − ∞ as x tends to ∞ .□

In fact, the coefficient e − b x , included in the model, affects on the early life and makes the decreasing FR function to increasing in a beginning period of life.

The proof of the following proposition could be trivially obtained by comparing the reliability functions of the assumed random variables and is omitted.

Proposition 3

Let X 1 ∼ GP ( a , b 1 , c ) and X 2 ∼ GP ( a , b 2 , c ) and b 1 < b 2 . Then, X 1 < X 2 is in stochastic order.

Figure 1 draws the FR function for some parameters and shows that under the conditions of Proposition 3, X 1 is not smaller than X 2 in FR ordering. Moreover, since likelihood ratio ordering is stronger than FR ordering, X 1 would not be smaller than X 2 in likelihood ratio ordering.

The MRL function of GP ( a , b , c ) is finite for a < 1 and is given by

m ( x ) = ∫ x ∞ R ( t ) d t R ( x ) = 1 + a c x e − b x 1 a ∫ x ∞ 1 + a c t e − b t − 1 a d t .

By Proposition 2, it results that the MRL function has an increasing or bathtub form, see Lai and Xie [31]. Figure 2 shows the MRL for some parameter values.

Figure 2

The MRL and median residual life functions of GP for some parameter values.

Another prominent dynamic tendency measure is the α -quantile residual life ( α -QRL) function, which is given by

q α ( x ) = F − 1 ( 1 − α ¯ ( 1 − F ( x ) ) ) − x .

The special case q 0.5 ( x ) is referred to the median residual life function. The α -QRL and specially the median residual life function are better than MRL for models with heavy right tail, specially the Pareto distribution since the MRL may be infinite. Figure 2 exhibits the median residual life for some values of the parameters. Note that when the FR function is UBT, the α -QRL will be increasing or bathtub shaped, see Kayid [32] for detailed information.

3 Inference

Assume x 1 ≤ x 2 ≤ … ≤ x n represents an ordered, independent and identically distributed sample of GP ( a , b , c ) . In this section, the maximum likelihood (ML), least squared errors (LSE), and Anderson–Darling (AD) methods are discussed for estimating the parameters of the model.

4 Simulation study

To simulate one random variable X from GP ( a , b , c ) , first we generate one random instance U from standard uniform distribution and solve the equation F ( X ) = U in terms of X , where F is the distribution function of X .

In this simulation study, some values for parameters are selected. Then, in every run, r = 1,000 repetitions of size n = 80 or 150 are simulated and the parameters are estimated by one of the ML, LSE, or AD approaches. The bias (B) and mean squared error (MSE) of the estimators are computed by the following relations, respectively.

and

and similarly for b and c . All computations are performed in R environment. The optimization problems are solved by the built-in “optim” function of R . The initial values needed in this function are selected randomly from uniform distribution, e.g., for a , from uniform distribution on the interval ( 0.9 a , 1.1 a ) . Table 1 abstracts the results of the simulation study. The results show that all studied estimators are consistent and efficient but the MLE outperforms the LSE and AD. On the other hand, AD estimator gives smaller MSE than LSE.

5 Application

Table 2 shows the consecutive flood discharges in terms of f t 3 / s for the Floyd river located in James, Iowa, USA, during 1935 to 1973, see Mudholkar and Hutson [33]. The flood discharges are extreme values and analyzing them could completely helpful in predicting the extreme flood occurrences. This data set was analyzed by Mudholkar and Hutson [33] and Merovcia and Puka [34]. The box plot of the data presented in Figure 3 shows one extra ordinary large value 71,500. Figure 3 shows the total time on test (TTT) plot and verifies a UBT FR function. In a comparative analysis, the proposed GP along with some other UBT FR models and some models with other FR forms are participated, and the results are abstracted in Table 3.

Figure 3

The box plot (left panel) and TTT plot (right panel) for the flood discharge data.

The alternative models are Pareto; exponentiated Pareto (EP); Marshal–Olkin Pareto (MOP); Dimitrakopoulou, Adamidis, and Loukas (DAL) modified Weibull model proposed by Dimitrakopoulou et al. [35]; inverse Weibull (IW); Marshal–Olkin inverse Weibull (MOIW); gamma, Marshal–Olkin gamma (MOG); Weibull and Pareto exponential competing risk (PECR).

The parameters of the mentioned models are estimated by the ML method. The R programming language was used for computations, and all optimizations were done by the built-in function “optim” of R. The Akaike information criterion (AIC), Cramer–von Mises (CVM) statistics, AD and Kolmogorov–Smirnov (KS) statistics are reported for every model. Clearly, the proposed GP and MOIW show a close-run. However, the GP outperforms other models and provides a good description of the data. Figure 4 draws the empirical and fitted distribution function for GP and some of the alternatives which show better fits. The estimated FR function is plotted in Figure 5 and confirms a UBT form for the FR function. Also, histogram of the data and estimated PDF are plotted in the right side of Figure 5.

Figure 4

The empirical distribution function along with estimated distribution function for some alternatives.

Figure 5

The estimated FR function of GP (left panel) and the histogram of the data with the estimated PDF of GP for flood discharge data.

6 Conclusion

One new flexible GP model which preserves the heavy right tail attribute but exhibits an early increasing FR function is introduced. The limiting behavior of the proposed model is similar to the baseline Pareto, but the attributes differ at beginning of the support. The proposed GP model has a UBT FR function. The simulation results show that the ML estimator is efficient and consistent. Applying the model on one flood discharge data of the Floyd river shows that the proposed GP model could be useful in describing many data sets which occur in a wide variety of physical phenomena. There are many future related topics. For example, studying a mixture of the proposed GP model or introducing proper extensions of the GP model based on the underlying physical justifications.