A New Burr XII-Weibull-Logarithmic Distribution for Survival and Lifetime Data Analysis : Model , Theory and Applications

A new compound distribution called Burr XII-Weibull-Logarithmic (BWL) distribution is introduced and its properties are explored. This new distribution contains several new and well known sub-models, including Burr XII-Exponential-Logarithmic, Burr XII-Rayleigh-Logarithmic, Burr XII-Logarithmic, Lomax-Exponential-Logarithmic, Lomax–Rayleigh-Logarithmic, Weibull, Rayleigh, Lomax, Lomax-Logarithmic, Weibull-Logarithmic, Rayleigh-Logarithmic, and Exponential-Logarithmic distributions. Some statistical properties of the proposed distribution including moments and conditional moments are presented. Maximum likelihood estimation technique is used to estimate the model parameters. Finally, applications of the model to real data sets are presented to illustrate the usefulness of the proposed distribution.


Introduction
Compound distributions have applications in various fields of study such as economics, engineering, public health, industrial reliability and medicine.New distributions have been developed by compounding well-known continuous distributions such as the exponential, Weibull, and exponentiated exponential distributions with the power series distribution that includes the Poisson, logarithmic, geometric and binomial distributions as particular cases [1,2].In recent years, compound distributions have been proposed because of their flexibility as they provide both monotonic and non-monotonic hazard rate functions that are encountered in real life.Amongst these is the Weibull-power series (WPS) distributions by Morais and Barreto-Souza [3].Silva et al. studied the extended Weibull power series family, which includes as special models the exponential power series and Weibull power series distributions [4].Silva and Cordeiro introduced a new family of Burr XII power series models [5].Oluyede et al. recently proposed a Log-logistic Weibull Poisson distribution that has applications in several areas including lifetime data, reliability and economics [6].
The primary motivation for this study is the development of a model that generalizes both the Burr XII and Weibull distributions and the modeling of lifetime data and other data types with a diverse model that takes into consideration not only shape, and scale but also skewness, kurtosis and tail variation.Hitherto, the nested and non-nested distributions such as gamma log-logistic Weibull (GLLoGW), beta modified Weibull (BetaMW),beta Weibull Poisson (BWP), gamma-Dagum (GD) and exponentiated Kumaraswamy Dagum (EKD) have lower precisions with high Sum of Squares (SS), Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) when compared to the new Burr XII-Weibull Logarithimic (BWL) distribution.
Weibull distribution is one of the most popular and well reputed to model failure time in life-testing and reliability theory.Nevertheless, hazard rate function (hrf) of Weibull distribution in modelling lifetime analysis has been reported to have monotonic behaviour and this is considered as a major shortcoming of the distribution [7].Tahir et al. [7] had suggested the need to search for some generalizations or modifications of Weibull distribution that can provide more flexibility in lifetime modeling since empirical hazard rate curves often exhibit non-monotonic shapes including bathtub, upside-down bathtub (unimodal) and others in real-life applications.
In addition, motivated by various applications of Burr XII, Weibull and logarithmic distributions in several areas including reliability, exponential tilting (weighting) in finance and actuarial sciences, as well as economics, where Burr XII distribution plays an important role in income, we construct and develop the statistical properties of this new class of generalized Burr XII-Weibull-type distribution called the Burr XII-Weibull-Logarithmic distribution and apply it to real lifetime data in order to demonstrate the usefulness of the proposed distribution.In this regard, we propose a new distribution, called the Burr XII-Weibull-Logarithmic (BWL) distribution.
Let X i , i = 1, . . ., N, be independent and identically distributed random variables from the Burr XII-Weibull distribution [8] whose cumulative distribution function (cdf) and probability density function (pdf) are given by respectively, for c, k, α, β > 0, and x ≥ 0. Now, let N be a discrete random following a power series distribution assumed to be truncated at zero, whose probability mass function is given by where C(θ) = ∑ ∞ n=1 a n θ n is finite, θ > 0, and {a n } n≥1 a sequence of positive real numbers.Let X (1) = min(X 1 , . . ., X N ).The conditional cdf of X (1) given N = n is given by The cdf of the Burr XII-Weibull Power Series (BWPS) class of distributions is the marginal cdf of X (1) , which is given by where x > 0, c > 0, k > 0, α > 0, β > 0 and θ > 0.
In this paper, we present the BWL distribution and our results are organized as followed.The BWL distribution and its quantile function functions are given in Section 2. In Section 3, moments and conditional moments are presented.The maximum likelihood estimates of the model parameters are given in Section 4. A Monte Carlo simulation study to examine the bias and mean square error of the maximum likelihood estimates are presented in Section 5. Section 6 contains applications of the new model to real data sets.A short conclusion is given in Section 7.

Burr XII-Weibull-Logarithmic Distribution
The BWL distribution is a special case of the BWPS distribution with C(θ) = − log(1 − θ) and a n = 1 n .From Equation (3), the cdf of the BWL distribution is given by for x > 0, c > 0, k > 0, α > 0, β > 0 and 0 < θ < 1.The corresponding pdf of the BWL distribution is given by The plots show that the BWL pdf can be decreasing or concave down with positive skewness as shown in Figure 1.

Quantile Function
In this sub-section, the quantile function of the BWL distribution is presented.The quantile function of the BWL distribution is obtained by solving the nonlinear equation by using numerical methods.Consequently, random numbers can be generated based on Equation (5).

Moments, Conditional Moments and Mean Deviations
In this section, moments, conditional moments, mean deviations, Lorenz and Bonferroni curves are given for the BWL distribution.Moments are necessary and crucial in any statistical analysis, especially in applications.Moments can be used to study the most important features and characteristics of a distribution (e.g., central tendency, dispersion, skewness and kurtosis).

Using the series expansions e
x k , the rth moment of the BWL distribution can be written as Let y = (1 + x c ) −1 , then the rth moment of the BWL distribution is given by:

Conditional Moments
For lifetime models, it may be useful to know about the conditional moments that is defined as E(X r | X > t).The rth conditional moment is given by where is the incomplete beta function.The mean residual lifetime function of the BWL distribution is given by E(X|X > t) − t.

Maximum Likelihood Estimation
Let X ∼ BWL(c, k, α, β, θ) and ∆ = (c, k, α, β, θ) T be the parameter vector.The log-likelihood = (∆) for a single observation of x of X is given by The first derivative of the log-likelihood function with respect to ∆ = (c, k, α, β, θ) T are given by .
The total log-likelihood function based on a random sample of n observations: x 1 , x 2 , ...., x n drawn from the BWL distribution is given by * n = ∑ n i=1 i (∆), where i (∆), i = 1, 2, ....., n is given by Equation (7).The equations obtained by setting the above partial derivatives to zero are not in closed form and the values of the parameters c, k, α, β, λ must be found by using iterative methods.The maximum likelihood estimates of the parameters, denoted by ∆ is obtained by solving the nonlinear equation ) T = 0, using a numerical method such as Newton-Raphson procedure.The Fisher information matrix is given by I ), i, j = 1, 2, 3, 4, 5 can be numerically obtained by MATLAB (version 9.1, Math Works Inc, Natick, Massachusetts, USA) or R Project for Statistical Computing (version 3.4.4,R Core Team, Vienna, Austria).The total Fisher information matrix nI(∆) can be approximated by For a given set of observations, the matrix given in Equation ( 8) is obtained after the convergence of the Newton-Raphson procedure in MATLAB or R software.

Asymptotic Confidence Intervals
In this sub-section, we present the asymptotic confidence intervals for the parameters of the BWL distribution.The expectations in the Fisher Information Matrix (FIM) can be obtained numerically.
Let ∆ = ( ĉ, k, α, β, θ) be the maximum likelihood estimate of ∆ = (c, k, α, β, θ).Under the usual regularity conditions and that the parameters are in the interior of the parameter space, but not on the boundary, we have: , where I(∆) is the expected Fisher information matrix.The asymptotic behavior is still valid if I(∆) is replaced by the observed information matrix evaluated at ∆, which is J( ∆).The multivariate normal distribution N 5 (0, J( ∆) −1 ), where the mean vector 0 = (0, 0, 0, 0, 0) T , can be used to construct confidence intervals and confidence regions for the individual model parameters and for the survival and hazard rate functions.That is, the approximate 100(1 − η)% two-sided confidence intervals for s, k, α, β and λ are given by:

Simulation Study
We study the performance and accuracy of maximum likelihood estimates of the BWL model parameters by conducting various simulations for different sample sizes and different parameter values.Equation (5) is used to generate random data from the BWL distribution.The simulation study is repeated N = 1000 times each with sample size n = 25, 50, 75, 100, 200, 400, 800, 1000 and parameter values I : c = 5.9, k = 0.9, α = 0.5, β = 0.6, θ = 0.7 and I I : c = 9.8, k = 2.5, α = 0.5, β = 0.4, θ = 0.7.The choice of N=1000 is based on the need to have a reasonably large N to yield a true sampling distribution of our data on which our estimates of the distribution are based.Four quantities are computed in this simulation study.
Table 1 presents the Average Bias, RMSE, CP and AW values of the parameters c, k, α, β, θ for different sample sizes.From the results, we can verify that as the sample size n increases, the RMSEs decay toward zero.We also observe that the bias decreases in general as the sample size n increases.
In addition, the average confidence widths decrease as the sample size increases.
The beta Weibull-Poisson (BWP) pdf [11] is given by for a, b, α, β, λ > 0, and x > 0. The GD pdf (Oluyede et al. [12]) is given by In addition, the pdf of the gamma log-logistic Weibull (GLLoGW) distribution (Oluyede et al. [13]) is given by The maximum likelihood estimates (MLEs) of the BWL model parameters ∆ = (c, k, α, β, θ) are computed by maximizing the objective function via the subroutine mle2 in R [14].The subroutine mle2 in R was applied and executed for wide range of initial values.The issues of existence and uniqueness of the MLEs are of theoretical interest and has been studied by several authors for different distributions including [15][16][17][18].
The estimated values of the model parameters (standard error in parenthesis), -2log-likelihood statistic, Akaike Information Criterion, AIC = 2p − 2 ln(L), Bayesian Information Criterion, BIC = p ln(n) − 2 ln(L), and Consistent Akaike Information Criterion, AICC = AIC + 2 p(p+1) where L = L( ∆) is the value of the likelihood function evaluated at the parameter estimates, n is the number of observations, and p is the number of estimated parameters are presented in each table of estimates.The goodness-of-fit statistics W * and A * , [19] are also presented in the table.These statistics can be used to verify which distribution fits better to the data.In general, the smaller the values of W * and A * , the better the fit.BWL distribution is fitted to the data set and these fits are compared to the fits using some of the submodels and several non-nested distributions given above.We can use the likelihood ratio test (LRT) to compare the fit of the BWL distribution with its sub-models for a given data set.For example, to test β = 1, the LR statistic is ω = 2[ln(L( ĉ, k, α, β, θ)) − ln(L( c, k, α, 1, θ))], where ĉ, k, α, β and θ are the unrestricted estimates, and c, k, α and θ are the restricted estimates.The LR test rejects the null hypothesis if ω > χ 2 , where χ 2 denote the upper 100 % point of the χ 2 distribution with one degree of freedom.
To obtain the probability plot, we plotted F(y (j) ; ĉ, k, α, β, θ) against j − 0.375 where y (j) are the ordered values of the observed data.The measures of closeness are given by the sum of squares

Waiting Time between Eruptions (Seconds)
The ocean swell produces spectacular eruptions of water through a hole in the cliff at Kiama, about 120 km south of Sydney, known as the Blowhole.The times at which 65 successive eruptions occurred from 1340 h on 12 July 1998 were observed using a digital watch.The data was collected and contributed by Jim Irish [20].Initial values for the BWL model in R are c = 0.1, k = 0.5, α = 0.5, β = 0.8, θ = 0.8.The parameter estimates, goodness-of-fit statistics and results for this data are given in Table 2.
The estimated variance-covariance matrix for the BWL distribution is given by: .
Plots of the fitted densities and histogram, observed probability versus predicted probability for waiting time between eruptions data are given in Figures 2 and 3, respectively.
The LRT statistics for testing H 0 : BE against H a : BWL and H 0 : LR against H a : BWL are 13.43 (p-value = 0.0012) and 24.22 (p-value < 0.0001), respectively.We conclude that there is a significant difference between the BE and the BWL distributions.There is also a significant difference between the LR and the BWL distributions.The LRT statistic for testing H 0 : BW against H a : BWL is 1.03 (p-value = 0.3102).We conclude that there is no significant difference between BW and BWL distributions; however, there is indeed clear and convincing evidence based on the goodness-of-fit statistics W * , A * and KS and its p-value that the BWL distribution is far better than the sub-models, and the non-nested models.In addition, the values of AIC and BIC shows that the BWL distribution is better than the non-nested GLLoGW, BetaMW, BWP, GD and EKD distributions.The values of SS from the probability plots is smallest for the BWL distributions when compared to the nested and non-nested models.

Time to Failure of Kevlar 49/Epoxy Strands Tested at Various Stress Levels
This real life example is taken from Cooray and Ananda [21], where 101 data points represent the stress-rupture life of kevlar 49/epoxy strands that are subjected to constant sustained pressure at the 90% stress level until all have failed, so that the complete data set with the exact times of failure is recorded.These failure times in hours are originally given in Andrews and Herzberg [22] and Barlow et al. [23].Initial values for the BWL model in R are c = 0.1, k = 0.1, α = 0.1, β = 4, and θ = 0.8.The parameter estimates, goodness-of-fit statistics and results for this data are given in Table 3.The estimated variance-covariance matrix for the BWL distribution is given by:  The LRT statistics for testing H 0 : BE against H a : BWL and H 0 : LEL against H a : BWL are 9.7 (p-value = 0.0078) and 24.16 (p-value < 0.0001), respectively.We conclude that there is a significant difference between the BE and the BWL distributions.There is also a significant difference between the LEL and the BWL distributions.The LRT statistic for testing H 0 : BW against H a : BWL is 9.6 (p-value = 0.0019), hence we conclude that there is a significant difference between BW and BWL distributions.There is clear and convincing evidence based on the goodness-of-fit statistics W * , A * , KS and its p-value, that the BWL distribution is far better than the sub-models and the non-nested models.In addition, the values of AIC and BIC show that the BWL distribution is better than the non-nested GLLoGW, BetaMW, BWP, GD and EKD distributions.The values of SS from the probability plots is the smallest for the BWL distributions when compared to the nested and non-nested models GLLoGW, BetaMW, BWP and EKD distributions.

Conclusions
A new generalized distribution called the Burr XII-Weibull-Logarithmic (BWL) distribution has been proposed and studied.Statistical properties including the moments and conditional moments were presented.A maximum likelihood estimation technique is used to estimate the model parameters.Finally, the BWL distribution is fitted to real data sets in order to illustrate its applicability and usefulness.We found BWL to be more precise than some other nested and non-nested models.

Figure 2 .
Figure 2. Fitted densities for waiting time between eruptions data.

Figure 3 .
Figure 3. Probability plots for waiting time between eruptions data.
fitted densities and histogram, observed probability versus predicted probability for waiting time between eruptions data are given in Figures4 and 5, respectively.

Figure 4 .
Figure 4. Fitted densities for time to failure of kevlar 49/epoxy strands data.

Figure 5 .
Figure 5. Probability plots for time to failure of kevlar 49/epoxy strands data.

Table 1 .
Monte Carlo simulation results: Average bias, Root Mean Square, Coverage Probability and Average Width.

Table 2 .
Model parameter estimates for the waiting time between eruptions data.

Table 3 .
Model parameter estimates for time to failure of kevlar 49/epoxy strands data.