Research Article  Open Access
Fan Yang, Hu Ren, Zhili Hu, "Maximum Likelihood Estimation for ThreeParameter Weibull Distribution Using Evolutionary Strategy", Mathematical Problems in Engineering, vol. 2019, Article ID 6281781, 8 pages, 2019. https://doi.org/10.1155/2019/6281781
Maximum Likelihood Estimation for ThreeParameter Weibull Distribution Using Evolutionary Strategy
Abstract
The maximum likelihood estimation is a widely used approach to the parameter estimation. However, the conventional algorithm makes the estimation procedure of threeparameter 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 threeparameter 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 threeparameter 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 threeparameter Weibull distribution parameters estimation is an extremely complicated procedure. To date, several estimation approaches have been conducted, such as graphic method [8–10], moment estimation [11–14], maximum likelihood estimation (ML) [15–17], 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 socalled 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 threeparameter 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 neurofuzzy 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 (CMAES) algorithm [31–33] 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 CMAES 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 ThreeParameter Weibull Distribution
The cumulative distribution function (CDF) and probability density function (PDF) of the threeparameter 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 probabilityweighted 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 threeparameter 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 loglikelihood 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 (CMAES) is used because it is robust and effective for global optimization.
3. Evolution Optimization
3.1. CMAES Algorithm
In this section, the Covariance Matrix Adaptation Evolution Strategy is briefly introduced for threeparameter Weibull distribution. The CMAES 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 CMAES is a secondorder 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 CMAES 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 CMAES 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 CMAES 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 CMAES, it is easy to obtain the maximum of loglikelihood function. Therefore, it is changed to an optimization problem with three parameters and an optimization goal based on the loglikelihood function, defined asHere, is loglikelihood 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 CMAES algorithm is described below.(1)Initialization. Give a vector, .(2)Set the parameters of CMAES, 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 stepsize 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 threeparameter 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 threeparameter Weibull distribution parameters in vector with different sample size . The parameters are considered here. The sample data are generated by the Monte Carlo acceptancerejection 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/X_{sup}, 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 CMAES is set as the mean value of the search space, namely, the boundaries of the parameters. The step size determines the coordinatewise 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 CMAES, GA, PSO, and SA are shown in Table 3. It can be seen from Table 3 that the result of CMAES is very close to the literature result. The result of CMAES 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.
 
Note. LR, least square regression; WLRF&T, weighted least squares regression with Faucher & Tyson’s equation [35]; WLRB, weighted least squares regression with Bergman’s equation [36]; CMAES(ML), the proposed method ML using CMAES; ks_stat/p, value of KS test/the probability of the CDF data similarity. 
It should be noted that the twoparameter 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 goodnessoffit 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, WLRF&T, WLRB, 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, KolmogorovSmirnov (KS) 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 loglikelihood, 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 KS 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 (CMAES). The fitness function of the CMAES is the likelihood function of the threeparameter 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 CMAES 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 goodnessoffit for empirical distribution functions.
For the further research, the proposed method is aimed to be applied to parameter estimation of other distributions, such as lognormal 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.
References
 W. Weibull, “Wide applicability,” Journal of Applied Mechanics, 1951. View at: Google Scholar
 E. S. Lindquist, “Strength of materials and the Weibull distribution,” Probabilistic Engineering Mechanics, vol. 9, no. 3, pp. 191–194, 1994. View at: Publisher Site  Google Scholar
 K. C. Datsiou and M. Overend, “Weibull parameter estimation and goodnessoffit for glass strength data,” Structural Safety, vol. 73, pp. 29–41, 2018. View at: Publisher Site  Google Scholar
 R. Danzer, P. Supancic, J. Pascual, and T. Lube, “Fracture statistics of ceramics—weibull statistics and deviations from Weibull statistics,” Engineering Fracture Mechanics, vol. 74, no. 18, pp. 2919–2932, 2007. View at: Publisher Site  Google Scholar
 L. Afferrante, M. Ciavarella, and E. Valenza, “Is Weibull's modulus really a material constant? Example case with interacting collinear cracks,” International Journal of Solids and Structures, vol. 43, no. 17, pp. 5147–5157, 2006. View at: Publisher Site  Google Scholar
 P. K. Chaurasiya, S. Ahmed, and V. Warudkar, “Comparative analysis of Weibull parameters for wind data measured from metmast and remote sensing techniques,” Journal of Renewable Energy, vol. 115, pp. 1153–1165, 2018. View at: Publisher Site  Google Scholar
 P. Wais, “Two and threeparameter Weibull distribution in available wind power analysis,” Journal of Renewable Energy, vol. 103, pp. 15–29, 2017. View at: Publisher Site  Google Scholar
 J. H. K. Kao, “A graphical estimation of mixed weibull parameters in lifetesting of electron tubes,” Technometrics, vol. 1, no. 4, pp. 389–407, 1959. View at: Publisher Site  Google Scholar
 R. Jiang and D. N. P. Murthy, “Modeling failuredata by mixture of 2 weibull distributions: a graphical approach,” IEEE Transactions on Reliability, vol. 44, no. 3, pp. 477–488, 1995. View at: Publisher Site  Google Scholar
 R. Jiang and D. N. P. Murthy, “The exponentiated weibull family: a graphical approach,” IEEE Transactions on Reliability, vol. 48, no. 1, pp. 68–72, 1999. View at: Publisher Site  Google Scholar
 J. A. Greenwood, J. M. Landwehr, N. C. Matalas, and J. R. Wallis, “Probability weighted moments: definition and relation to parameters of several distributions expressable in inverse form,” Water Resources Research, vol. 15, no. 5, pp. 1049–1054, 1979. View at: Publisher Site  Google Scholar
 J. R. M. Hosking, J. R. Wallis, and E. F. Wood, “Estimation of the generalized extremevalue distribution by the method of probabilityweighted moments,” Technometrics, vol. 27, no. 3, pp. 251–261, 1985. View at: Publisher Site  Google Scholar
 A. C. Cohen, B. J. Whitten, and Y. Ding, “Modified moment estimation for the threeparameter weibull distribution,” Journal of Quality Technology, vol. 16, no. 3, pp. 159–167, 1984. View at: Publisher Site  Google Scholar
 S. A. Akdağ and Ö. Güler, “Alternative Moment Method for wind energy potential and turbine energy output estimation,” Journal of Renewable Energy, vol. 120, pp. 69–77, 2018. View at: Publisher Site  Google Scholar
 N. Balakrishnan and M. Kateri, “On the maximum likelihood estimation of parameters of Weibull distribution based on complete and censored data,” Statistics & Probability Letters, vol. 78, no. 17, pp. 2971–2975, 2008. View at: Publisher Site  Google Scholar
 P. Groeneboom and J. A. Wellner, Information Bounds and Nonparametric Maximum Likelihood Estimation, vol. 19, Birkhäuser, 2012.
 W. Chen, M. Xie, and M. Wu, “Modified maximum likelihood estimator of scale parameter using moving extremes ranked set sampling,” Communications in Statistics—Simulation and Computation, vol. 45, no. 6, pp. 2232–2240, 2016. View at: Publisher Site  Google Scholar  MathSciNet
 D. Marković, D. Jukić, and M. Benšić, “Nonlinear weighted least squares estimation of a threeparameter Weibull density with a nonparametric start,” Journal of Computational and Applied Mathematics, vol. 228, no. 1, pp. 304–312, 2009. View at: Publisher Site  Google Scholar  MathSciNet
 Y.J. Kang, Y. Noh, and O.K. Lim, “Kernel density estimation with bounded data,” Structural and Multidisciplinary Optimization, vol. 57, no. 1, pp. 95–113, 2018. View at: Publisher Site  Google Scholar  MathSciNet
 T. C. Carneiro, S. P. Melo, P. C. M. Carvalho, and A. P. D. S. Braga, “Particle swarm optimization method for estimation of weibull parameters: a case study for the Brazilian northeast region,” Journal of Renewable Energy, vol. 86, pp. 751–759, 2016. View at: Publisher Site  Google Scholar
 A. Yalçınkaya, B. Şenoğlu, and U. Yolcu, “Maximum likelihood estimation for the parameters of skew normal distribution using genetic algorithm,” Swarm and Evolutionary Computation, vol. 38, pp. 127–138, 2018. View at: Publisher Site  Google Scholar
 C. Leahy, E. J. O'Brien, B. Enright, and R. T. Power, “Estimating characteristic bridge traffic load effects using Bayesian statistics,” in Proceedings of the 12th International Conference on Applications of Statistics and Probability in Civil Engineering (ICASP12), Vancouver, Canada, July, 2015. View at: Google Scholar
 M. Keller, A.L. Popelin, N. Bousquet, and E. Remy, “Nonparametric estimation of the probability of detection of flaws in an industrial component, from destructive and nondestructive testing data, using approximate bayesian computation,” Risk Analysis, vol. 35, no. 9, pp. 1595–1610, 2015. View at: Publisher Site  Google Scholar
 R. Aghmasheh, V. Rashtchi, and E. Rahimpour, “Gray box modeling of power transformer windings for transient studies,” IEEE Transactions on Power Delivery, vol. 32, no. 5, pp. 2350–2359, 2017. View at: Publisher Site  Google Scholar
 T. P. Talafuse and E. A. Pohl, “Small sample reliability growth modeling using a grey systems model,” Quality Engineering, vol. 29, no. 3, pp. 455–467, 2017. View at: Publisher Site  Google Scholar
 H. H. Örkcü, E. Aksoy, and M. I. Dogan, “Estimating the parameters of 3p Weibull distribution through differential evolution,” Applied Mathematics and Computation, vol. 251, pp. 211–224, 2015. View at: Publisher Site  Google Scholar  MathSciNet
 I. Usta, “An innovative estimation method regarding Weibull parameters for wind energy applications,” Energy, vol. 106, pp. 301–314, 2016. View at: Publisher Site  Google Scholar
 H. H. Örkcü, V. S. Özsoy, E. Aksoy, and M. I. Dogan, “Estimating the parameters of  Weibull distribution using particle swarm optimization: a comprehensive experimental comparison,” Applied Mathematics and Computation, vol. 268, pp. 201–226, 2015. View at: Publisher Site  Google Scholar  MathSciNet
 D. Petković, “Adaptive neurofuzzy approach for estimation of wind speed distribution,” International Journal of Electrical Power & Energy Systems, vol. 73, pp. 389–392, 2015. View at: Publisher Site  Google Scholar
 M. Nassar, A. Z. Afify, S. Dey, and D. Kumar, “A new extension of Weibull distribution: properties and different methods of estimation,” Journal of Computational and Applied Mathematics, vol. 336, pp. 439–457, 2018. View at: Publisher Site  Google Scholar  MathSciNet
 N. Hansen and A. Ostermeier, “Completely derandomized selfadaptation in evolution strategies,” Evolutionary Computation, vol. 9, no. 2, pp. 159–195, 2001. View at: Publisher Site  Google Scholar
 N. Hansen, S. D. Müller, and P. Koumoutsakos, “Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMAES),” Evolutionary Computation, vol. 11, no. 1, pp. 1–18, 2003. View at: Publisher Site  Google Scholar
 N. Hansen and S. Kern, “Evaluating the CMA evolution strategy on multimodal test functions,” in Parallel Problem Solving from Nature—PPSN VIII, vol. 3242 of Lecture Notes in Computer Science, pp. 282–291, Springer, Berlin, Germany, 2004. View at: Publisher Site  Google Scholar
 S. F. Duffy, J. L. Palko, and J. P. Gyekenyesi, “Structural reliability analysis of laminated CMC components,” Journal of Engineering for Gas Turbines and Power, vol. 115, no. 1, pp. 103–108, 1993. View at: Publisher Site  Google Scholar
 B. Faucher and W. R. Tyson, “On the determination of Weibull parameters,” Journal of Materials Science Letters, vol. 7, no. 11, pp. 1199–1203, 1988. View at: Publisher Site  Google Scholar
 B. Bergman, “Estimation of Weibull parameters using a weight function,” Journal of Materials Science Letters, vol. 5, no. 6, pp. 611–614, 1986. View at: Publisher Site  Google Scholar
 J. E. Gentle, Random Number Generation and Monte Carlo Methods, Springer Science & Business Media, 2006.
 S. Acitas, C. H. Aladag, and B. Senoglu, “A new approach for estimating the parameters of Weibull distribution via particle swarm optimization: an application to the strengths of glass fibre data,” Reliability Engineering & System Safety, vol. 183, pp. 116–127, 2019. View at: Google Scholar
Copyright
Copyright © 2019 Fan Yang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.