On the Bias of the Maximum Likelihood Estimators of Parameters of the Weibull Distribution

Usually, the parameters of a Weibull distribution are estimated by maximum likelihood estimation. To reduce the biases of the maximum likelihood estimators (MLEs) of two-parameter Weibull distributions, we propose analytic bias-corrected MLEs. Two other common estimators of Weibull distributions, least-squares estimators and percentiles estimators, are also introduced. Based on a comparison of their performances in the simulation study, we strongly recommend the analytic bias-corrected MLEs for the parameters of Weibull distributions, especially when the sample size is small.


Introduction
We investigate several different estimations of the parameters of the two-parameter Weibull distribution, and analyze two different methods to reduce the biases of maximum likelihood estimators (MLEs).
The following formula is the cumulative distribution function of a Weibull distribution: and its probability density function is: The Weibull distribution was first identified by Fréchet in 1927, and it was named after Waloddi Weibull, a Swedish mathematician who described this distribution in detail in 1951.Rammler [1] used it to explain a particle size distribution.
The Weibull distribution is not a single distribution, but is in fact a family of distributions, because of the distribution having the shape parameter.The shape parameter enables Weibull distributions to take on various shapes, contingent on the value of the shape parameter.Those distributions are particularly useful in modeling applications since they are flexible enough to model various datasets.
The Weibull distribution has a relatively simple form.Nevertheless, the shape parameter enables the Weibull to assume a variety of shapes.This combination of simplicity and flexibility in the shape of the Weibull distribution makes it an effective distributional model in reliability applications.The ability to model various distributional shapes by using a relatively simple distributional form is possible with many other distributional families, as well.
The Weibull distribution is associated with various other distributions.For instance, it is well known that a Weibull distribution contains the exponential distribution (when k = 1) and the Rayleigh distribution (when k = 2).Recently, Ling and Giles [2] studied the Rayleigh distribution and the bias adjustment of the Rayleigh distribution.The Weibull distribution is a special condition of the generalized extreme value distribution, which was proposed by Fréchet [3].This closely-related Fréchet distribution has the probability density function: Actually, f Frechet (x; k, λ) = − f Weibull (x; −k, λ).A Weibull distribution is also characterized in terms of a uniform distribution; when U is uniformly distributed on (0, 1), then the random variable λ(− ln(U)) 1/k is Weibull distributed with parameters k and λ.
The Weibull distribution has been used in numerous fields, such as production, industry, data analysis, reliability in devices, materials testing and environmental assessment.This distribution was studied and used in an aircraft system to investigate operational reliability by Kaltschmidt et al. [4].Nikolaj [5] showed that this distribution can be used in oil pollution examination.Singh [6] studied the Weibull distribution and worked in hydrology by using this distribution.Aarset [7] applied a Weibull distribution to test failure times of devices, and the data were further analyzed and modified by Lai et al. [8].It was used in breakdown voltage estimation by Hirose [9], and Fabiani [10] used it to test electrical breakdown of insulating materials.It was applied to censored data by Ghitany et al. [11].In the study of wind energy, the Weibull distribution was also applied by Akda ǧ et al. [12].
Additionally, we have also considered the exponential distribution.When the value of k is one, the Weibull distribution becomes the exponential distribution, which is studied widely.The cumulative distribution function of the exponential distribution is: and its probability density is: Maximum-likelihood estimators of the Weibull distribution are commonly used.However, they are biased, and thus, in this paper, we propose analytic bias-corrected MLEs for the two-parameter Weibull distribution.The analytic bias-corrected MLEs are considered in Section 2 and compared with the other bias-corrected MLEs, which are also known as bootstrap bias-corrected MLEs.In Sections 3 and 4, we introduce two other commonly-used estimators.Monte Carlo experiments are conducted for the simulation study.Tabular and graphical representations of the results are presented in this section.The results indicate that the analytic bias-corrected estimators are superior to not only traditional bootstrap bias-corrected estimators, but also the other two commonly-used estimators.Thus, the analytic bias-corrected estimators should be highly recommended.Finally, we apply these estimators to an actual situation and draw conclusions.

Maximum Likelihood Estimators
Suppose the sample (x 1 , x 2 , ...x n ) is generated from the log-normal distribution.The likelihood function of the Weibull distribution is: and the formula for the log-likelihood function is: The MLEs can be obtained by solving the following equations: The Newton-Raphson method is used to solve these two equations, and solutions k and λ were obtained.

Analytic Bias-Corrected Maximum Likelihood Estimators
Suppose l(σ) is the log-likelihood function of the Weibull distribution based on a sample size n.The joint cumulants of the derivatives of the log-likelihood function are: In this part, the function l(σ) is regular with respect to all derivatives up to the third order.We assume that Equations ( 10)-( 13) are O(n).Cox and Snell [13] showed that for sample data that are independent, but not always identically distributed, the bias of the s-th element of the MLE of σ is: where h ij is the (i, j)-th element of the inverse of Fisher's information matrix, H = {−h ij }.Cordeiro and Klein [14] observed that the bias expression still holds even if the observations are non-independent, proved that all of the h terms are O(n) and recommended another form for the expression: In Equation (15), we define: and ( 15) can be written as: where vec(H −1 ) is a transformation of the matrix H −1 into a vector form, and the matrix A is: A = [A (1) Then, the bias-corrected MLEs are: where H = (H)| σ and A = (A)| σ .The bias of σ can be proven to be O(n −2 ).The values of the elements of σ can be obtained based on the maximum likelihood equations.This method was applied to the generalized Pareto distribution by Giles et al. [15] successfully, and Min Wang [16] used it for the weighted Lindley distribution with Wentao Wang.
Specific calculation procedures of the bias-corrected MLEs of the parameters of the Weibull distribution are listed in Appendix A.

Bootstrap Bias-Corrected Estimators
In this section, the bootstrap bias-corrected estimators of the Weibull distribution are considered, which were suggested by Mackinnon and Smith [17].
The solutions k and λ are the particular values of estimators based on MLEs, and we regard k and λ as the particular parameters of the log-normal distribution, from which the N B replicative samples are generated.Thus, the corresponding estimates j=1 λ j (the parameters k j and λ j are the MLEs of k and λ from the j-th sample of the N B replication samples).
This method is fairly effective and is worth applying to various situations, which causes that many experts have investigated it.Ferrari and Cribari-Neto [18] explored the relationship between bias correction of maximum likelihood estimators through the bootstrap method, and the bootstrap method was presented by Ho and Fernandes da Silva [19] to make improvements in the MLE of Mean Time to Failure (MTTF) and p-quantiles.Focarelli [20] used the bootstrap bias-correction procedure to estimate long-run relationships from dynamic panels, with an application to money demand in the euro area.

Percentile Estimators
The estimators based on percentiles are also investigated.This approach was first studied by Kao [21,22] and then was explored by Mann, Schafer and Singpurwalla [23].This method is suitable for the Weibull distribution.
Let us consider the situation where both the parameters λ and k are unknown.From the cumulative distribution function (1), we can obtain: Suppose that (x 1 , x 2 , ...x n ) is a random sample of data, and the (x (1) , x (2) , ...x (n) ) is the ordered observation.Let p i be the estimate of F(x (i) ; k, λ), then we can obtain estimates of λ and k by minimizing the equation: Generally, some p i 's are used to be the estimates of the F(x i ; k, λ), and the p i = (i/(n + 1)) is most frequently used, because (i/(n + 1)) is the expected value of F(x (i) ; k, λ).In this article, we also use this p i .
Equation ( 22) is a non-linear function.To obtain estimates of the two parameters, we can solve ( 22) by using some non-linear regression techniques.These estimators are called percentile estimators, which were also applied and introduced by Gupta and Kundu [24].

Least Squares Estimators
In this section, we introduce the expressions of another estimation.This method was first applied by Swain, Venkatraman and Wilson [25] to estimate the parameters of the beta distribution, and it can be utilized in other cases, as well.Suppose that (x 1 , x 2 , ..., x n ) is a random sample of size n from a distribution whose distribution function is F(x) and that (x (1) , x (2) , ..., x (n) ) denotes the ordered set of sample observations.Jonson, Kotz and Balakrishnan [26] proposed expressions for the expectations and variances, The least squares estimators can be obtained by minimizing: We call this estimator as least squares estimator.Therefore, in the case of the Weibull distribution, the least squares estimators of k and λ are denoted by k Lse and λ Lse , respectively.These values are obtained by minimizing:

Procedures
In this part, the Monte Carlo experiment is described, and the simulation results are analyzed.The R 3.2.2software [27] routines for log-normal random variates are used to generate the sample from the Weibull distribution.Then, the MLE, percentile estimation and least squares estimation can be realized by the non-linear minimization function.In the R software, the non-linear minimization function can solve a minimization problem using a Newton-type algorithm.After substituting the MLEs into the corrected formula mentioned previously, we can obtain the analytic bias-corrected estimators.The bootstrap bias-corrected estimators will be acquired after conducting the replicative MLE experiments.

Comparison of Maximum Likelihood Estimators and Bias-Corrected Maximum Likelihood Estimators
Comparing the theoretical performances of the different estimators mentioned above is not convincing.Therefore, we performed extensive simulations to compare the performances of the different methods, mainly with respect to their real absolute percentage biases and percentage mean squared errors (MSEs) for different sample sizes.The percentage MSE is defined as 100 × (MSE( k)/k 2 ), and the absolute percentage bias is defined as 100 × |(mean( k) − k)/k|.
In each section of this experiment, 20,000 Monte Carlo replications are performed, and 10,000 bootstrap samples are used for every replication.
We also investigate general Weibull distributions with some common parameter values.The exponential distribution with particular parameter value (k = 1) is also investigated here.We select parameters values for the Weibull distribution (k = 1) that are consistent with the data generating process being an exponential distribution.
Absolute percentage biases and MSEs of those estimators for the exponential distribution are presented in Table 1, and those for the Weibull distribution are presented in Table 2.
These tables show that the percentage biases of the analytic bias-corrected MLEs are much smaller than those of the MLEs.Meanwhile, both the analytic bias correction and bootstrap bias correction perform well.However, the Cox-Snell bias correction is more effective than the bootstrap reduction, especially with the modest sample size.
The absolute values of the percentage biases tend to decline as the sample size increases, and the percentage MSEs of the Cox-Snell estimators are smaller than those of the bootstrap estimators, except in a few cases.The MSEs and biases of the analytic bias-correction estimators are smaller than those of the other estimators.

Analysis for All Estimators
The absolute values of the percentage biases of all five estimators are shown in the Figures 1-3.Some of the points are clear based on the experimental results.For all methods, as the sample size increases, the biases decrease.For the most part, the biases of the percentile estimators are larger than those of other estimators, except in specific cases, as shown in Figures 1a and 2a.The bootstrap bias-corrected MLEs and analytic bias-corrected MLEs are generally both superior to MLEs.Furthermore, the analytic bias-corrected MLEs are superior to the other estimators.
The performance of the least squares estimators is uncertain.Sometimes, the least squares estimators perform better than the MLEs, whereas at other times, the latter's performance stands out.Clearly, the biases of the least squares estimators become extremely similar to those of the bias-corrected MLEs as the sample size increases.
However, the MSEs of the least squares estimators are often larger than those of the other estimators.Meanwhile, the MSEs of the MLEs, analytic bias-corrected MLEs, bootstrap bias-corrected MLEs and percentile estimators are close, as evidenced by the results shown in Figures 4-6.These results indicate that the analytic bias-corrected MLEs can effectively reduce the biases of MLEs.Moreover, to obtain the least biased estimators of the parameters of the Weibull distribution, the analytic bias-corrected MLEs are much more appropriate than the other common estimators mentioned above.

Real Illustrative Example
In this section, the practical performance of the Cox-Snell bias correction is explored, and we apply this estimation to real data (i.e., the data provided by Aarest [7]).Aarest tested the times to failure of 50 devices.
We fit the Weibull distribution to these data using MLEs, percentile estimators and least squares estimators.Then, we obtain the bias-corrected parameters by solving the analytic expressions mentioned previously.Table 3 reports the failure time of 50 devices.Table 4 reports the maximum likelihood estimates, bootstrap maximum likelihood estimates, analytic bias-corrected maximum likelihood estimates, percentile estimates and least squares estimates after fitting the Weibull distribution to these data.Based on the MLE of the data shown in Table 3, Figure 7 presents the fitted probability density function curve of Weibull distribution.For these data, the analytic bias correction decreases the height of the Weibull distribution density, and the shape of the density becomes slightly "fatterr" and "shorter".Thus, the distribution is more uniform.

Conclusions
Usually, the two parameters of a Weibull distribution are estimated using maximum likelihood estimation, which is realized under non-linear and first-order conditions.Nevertheless, after comparing all of the methods, it is clear that in terms of the minimization, analytic bias-corrected MLEs exhibited the best performance in almost all cases.Therefore, the analytic bias-corrected MLEs of the two-parameter Weibull distribution are proposed to more effectively reduce the biases of the MLEs.
We conduct a Monte Carlo experiment, which is described in the simulation study section.According to the MSEs and absolute percentage biases of the estimates, we explore the effectiveness of these estimations.Clearly, the biases of the analytic bias-corrected MLEs are smaller than those of the other estimators.Consequently, we can conclude that the analytic bias-corrected MLEs are superior to the others tested here.
All in all, the analytic bias-corrected MLEs discussed in this paper are worth popularizing for utilization on actual data.They should be highly recommended, especially for small sample sizes.

= − kn
x k i log λ (A10) and then, we can obtain the following formulas: The following equations are also easy to get: and the following equation is required for the next calculation part: Where γ is a constant called eulergamma and its value is about 0.57721; and ζ(3) is another constant, and its value is about 1.20206; then, these results can be obtained: Then, based on Equation ( 16), a series of results is obtained: From the formula H ij = {−h ij }, we can know that the information matrix is:

Table 2 .
Percentage bias and MSEs for the Weibull distribution.

Table 3 .
Calculated failure time of 50 devices.

Table 4 .
Five estimates of Parameters.