Mathematical Problems in Engineering
Volume 2019 (2019), Article ID 6281781, 8 pages
https://doi.org/10.1155/2019/6281781

Research Article

Maximum Likelihood Estimation for Three-Parameter Weibull Distribution Using Evolutionary Strategy

Fan Yang,1,2 Hu Ren,³ and Zhili Hu1,2

¹State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China
²College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China
³Wuxi Hengding Supercomputing Center Ltd., China

Correspondence should be addressed to Fan Yang; fanyang@nuaa.edu.cn and Zhili Hu; zhili.hu@nuaa.edu.cn

Received 12 November 2018; Revised 18 March 2019; Accepted 2 May 2019; Published 23 May 2019

Academic Editor: Carmen Castillo

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.

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 [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 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 [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 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.

Figure 1: Weibull PDF with various values of when assuming and .

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 as

The 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 as

Then, 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, as

We 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 as

where 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 by

where 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 as

Here, 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⁻¹(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.

Table 1: Estimation results of different methods for .

Sample size

Time/s

Exact value
0
2
1.5

100
0.0123
1.9506
1.4804
8.12

500
0.0030
1.9677
1.4901
45.36

1000
3.31e-4
1.9718
1.4928
93.19

2500
6.306e-5
2.0184
1.4981
248.48

Figure 2: History of the optimization parameter.

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].

Table 2: Failure data of ceramic material/MPa.

Raw failure data (MPa)

307
308
322
328
328
329
331
332
335

337
343
345
347
350
352
353
355
356

357
364
371
373
374
375
376
376
381

385
388
395
402
411
413
415
456

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.

Table 3: Parameters obtained by optimization method.

Methods

Ref [34]
298
2.0
73

CMA-ES
298.29
1.9973
72.92

GA
305.41
1.78
64.58

SA
300.29
1.72
63.83

PSO
298.63
2.02
71.66

Table 4: Results for different search space.

Search space

298.29
1.9973
72.92

302.76
2.38
77.35

350.21
10.17
144.68

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.

Table 5: Raw strength data.

Failure strength, (MPa)

24.12
24.13
28.52
29.18
29.67
30.48

32.98
35.91
35.92
36.38
37.60
37.70

39.71
49.10
52.43
52.46
52.61
61.72

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.

Table 6: Results of different estimation methods [3].

Method

K-S(ks_stat/p)
AIC

LR (Ref [3])
∖
4.2
42.2
(0.2222/0.7088)
7.9410

WLR-F&T (Ref [3])
∖
3.4
41.7
(0.1667/0.9448)
7.9438

WLR-B (Ref [3])
∖
3.2
41.3
(0.1667/0.9448)
7.9727

CMA-ES(ML)
22.4101
1.3954
17.4991
(0.1111/0.9997)
7.6697

Note. LR, least square regression; WLR-F&T, weighted least squares regression with Faucher & Tyson’s equation [35]; WLR-B, weighted least squares regression with Bergman’s equation [36]; CMA-ES(ML), the proposed method ML using CMA-ES; ks_stat/p, value of K-S test/the probability of the CDF data similarity.

Figure 3: Iteration history of objective function.

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.

Figure 4: Empirical distribution function and assumed cumulative distribution function.

Figure 5: Empirical density function and assumed density distribution function.

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 as

where 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.

References

W. Weibull, “Wide applicability,” Journal of Applied Mechanics, 1951.

E. S. Lindquist, “Strength of materials and the Weibull distribution,” Probabilistic Engineering Mechanics, vol. 9, no. 3, pp. 191–194, 1994.

K. C. Datsiou and M. Overend, “Weibull parameter estimation and goodness-of-fit for glass strength data,” Structural Safety, vol. 73, pp. 29–41, 2018.

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.

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.

P. K. Chaurasiya, S. Ahmed, and V. Warudkar, “Comparative analysis of Weibull parameters for wind data measured from met-mast and remote sensing techniques,” Journal of Renewable Energy, vol. 115, pp. 1153–1165, 2018.

P. Wais, “Two and three-parameter Weibull distribution in available wind power analysis,” Journal of Renewable Energy, vol. 103, pp. 15–29, 2017.

J. H. K. Kao, “A graphical estimation of mixed weibull parameters in life-testing of electron tubes,” Technometrics, vol. 1, no. 4, pp. 389–407, 1959.

R. Jiang and D. N. P. Murthy, “Modeling failure-data by mixture of 2 weibull distributions: a graphical approach,” IEEE Transactions on Reliability, vol. 44, no. 3, pp. 477–488, 1995.

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.

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.

J. R. M. Hosking, J. R. Wallis, and E. F. Wood, “Estimation of the generalized extreme-value distribution by the method of probability-weighted moments,” Technometrics, vol. 27, no. 3, pp. 251–261, 1985.

A. C. Cohen, B. J. Whitten, and Y. Ding, “Modified moment estimation for the three-parameter weibull distribution,” Journal of Quality Technology, vol. 16, no. 3, pp. 159–167, 1984.

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.

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.

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.

D. Marković, D. Jukić, and M. Benšić, “Nonlinear weighted least squares estimation of a three-parameter Weibull density with a nonparametric start,” Journal of Computational and Applied Mathematics, vol. 228, no. 1, pp. 304–312, 2009.

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.

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.

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.

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.

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.

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.

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.

H. H. Örkcü, E. Aksoy, and M. I. Dogan, “Estimating the parameters of 3-p Weibull distribution through differential evolution,” Applied Mathematics and Computation, vol. 251, pp. 211–224, 2015.

I. Usta, “An innovative estimation method regarding Weibull parameters for wind energy applications,” Energy, vol. 106, pp. 301–314, 2016.

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.

D. Petković, “Adaptive neuro-fuzzy approach for estimation of wind speed distribution,” International Journal of Electrical Power & Energy Systems, vol. 73, pp. 389–392, 2015.

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.

N. Hansen and A. Ostermeier, “Completely derandomized self-adaptation in evolution strategies,” Evolutionary Computation, vol. 9, no. 2, pp. 159–195, 2001.

N. Hansen, S. D. Müller, and P. Koumoutsakos, “Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES),” Evolutionary Computation, vol. 11, no. 1, pp. 1–18, 2003.

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.

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.

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.

B. Bergman, “Estimation of Weibull parameters using a weight function,” Journal of Materials Science Letters, vol. 5, no. 6, pp. 611–614, 1986.

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.