New generalized-X family: Modeling the reliability engineering applications

As is already known, statistical models are very important for modeling data in applied fields, particularly in engineering, medicine, and many other disciplines. In this paper, we propose a new family to introduce new distributions suitable for modeling reliability engineering data. We called our proposed family a new generalized-X family of distributions. For the practical illustration, we introduced a new special sub-model, called the new generalized-Weibull distribution, to describe the new family’s significance. For the proposed family, we introduced some mathematical reliability properties. The maximum likelihood estimators for the parameters of the new generalized-X distributions are derived. For assessing the performance of these estimators, a comprehensive Monte Carlo simulation study is carried out. To assess the efficiency of the proposed model, the new generalized-Weibull model is applied to the coating machine failure time data. Finally, Bayesian analysis and performance of Gibbs sampling for the coating machine failure time data are also carried out. Furthermore, the measures such as Gelman-Rubin, Geweke and Raftery-Lewis are used to track algorithm convergence.


Introduction
Within the area of reliability engineering and other related fields, modeling of data related to lifetime events is very crucial. A range of probability models, such as Weibull, gamma, exponential, etc., are available for modeling lifetime data. However, in many cases, these classical models are not suitable for modeling lifetime data, and there is always a clear need for modified forms of these existing distributions. Therefore, the researchers have introduced new families that are more advantageous than the old ones. For more reading about statistical distributions using different approaches; see [1].
In the literature, the Weibull distribution is the most outstanding one that has broadly been used in reliability engineering and in other various areas of research; see, for instance, [2]. Although the Weibull distribution is often used, the confined structure of its hazard function (hf) can only be decreasing, increasing or constant. Generally, many practical problems require a flexible range of hf, for example, the lifetime events that exhibit a bathtub-shaped hf such as human mortality and life cycles of electronic machines and components. Researchers in the past couple of years developed more flexible extensions of the Weibull model to model reliability data adequately. For the recent survey about such distributions, we refer to [3]. The foremost goal of this research is to introduce a new flexible modification of the Weibull model by means of inducting one additional parameter. The induction of the additional parameter leads to the greater flexibility to enhance goodness-of-fit to reliability data. In fact, we show empirically that the new extension of the Weibull distribution offers the best fit to the coating machine failure time data than the two-parameter, three-parameter, and four-parameter competitive distributions (see section 5). The practical example really shows that the proposed distribution is a good alternative candidate for modeling reliability data. Now we introduce the proposed family of distributions called a new generalized-X (NG-X) family. A random variable X is said to follow a NG-X family, if its cumulative distribution function (cdf) is given by where θ is an extra shape parameter, and F(x;ξ) is the baseline cdf which may depend on the parameter vector x 2 R. By adding the additional shape parameter, the NG-X distributions can provide best fit to reliability engineering data. The corresponding density function is g x; y; x ð Þ ¼ f x; x ð Þ½1 À Fðx; xÞ� yÀ 1 ð1 þ yÞ À Fðx; xÞ e Fðx;xÞ The reliability function S(x;θ, ξ) and the hf h(x;θ, ξ) of NG-X distributions are given by S x; y; x ð Þ ¼ ½1 À Fðx; xÞ� y e Fðx;xÞ ; x > 0; and h x; y; x ð Þ ¼ f ðx; xÞ 1 À Fðx; xÞ ð1 þ yÞ À F x; x ð Þ ½ �; x > 0; respectively. Using the NG-X distributions approach, we introduce a new form of the Weibull model called, a new generalized Weibull (NG-Weibull) distribution. Furthermore, we consider maximum likelihood (Non-Bayesian) and Bayesian procedures in order to estimate the unknown parameters of the NG-Weibull model. In the Bayesian discussion, we consider different types of symmetric and asymmetric loss functions such as squared error loss, weighted squared error, precautionary, K-loss, and modified squared error loss function to estimate the unknown parameters of the NG-Weibull model. Since all the parameters are positive, we use gamma prior distributions. Bayesian 95% credible and highest posterior density (HPD) intervals (see [4]) are obtained for every parameter of the NG-Weibull model. We used the Gibbs sampling technique to get posterior samples. From a graphical point of view, we graphed the posterior density function plots. Next, for evaluating the MCMC procedure in Bayesian analysis, we reported diagnostics measures such as Gelman-Rubin, Geweke, and Raftery-Lewis for checking the convergence of the algorithm. This paper is outlined in the following manner: In Section 2, we define a NG-X family. Some mathematical properties of NG-X distributions are derived in Section 3. Section 4 is specified for obtaining the estimates using the maximum likelihood estimation, and the Monte Carlo simulation study is also provided in the same section. Section 5 is concerned with the goodness of fit of the proposed distribution. In this section, we showed that NG-Weibull model provides fit to reliability engineering data. Section 6 offers the Bayesian analysis. The future research directions are provided in Section 7. Some concluding comments are presented in Section 8.

A new generalized Weibull model
Let Fðx; xÞ ¼ 1 À e À gx a ; x � 0; be the cdf of the Weibull distribution, where ξ = (α, γ). Then, a random variable say X is said to follow the NG-Weibull distribution, if its cdf is given by The density function corresponding to Eq (3) is The reliability function and hf of the NG-Weibull distribution are given by x > 0; and hðx; y; xÞ ¼ agx aÀ 1 ½y þ e À gx a �; x > 0; respectively. In the Fig 1, we have sketched the density function plots of the NG-Weibull distribution. Fig 1 shows that the NG-Weibull density can be reverse J-shape, symmetric, positively skewed, negatively skewed, and bi-model. The hf plots of the NG-Weibull model are presented in Fig  2. The NG-Weibull hf can be monotonically decreasing, increasing, uni-modal, and modified uni-modal shaped.

Mathematical properties
In this part of the paper, we derived the mathematical properties associated with the NG-X distributions, which include identifiability, quantile function, random number generation, r th noncentral moments, and the Renyi entropy with numerical illustrations. A characterization theorem extending the NG-X class of distributions in terms of the hf is also provided.

Identifiability
The identifiability is an important statistical property that a model must obey to make sure that the inference should be precise. In this subsection, we prove the identifiability property of the NG-X distributions. To prove the identifiability property of the NG-X distributions, we have to show that θ 1 = θ 2 . Let θ 1 and θ 2 be the two parameters having the NG-X distributions with cdfs given by G(x;θ 1 , ξ) and G(x;θ 2 , ξ), respectively. From the definition of identifiability,

Random number generation
If U � Uniform(0, 1), then random numbers from NG-X distributions can be obtained from where F −1 is the quantile of the distribution with cdf F(x), and W(z) gives the principal solution for w in z = we w .

The r th non central moments Theorem 2. The r th Non Central Moments of NG-X distributions can be expressed as
where O i,k, n, q is defined as in the proof of the theorem, E[�] denotes an expectation, and U � Uniform(0, 1). Proof. From subsection 3.3, we know if U � Uniform(0, 1), then the following random variable where F −1 is the quantile of the distribution with cdf F(x), and W(z) gives the principal solution for w in z = we w , follows the NG-X family of distributions. According to [5], we can write where the coefficients are suitably chosen real numbers that depend on the parameters of the F(x) distribution. For a power series raised to a positive integer r � 1, we have where δ r,i are obtained from d r;i ¼ ðih 0 Þ À 1 P i s¼1 ½sðr þ 1Þ À i�h s d r;iÀ s with d r;0 ¼ h r 0 for i = 1, 2, . . .; see [6]. Thus we have the following where E(�) is an expectation. By the binomial series, we can write By the binomial series we can write It follows that Some numerical description of the ordinary moments are presented in Table 2.

Renyi entropy
Theorem 3. The Renyi entropy of the NG-X distributions, for δ 6 ¼ 1, δ > 0, can be expressed as where X is a random variable with cdf F(x) and pdf f(x), and O k,q, r is defined as in the proof of the theorem. Proof. Recall the pdf of NG-X distributions is given by We first find an expansion for g(x) δ where δ 6 ¼ 1 and δ > 0. By the binomial series we can write By the binomial theorem we can write By the power series representation for the exponential function, we can write It follows that Therefore, the Renyi entropy is

Characterization theorem
It is clear that hf, of a function, F, that can be differentiated twice satisfies the following differential equation In this section, we present a Kumaraswamy NG-X type distribution. The result here is inspired by [7]. First, let us introduce the following.
Definition: We say a random variable X follows a Kumaraswamy-G type distribution if its cdf is given by where G is some baseline distribution, x 2 Supp(G), and ξ is a vector of parameters in the baseline distribution whose support depends on G. Remark: Note that if we take λ = 1 and φ = 2 in Eq (1) of [8], then we get the cdf in the above definition.
The pdf of the Kumaraswamy-G type distribution is given by where g is the pdf of the baseline distribution. Clearly the hf of the Kumaraswamy-G type distribution is given by with boundary condition h(0) = 2g(0). Proof. If X has pdf as stated in the theorem, then the differential equation as stated holds. Now if the stated differential equation holds, then d dx which is the hf of the Kumaraswamy-G type distribution. Clearly, a characterization of the Kumaraswamy NG-X type distribution is obtained from the above theorem by letting the baseline pdf and cdf given in section 1.

Maximum likelihood estimation and Monte Carlo simulation
This section is devoted to estimating the NG-X parameters using the maximum likelihood estimation approach and providing a comprehensive Monte Carlo (MC) simulation study to assess the maximum likelihood estimators (MLEs) performance.

Maximum likelihood estimation
Let x 1 , x 2 , . . ., x n be the observations of a random sample of size n taken from the NG-X distribution with the parameter vector Θ = (θ, ξ) T . For Θ, the log-likelihood function (LLF) is given by Fðx; xÞ: ð5Þ The partial derivatives of the LLF are given by respectively. The MLEs of the unknown parameters θ and ξ of the NG-X distributions can be obtained by maximizing @ @y 'ðYÞ ¼ 0 and @ @x 'ðYÞ ¼ 0, respectively.

Monte Carlo simulation study
In the following sub-section, we assess the behavior of the MLEs of NG-Weibull distribution by means of the MC simulation study. The process is carried out by maximizing the LLF using the optim() R-function with the argument method = "L-BFGS-B". We made 1000 MC-iterations using different sizes of the samples as follows, n = 25, 50,. . .,1000. We computed the average MLEs, the associated mean square errors (MSE), biases and absolute biases.

Application on our model from reliability engineering
In this section, we evaluate the usefulness of the NG-Weibull distribution by means of analyzing reliability engineering data taken from [9]. The data represents the failure time of the coating machine given by:  [10], Kumaraswamy Weibull (Ku-W) distribution [11], Extended Alpha Power Transformed Weibull (Ex-APTW) [12] and type-I heavy-tailed Weibull (TI-HTW) distribution [13].
To figure out about the goodness of fit amongst the competitive distributions, we consider certain goodness of fit measures such as Cramer-Von-Mises (CM) statistic, Anderson Darling (AD) statistic, and Kolmogorov-Smirnov (KS)statistics alongside with its p-values.
Corresponding to the coating machine failure time data, the model parameters' estimated values with standard errors in the parenthesis are presented in Table 3. The goodness of fit measures of the competitive models are provided in Table 4. From the results reported in Table 4, we can see that the proposed model has lower values of the goodness of fit measures and a high p-value indicating the best fit for the reliability data.
In support of the goodness of fit measures given in Table 4, the estimated cdf and Kaplan-Meier survival plot of the NG-Weibull distribution are plotted in Fig 7. The probability-probability (PP) and quantile-quantile (QQ) plots are presented in Fig 8. These figures confirm the best fitting of the NG-Weibull to the coating machine failure time data.
Furthermore, for the coating machine failure time data, we calculated the KS statistic values of the NG-Weibull and other considered models. Then, we utilized the parametric bootstrap

PLOS ONE
New generalized-x family approach [14], and bootstrapped the p-value for all models. The KS statistic and the corresponding bootstrapped p-value are reported in Table 5. According to the results provided in Table 5, we observe that the NG-Weibull is a good candidate model amongst the competing distributions for modeling engineering reliability data.

Bayesian estimation
In this section, we consider different types of symmetric and asymmetric loss functions such as squared error loss function (SELF), weighted squared error loss function (WSELF),   Table 6. For more detail, we refer to [15][16][17][18]. Next, we provide a Bayesian estimation approach for estimating the parameters of NG-Weibull distribution via analyzing complete sample data.

Bayesian point estimation
Under the marginal posterior density function as in Eq (10) and the loss functions which are given in Table 6. The Bayesian point estimation for the parameter vector C = (C 1 , C 2 , C 3 ) = (α, γ, θ) is obtained via minimizing the expectation of loss function under the marginal posterior density as follows: Eðc À 1 jxÞ Eðc À 2 jxÞ 1 À Eðc À 1 jxÞ 2 However, in practice, because of the intractable integral in relation Eq (11), it is suggested to use the well-known Gibbs sampler [19], or Metropolis Hastings algorithms to generate posterior samples [20]. We will argue about this issue more precisely in subsection 6.5.

Credibility interval
In the Bayesian framework, interval estimation can be done via credibility interval conception. Consider the parameter vector C = (C 1 , C 2 , C 3 ) = (α, γ, θ), which is associated with the NG-Weibull distribution and pðC i jxÞ denote the marginal posterior pdf of the parameter C j ; (j = 1, 2, 3) as in Eq (10). For a given value of η 2 (0, 1), the (1 − η)100% credibility interval By considering the relation Eq (12), it is very difficult to obtain the marginal density from the joint posterior density. We use the Gibbs sampler to generate posterior samples. Let C 1 , be a posterior random sample of size k, which is extracted from the joint posterior density as in Eq (8). Using these generated posterior samples, the marginal posteriors densities of C j given x can be given by where the C i À j shows the vector of posterior samples when jth component is removed. Using Eq (13) in Eq (12), one can be able to compute the credibility intervals for C j , j = 1, 2, 3 as follows

Highest posterior density interval
The HPD interval is a kind of credibility interval with a specific restriction. The (1 − η)100% (i = 1, . . ., p) HPD interval for C j , j = 1, 2, 3 is the simultaneous solution of the following integral equations

Generating posterior samples
It is clear from Eqs (8) and (12) that there is no closed-form for point estimators using different loss functions, because of intractable integrals. So we will try to solve those integrals numerically using MCMC methods. There are a lot of possible methods. One of these methods is known as the Metropolis-Hastings algorithm. Another method for approximation of unsolvable integrals is known asthe Gibbs sampling. Suppose that the general model f ðxjψÞ is associated with the parameter vector ψ = (ψ 1 , ψ 2 , . . ., ψ p ) and observed data x. Thus, the joint posterior distribution is pðc 1 ; c 2 ; . . . ; c p jxÞ. We also assume that c 0 ¼ ðc ð0Þ 1 ; c ð0Þ 2 ; . . . ; c ð0Þ p Þ is the initial vector to start the Gibbs sampler. The steps for any iteration, say iteration k, are as follows: • Starting with an initial estimate ðc ð0Þ 1 ; c ð0Þ 2 ; . . . ; c ð0Þ p Þ • Draw c k 1 from pðc 1 jc kÀ 1 2 ; c kÀ 1 3 ; . . . ; c kÀ 1 p ; xÞ • Draw c k 2 from pðc 2 jc k 1 ; c kÀ 1 3 ; . . . ; c kÀ 1 p ; xÞ; and so on down to • Draw c k p from pðc p jc k 1 ; c k 2 ; . . . ; c k pÀ 1 ; xÞ: In case of the NG-Weibull distribution, by considering the parameter vector C = (α, γ, θ) and initial parameter vector C 0 = c(α 0 , γ 0 , θ 0 ), the posterior samples are extracted by above Gibbs sampler where the full conditional distributions are given as and pðyja kÀ 1 ; g kÀ 1 ; xÞ / y y 0 À n e À y 1 y Y n i¼1 e À ygx a i fy þ e À gx a i g: Gibbs sampling processes can be carried out via OpenBUGS software, which is an available version of WinBUGS. Here, since there aren't any prior information about hyper-parameters in Eq (6), we implement the idea of [21], and the hyper-parameters values are setted as α 0 = α 1 = γ 0 = γ 1 = θ 0 = θ 1 = 0.0001. We can use the MCMC procedure to extract posterior samples of Eq (8) by means of the Gibbs sampling process in OpenBUGS software.
Next, we provide Bayesian estimation results. It is evident from equation Eq (10), there are no closed-form expressions for Bayesian estimators, which are extracted based on the loss functions in Table 6. Therefore, an MCMC procedure via the Gibbs sampler process is designed using the expressions provided in Eqs (15), (16) and (17), with 10,000 replicates to obtain the Bayesian estimators. In Table 7, we provide the corresponding point and posterior risk estimations. Furtherer, 95% credible, and HPD intervals are provided in Table 8. In order to provide a visual inspection, the posterior plots such as trace plots are provided in Fig 9, autocorrelation plots are presented in Fig 10, and the histogram plots are sketched in Fig 11. These plots verify that the convergence of Gibbs sampling process has occurred. Next, for evaluation the MCMC procedure in Bayesian analysis, we report some diagnostics measures such as Gelman-Rubin (GR), Geweke (G) and Raftery-Lewis (RL) for checking the convergence of the Gibbs algorithm are provided in Table 9. For more details about these indexes; see [22]. The GR diagnostic for parameters α, γ, and θ is equal to 1. Hence, based on the GR diagnostic measure, the chains are acceptable. Fig 12 shows that the estimates come from state spaces of the corresponding parameters. From Table 9, Geweke's test statistics for parameters α, γ and θ, are −0.0795, 0.0807 and 0.59370, respectively. Hence, the G diagnostic measure also confirms the acceptance of chains as shown in Figs 13 and 14. Moreover, the reported diagnostic statistics for parameters α, γ and θ based on the RL method do not show a significant degree of dependence between estimates.

Discussion and future framework
The statistical decision theory helps to address the state of uncertaintyand provides potentially a sound framework for dealing with problems of bio-medical, reliability, actuarial, economic and financial decision-making. Among the applied fields, the reliability engineering has received a serious consideration. In the field of reliability theory, modeling of lifetime data is very crucial. The data related to the lifetime of electronic product or any entity, etc., are usually positive and skewed to the right. In lifetime analysis and reliability theory, the failure rate function (also known as hazard rate function) is a prominent reliability characteristic. Among the possible failure rate functions, the unimodal, modified unimodal or bathtub-shaped failure rate curves are well-known in reliability literature. The classical models are not so flexible to model such complex forms of data. Due to the importance and application of the statistical models in reliability, medical and financial sciences, a reasonable work has been done in the literature aiming to improve the characteristics of the classical models. Although the new improvement has achieved the respective goal, unfortunately, the numbers of parameters have been increased, and the estimation of parameters, statistical inference and derivation of mathematical problems become complicated.
To provide a better description and best fitting to the reliability data, therefore, in the present study, a new family of distributions has been studied. The key goal of introducing the class of distributions is to improve the characteristics of the classical distributions. From the above theory and discussion, it is quite clear that the researchers are always in search of new flexible distributions. Therefore, to introduce new flexible distributions and bring further flexibility in the model proposed in this paper, we suggest to introduce modified forms of the proposed model. As a future research direction, we suggest to introduce new models useful for modeling data in reliability engineering and financial sciences.
As we stated above that the statistical models with bathtub-shaped failure rate function are very useful in reliability engineering. Here, we suggest a new modification of the proposed distribution to model lifetime data with a bathtub-shaped failure function.
Furthermore, in the practice of actuarial and financial sciences, the heavy-tailed distributions are useful for modeling heavy-tailed financial data sets. Heavy-tailed distributions are those, whose right tail probabilities are heavier than the exponential distribution, and satisfies lim x!1 e px ½1 À FðxÞ� ¼ 1; p > 0: Here, we suggest a new heavy-tailed distribution to provide the best description of the heavy-tailed financial data set. The new heavy-tailed distribution can obtained by using the Burr-XII (B-XII) distribution as a special case of the NG-X family.
Let F(x;ξ) = 1 − (1 + x c ) −k , be the cdf of the two-parameter B-XII distribution, where ξ = (c, k). Then, a random variable say X is said to follow the new generalized B-XII (NGB-XII) distribution, if its cdf is given by x � 0: The corresponding pdf is given by G x; y; x ð Þ ¼ ckx cÀ 1 ð1 þ x c Þ À ky y þ ð1 þ x c Þ À k e 1À ð1þx c Þ À k ( ) x > 0: From the plots provided in Fig 16, we can see that increasing the value of c and θ, the NGB-XII distribution tends to a heavy-tailed distribution. Henceforth, the NGB-XII distribution can be a good candidate distribution to model heavy-tailed financial and other related data sets. In the future, therefore, we are intended to use the proposed method to introduce new heavy-tailed distributions.

Concluding remarks
Due to the great importance of statistical distributions in modeling data in reliability engineering, we introduced and studied a new family of distributions, called a NG-X family. The MLEs of the model parameters along with some mathematical properties, are derived. Based on the proposed approach, a new modified form of the Weibull model called an NG-Weibull distribution is introduced and studied in detail. The NG-Weibull model is very versatile and is able to cater to the different patterns of failure rates. Due to the flexible behavior of the hf, the proposed model is capable to describe adequately the failure behavior of several lifetime datasets, particularly reliability engineering data. The usefulness of the NG-Weibull distribution is proved by analyzing the coating machine failure time data. We performed Bayesian estimation and estimated the model parameters using five different loss functions. The diagnostics measures such as the Gelman-Rubin, Geweke, and Raftery-Lewis are discussed to evaluate the MCMC procedure in the Bayesian analysis. As a future work, new models based on the proposed approach will be introduced to model complex form of reliability engineering and heavy-tailed financial data sets.