The Log-Logistic Weibull Distribution with Applications to Lifetime Data

In this paper, a new generalized distribution called the log-logistic Weibull (LLoGW) distribution is developed and presented. This distribution contain the log-logistic Rayleigh (LLoGR), log-logistic exponential (LLoGE) and log-logistic (LLoG) distributions as special cases. The structural properties of the distribution including the hazard function, reverse hazard function, quantile function, probability weighted moments, moments, conditional moments, mean deviations, Bonferroni and Lorenz curves, distribution of order statistics, Lmoments and Rényi entropy are derived. Method of maximum likelihood is used to estimate the parameters of this new distribution. A simulation study to examine the bias, mean square error of the maximum likelihood estimators and width of the confidence intervals for each parameter is presented. Finally, real data examples are presented to illustrate the usefulness and applicability of the model.


Introduction
There are several generalizations of univariate distributions including those of (Eugene, Lee, and Famoye 2002) dealing with the beta-normal distribution, as well general family of univariate distributions generated from the Weibull distribution that was introduced by Gurvich, Dibenedetto, and Ranade (1997).The cumulative distribution function (cdf) given by (Gurvich et al. 1997) where C is a subset of R, and H(x; Θ) is a non-negative monotonically increasing function that depends on the vector of parameters Θ.The corresponding probability density function (pdf) is given by g(x; α, Θ) = α exp[−αH(x; Θ)]h(x; Θ), where h(x; Θ) is the derivative of H(x; Θ).The choice of the function H(x; Θ) can lead to The Log-logistic Weibull Distribution with Applications to Lifetime Data different models including for example, exponential distribution with H(x; Θ) = x, Rayleigh distribution is obtained from H(x; Θ) = x 2 and Pareto distribution from setting H(x; Θ) = log(x/k).
There are several ways of generating new probability distributions from classic ones to relative new distributions in the literature.Nelson mentioned in (Nelson 1982) that distributions with bathtub-shaped failure rate are sufficiently complex and, therefore, difficult to model.The distribution proposed by (Hjorth 1980) is such an example.Later on, (Rajarshi and Rajarshi 1988) presented a revision of these distributions, and (Haupt and Schäbe 1992) put forward a new lifetime model with bathtub-shaped failure rates.Unfortunately, these models are not sufficient to address various practical situations, so new classes of distributions were presented based on the modifications of the Weibull distribution to satisfy non-monotonic failure rate.For a review of these models, the reader can refer to (Mudholkar and Srivastava 1993), and (Pham and Lai 2007), where the authors summarized some generalizations of Weibull distribution in their papers.Other generalizations include the exponentiated Weibull (EW) (Gupta and Kundu 2001), the modified Weibull (MW) (Lai, Xie, and Murthy 2003), and the beta exponential (BE) (Nadarajah and Kotz 2006).Some more recent extensions are the generalized modified Weibull (GMW) (Carrasco, Ortega, and Cordeiro 2008), the beta modified Weibull (BMW) (Silva, Ortega, and Cordeiro 2010), (Nadarajah, Cordeiro, and Ortega 2011), the Weibull-G family (Bourguignon, Silva, and Cordeiro 2014) and the Gamma-exponentiated Weibull distributions (GEW) (Pinho, Cordeiro, and Nobre 2012).(Gurvich et al. 1997) developed a new statistical distribution for characterizing the random strength of brittle materials.
To motivate the model under study, consider a series system and assume that the lifetime of the components follow the log-logistic and Weibull distributions with with reliability functions R 1 (t) = (1 + ( t s ) c ) −1 and R 2 (t) = e −αt β , respectively.The reliability R(t) = P (T > t) of the system is given by (3) In some context, a series model is referred to as a competing risk model.
Also, a primary motivation for developing this model is the advantages presented by this generalized distribution with respect to having a hazard function that exhibits increasing, decreasing and bathtub shapes, as well as the versatility and flexibility of the log-logistic and Weibull distributions in modeling lifetime data.We propose and study this new distribution called the log-logistic Weibull distribution which inherits these desirable properties and also covers quite a variety of shapes.
There is an added advantage to this model, in that it has an additional dispersion parameter, depending on the overall form that accounts for the scale of the underlying random variable.
The distribution also has exponential dumping in the upper tail making the distribution suitable for modeling samples that display power behavior for intermediate observations and decrease in tail probability for large observations or beyond a certain threshold or specified value.
The proposed new distribution generalizes the log-logistic and Weibull distributions.Some structural properties of this distribution are obtained and estimation the parameters via the method of maximum likelihood presented.
This paper is organized as follows.In section 2, we present the generalized distribution including the corresponding probability density functions (pdf), hazard and reverse hazard functions, quantile function and various sub-models.In section 3, the probability weighted moments, moments and conditional moments are presented.Section 4 contain the derivation of the mean deviations, Bonferroni and Lorenz curves.Section 5 is concerned with Rényi entropy, distribution of order statistics and L-moments.Estimation of model parameters is presented in section 6. Monte Carlo simulation study is conducted in section 7 to examine the bias and mean square error of the maximum likelihood estimators and the width of the confidence intervals for each parameter.Applications of the proposed model to real data are given in section 8, followed by concluding remarks.

The log-logistic Weibull distribution
In this section, we present some statistical properties of the new log-logistic Weibull (LLoGW) distribution, including pdf, cdf, quantile function, hazard and reverse hazard functions.Plots of the hazard rate function for selected values of the model parameters are also given.We first of all present the Burr-XII, log-logistic and Weibull distributions.The very popular Burr Type III and Type XII distributions attract special attention because they include several families of distributions (e.g., the gamma distribution) with varying degrees of skewness and kurtosis.Further, these distributions have applications in a wide variety of areas in statistics and applied mathematics including modeling events associated with fracture roughness, life testing, operational risk, option market price distributions, forestry, meteorology, modeling crop prices, software reliability growth, and reliability analysis.See (Burr 1942), (Burr 1973) for additional details.The cdf and pdf of Burr XII distribution are given by and , for s, c, k, and x ≥ 0, respectively.The reliability and hazard rate functions are given by , respectively.Note that the pdf is unimodal with mode at x 0 = ((c − 1)/(ck + 1)) 1/c when c > 1, and L−shaped when c = 1.The r th non-central moment is given by Note that k and c are shape parameters and s is a scale parameter.When k = 1 we obtain the log-logistic distribution.The cdf of the well known Weibull (W) distribution is given by where α and β are scale and shape parameter, respectively.Now, consider the log-logistic Weibull (LLoGW) distribution obtained by taking for s, c, α, β > 0 and x ≥ 0. If a random variable X has the LLoGW cdf, we write X ∼ LLoGW (s, c, α, β).The corresponding LLoGW pdf is given by s, c, α, β > 0, and x ≥ 0. Plots of the pdf for selected values of the model parameters are given in Figure 1.The plots suggests that the LLoGW pdf can be right skewed or decreasing for the selected values of the parameters.
The quantile function of the LLoGW distribution is obtained by solving the equation using numerical methods.Consequently, random number can be generated based on equation (8).Table 1 lists the quantile for selected values of the parameters of the LLoGW distribution.

Some new and known sub-models
There are several new as well as well known distributions that can be obtained from the LLoGW distribution.Note that when s = m −1 , we have the log-logistic Weibull (LLoGW) distribution with the survival function S(x; m, c, α, β) = [1 + (xm) c ] −1 e −αx β .When c = 1 it reduces to the generalized Pareto type II Weibull (GP-II-W) distribution.
• If c = 1, then the LLoGW cdf reduces to the three parameter distribution with cdf given by for s, α, β > 0, and x ≥ 0.
• If c = β = 1 then the LLoGW cdf reduces to to the two parameter distribution given by for s, α > 0, and x ≥ 0.
• If c = 1 and β = 2, then the LLoGW cdf reduces to the two parameter model for s, α > 0, and x ≥ 0.

Hazard and reverse hazard functions
In general, if X is a continuous random variable with cdf F, and pdf f, then the hazard function, reverse hazard function and mean residual life function are given by respectively.The functions λ F (x), δ F (x), and F (x) are equivalent.See (Shaked and Shanthikumar 1994) and references therein.In this subsection, the hazard and reverse hazard functions of the LLoGW distribution are presented.The hazard and reverse hazard functions of the LLoGW distribution are for x ≥ 0, s, c, α, β > 0, respectively.The limiting behavior of the hazard function of the LLoGW distribution, which can be readily established is as follows: The Log-logistic Weibull Distribution with Applications to Lifetime Data • For β = 1, and for each c > 0, lim x→∞ h G (x) = α.
Plots of the hazard function are given in Figure 2. The graphs exhibit increasing, decreasing, bathtub, upside down bathtub, and upside down bathtub followed by bathtub shapes for the selected values of the model parameters.This very attractive flexibility makes the LLoGW hazard function useful and suitable for non-monotonic empirical hazard behaviors which are more likely to be encountered in practice or real life situations.

Probability weighted moments, moments and conditional moments
In this section, we obtain probability weighted moments (PWMs) (Greenwood, Landwehr, Matalas, and Wallis 1979), moments and conditional moments for the LLoGW distribution.The probability weighted moments (PWMs) of the LLoGW distribution is given by Now, by setting y = (1 + (x/s) c ) −1 , the PWMs of the LLoGW distribution can be written as: −1 dy .
Consequently, the PWMs of the LLoGW distribution is given by Remarks: special cases • When l = m = 0, we obtain the r th non-central moment µ r given by • When l = 0, the LLoGW PWMs reduces to • When m = 0, the LLoGW PWMs reduces to • When r = m = 0, we have • When r = 0, we have Note that the r th raw moment µ r can also be obtained as follows: We apply the following results which holds for m and k positive integers, w > −1 and p > 0 ( (Prudnikov, Brychkov, and Marichev 1992), page 21): where ∆(k, a) = a k , a+1 k , ...., a+k k , and the Meijer G function is defined by x −t dt, where i = √ −1 is the complex unit and L is the integration path (see (Gradshtein and Ryzhik 2000), sec.9.3 for a description of this path).
Consequently, the r th moment of the LLoGW distribution is given by .
Alternatively, the following direct computation also gives the r th moment of the LLoGW distribution when we use the fact that e for c > kβ+β+r.To obtain the moment generating function (MGF) of the LLoGW distribution, recall the Taylor's series expansion of the function e tx = ∞ j=0 (tx) j j! , so that, we have Table 2 lists the first six moments together with the standard deviation (SD or σ), coefficient of variation (CV), coefficient of skewness (CS) and coefficient of kurtosis (CK) of the LLoGW distribution for selected values of the parameters, by fixing α = 1.5 and β = 1.5.And Table 3 lists the first six moments, SD, CV, CS and CK of the LLoGW distribution for selected values of the parameters, by fixing s = 1.5 and c = 1.0.These values can be determined numerically using R and MATLAB.The SD, CV, CS and CK are given by σ Table 2: Moments of the LLoGW distribution for some parameter values; α = 1.5 and β = 1.5.

Conditional moments
For lifetime models, it is also of interest to find the conditional moments and the mean residual lifetime function.The r th conditional moments for LLoGW distribution is given by where B x (a, b) = x 0 y a−1 (1 − y) b−1 dy is the incomplete beta function, and c > kβ + β + r.Alternatively, consider the following integral for q > 0 and b > 0, (Paranaíba, Ortega, Cordeiro, and Pescim 2011), where 2 F 1 is the hypergeometric function given by and (a) k = a(a + 1)....(a + k − 1) is the ascending factorial.Now, consider the integral Consequently, the r th conditional moment of the LLoGW distribution is given by The mean residual lifetime function of the LLoGW distribution is E(X|X > t) − t.

Mean deviations
The amount of scatter in a population is evidently measured to some extent by the totality of deviations from the mean and median.These are known as the mean deviation about the mean and the mean deviation about the median, and defined by respectively, where µ = E(X) is the mean and M =Median (X) is the median.The measures δ 1 (x) and δ 2 (x) can be calculated using the relationships and respectively.When r = 1, we get the mean µ = E(X).Note that and respectively.Alternatively, Consequently, the mean deviation about the mean is and the mean deviation about the median is

Bonferroni and Lorenz curves
In this subsection, we present Bonferroni and Lorenz Curves.Bonferroni and Lorenz curves have applications not only in economics for the study income and poverty, but also in other fields such as reliability, demography, insurance and medicine.Bonferroni and Lorenz curves for the LLoGW distribution are given by and respectively, where T (q) = ∞ q xg(x)dx, and q = G −1 (p), 0 ≤ p ≤ 1.

Order statistics, L-moments and Rényi entropy
The concept of entropy plays a vital role in information theory.The entropy of a random variable is defined in terms of its probability distribution and can be shown to be a good measure of randomness or uncertainty.In this section, we present the distribution of the order statistic, L-moments and Rényi entropy for the LLoGW distribution.

Order statistics
Order statistics play an important role in probability and statistics.In this subsection, we present the distribution of the i th order statistic from the LLoGW distribution.The pdf of the i th order statistic from the LLoGW pdf g(x) is given by where g m+i (x) is the pdf of the exponentiated LLoGW distribution with parameters s, c, α, β and (m + i), B(., .) is the beta function and the weights w i,m are given by The t th moment of the i th order statistics from the LLoGW distribution can be derived via a result of (Barakat and Abdelkader 2004) as follows: Note that where B(a, b) = 1 0 t a−1 (1 − t) b−1 dt is the beta function.

Rényi entropy
In this subsection, Rényi entropy of the LLoGW distribution is derived.An entropy is a measure of uncertainty or variation of a random variable.Rényi entropy is an extension of Shannon entropy and is defined to be Rényi entropy tends to Shannon entropy as v → 1.Note that [g(x; s, c, α, β)] v = g v (x) can be written as Consequently, Rényi entropy is given by for v = 1, and v > 0.

Maximum likelihood estimation
Let X ∼ LLoGW (s, c, α, β) and ∆ = (s, c, α, β) T be the parameter vector.The log-likelihood The Log-logistic Weibull Distribution with Applications to Lifetime Data = (∆) for a single observation x of X is given by The first derivative of the log-likelihood function with respect to ∆ = (s, c, α, β) T are given by and The total log-likelihood function based on a random sample of n observations: x 1 , x 2 , ...., x n drawn from the LLoGW distribution is given by * , where i (∆), i = 1, 2, ....., n is given by equation ( 25).The equations obtained by setting the above partial derivatives to zero are not in closed form and the values of the parameters s, c, α, β must be found by using iterative methods.The maximum likelihood estimates of the parameters, denoted by ∆ is obtained by solving the nonlinear equation ( ∂ ∂s , ∂ ∂c , ∂ ∂α , ∂ ∂β ) T = 0, using a numerical method such as Newton-Raphson procedure.The Fisher information matrix is given by ∂θ i ∂θ j ), i, j = 1, 2, 3, 4, can be numerically obtained by MATLAB, SAS or R software, where (θ 1 , θ 2 , θ 3 , θ 4 ) = (s, c, α, β).The total Fisher information matrix nI(∆) can be approximated by For a given set of observations, the matrix given in equation ( 26) is obtained after the convergence of the Newton-Raphson procedure.

Asymptotic confidence intervals
In this subsection, we present the asymptotic confidence intervals for the parameters of the LLoGW distribution.The expectations in the Fisher Information Matrix (FIM) can be obtained numerically.Let ∆ = (ŝ, ĉ, α, β) be the maximum likelihood estimate of ∆ = (s, c, α, β).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 ∆, that is J( ∆).The multivariate normal distribution N 4 (0, J( ∆) −1 ), where the mean vector 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, c, α, and β are given by: respectively, where I −1 ss ( ∆), I −1 cc ( ∆), I −1 αα ( ∆), and I −1 ββ ( ∆), are the diagonal elements of I −1 n ( ∆) = (nI( ∆)) −1 , and Z η 2 is the upper η 2 th percentile of a standard normal distribution.

Simulation study
In this section, we study the performance and accuracy of maximum likelihood estimates of the LLoGW model parameters by conducting various simulations for different sample sizes and different parameter values.Equation ( 8) is used to generate random data from the LLoGW distribution.The simulation study is repeated N = 5, 000 times, each with sample size n = 25, 50, 75, 100, 200, 400 and parameter values I : s = 5.5, c = 2.5, α = 0.8, β = 0.2 and II : s = 5.5, c = 8.5, α = 0.5, β = 0.5.Four quantities are computed in this simulation study.
(b) Average bias of the MLE θ of the parameter ϑ = s, c, α, β : (b) Root mean squared error (RMSE) of the MLE θ of the parameter ϑ = s, c, α, β : (c) Coverage probability (CP) of 95% confidence intervals of the parameter ϑ = s, c, α, β, i.e., the percentage of intervals that contain the true value of parameter ϑ.
Table 4 presents the Average Bias, RMSE, CP and AW values of the parameters s, c, α and β 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 for all the parametric values, the biases decrease as the sample size n increases.Also, the table shows that the coverage probabilities of the confidence intervals are quite close to the nominal level of 95% and that the average confidence widths decrease as the sample size increases.Consequently, the MLE's and their asymptotic results can be used for estimating and constructing confidence intervals even for reasonably small sample sizes.

Applications
In this section, we present examples to illustrate the flexibility of the LLoGW distribution and its sub-models for data modeling.Estimates of the parameters of LLoGW distribution (standard error in parentheses), Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), Cramer von Mises (W * ), Andersen-Darling (A * ), and sum of squares (SS) from the probability plots are presented for each data set.We also compare the LLoGW distribution with the gamma-Dagum (GD) (Oluyede, Huang, and Pararai 2014) and beta Weibull (BW) (Lee, Famoye, and Olumolade 2007), (Famoye, Lee, and Olumolade 2005) distributions.The GD and BW pdfs are given by and respectively.
The maximum likelihood estimates (MLEs) of the LLoGW parameters ∆ = (s, c, α, β) are computed by maximizing the objective function via the subroutine NLMIXED in SAS as well as the function nlm in R (R Development Core Team 2011).The estimated values of the parameters (standard error in parenthesis), -2log-likelihood statistic, Akaike Information Criterion, AIC = 2p − 2 ln(L), and Bayesian Information Criterion, BIC = p ln(n) − 2 ln(L), 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 Tables 5 and 6.The LLoGW distribution is fitted to the data sets and these fits are compared to the fits using LLoGR, LLoGE and LLoG distributions.
As stated earlier, we maximize the likelihood function using NLmixed in SAS as well as the function nlm in R (R Development Core Team 2011).These functions were applied and executed for wide range of initial values.This process often results or lead to more than one maximum, however, in these cases, we take the MLEs corresponding to the largest value of the maxima.In a few cases, no maximum was identified for the selected initial values.In these cases, a new initial value was tried in order to obtain a maximum.The issues of existence and uniqueness of the MLEs are theoretical interest and has been studied by several authors for different distributions including (Seregin 2010), (Silva and Tenreyro 2010), (Zhou 2009), and (Xia, Mi, and Zhou 2009).
We can use the likelihood ratio (LR) test to compare the fit of the LLoGW distribution with its sub-models for a given data set.For example, to test β = 1, the LR statistic is ω = 2[ln(L(α, β, ŝ, ĉ)) − ln(L( α, 1, s, c))], where α, β, ŝ, and ĉ are the unrestricted estimates, and α, s, and c 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 1 degrees of freedom.
Plots of the fitted densities, the histogram of the data and probability plots (Chambers, Cleveland, Kleiner, and Tukey 1983) are given in Figures 3 and 4 for the first data set and Figures 5 and 6 for the second dataset.For the probability plot, we plotted G(x (j) ; ŝ, ĉ, α, β) against j − 0.375 n + 0.25 , j = 1, 2, • • • , n, where x (j) are the ordered values of the observed data.The measure of closeness given by the sum of squares SS = n j=1 G(x (j) )− j − 0.375 n + 0.25 2 , was computed for each fitted model.
The goodness-of-fit statistics W * and A * , described by (Chen and Balakrishnan 1995) are also presented in the tables.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.Let G(x; ∆) be the cdf, where the form of G is known but the k-dimensional parameter vector, say ∆ is unknown.We can obtain the statistics W * and A * as follows: (i) Compute u i = G(x i ; ∆), where the x i 's are in ascending order; (ii) Compute y i = Φ −1 (u i ), where Φ(.) is the standard normal cdf and Φ −1 (.) its inverse; (iii) Compute v i = Φ((y i − y)/s y ), where y = n −1 n i=1 y i and s 8.1.Fracture toughness of alumina (Al 2 O 3 )(in the units of MPa m 1/2 ) This data set consists of 119 observations on fracture toughness of Alumina (Al 2 O 3 )(in the units of MPa m 1/2 .The data are taken from the web-site: http://www.ceramics.nist.gov/srd/summary/ftmain.htm.The same data set has also been studied by (Nadarajah and Kotz 2007).Initial values for the LLoGW model in R code are s = 2.34, c = 3.42, α = 0.155, β = 1.Estimates of the parameters of LLoGW distribution and its related sub-models (standard error in parentheses), AIC, BIC, W * , A * and SS are give in Table 5. Plots of the fitted densities and the histogram are given in Figure 3, and plots of the observed probability vs predicted probability are given in Figure 4.The LR test statistic of the hypothesis H 0 : LLoGE against H a : LLoGW, H 0 :LLoG against H a : LLoGW, and H 0 :LLoGR against H a : LLoGW are 15.3 (p-value < 0.0001), 22.0 (p-value < 0.0001), and 4.9 (p-value=0.02686< We can conclude that there are significant differences between LLoGW and LLoGE, LLoG, LLoGR distributions.The values of the statistics: AIC and BIC also shows that the LLoGW distribution is a better fit than the non-nested GD and BW distributions for the fracture toughness of alumina data.There is also clear evidence based on the goodness-of-fit statistics W * and A * that the LLoGW distribution is by far the better fit for the fracture toughness of alumina data.

Breaking stress of carbon fibres (in Gba)
This data set consists of 100 uncensored data on breaking stress of carbon fibres (in Gba), (Nichols and Padgett 2006).Initial values for the LLoGW model in R code are s = 1, c = 1, α = 1, β = 3.Estimates of the parameters of LLoGW distribution and its related sub-models (standard error in parentheses), AIC, BIC, W * , A * and SS are give in Table 6.Plots of the fitted densities and the histogram, as well as observed probability vs predicted probability are given in Figures 5 and 6 The LR test statistic of the hypothesis H 0 : LLoGE against H a : LLoGW and H 0 :LLoG against H a : LLoGW are 7.5791 (p-value = 0.0059) and 10.1860 (p-value = 0.0014).We can conclude that there are significant difference between LLoGW and LLoGE distributions as well between LLoGW and LLoG distributions.There is also very clear and convincing evidence based on the goodness-of-fit statistics W * and A * that the LLoGW distribution is by far the better fit than the sub-models.The SS value of 0.0584 for the LLoGW distribution is smaller than the values for the non-nested GD and BW distributions.The values of AIC and BIC also shows that the LLoGW distribution is a better fit than the non-nested GD and BW distributions for the breaking stress of carbon fibres data.

Concluding remarks
We have presented a new distribution called the log-logistic Weibull (LLoGW) distribution that is suitable for applications in various areas including reliability, survival analysis and actuarial

Figure 3 :
Figure 3: Fitted densities for fracture toughness of alumina data

Figure 5 :
Figure 5: Fitted densities for breaking stress of carbon fibres data

Table 4 :
Monte Carlo simulation results: average bias, RMSE, CP and AW

Table 5 :
Estimates of models for fracture toughness of alumina dataNote.Standard errors are in parentheses.

Table 6 :
Estimates of models for breaking stress of carbon fibres dataNote.Standard errors are in parentheses.

The Graph of Observed vs Expected Probability
Figure 6: Probability plots for breaking stress of carbon fibres data Barakat HM, Abdelkader YH (2004)."Computing the Moments of Order Statistics from Nonidentical Random Variables."Statistical Methods and Applications, 13(1), 15-26.