Abstract

The maximum likelihood estimation is a widely used approach to the parameter estimation. However, the conventional algorithm makes the estimation procedure of three-parameter Weibull distribution difficult. Therefore, this paper proposes an evolutionary strategy to explore the good solutions based on the maximum likelihood method. The maximizing process of likelihood function is converted to an optimization problem. The evolutionary algorithm is employed to obtain the optimal parameters for the likelihood function. Examples are presented to demonstrate the proposed method. The results show that the proposed method is suitable for the parameter estimation of the three-parameter Weibull distribution.

1. Introduction

Due to its extremely high flexibility, the Weibull distribution [1] is widely used for fitting engineering data, such as the strength of materials [2, 3], fracture of brittle materials [4, 5], and wind speed [6, 7]. The three-parameter Weibull distribution is very flexible for random data fitting so that it has a strong adaptability for different types of probability distribution. When the three parameters are well chosen, it can be equal or approximate to some other statistical distributions. However, it is vitally important to estimate the three parameters precisely in order to utilize the Weibull model successfully.

It is a known fact that the three-parameter Weibull distribution parameters estimation is an extremely complicated procedure. To date, several estimation approaches have been conducted, such as graphic method [810], moment estimation [1114], maximum likelihood estimation (ML) [1517], kernel density estimation [18, 19], least squares estimation with particle swarm optimization [20], ML with genetic algorithm [21], Bayesian approach [22, 23], grey model [24, 25], etc.

The maximum likelihood method is an effective and important approach for parameter estimation. The estimation procedure is converted to maximize the so-called likelihood function with respect to three undetermined parameters. However, it is complicated to solve the maximum likelihood equations by conventional numerical methods. Therefore, various methods have been studied for parameter estimation of the three-parameter Weibull distribution. Örkcü et al. [26] proposed an approach based on the differential evolution algorithm to search the maximum value of the likelihood function. Usta et al. [27] combined the probability weighted moments and the power density method for estimating the Weibull parameters. In [28], the particle swarm optimization (PSO) was adopted to provide accurate estimations of the Weibull parameters. Petkovic et al. [29] proposed an adaptive neuro-fuzzy inference system to predict the parameters of the Weibull distribution. Mazen et al. [30] discussed various approaches including maximum likelihood method, percentile based method, least squares method, weighted least squares method, and maximum product of spacing estimators by numerical simulations. However, all these methods depend heavily on the sample data.

The Covariance Matrix Adaptation Evolution Strategy (CMA-ES) algorithm [3133] is a new global optimization algorithm, which has a powerful searching ability and high efficiency. In this article, the maximum likelihood estimation combined with evolutionary algorithm is proposed to obtain the estimates of the Weibull parameters. The parameter estimation problem is converted to an optimization problem and is solved by the CMA-ES method. The three unknown parameters are taken as variables in the optimization problem and the fitness for optimization is the likelihood function. The advantage of this method is that it does not need to solve the maximum likelihood equations, and we can obtain the estimates of the three parameters directly by the optimization process. The results are compared with other estimation methods and literatures.

The rest of this article is organized as follows. Section 2 gives a short review of maximum likelihood method for Weibull distribution. In Section 3, the maximum likelihood estimation and the evolution optimization are presented. Section 4 presents discussions and comparisons with other estimation methods of two cases. Finally, conclusions are given in Section 5.

2. Maximum Likelihood Estimation for Three-Parameter Weibull Distribution

The cumulative distribution function (CDF) and probability density function (PDF) of the three-parameter Weibull distribution are given by

Here, , , and are location, shape, and scale parameters, respectively. Figure 1 illustrates shapes of PDF for different parameters. As can be seen, the shape of Weibull PDF is very flexible so that it can fit into a wide range of experiment data.

In general, there are many methods to estimate the parameters of a distribution, such as probability-weighted moment, maximum likelihood method, and least square method. Among them, the ML estimators are asymptotically unbiased with the minimum variance under regularity conditions. Then, the MLE for three-parameter Weibull distribution is described briefly.

Let be a random sample of size ; is noted as the Weibull model parameters which are to be estimated; namely, . The likelihood function is written asThe aim of estimation is to determine the unknown vector and the three unknown parameters , , and by maximizing the likelihood function. Since it contains exponential term, it is easier to obtain the maximum by its logarithm. By this way, the complexity of calculations is reduced. The logarithm of the likelihood function is shown asThen, the vector can be obtained by maximizing of the likelihood function. To achieve this, the conventional approach is to take the partial derivation of the likelihood function in terms of vector and set the partial equations to zero, asWe substitute the log-likelihood function into the above equations. The following equations are obtained:

It is well known that it is difficult to obtain the estimates of the unknown parameters by solving the above equations with numerical methods. Inspired by the ideas of heuristic algorithm, the direct optimization searching methods for unknown parameters based only on the evaluation of the objective function are proposed. The optimization strategy is used to overcome the difficulties of solving the above equations. In this study, the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) is used because it is robust and effective for global optimization.

3. Evolution Optimization

3.1. CMA-ES Algorithm

In this section, the Covariance Matrix Adaptation Evolution Strategy is briefly introduced for three-parameter Weibull distribution. The CMA-ES searches for a minimization of an objective function without any requirement of objective function derivatives. It has well performance on nonlinear and badly conditioned problems because the CMA-ES is a second-order approach estimating a positive definite matrix within an iterative procedure. The searching is taken by controlling the step and direction of evolution to generate new individuals in iteration. The CMA-ES starts with a given or random initial search points. The next generation is got by mutation, competition, selection, and recombination. Such repetition goes on until the best solution is achieved.

The mutation direction of group depends on the step size and covariance matrix . The basic function of CMA-ES algorithm is shown aswhere is the individual in the (g+1)th generation, is the distribution mean of the (g)th generation, is the step size, is the normal distribution function with mean 0 and standard deviation , is positive definite covariance matrix in (g)th generation, and is the number of the solutions of on (g+1)th generation. The principle of updating is to increase the variance along the successful search direction. The covariance matrix of (g+1)th generation is calculated by the following written expression: and is the evolution path and is calculated bywhere is the learning rate of the covariance matrix, is the learning rate of the evolution path, which is inversely proportional to the number of parameters, is the weight coefficient, and is the mean of individuals in (g)th generation. The CMA-ES needs fewer populations than other evolution strategies. The history path is expressed by the vector sum of mean of each population.

3.2. Application in Weibull Parameter Estimation

For the random samples of the Weibull distribution, , the parameter estimation is done by maximizing (4). It is hard to solve the equations by gradient method due to the difficulties of evaluating the gradient terms. Owing to the advantages of the CMA-ES, it is easy to obtain the maximum of log-likelihood function. Therefore, it is changed to an optimization problem with three parameters and an optimization goal based on the log-likelihood function, defined asHere, is log-likelihood function. are the upper and lower boundaries of , respectively. For an optimization problem, it is important to choose the boundaries. In this article, firstly, we obtain the initial values of the unknown parameters by rank regression method (RRM) [18]; then the left boundaries of the parameters are defined as and the right boundaries are set as five times of the initial values. To estimate the three parameters of the Weibull distribution, the pseudocode of the CMA-ES algorithm is described below.(1)Initialization. Give a vector, .(2)Set the parameters of CMA-ES, such as upper and lower boundary of independent variable.(3)Set number of samples per iteration (4)Initialize state variables, (5)While (true), start iterate(5.1)For in , evaluate new solutions for each (5.2): sample multivariate normal. Here, m: mean, : covariance matrix(5.3): fitness ()(5.4) with sort , sort solutions.(5.5), replace m.(5.6)Move mean to better solutions, update m  (5.7)Update isotropic evolution path. Update ps ()(5.8)Update anisotropic evolution path. Update pc ()(5.9)Update the covariance matrix. Update C (C, , )(5.10)Update step-size using isotropic path length. Update ()(5.11)Return (5.12)Calculate the objective function, (6)Stop when solution is converged, output the parameters .

4. Results and Discussion

To validate the proposed approach, several examples of three-parameter Weibull distribution estimation are considered to reveal the effect of sample size and compare with other intelligent global optimization algorithms. All algorithm evaluations are performed on the processing unit of 2.2GHz CPU and 8 GB RAM with MATLAB 2017B.

4.1. Example 1

The first example focuses on the three-parameter Weibull distribution parameters in vector with different sample size . The parameters are considered here. The sample data are generated by the Monte Carlo acceptance-rejection method [37]. It is a widely used method for the generation of pseudorandom numbers according to their probability density function. The algorithm is briefly introduced as follows:(1)Generate pseudorandom numbers u, v from the uniform distribution on (a, b) and (0,1), respectively.(2)Calculate by the inverse of the cumulative distribution function, X=F-1(u).(3)If v<X/Xsup, accept X; otherwise reject X.(4)Repeat Steps 1 to 3, until the necessary number of samples has been accepted.

Here, (a, b) is the our interval of interest and is the supremum of probability density function in the interval.

The sample sizes of 100, 500, 1000, and 2500 are chosen for each example. The initial value of the CMA-ES is set as the mean value of the search space, namely, the boundaries of the parameters. The step size determines the coordinate-wise standard deviations for the search. Here, we set the initial step size as the standard deviation of the initial value.

The optimization history of three parameters for the case of sample size 100 is shown in Figure 2. The algorithm converges in about 632 iterations. Once the objective function value is converged, the estimates of the three parameters are obtained. Then, to obtain the statistical results, the program has been performed ten times for averaging the results. The results are tabulated in Table 1 for estimation scenario. As can be seen, the estimation results become better as the sample size increases. Meanwhile, the running time increases as the sample size increases, and for larger sample size, it takes a long time to converge.

4.2. Example 2

In this example, the method is applied to experimentally obtained data—failure data of ceramic material. The failure data is shown in Table 2 [34].

The parameters from literature [34] are . In this example, the comparison with intelligent global optimization algorithms was carried out. Firstly, the initial values of the unknown parameters are calculated by RRM, with a searching space of . These algorithms include genetic algorithm (GA), particle swarm optimization (PSO) algorithm, and simulated annealing (SA) algorithm. The parameters of GA are set as maximum generation 100, population size 20, crossover probability 0.4, and mutation probability 0.2. The parameters of PSO are set as cognitive learning parameter 2, social learning parameter 2, inertia weight 0.8, and individuals 20. The cooling schedule parameters of SA are set as Markov chain length 10, decay scale 0.95, step factor 0.02, initial temperature 30, and accept points 0. The Weibull parameters calculated from CMA-ES, GA, PSO, and SA are shown in Table 3. It can be seen from Table 3 that the result of CMA-ES is very close to the literature result. The result of CMA-ES is better than the other methods. To further test the proposed method, the calculation is carried out for different search spaces [38]. The results are listed in Table 4, which is average of ten calculations. As can be seen, the results of small search space are better than large search space and the large search space can lead to unreliable results.

4.3. Example 3

In this example, a glass strength data set is adopted to test the proposed method. The data are from literature [3], where only the data of the naturally aged glass is utilized, listed in Table 5.

There are 18 raw data and we use the proposed method to estimate the three unknown parameters. The history curve of objective function is plotted in Figure 3. As can be seen, the optimal values converged after 294 iterations. The simulation results are listed in Table 6.

It should be noted that the two-parameter Weibull distribution is used to fit the raw strength data. The shape and scale parameters by the proposed method are close to the literature. It can be seen that the proposed method is valid for this case. Furthermore, to compare the goodness-of-fit for glass strength data, the cumulative distribution functions (CDF) and empirical distribution function (EDF) curves are plotted in Figure 4. The mean square values of empirical distribution function and cumulative distribution functions of the sample data are 0.6035, 0.0191, 0.2317, and 0.0025 for LR, WLR-F&T, WLR-B, and proposed method, respectively. The frequency of raw experimental data and the estimated probability density function are shown in Figure 5. As can be seen, the probability density function curve of the proposed method is closest to the frequency distribution curves.

To further test the effectiveness of the proposed method, Kolmogorov-Smirnov (K-S) tests and the Akaike Information Criterion (AIC) are adopted. In this paper, AIC is calculated aswhere p is the number of parameters, L is the sum of log-likelihood, and n is the sample size. The smaller value of AIC means the model is better. The significance level is set as 0.05 for K-S test. As can be seen, the CDF data by our proposed method has the largest probability, same as the EDF. The AIC value of the proposed method is smallest, which is 7.6697. It indicates that the result of the proposed method outperforms the others.

5. Conclusions

In this article, the estimation of Weibull distribution parameters is converted to an optimization problem and solved by the covariance matrix adaptation evolution strategy (CMA-ES). The fitness function of the CMA-ES is the likelihood function of the three-parameter Weibull model. The optimization variables are the unknown parameters, namely, location, scale, and shape parameters. The results show that better estimation results are obtained with larger samples. Compared with other heuristic approaches, the CMA-ES has a good performance for parameter estimation. The proposed method is also suitable for real data sets and performs better than other methods in terms of goodness-of-fit for empirical distribution functions.

For the further research, the proposed method is aimed to be applied to parameter estimation of other distributions, such as log-normal distribution, skew normal distribution, and some multiparameter distributions.

Data Availability

The digital data supporting this article are from previously reported studies and datasets, which have been cited.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work is supported by the Fundamental Research Funds for the Central Universities, NO. NS2018004.