The Weibull Birnbaum-Saunders distribution and its applications

A new lifetime model, with four positive parameters, called the Weibull Birnbaum-Saunders distribution is proposed. The proposed model extends the Birnbaum-Saunders distribution and provides great flexibility in modeling data in practice. Some mathematical properties of the new distribution are obtained including expansions for the cumulative and density functions, moments, generating function, mean deviations, order statistics and reliability. Estimation of the model parameters is carried out by the maximum likelihood estimation method. A simulation study is presented to show the performance of the maximum likelihood estimates of the model parameters. The flexibility of the new model is examined by applying it to two real data sets.


Introduction
For modeling the material fatigue failure time process, Birnbaum and Saunders [5] developed a two-parameter family of continuous probability models, known as the Birnbaum-Saunders (BS) distribution or as the fatigue life distribution. Since the BS distribution was proposed, it has received much attention in the literature and has been used in several fields such as engineering [19], insurance [26], medical sciences [2], wind analysis [21] and econometrics [4], among others. For more details on the BS distribution, we refer to [3], [17] and [18].
A random variable T follows the BS distribution with shape parameter α > 0 and scale parameter β > 0, if it can be written as follows where Z is a standard normal random variable. The cumulative distribution function (cdf) of T is where Φ(·) is the standard normal distribution function, v = α −1 ρ(t/β) and ρ(z) = z 1/2 − z −1/2 . The probability density function (pdf) corresponding to (1) is given by and the function K ϵ (s) denotes the modified Bessel function of the third kind with ϵ representing its order and s the argument. The parameter β is the median of the BS distribution, because G(β) = Φ(0) = 1/2.
The BS model is not appropriate for modeling highly skewed and heavy-tailed data sets. Also, the BS distribution can exhibit unimodal (upside-down bathtub) shaped but not more other shapes like the monotone (decreasing or increasing) and bathtub-shaped hazard rate functions. For these reasons, using different methods, several authors have proposed various generalizations and extensions of the BS distribution. For example, Cordeiro and Lemonte [8] proposed the beta BS distribution. Saulo et al. [29] introduced the Kumaraswamy BS distribution. Cordeiro et al. [11] defined the McDonald BS distribution. Cordeiro et al. [10] introduced the gamma BS distribution. Hashemi et al. [15] proposed the normal mean-variance Lindley BS distribution. Naderi et al. [24] proposed the normal mean-variance generalized BS distribution. Naderi et al. [25] presented the skew-Laplace BS distribution.
Based on the Weibull cdf F W (w) = 1 − e −aw b for w > 0, a > 0 and b > 0, Bourguignon et al. [7] proposed a more general family with extra parameters a and b, called the Weibull-G family of distributions. The cdf of the Weibull-G family is given by where G(t, θ) is the baseline cdf depending on a parameter vector θ andḠ(t, θ) = 1 − G(t, θ). Its corresponding pdf is f (t; a, b, θ) = abg(t, θ) where g(t, θ) is the baseline pdf. Many authors have used this family to propose an extension of some classical distributions. Tahir et al. ([30], [31], [32]) defined the Weibull-Pareto, Weibull-Lomax and Weibull-Dagum distributions by taking G(t, θ) to be the cdf of the Pareto, Lomax and Dagum distributions, respectively. Afify et al. [1] defined and studied the Weibull Burr XII distribution. In a similar way, we propose a new extension of the BS distribution called the Weibull BS (WBS) distribution.
The rest of the paper is organized as follows. In Section 2, we introduce the WBS distribution and plot its density and failure rate functions. Mathematical properties of the new distribution are derived in Section 3. In Section 4, we discuss maximum likelihood estimation of the WBS parameters and derive the observed information matrix. In Section 5, we carry out a simulation study to examine the performance of the maximum likelihood estimators for WBS model. Two applications are presented in Section 6 to show the flexibility of the new distribution. Conclusions are given in Section 7.

The WBS distribution
The cdf of the WBS distribution is obtained by inserting (1) in (4) as follows (we have omitted the dependence on the parameters α, β, a and b)

63
The pdf corresponding to (6) is given by where β is a scale parameter and α, a and b are positive shape parameters. It is clear that the BS distribution is not a special case of WBS distribution. Henceforth, a random variable T with pdf (7) is denoted by T ∼ WBS (α, β, a, b). The reliability and the failure rate function of T are, respectively, given by The plots of pdf and failure rate function for the WBS distribution are represented in Figures 1 and 2 respectively for different values of parameters α, β, a and b. It is clear that, from Figure 1, the WBS pdf is symmetric, rightskewed and left-skewed depending on the parameter values. From Figure 2, we observe that the failure rate function of the WBS distribution can be increasing, decreasing, upside-down bathtub (unimodal) shaped, bathtub-shaped or modified bathtub shaped (unimodal shape followed by increasing) depending on the parameter values. So, the WBS distribution is quite flexible and can be used to fit various types of lifetime data sets in different fields.

Mathematical properties
Some characteristics of the WBS distribution can be found through an algebraic expansion which is more efficient than computing those directly by numerical integration of its pdf (7). In this section, using algebraic expansions, we give some mathematical properties of the WBS distribution.

Expansions for pdf and cdf
In this subsection, the expansions for the pdf and cdf of the WBS distribution are derived. We can write the pdf and cdf of the WBS distribution as a linear combination of the pdf and cdf of exponentiated BS (EBS) distribution, respectively. A random variable X is said to have a EBS distribution with parameters α, β and θ > 0, denoted by X ∼EBS(α, β, θ), if its cdf and pdf, respectively, are given by where g is given in (2). Several authors have studied the properties of exponentiated distributions. For example, Mudholkar and Srivastava [23] studied the exponentiated Weibull distribution, Gupta and Kundu [14] studied the exponentiated exponential distribution and Sarhan and Apaloo [28] proposed the exponentiated modified Weibull extension distribution. The pdf of the WBS distribution (7) can be written as follows Using the series expansion of the exponential function, we obtain Inserting (9) in (8), we obtain Since 0 < Φ(v) < 1, for t > 0 and (bk + b + 1) > 0, then by using the binomial series expansion of [ where Γ(·) is the complete gamma function. So, the pdf of the WBS distribution can be expressed as If b is a positive real non-integer, we can expand where Inserting (11) in (10), we get where and h r+1 (t) is the EBS(α, β, r + 1) density function. By integrating (12), we get Equation (12) means that the pdf of the WBS distribution is a linear combination of the pdf of EBS distribution. Based on this equation, several properties of the WBS distribution can be obtained directly from the EBS distribution.

Moments
We can use the moments to study the most important features and characteristics of a distribution such as dispersion, skewness and kurtosis. In this subsection, we derive the sth moment of the WBS random variable T from the probability weighted moments of the BS distribution. The probability weighted moments of the BS distribution are 66 THE WEIBULL BIRNBAUM-SAUNDERS DISTRIBUTION AND ITS APPLICATIONS defined, for p and r non-negative integers, by To compute the integral (15) numerically, we can use any softwares such as MAPLE, MATLAB and R. From [8], we have an alternative representation to compute τ p,r that is where (3) in terms of the modified Bessel function of the third kind. Therefore, the sth moment of T can be written from (12) as where τ s,r is obtained from (16) and d r is given by (13). We can compute numerically the sth moment in any symbolic software by taking in the sum a large number of summands instead of infinity.
In Figures 3 and 4, we plot the skewness and kurtosis, for selected values of a and b where α = 5 and β = 2.

Moment generating function
, is an alternative specification of its probability distribution. From (12), we obtain Using the series expansion for the exponential function, we obtain where τ s,r is obtained from (16) and d r is given by (13).

Quantile function
The WBS quantile function is obtained, by inverting F (·) given in (6), as follows where, Φ −1 (·) is the standard normal quantile function and

THE WEIBULL BIRNBAUM-SAUNDERS DISTRIBUTION AND ITS APPLICATIONS
Therefore, it is easy to simulate the WBS distribution. Let U be a continuous uniform variable on the unit interval (0, 1]. So, using the inverse transformation method, the random variable T given by has the WBS distribution. Equation (18) may be used to generate random numbers from the WBS distribution when the parameters are known. β, a, b). The mean deviations of T about the mean and about the median can be used as measures of spread in a population. They are given by

Mean deviations
respectively, where the mean µ 1 is calculated from (17) and the median m is given by m = Q W BS (1/2). The measures δ 1 and δ 2 can be expressed as (12) and (15), J(q) can be written as where s j and A(k 1 , ..., k j ) are defined in (16). Consider From [33], we can write where, K ω (u 1 , u 2 ) is the incomplete Bessel function with order ω and arguments u 1 and u 2 . Then, we obtain which can be calculated from the function D(p, q). So, we can use this expression for ϕ(q, r) to compute J(q). From (19), we obtain the Bonferroni and Lorenz curves defined by B(p) = J(q)/pµ 1 and L(t) = J(q)/µ 1 , respectively.

Reliability
In the stress-strength modelling, R = (T 2 < T 1 ) is a measure of component reliability when it is subjected to random stress T 2 and has strength T 1 . The component fails when T 2 < T 1 . The parameter R is referred to as the reliability parameter. In this subsection, we derive the reliability R when T 1 and T 2 have independent WBS(α, β, a 1 , b 1 ) and WBS(α, β, a 2 , b 2 ) distributions. The pdf of T 1 and the cdf of T 2 can be obtained from (12) and (14) as From (15), we can write where τ 0,2r is obtained from (16).

Order statistics
In this section, the distribution of the ith order statistic for the WBS distribution are presented. If T 1 , ..., T n are a random sample from WBS distribution with cdf (6) and pdf (7), and T 1,n ≤ · · · ≤ T n,n are the order statistics obtained from this sample, the pdf of T i,n is given by From (6), we have

THE WEIBULL BIRNBAUM-SAUNDERS DISTRIBUTION AND ITS APPLICATIONS
Using the binomial series expansion, we get Inserting (7) and (21) in (20), we obtain Using the power series for the exponential function, we get Therefore, the pdf of the ith order statistic for WBS distribution is where and h bl+b+m is the EBS density function with power parameter bl + b + m. Equation (22) means that the density function of the WBS order statistics is a linear mixture of the EBS densities. Then, we can easily obtain the mathematical properties for T i,n .
Using the normal approximation of the MLE of ξ, we can construct approximate confidence intervals for the parameters. Under some regular conditions (see [12]) that are fulfilled for parameters in the interior of the parameter space but not on the boundary, the asymptotic distribution of √ n , where J −1 (ξ) is the observed information matrix evaluated at ξ. The observed information matrix is given by where the elements are given in the Appendix. The approximate 100(1 − η)% two-sided confidence intervals for α, β, a and b are given by α ± z η/2 √ var(α), β ± z η/2 √ var(β), a ± z η/2 √ var(â) and b ± z η/2 √ var(b) respectively, where z η/2 is the quantile (1 − η/2) of the standard normal distribution and var(·) is the diagonal element of J −1 (ξ) corresponding to each parameter.

Simulation study
In this section, a simulation study is performed (by means of the statistical software R) to assess the finite sample behavior of the MLEs of the parameters of the WBS distribution. The inversion method is used to generate random samples from the WBS distribution using equation (18). The simulation study is repeated for N = 1000 times each with different sample sizes n = 20, 50, 100, 200, 500, and with the following cases for the true parameters The evaluation of the performance is based on the bias and the root mean squared errors (RMSE) defined as follows: and RM SE = 1 1000 The results of our simulation study are summarized in Table 1. We can see that the bias and RMSE of the MLEs converge to zero when the sample size is increased, as expected.

Applications
In this section, we show the flexibility of the WBS distribution by means of two well-known real data sets with different shapes. The first data set is given in [22] while the second data set is given in [34]. The first data set has a bathtub shaped failure rate function whereas the second data set has an increasing failure rate function. For these data sets, we compare the WBS distribution with the beta BS (BBS) [8], Kumaraswamy BS (KBS) [29], McDonald BS (McBS) [10], Marshall-Olkin extended BS (MOEBS) [20], gamma BS (GBS) [9], exponentiated generalized BS (EGBS) [9], transmuted BS (TBS) [6], EBS and BS distributions using the graphical method, minus twice the maximized log-likelihood (−2l), Akaike information criterion (AIC), Bayesian information criterion (BIC), Consistent akaike information criterion (CAIC) and Kolmogorov-Smirnov (K-S) test. The pdfs of the BBS, KBS, McBS, MOEBS, GBS, EGBS and TBS distributions (for t > 0) are, respectively, given by is the BS(α, β) pdf (2) and α, β, a, b, c > 0.

Meeker and Escobar data
This data set represents failure and running times of 30 devices provided by Meeker and Escobar [22]. The data set is: 2, 10, 13, 23, 23, 28, 30, 65 ,80, 88, 106, 143, 147, 173, 181, 212, 245, 247, 261, 266, 275, 293, 300, 300, 300, 300, 300, 300, 300, 300. The total time on test (TTT) plot for the Meeker and Escobar data in Figure 5(a) shows a convex shape followed by a concave shape. This corresponds to a bathtub-shaped failure rate. So, the WBS distribution is suitable for this data set. Table 2 Table 3, we observe that the WBS distribution has the largest p-value and the smallest −2l, AIC, BIC, CAIC and K-S values. So, the WBS distribution gives an excellent fit than the others models for Meeker and Escobar data set. Figure 6 shows of the data set and the estimated densities of all models. From this Figure, we observe that the WBS model is the closest to the empirical histogram than the other fitted models. Therefore, the WBS model could be chosen as the best model for Meeker and Escobar data. The probability-probability (P-P) plots of the fitted distributions, in Figure 7, confirm this result.

Turbochargers failure data
The data set represents the time to failure(10 3 h) of turbocharger of one type of engine given in [34]. The data set is: 1.6, 2.0, 2.6, 3.0, 3.5, 3.9, 4.5, 4.6, 4.8, 5.0, 5.1, 5.3, 5.4, 5.6, 5.8, 6.0, 6.0, 6.1, 6.3, 6.5, 6.5, 6.7, 7.0, 7.  Figure 5(b) shows a concave shape then the data has an increasing failure rate function. Table 4 gives the MLEs of the parameters of all models used here and their corresponding standard errors in parentheses. The statistics −2l, AIC, BIC, CAIC, K-S and its p-value are listed in Table 5. From this Table, we can see the WBS distribution as the best fit for the data set among all the seven models. The histogram of this data set and the plots of the estimated densities of all models are shown in Figure 8. So, the WBS model provides a better fits. This result is confirmed by the P-P plots in Figure 9.

Conclusion
In this paper, we have proposed a new four-parameter model, referred to as the Weibull Birnbaum-Saunders (WBS) distribution. The proposed distribution gives a more flexible model since it appropriate for modeling highly skewed and heavy-tailed data sets and and has the capability to capture decreasing, increasing, bathtub, unimodal (upside-down bathtub) and modified unimodal shaped hazard rates. The properties of the new distribution including expansions for the density function, moments, generating function, order statistics, quantile function, mean deviations and reliability are provided. We have estimated the model parameters by maximum likelihood and determined the observed information matrix. We have shown, by simulation study, that the maximum likelihood method performs well for estimating the parameters. We have also shown that the new model fits two well-known data sets better than existing other extensions of the BS distribution. A possible future works are to: ,