Perturbative method for maximum likelihood estimation of the Weibull distribution parameters

The two-parameter Weibull distribution is the predominant distribution in reliability and lifetime data analysis. The classical approach for estimating the scale \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\alpha )$$\end{document}(α) and shape \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\beta )$$\end{document}(β) parameters employs the maximum likelihood estimation (MLE) method. However, most MLE based-methods resort to numerical or graphical techniques due to the lack of closed-form expressions for the Weibull \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}β parameter. A Weibull \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}β parameter estimator based on perturbation theory is proposed in this work. An explicit expression for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}β is obtained, making the estimation of both parameters straightforward. Several right-censored lifetime data sets with different sample sizes and censoring percentages were analyzed in order to assess the performance of the proposed estimator. Study case results show that our parameter estimator provides solutions of high accuracy, overpassing limitations of other parameter estimators.

data, which is optimized to find its extremum with respect to the distribution parameters. The MLE method can handle survival and interval data better than rank regression approaches, particularly when dealing with heavily censored data sets that contain few points of highly accurate observed data. Teimouri et al. (2013) compares the MLE method with other four methods [the Method of Logarithm Moment (MLM), the Percentile Method (PM), the L-Moments Method (LM), and the Method of Moments (MM)] to determine Weibull parameters. One of the main findings of this work is that estimation of parameters is better performed using MLE and LM estimators. However, MLE leads to likelihood equations that need to be solved numerically. Therefore, low convergence rates and efficient iterative methods must be properly addressed, which can be particularly difficult with censored data (Balakrishnan and Kateri 2008).
Recent research has been focused on obtaining new efficient numerical and statistical inference methods in order to deal with this problem. Joarder et al. (2011) consider statistical inferences of the unknown parameters of the Weibull distribution with rightcensored data samples, stating that the MLE cannot lead to explicit forms of the Weibull distribution. Therefore, they propose approximate maximum likelihood estimators (AMLE), which are obtained by expanding the MLE equations in Taylor series. Also, the authors propose a fixed-point algorithm to compute the maximum likelihood estimators. Balakrishnan and Mitra (2012) use an expectation-maximization (EM) algorithm to estimate the model parameters of the Weibull distribution of left-truncation and rightcensored data. The algorithm consists of two steps: expectation step (E-step) and maximization step (M-step). The conditional expectation of the complete data likelihood is obtained with the E-step, using the incomplete observed data and current estimated value of the parameter. This expected likelihood is essentially a function of the involved parameter and its current value under which the expectation has been calculated. The expected likelihood is then maximized with respect to the parameter using the EM gradient algorithm. The E-and M-steps are then iterated until convergence. MLE and Bayes estimators are applied to calculate the survival function and the failure rate of the Weibull distribution for censored data in Guure and Ibrahim (2012). In order to estimate the survival and the failure rate functions under the MLE, the authors applied the Newton-Raphson method. Bayes estimators are obtained using a linear exponential, general entropy and squared error loss functions while a prior noninformative Bayesian approach is employed to estimate the survival function and failure rate. However, the aposteriori distribution function cannot be reduced to a closed form because it involves a ratio of complicated integrals. More work concerning Weibull parameter estimation can be found in Jabeen et al. (2013), Yang and Scott (2013), Guure and Ibrahim (2014), Mohammed Ahmed (2014) and Wang and Ye (2015).
Most parameter estimation methods presented in the literature are useful tools for solving practical problems, showing that the Weibull parameter estimation problem continues to be important in the research field of data analysis. Hence, it is clear that the development of general and new methods for a wider range of applications is desirable.
In this paper, an approximate analytical method to estimate the β Weibull parameter for complete and right-censored data is proposed using perturbation theory. The method involves a systematic construction of an analytical solution to the likelihood equation for β, taking advantage of the presence of a small parameter. The solution is developed as a power series with respect to this parameter. As a result, the likelihood equation for β is replaced by a set of simple solvable algebraic equations. These equations are explicitly solved one by one in order to obtain an increasingly accurate approximation to the true solution.

Problem statement
Let t k |k = 1, N , where 0 < t 1 ≤ t 2 ≤ · · · ≤ t N , be a set of lifetime data, collected from N products or components, consisting of n observed ages of failed components and N − n ages of surviving components, i.e., the so-called right-censored data. Let also δ k = 1 if t k is the age of a failed component and δ k = 0 if t k is the age of a surviving component, so the number of failed components is: Let us assume that the lifetime data set follows the Weibull distribution W (α, β) with a probability density function where α and β are the scale and shape parameters, respectively. Therefore, the corresponding log-likelihood function can be recast as: where The MLE method states that the most probable values of α and β correspond to the extremum of (3), or equivalently to the existence and uniqueness of the solution (α * , β * ) of the following system equations: The system of Eq. (4) can be expressed as: ( β * is then estimated solving Eq. (6). Equation (5) is the closed-form expression for the MLE of α. Thus, α * can be calculated by substituting the estimated value of β * into Eq. (5). β is usually solved for employing numerical or graphical methods, which involve inaccuracies and numerical problems. So, the main aim of this research is to obtain an analytical solution to Eq. (6), allowing the shape parameter β to be explicitly found.

Existence and uniqueness of the likelihood estimate
The existence and uniqueness of the solution of the MLE equation have already been proved in Balakrishnan and Kateri (2008) and Farnum and Booth (1997) using Cauchy-Schwarz inequality. A different proof is presented here, leading to our proposed analytical solution for the β parameter. Let us denote: where 0 < x k ≤ 1 and x N = 1. Hence, it can be seen that after multiplying Eq. (6) by β 1 . Here, In order to prove the existence and uniqueness of the solution of Eq. (9), the global monotonicity of Z(ζ ) and its asymptotic behavior, in the limits ζ → 0+ and ζ → +∞, is developed. Let us consider that 0 < x 1 ≤ x 2 ≤ · · · ≤ x N = 1, where N ≥ 2, is a monotonous nondecreasing data sequence. Also, let us suppose that there must be at least two different statistical data sets, i.e., It can be proved that r(ζ ) is continuous and monotonously increasing on ζ > 0, since the derivative of r(ζ ) is positive for all ζ > 0: This implies that Z(ζ ) is also continuous and monotonously increasing on ζ > 0. On the other hand, r(ζ ) is bounded and its boundaries can be obtained from the asymptotic behavior of r(ζ ) in the limits ζ → 0+ and ζ → +∞. In the limit ζ → 0+, the asymptotic series expansions of ζ is given by: and where O(ζ ) is the big Landau O notation. Substitution of these results into (10) yields the following asymptotic equation for r(ζ ) in the limit ζ → 0+: due to (11) and noticing that x k = x N = 1 for all m < k ≤ N and 0 < x m < 1, x ζ m → 0 as ζ → +∞. Therefore, Eq. (10) and this last result, along with the fact that ln x k = 0 for k > m, show that the asymptotic expansion for r(ζ ) in the limit ζ → +∞ is given by: i.e., r(ζ ) → 0 as ζ → +∞.
It follows from the asymptotic equations (12) and (13) that Z(ζ ) is a continuous and monotonously increasing function for ζ > 0 with the following asymptotic behavior: As a result, there exists a unique value ζ * such that Z(ζ * ) = 0 since (14) and (15).

A perturbative approach to estimate the shape parameter
Perturbation theory is employed in this section to solve Eq. (20). It allows the representation of ζ * to be asymptotically expanded, which in turn can be conveniently truncated to obtain an analytical solution to Eq. (20).
After expanding each term of Eq. (20) in Taylor series and dividing the entire equation by the expression N k=1 x k (ln x k ) 2 + ln x k + 1 , it can be written as: where The simplest analytical solution to Eq. (21) can be found by truncating the power series in (21) at m = 1. As a result, the following expression can be obtained: The solution to this equation is given by: and subsequently denote where 0 < −σ 0 < 1 (see "Appendix"). Therefore, the exact solution to Eq. (21) can be expanded in a power series with respect to ε, so that a solution to z can be written as: Noticing that z m can be recast as: where q 1 ···q m m k=1 q k =Q means summation over all q 1 , q 2 , . . . , q m such that m k=1 q k = Q . Substitution of (25) into (21), and taking on account of (23), Eq. (21) has the following expression: Each coefficient of the power series in the left-hand-side of Eq. (26) is zero according to the perturbation theory. As a result, the following infinite system of algebraic equations can be obtained: All coefficients of the asymptotic series (24) can be determined by iteratively solving the system of equations (27). Thus, the complete asymptotic series (24) is fully established. Generally, higher-order terms in the series (24) become successively smaller for ε small. Therefore a good approximation is obtained when the power series is truncated using a few terms. For example, if series (27) is truncated at term Q = 4, the system becomes: A consistent solution to (28) can be obtained by successively solving each of its equations: Hence, substitution of these coefficients into Eq. (24) yields the asymptotic solution to Eq. (26): The asymptotic solution to ζ * in Eq. (9) is obtained by substituting (30) into (19): Finally, it can be seen from (8) that β = β 1 · ζ. The approximate analytical solution for β * is then obtained as: In turn, the scale parameter α * can be determined by substituting the result of (32) into Eq. (5).

Cases of study
Three study cases are shown in this section to illustrate the application of our proposed analytical method for the estimation of the Weibull β parameter. The first study considers right-censored data set found in Balakrishnan and Kateri (2008), where a graphical solution for the determination of the MLE shape parameter is employed. In a second study, the proposed method is applied to right-censored data used in Balakrishnan and Mitra (2012). Finally, sets of lifetime data are randomly generated combining different censoring rates and sample sizes, in order to cover a wider range of data sampling scenarios that might be encountered in practical applications. Corresponding Weibull parameters for each data set are accordingly estimated. For the first two cases, the β parameter was also estimated using a Newton-Rapshon algorithm with the purpose of illustrating the advantage of our proposed analytical method, where β is obtained by a single equation.

Case 1
The censored data set provided by Dodson (2006) and analyzed by Balakrishnan and Kateri (2008) is shown in Table 1. It provides twenty identical grinders which were tested with a ending time t = 152.7. Twelve grinders failed in this period of time. Values of σ m , m = 0, 4, are obtained from Eq. (22): σ 0 = −0.2244, σ 1 = 1.0, σ 2 = −0.3397, σ 3 = 0.3617 , σ 4 = −0.3320. The Weibull β parameter is estimated from Eq. (32) and these values of σ m . The result is substituted into Eq. (5) to obtain α. Results are presented in Table 2 along with the parameters estimated by a graphical method in Balakrishnan and Kateri (2008). The MLE was combined with the Newton-Raphson method using a convergence tolerance set to 0.001 and an initial guess β 0 = 2, which corresponds to our β solution when it is rounded up to the next integer. This last criterion is adopted from the wellknown fact that an appropriate initial value (close to the desired solution) ensures the convergence of the NR method within few iterations.
It can be observed from Table 2 that the parameters obtained from our analytical method match closely the estimates of the NR method and graphical method proposed in Balakrishnan and Kateri (2008). Therefore, it can be stated that our proposed parameter estimation methodology not only works well for this case, but it also directly provides β from an explicit expression [Eq. (32)].

Case 2
In this case, the data set is provided by Balakrishnan and Mitra (2012) and reproduced in Table 3. It is used to assess again the efficacy of our proposed parameter estimation method. This data set was numerically produced using the Weibull distribution with parameter values α true = 35 and β true = 3. It can be considered as lifetime data of power transformers in the electrical industry, with information of installation and failure dates. The data considers units that failed before year 2008, which was set as the year of censoring.
The parameter values obtained for this case are presented in Table 4. Similarly to Case 1, selection for the initial value of the NR method is the round up integer of the β solution obtained using our analytical method, with a tolerance of 0.001.
It can be observed again that the parameters obtained with our analytical method match quite well the estimates of the NR program. This case is another example that shows the efficacy of our proposed analytical method for the determination of Weibull parameters.

Case 3
The simulation approach of Zhou et al. (2013) was adopted to generate different sets of lifetime data at prespecified points of time. The generated data mimic a time-censored sampling scenario for a hypothetical number of transformer units, which operate at the same time. In addition, it is assumed that a population of units is homogenous with a fixed censoring time C, and each individual unit has a lifetime t k , k = 1, M, where M denotes the total number of units. Each t k is identically considered an independent random variable that follows a specific probability distribution. Finally, lifetime data is characterized by the censoring rate CR, defined as the proportion of censored data and calculated as the number of suspensions divided by the sample size.
The sample sizes employed in this study were 10, 20, 50, 100, 500 and 1000. Censoring rates were fixed at 0, 20 and 80 %. Each group of simulated lifetime data required M  numbers of t k that were randomly generated from a two-parameter Weibull distribution with prespecified values of α true = 3.0 and β true = 1.5. Censoring times were chosen to have a common value, which is calculated as F −1 x (p; α, β), where p is the probability of a unit, starting at time 0, fails before reaching censoring time C. p was fixed for each CR at 1.0, 0.8 and 0.2. Then, a lifetime data set is generated through the comparison of lifetime units and a selected censoring time: If t k is less than or equal to C, the unit is failed. Otherwise, the unit is in suspension with lifetime data censored at time C.
This study was specially designed to bring about the effectiveness of our proposed analytical MLE method for Weibull parameters. Our proposed method was also compared in this work with the L-Moments estimation method presented in Teimouri et al. (2013), which is based on linear combination of order statistics and provides closed-form expressions for Weibull parameters.
Weibull parameters were obtained for each simulated data set using our analytical MLE method and the L-Moments method (see Table 5). It can be observed from these results that both methods provide a close match to α true and β true . Our analytical MLE method provides estimates that are closer to α true and β true for censored data sets. Therefore, we can state from these results that the L-Moments method is effective for  complete data sets, whereas our analytical MLE method is effective for complete data sets as well as right-censored data sets. An additional analysis was carried out to compare both methods in terms of the mean-squared error and bias. For this purpose, the bias and the mean-squared error were computed for data sets of size M = 10, 20, . . . , 1000, which were randomly generated from a two-parameter Weibull distribution with prespecified values of (α true , β true ) = (0.5, 0.5), (1.0, 1.0), (3.0, 1.5), and the number of simulations for each data set was 100. The mean-squared error and bias are defined as follows (Teimouri et al. 2013): and where η is the number of simulations and α * i and β * i are estimators of α and β for the i-th simulation, respectively. Figures 2, 3 and 4 illustrate the behavior of the mean-squared error and bias for α * and β * . It can be concluded from these results that the L-Moments method has more bias in the estimates of β * , whereas our analytical MLE method has more bias in the estimates of α * . However, the bias and the mean-squared error in both methods tend to consistently decrease with the increase of the size of the data sets. The difference between both methods, with respect to the mean-squared error criterion, is very small, mainly for large M where the difference is practically indistinguishable. In addition, the amount of bias and mean-squared error behaves in a consistent manner over different values of the parameters.

Conclusions
An analytical approach is developed in this work to estimate the Weibull parameter β for complete and right-censored data using perturbation theory. The idea behind this method is to formally expand the β solution to its likelihood equation around point 1.0 as a power series in ε, which turns out to be a small parameter. In fact, if ε is zero, the equation is exactly solvable. Therefore, the problem is reduced to find the asymptotic behavior of the best approximation to the true solution within ε, ε 2 , . . . Thus, perturbation theory leads to an expression for the desired solution in terms of a formal power series in a "small" parameter that quantifies the deviation from the exactly solvable problem. Hence, an approximate analytical solution for β parameter is obtained by truncating the series at a prespecified order.
Our analytical method for estimation of the Weibull parameter was tested on several lifetime data sets. This way, it was concluded that the performance of our proposed method was satisfactory for all lifetime data sets with different combinations of sample sizes with small and heavy censoring. These data sets cover a wide range of practical scenarios that our method can easily deal with.
The main conclusion that can be drawn from this work is that the use of the formulations described in "Existence and uniqueness of the likelihood estimate" and "A perturbative approach to estimate the shape parameter" sections allows the analytical obtention of β. This method was not only numerically tested using common and practical data sets, but it was also theoretically and mathematically proved. Our approach Finally, it is worth mentioning that although the estimation of Weibull parameters under right-censored scheme was considered, the proposed method can be extended to other censoring schemes such as left-truncation and hybrid. Additional work is required in this direction, which is currently considered by the authors.