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

Authors: Fan Yang; Hu Ren; Zhili Hu
Journal: Mathematical Problems in Engineering
 (2019)

---

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


---

## Body

## 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. Section2 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(1)Ft,α,β,η=1-e-t-α/ηβ,t>α0,t≤α(2)ft,α,β,η=βηt-αηβ-1e-t-α/ηβ,t>α0,t≤αHere,α≥0, β≥0, and η≥0 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 α=0 and η=1.5.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.Lett1,t2,⋯,tn be a random sample of size n; θ⃑ is noted as the Weibull model parameters which are to be estimated; namely, θ⃑=(α,β,η). The likelihood function is written as(3)L=∏i=1nftiti;θ⃑=∏i=1nβηti-αηβ-1e-ti-α/ηβ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(4)ln⁡Lti;θ⃑=∑i=1nln⁡β+β-1ln⁡ti-α-βln⁡η-ti-αηβ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(5)∂ln⁡Lti;θ∂α=0∂ln⁡Lti;θ∂β=0∂ln⁡Lti;θ∂η=0We substitute the log-likelihood function into the above equations. The following equations are obtained:(6)L1=∑i=1n1β+ln⁡ti-α-ln⁡η-ti-αηβln⁡ti-αη=0L2=∑i=1n-βη+βηti-αηβ=0L3=∑i=1n-β-1ti-α+βηti-αηβ-1=0It 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 C. The basic function of CMA-ES algorithm is shown as(7)xkg+1=mg+σgN0,Cg,k=1,2,…,λwhere xk(g+1)∈Rn is the individual in the (g+1)th generation, m(g) is the distribution mean of the (g)th generation, σ(g) is the step size, N(0,C(g)) is the normal distribution function with mean 0 and standard deviation C(g), C(g) is positive definite covariance matrix in (g)th generation, and λ is the number of the solutions of x on (g+1)th generation. The principle of updating C is to increase the variance along the successful search direction. The covariance matrix C(g+1) of (g+1)th generation is calculated by the following written expression: (8)Cg+1=1-ccov·Cg+ccov·pcg+1pcg+1Tand pc(g+1) is the evolution path and is calculated by(9)pcg+1=1-cc·pcg+cc·2-cc·μσgxμg+1-xμgwhere ccov is the learning rate of the covariance matrix, cc is the learning rate of the evolution path, which is inversely proportional to the number of parameters, μ is the weight coefficient, and xμ(g) 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,t1,t2,⋯,tn, 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(10)minf=∑i=13Li2θ⃑,θ⃑=α,β,ηs.tαl≤α≤αuβl≤β≤βuηl≤η≤ηuHere, Li(θ⃑) is log-likelihood function. αl,αu,βl,βu,ηl,ηu 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 [0,0,0] 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,x⃑=θ⃑(α0,β0,η0).(2)
Set the parameters of CMA-ES, such as upper and lower boundary of independent variable.(3)
Set number of samples per iterationλ=8(4)
Initialize state variables,m,σ,C=I,pσ=0,pc=0(5)
While (true), start iterate
(5.1)
Fori in {1…λ}, evaluate new solutions for each λ(5.2)
xi: sample multivariate normal. Here,m: mean, σ2C: covariance matrix(5.3)
fi: fitness (xi)(5.4)
x1…λ←xs(1)…s(λ) with s(i) sort (f1…λ,i), sort solutions.(5.5)
m′=m, replacem.(5.6)
Move mean to better solutions, updatem(x1,…,xλ)(5.7)
Update isotropic evolution path. Updateps (pσ,σ-1C-1/2(m-m′))(5.8)
Update anisotropic evolution path. Updatepc (pc,σ-1(m-m′),pσ)(5.9)
Update the covariance matrix. UpdateC (C, pc, (x1-m′)/σ,…,(xλ-m′)/σ)(5.10)
Update step-size using isotropic path length. Updateσ (σ,pσ)(5.11)
Returnx(5.12)
Calculate the objective function,∑i=13Li2(θ⃑)(6)
Stop when solution is converged, output the parametersx⃑=θ⃑(α,β,η).


## 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 C. The basic function of CMA-ES algorithm is shown as(7)xkg+1=mg+σgN0,Cg,k=1,2,…,λwhere xk(g+1)∈Rn is the individual in the (g+1)th generation, m(g) is the distribution mean of the (g)th generation, σ(g) is the step size, N(0,C(g)) is the normal distribution function with mean 0 and standard deviation C(g), C(g) is positive definite covariance matrix in (g)th generation, and λ is the number of the solutions of x on (g+1)th generation. The principle of updating C is to increase the variance along the successful search direction. The covariance matrix C(g+1) of (g+1)th generation is calculated by the following written expression: (8)Cg+1=1-ccov·Cg+ccov·pcg+1pcg+1Tand pc(g+1) is the evolution path and is calculated by(9)pcg+1=1-cc·pcg+cc·2-cc·μσgxμg+1-xμgwhere ccov is the learning rate of the covariance matrix, cc is the learning rate of the evolution path, which is inversely proportional to the number of parameters, μ is the weight coefficient, and xμ(g) 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,t1,t2,⋯,tn, 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(10)minf=∑i=13Li2θ⃑,θ⃑=α,β,ηs.tαl≤α≤αuβl≤β≤βuηl≤η≤ηuHere, Li(θ⃑) is log-likelihood function. αl,αu,βl,βu,ηl,ηu 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 [0,0,0] 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,x⃑=θ⃑(α0,β0,η0).(2)
Set the parameters of CMA-ES, such as upper and lower boundary of independent variable.(3)
Set number of samples per iterationλ=8(4)
Initialize state variables,m,σ,C=I,pσ=0,pc=0(5)
While (true), start iterate
(5.1)
Fori in {1…λ}, evaluate new solutions for each λ(5.2)
xi: sample multivariate normal. Here,m: mean, σ2C: covariance matrix(5.3)
fi: fitness (xi)(5.4)
x1…λ←xs(1)…s(λ) with s(i) sort (f1…λ,i), sort solutions.(5.5)
m′=m, replacem.(5.6)
Move mean to better solutions, updatem(x1,…,xλ)(5.7)
Update isotropic evolution path. Updateps (pσ,σ-1C-1/2(m-m′))(5.8)
Update anisotropic evolution path. Updatepc (pc,σ-1(m-m′),pσ)(5.9)
Update the covariance matrix. UpdateC (C, pc, (x1-m′)/σ,…,(xλ-m′)/σ)(5.10)
Update step-size using isotropic path length. Updateσ (σ,pσ)(5.11)
Returnx(5.12)
Calculate the objective function,∑i=13Li2(θ⃑)(6)
Stop when solution is converged, output the parametersx⃑=θ⃑(α,β,η).


## 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 n. The parameters θ⃑=(0,2,1.5) 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 numbersu,v from the uniform distribution on (a,b) and (0,1), respectively.(2)
CalculateX by the inverse of the cumulative distribution function,X=F-1(u).(3)
Ifv<X/Xsup, acceptX; otherwise rejectX.(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 Xsup 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 Figure2. 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 θ⃑=(0,2,1.5) 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θ⃑=(0,2,1.5).
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.48Figure 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 Table2 [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 456The parameters from literature [34] are θ⃑=(298,2,73). In this example, the comparison with intelligent global optimization algorithms was carried out. Firstly, the initial values θ⃑0 of the unknown parameters are calculated by RRM, with a searching space of [0,5θ⃑0]. 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.66Table 4
Results for different search space.
Search space α β η [ 0 , 5 θ ⃑ 0 ] 298.29 1.9973 72.92 [ 0 , 10 θ ⃑ 0 ] 302.76 2.38 77.35 [ 0 , 50 θ ⃑ 0 ] 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,σf(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.72There 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 Figure3. 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 Figure4. 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(11)AIC=2p-2Lnwherep is the number of parameters,L is the sum of log-likelihood, andn 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.


## 4.1. Example 1

The first example focuses on the three-parameter Weibull distribution parameters in vectorθ⃑=(α,β,η) with different sample size n. The parameters θ⃑=(0,2,1.5) 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 numbersu,v from the uniform distribution on (a,b) and (0,1), respectively.(2)
CalculateX by the inverse of the cumulative distribution function,X=F-1(u).(3)
Ifv<X/Xsup, acceptX; otherwise rejectX.(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 Xsup 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 Figure2. 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 θ⃑=(0,2,1.5) 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θ⃑=(0,2,1.5).
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.48Figure 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 Table2 [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 456The parameters from literature [34] are θ⃑=(298,2,73). In this example, the comparison with intelligent global optimization algorithms was carried out. Firstly, the initial values θ⃑0 of the unknown parameters are calculated by RRM, with a searching space of [0,5θ⃑0]. 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.66Table 4
Results for different search space.
Search space α β η [ 0 , 5 θ ⃑ 0 ] 298.29 1.9973 72.92 [ 0 , 10 θ ⃑ 0 ] 302.76 2.38 77.35 [ 0 , 50 θ ⃑ 0 ] 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,σf(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.72There 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 Figure3. 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 Figure4. 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(11)AIC=2p-2Lnwherep is the number of parameters,L is the sum of log-likelihood, andn 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.



---

*Source: 6281781-2019-05-23.xml*