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 [Formula: see text] and shape [Formula: see text] 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 [Formula: see text] parameter. A Weibull [Formula: see text] parameter estimator based on perturbation theory is proposed in this work. An explicit expression for [Formula: see text] 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.

Background

The two-parameter Weibull distribution is widely used in reliability engineering and lifetime data analysis because of its flexibility to properly model increasing and decreasing failure rates. It has gained the interest of researchers who have worked on its various aspects, such as inference, application and parameter estimation (see Nelson 1982; Cohen 1991; Johnson et al. 1994; Meeker and Escobar 1998). Traditional parameter estimation methods call on probability plotting, least squares and maximum likelihood estimation (Lawless 1982).

A probability plotting approach is straightforward and it is best used for small size data samples. However, this estimation method has not been sufficiently accurate as reported in Mao and Li (2007). The least squares (or rank regression) method is essentially a probability plotting method that applies least squares to determine lines through points. The main disadvantage of this method is that it assigns a large weight for extreme observations, producing a large variance (Genschel and Meeker 2010).

Maximum likelihood estimation (MLE) is considered one of the most robust parameter estimation techniques. It constructs a likelihood function for a set of statistical 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 right-censored 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 right-censored 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 , where , be a set of lifetime data, collected from N products or components, consisting of n observed ages of failed components and ages of surviving components, i.e., the so-called right-censored data. Let also if is the age of a failed component and if 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 with a probability density functionwhere and are the scale and shape parameters, respectively.

Therefore, the corresponding log-likelihood function can be recast as:whereThe 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: where 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 and . Hence, it can be seen thatafter multiplying Eq. (6) by . Here,In order to prove the existence and uniqueness of the solution of Eq. (9), the global monotonicity of and its asymptotic behavior, in the limits and , is developed. Let us consider that , where , 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 is continuous and monotonously increasing on , since the derivative of is positive for all :This implies that is also continuous and monotonously increasing on . On the other hand, is bounded and its boundaries can be obtained from the asymptotic behavior of in the limits and . In the limit , the asymptotic series expansions of is given by:andwhere is the big Landau O notation.

Substitution of these results into (10) yields the following asymptotic equation for in the limit :which states that as .

Now, it can be shown that as . Hence,due to (11) and noticing that for all and , as . Therefore, Eq. (10) and this last result, along with the fact that for , show that the asymptotic expansion for in the limit is given by:i.e., as .

It follows from the asymptotic equations (12) and (13) that is a continuous and monotonously increasing function for with the following asymptotic behavior: As a result, there exists a unique value such that since (14) and (15).

Let us determine the solution interval of . It follows from Eq. (9) that since . As a result . On the other hand, it can be seen that for , i.e., since is monotonously increasing. At the same time, . Therefore, , whereSummarizing, it can be concluded thatThe asymptotical behavior of (in the limits and ) and are illustrated in Fig. 1. However, the expression can be negative for some statistical data sets. In this case, the boundary estimation for becomes:It can be observed that this interval is too large. Then, a better estimation of the right boundary is required. For this, let us denote:where in virtue of (17). Then, Eq. (9) can be written as:

Fig. 1

Plot of the function

The solution is found from Eq. (20), which can be substituted into according to Eqs. (8) and (19). This way, can be calculated by substituting the estimated value into Eq. (5).

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 , it can be written as:whereThe simplest analytical solution to Eq. (21) can be found by truncating the power series in (21) at . As a result, the following expression can be obtained:The solution to this equation is given by:and subsequently denotewhere (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 can be recast as:where means summation over all such that . 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 , 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 . 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 . Twelve grinders failed in this period of time. Values of , , are obtained from Eq. (22): , , , , . The Weibull parameter is estimated from Eq. (32) and these values of . 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).

Table 1

Lifetime data for Case 1

k	t_k	δ_k	k	t_k	δ_k	k	t_k	δ_k	k	t_k	δ_k
1	12.5	1	6	95.5	1	11	125.6	1	16	152.7	0
2	24.4	1	7	96.6	1	12	152.7	1	17	152.7	0
3	58.2	1	8	97.0	1	13	152.7	0	18	152.7	0
4	68.0	1	9	114.2	1	14	152.7	0	19	152.7	0
5	69.1	1	10	123.2	1	15	152.7	0	20	152.7	0

Table 2

Parameter estimates for Case 1

Parameters	Approximate analytical method	Approach in Balakrishnan and Kateri (2008)	NR method β_0=2
β	1.6466	1.647	1.6466 (6 iterations)
α	162.22306	162.223	162.2330

The MLE was combined with the Newton–Raphson method using a convergence tolerance set to 0.001 and an initial guess , which corresponds to our solution when it is rounded up to the next integer. This last criterion is adopted from the well-known 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 and . 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.

Table 3

Simulated data set provided by Balakrishnan and Mitra (2012)

Simulated data set provided by Balakrishnan and Mitra (2012) No.	Installation year	Failure year	No.	Installation year	Failure year	No.	Installation year	Failure year	No.	Installation year	Failure year
1	1984	–	26	1986	–	51	1982	–	76	1974	2006
2	1990	2001	27	1987	–	52	1981	–	77	1978	1995
3	1983	2002	28	1990	1997	53	1986	–	78	1962	1993
4	1981	2000	29	1980	1996	54	1980	1990	79	1963	–
5	1985	–	30	1980	–	55	1980	1994	80	1960	1998
6	1991	–	31	1981	–	56	1982	–	81	1962	2007
7	1982	–	32	1983	1997	57	1990	2008	82	1960	1990
8	1990	–	33	1980	–	58	1985	–	83	1962	1980
9	1983	1999	34	1984	–	59	1983	–	84	1961	1981
10	1992	–	35	1982	–	60	1982	–	85	1964	1989
11	1983	–	36	1980	–	61	1963	1996	86	1964	1987
12	1989	–	37	1985	2007	62	1963	2001	87	1960	2006
13	1985	–	38	1993	–	63	1961	1998	88	1961	1992
14	1982	–	39	1983	–	64	1961	1992	89	1964	–
15	1983	–	40	1980	–	65	1960	1984	90	1963	1991
16	1981	–	41	1981	2001	66	1964	2004	91	1973	–
17	1985	–	42	1989	–	67	1961	1994	92	1964	–
18	1981	–	43	1993	–	68	1977	1998	93	1972	1984
19	1988	2002	44	1983	–	69	1963	1987	94	1962	2007
20	1983	–	45	1993	–	70	1960	1991	95	1963	1997
21	1984	–	46	1987	–	71	1961	1983	96	1964	1987
22	1989	–	47	1994	–	72	1964	1995	97	1964	2002
23	1988	–	48	1985	2007	73	1963	1998	98	1971	–
24	1982	–	49	1981	–	74	1961	2001	99	1965	1990
25	1981	–	50	1983	2004	75	1960	1988	100	1962	1994

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.

Table 4

Parameter estimates for Case 2

Parameters	Approximate analytical method	NR method β_0=4
β	3.205	3.2506 (8 iterations)
α	35.245	35.2084

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 , , where M denotes the total number of units. Each 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 . Each group of simulated lifetime data required M numbers of that were randomly generated from a two-parameter Weibull distribution with prespecified values of and . Censoring times were chosen to have a common value, which is calculated as , 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 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 and . Our analytical MLE method provides estimates that are closer to and 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.

Table 5

Case 3. Parameter estimates for different simulated data sets

CR (%)	M	L-Moments method		Approximate analytical method
		α	β	α	β
0	10	3.0522	2.1521	3.0563	2.5286
	20	3.4026	1.3887	3.4289	1.4851
	50	2.7332	1.3915	2.7112	1.3497
	100	2.9970	1.5129	2.9860	1.4877
	500	3.1178	1.5677	3.0670	1.4772
	1000	3.0152	1.5289	2.9624	1.4342
20	10	2.6830	1.9116	2.9353	1.7227
	20	2.9035	1.8576	3.3553	1.4175
	50	2.5856	1.8822	2.8145	1.5906
	100	2.8119	2.1616	2.9557	1.9031
	500	2.7591	1.8930	3.0847	1.4296
	1000	2.7483	1.9360	3.0350	1.5349
80	10	1.0279	4.5619	1.4874	2.1709
	20	1.0871	12.9564	2.4157	2.2982
	50	1.0709	7.8137	3.0147	1.4884
	100	1.0553	6.1198	3.4455	1.2138
	500	1.0757	9.3024	2.8698	1.6711
	1000	1.0709	8.2974	2.9130	1.5363

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 , which were randomly generated from a two-parameter Weibull distribution with prespecified values of , 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):andwhere is the number of simulations and and 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.

Fig. 2

Biases of , and mean-squared errors of , for and

Fig. 3

Biases of , and mean-squared errors of , for and

Fig. 4

Biases of , and mean-squared errors of , for and

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 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 efficiently estimates employing a single equation, with no need of graphical or iterative procedures.

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.