Maximum Likelihood Estimation for Three-Parameter 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 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.

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.

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$$\begin{array}{}\text{(1)}& F\left(t,\alpha ,\beta ,\eta \right)=\left\{\begin{array}{cc}\mathrm{1}-e^{-{\left(\left(t-\alpha \right)/\eta \right)}^{\beta }},\hfill & t>\alpha \hfill \\ \mathrm{0},\hfill & t\le \alpha \hfill \end{array}\right.\text{(2)}& f\left(t,\alpha ,\beta ,\eta \right)=\left\{\begin{array}{cc}\frac{\beta }{\eta }{\left(\frac{t-\alpha }{\eta }\right)}^{\beta -\mathrm{1}}{e}^{-{\left(\left(t-\alpha \right)/\eta \right)}^{\beta }},\hfill & t>\alpha \hfill \\ \mathrm{0},\hfill & t\le \alpha \hfill \end{array}\right.\end{array}$$

Here, $\alpha \ge \mathrm{0}$, $\beta \ge \mathrm{0}$, and $\eta \ge \mathrm{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 omitted; refer to PDF]

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 $t_{\mathrm{1}},t_{\mathrm{2}},\cdots ,t_n$ be a random sample of size $n$; $\stackrel{⃑}{\theta }$ is noted as the Weibull model parameters which are to be estimated; namely, $\stackrel{⃑}{\theta }=(\alpha ,\beta ,\eta )$. The likelihood function is written as$$\begin{array}{}\text{(3)}& L={\displaystyle \prod }_{i=\mathrm{1}}^n{f}_{t_i}\left(t_i;\stackrel{⃑}{\theta }\right)={\displaystyle \prod }_{i=\mathrm{1}}^n\frac{\beta }{\eta }{\left(\frac{t_i-\alpha }{\eta }\right)}^{\beta -\mathrm{1}}{e}^{-{\left(\left(t_i-\alpha \right)/\eta \right)}^{\beta }}\end{array}$$The aim of estimation is to determine the unknown vector $\stackrel{⃑}{\theta }$ and the three unknown parameters $\alpha$, $\beta$, and $\eta$ 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$$\begin{array}{}\text{(4)}& \ln \left[L\left(t_i;\stackrel{⃑}{\theta }\right)\right]={\displaystyle \sum }_{i=\mathrm{1}}^n\left[\ln \left(\beta \right)+\left(\beta -\mathrm{1}\right)\ln \left(t_i-\alpha \right)-\beta \ln \left(\eta \right)-{\left(\frac{\left(t_i-\alpha \right)}{\eta }\right)}^{\beta }\right]\end{array}$$Then, the vector $\stackrel{⃑}{\theta }$ 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 $\stackrel{⃑}{\theta }$ and set the partial equations to zero, as$$\begin{array}{c}\text{(5)}& \frac{\partial \ln \left[L\left(t_i;\theta \right)\right]}{\partial \alpha }=\mathrm{0}\frac{\partial \ln \left[L\left(t_i;\theta \right)\right]}{\partial \beta }=\mathrm{0}\\ \frac{\partial \ln \left[L\left(t_i;\theta \right)\right]}{\partial \eta }=\mathrm{0}\end{array}$$We substitute the log-likelihood function into the above equations. The following equations are obtained:$$\begin{array}{c}\text{(6)}& L_{\mathrm{1}}={\displaystyle \sum }_{i=\mathrm{1}}^n\left[\frac{\mathrm{1}}{\beta }+\ln \left(t_i-\alpha \right)-\ln \left(\eta \right)-{\left(\frac{\left(t_i-\alpha \right)}{\eta }\right)}^{\beta }\ln \left(\frac{\left(t_i-\alpha \right)}{\eta }\right)\right]=\mathrm{0}{L}_{\mathrm{2}}={\displaystyle \sum }_{i=\mathrm{1}}^n\left[-\frac{\beta }{\eta }+\left(\frac{\beta }{\eta }\right){\left(\frac{\left(t_i-\alpha \right)}{\eta }\right)}^{\beta }\right]=\mathrm{0}\\ L_{\mathrm{3}}={\displaystyle \sum }_{i=\mathrm{1}}^n\left[-\frac{\left(\beta -\mathrm{1}\right)}{\left(t_i-\alpha \right)}+\left(\frac{\beta }{\eta }\right){\left(\frac{\left(t_i-\alpha \right)}{\eta }\right)}^{\beta -\mathrm{1}}\right]=\mathrm{0}\end{array}$$

It is well known that it is difficult to obtain the estimates of the unknown parameters $(\widehat{\alpha },\widehat{\beta },\widehat{\eta })$ 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 $\sigma$ and covariance matrix $C$. The basic function of CMA-ES algorithm is shown as$$\begin{array}{}\text{(7)}& x_k^{\left(g+\mathrm{1}\right)}=m^{\left(g\right)}+{\sigma }^{\left(g\right)}N\left(\mathrm{0},C^{\left(g\right)}\right),\hspace{1em}k=\mathrm{1,2},\dots ,\lambda \end{array}$$where $x_k^{(g+\mathrm{1})}\in R^n$ is the individual in the (g+1)th generation, $m^{(g)}$ is the distribution mean of the (g)th generation, ${\sigma }^{(g)}$ is the step size, $N(\mathrm{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 $\lambda$ 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: $$\begin{array}{}\text{(8)}& C^{\left(g+\mathrm{1}\right)}=\left(\mathrm{1}-c_{\mathrm{c}\mathrm{o}\mathrm{v}}\right)·C^{\left(g\right)}+c_{\mathrm{c}\mathrm{o}\mathrm{v}}·p_c^{\left(g+\mathrm{1}\right)}{\left(p_c^{\left(g+\mathrm{1}\right)}\right)}^T\end{array}$$and $p_c^{(g+\mathrm{1})}$ is the evolution path and is calculated by$$\begin{array}{}\text{(9)}& p_c^{\left(g+\mathrm{1}\right)}=\left(\mathrm{1}-c_c\right)·p_c^{\left(g\right)}+\sqrt{c_c·\left(\mathrm{2}-c_c\right)}·\frac{\sqrt{\mu }}{{\sigma }^{\left(g\right)}}\left({⟨x⟩}_{\mu }^{\left(g+\mathrm{1}\right)}-{⟨x⟩}_{\mu }^{\left(g\right)}\right)\end{array}$$where $c_{\mathrm{c}\mathrm{o}\mathrm{v}}$ is the learning rate of the covariance matrix, $c_c$ is the learning rate of the evolution path, which is inversely proportional to the number of parameters, $\mu$ is the weight coefficient, and ${⟨x⟩}_{\mu }^{(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, $t_{\mathrm{1}},t_{\mathrm{2}},\cdots ,t_n$, 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$$\begin{array}{c}\text{(10)}& \mathrm{m}\mathrm{i}\mathrm{n}\hspace{1em}f={\displaystyle \sum }_{i=\mathrm{1}}^{\mathrm{3}}{L}_i^{\mathrm{2}}\left(\stackrel{⃑}{\theta }\right),\hspace{1em}\stackrel{⃑}{\theta }=\left(\alpha ,\beta ,\eta \right)\mathrm{s}.\mathrm{t}\hspace{1em}{\alpha }^l\le \alpha \le {\alpha }^u\\ \hspace{1em}{\beta }^l\le \beta \le {\beta }^u\\ \hspace{1em}{\eta }^l\le \eta \le {\eta }^u\end{array}$$Here, $L_i(\stackrel{⃑}{\theta })$ is log-likelihood function. ${\alpha }^l,{\alpha }^u,{\beta }^l,{\beta }^u,{\eta }^l,{\eta }^u$ are the upper and lower boundaries of $\alpha ,\beta ,\eta$, 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 $[\mathrm{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, $\stackrel{⃑}{x}=\stackrel{⃑}{\theta }({\alpha }_{\mathrm{0}},{\beta }_{\mathrm{0}},{\eta }_{\mathrm{0}})$.

(2)

Set the parameters of CMA-ES, such as upper and lower boundary of independent variable.

(3)

Set number of samples per iteration $\lambda =8$

(4)

Initialize state variables, $m,\sigma ,C=I,p_{\sigma }=\mathrm{0},p_c=\mathrm{0}$

(5)

While (true), start iterate

(5.1)

For $i$ in $\{1\dots \lambda \}$, evaluate new solutions for each $\lambda$

(5.2)

$x_i$: sample multivariate normal. Here, m: mean, ${\sigma }^{\mathrm{2}}C$: covariance matrix

(5.3)

$f_i$: fitness ($x_i$)

(5.4)

$x_{\mathrm{1}\dots \lambda }\leftarrow x_{s(\mathrm{1})\dots s(\lambda )}$ with $s(i)$ sort $(f_{\mathrm{1}\dots \lambda },i)$, sort solutions.

(5.5)

$m^{\prime }=m$, replace m.

(5.6)

Move mean to better solutions, update m  $(x_1,\dots ,x_{\lambda })$

(5.7)

Update isotropic evolution path. Update ps ($p_{\sigma },{\sigma }^{-\mathrm{1}}{C}^{-\mathrm{1}/\mathrm{2}}(m-m^{\prime })$)

(5.8)

Update anisotropic evolution path. Update pc ($p_c,{\sigma }^{-\mathrm{1}}(m-m^{\prime }),‖p_{\sigma }‖$)

(5.9)

Update the covariance matrix. Update C (C, $p_c$, $(x_{\mathrm{1}}-m^{\prime })/\sigma ,\dots ,(x_{\lambda }-m^{\prime })/\sigma$)

(5.10)

Update step-size using isotropic path length. Update $\sigma$ ($\sigma ,‖p_{\sigma }‖$)

(5.11)

Return $x$

(5.12)

Calculate the objective function, ${\sum }_{i=\mathrm{1}}^{\mathrm{3}}{L}_i^{\mathrm{2}}(\stackrel{⃑}{\theta })$

(6)

Stop when solution is converged, output the parameters $\stackrel{⃑}{x}=\stackrel{⃑}{\theta }(\alpha ,\beta ,\eta )$.

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 $\stackrel{⃑}{\theta }=(\alpha ,\beta ,\eta )$ with different sample size $n$. The parameters $\stackrel{⃑}{\theta }=(\mathrm{0,2},\mathrm{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 numbers u, v from the uniform distribution on (a, b) and (0,1), respectively.

(2)

Calculate $X$ 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 $X_{\mathrm{s}\mathrm{u}\mathrm{p}}$ 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 $\sigma$ 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 $\stackrel{⃑}{\theta }=(\mathrm{0,2},\mathrm{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 $\stackrel{⃑}{\theta }=(\mathrm{0,2},\mathrm{1.5})$.

Sample size	$$\alpha$$	$$\beta$$	$$\eta$$	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 omitted; refer to PDF]

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 $\stackrel{⃑}{\theta }=(\mathrm{298,2},\mathrm{73})$. In this example, the comparison with intelligent global optimization algorithms was carried out. Firstly, the initial values ${\stackrel{⃑}{\theta }}_{\mathrm{0}}$ of the unknown parameters are calculated by RRM, with a searching space of $[\mathbf{0},\mathrm{5}{\stackrel{⃑}{\theta }}_{\mathrm{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	$$\alpha$$	$$\beta$$	$$\eta$$
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	$$\alpha$$	$$\beta$$	$$\eta$$
$$[\mathbf{0},\mathrm{5}{\stackrel{⃑}{\theta }}_{\mathrm{0}}]$$	298.29	1.9973	72.92
$$[\mathbf{0},\mathrm{10}{\stackrel{⃑}{\theta }}_{\mathrm{0}}]$$	302.76	2.38	77.35
$$[\mathbf{0},\mathrm{50}{\stackrel{⃑}{\theta }}_{\mathrm{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, ${\sigma }_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.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	$$\alpha$$	$$\beta$$	$$\eta$$	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 omitted; refer to PDF]

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 omitted; refer to PDF] [figure omitted; refer to PDF]

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$$\begin{array}{}\text{(11)}& \mathrm{A}\mathrm{I}\mathrm{C}=\frac{\mathrm{2}p-\mathrm{2}L}{n}\end{array}$$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.

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

[1] W. Weibull, "Wide applicability," Journal of Applied Mechanics, 1951.

[2] E. S. Lindquist, "Strength of materials and the Weibull distribution," Probabilistic Engineering Mechanics, vol. 9 no. 3, pp. 191-194, DOI: 10.1016/0266-8920(94)90004-3, 1994.

[3] K. C. Datsiou, M. Overend, "Weibull parameter estimation and goodness-of-fit for glass strength data," Structural Safety, vol. 73, pp. 29-41, DOI: 10.1016/j.strusafe.2018.02.002, 2018.

[4] R. Danzer, P. Supancic, J. Pascual, T. Lube, "Fracture statistics of ceramics—weibull statistics and deviations from Weibull statistics," Engineering Fracture Mechanics, vol. 74 no. 18, pp. 2919-2932, DOI: 10.1016/j.engfracmech.2006.05.028, 2007.

[5] L. Afferrante, M. Ciavarella, 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, DOI: 10.1016/j.ijsolstr.2005.08.002, 2006.

[6] P. K. Chaurasiya, S. Ahmed, 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, DOI: 10.1016/j.renene.2017.08.014, 2018.

[7] P. Wais, "Two and three-parameter Weibull distribution in available wind power analysis," Journal of Renewable Energy, vol. 103, pp. 15-29, DOI: 10.1016/j.renene.2016.10.041, 2017.

[8] 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, DOI: 10.2307/1266719, 1959.

[9] R. Jiang, 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, DOI: 10.1109/24.406588, 1995.

[10] R. Jiang, D. N. P. Murthy, "The exponentiated weibull family: a graphical approach," IEEE Transactions on Reliability, vol. 48 no. 1, pp. 68-72, DOI: 10.1109/24.765929, 1999.

[11] J. A. Greenwood, J. M. Landwehr, N. C. Matalas, 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, DOI: 10.1029/wr015i005p01049, 1979.

[12] J. R. M. Hosking, J. R. Wallis, 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, DOI: 10.1080/00401706.1985.10488049, 1985.

[13] A. C. Cohen, B. J. Whitten, Y. Ding, "Modified moment estimation for the three-parameter weibull distribution," Journal of Quality Technology, vol. 16 no. 3, pp. 159-167, DOI: 10.1080/00224065.1984.11978908, 1984.

[14] S. A. Akdağ, Ö. Güler, "Alternative Moment Method for wind energy potential and turbine energy output estimation," Journal of Renewable Energy, vol. 120, pp. 69-77, DOI: 10.1016/j.renene.2017.12.072, 2018.

[15] N. Balakrishnan, 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, DOI: 10.1016/j.spl.2008.05.019, 2008.

[16] P. Groeneboom, J. A. Wellner, Information Bounds and Nonparametric Maximum Likelihood Estimation, vol. 19, 2012.

[17] W. Chen, M. Xie, 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, DOI: 10.1080/03610918.2014.904520, 2016.

[18] D. Marković, D. Jukić, 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, DOI: 10.1016/j.cam.2008.09.025, 2009.

[19] Y.-J. Kang, Y. Noh, O.-K. Lim, "Kernel density estimation with bounded data," Structural and Multidisciplinary Optimization, vol. 57 no. 1, pp. 95-113, DOI: 10.1007/s00158-017-1873-3, 2018.

[20] T. C. Carneiro, S. P. Melo, P. C. M. Carvalho, 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, DOI: 10.1016/j.renene.2015.08.060, 2016.

[21] A. Yalçınkaya, B. Şenoğlu, U. Yolcu, "Maximum likelihood estimation for the parameters of skew normal distribution using genetic algorithm," Swarm and Evolutionary Computation, vol. 38, pp. 127-138, DOI: 10.1016/j.swevo.2017.07.007, 2018.

[22] C. Leahy, E. J. O'Brien, B. Enright, R. T. Power, "Estimating characteristic bridge traffic load effects using Bayesian statistics," Proceedings of the 12th International Conference on Applications of Statistics and Probability in Civil Engineering (ICASP12), .

[23] M. Keller, A.-L. Popelin, N. Bousquet, 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, DOI: 10.1111/risa.12484, 2015.

[24] R. Aghmasheh, V. Rashtchi, E. Rahimpour, "Gray box modeling of power transformer windings for transient studies," IEEE Transactions on Power Delivery, vol. 32 no. 5, pp. 2350-2359, DOI: 10.1109/TPWRD.2017.2649484, 2017.

[25] T. P. Talafuse, E. A. Pohl, "Small sample reliability growth modeling using a grey systems model," Quality Engineering, vol. 29 no. 3, pp. 455-467, DOI: 10.1080/08982112.2017.1318920, 2017.

[26] H. H. Örkcü, E. Aksoy, M. I. Dogan, "Estimating the parameters of 3-p Weibull distribution through differential evolution," Applied Mathematics and Computation, vol. 251, pp. 211-224, DOI: 10.1016/j.amc.2014.10.127, 2015.

[27] I. Usta, "An innovative estimation method regarding Weibull parameters for wind energy applications," Energy, vol. 106, pp. 301-314, DOI: 10.1016/j.energy.2016.03.068, 2016.

[28] H. H. Örkcü, V. S. Özsoy, E. Aksoy, M. I. Dogan, "Estimating the parameters of 3 - p Weibull distribution using particle swarm optimization: a comprehensive experimental comparison," Applied Mathematics and Computation, vol. 268, pp. 201-226, DOI: 10.1016/j.amc.2015.06.043, 2015.

[29] D. Petković, "Adaptive neuro-fuzzy approach for estimation of wind speed distribution," International Journal of Electrical Power & Energy Systems, vol. 73, pp. 389-392, DOI: 10.1016/j.ijepes.2015.05.039, 2015.

[30] M. Nassar, A. Z. Afify, S. Dey, 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, DOI: 10.1016/j.cam.2017.12.001, 2018.

[31] N. Hansen, A. Ostermeier, "Completely derandomized self-adaptation in evolution strategies," Evolutionary Computation, vol. 9 no. 2, pp. 159-195, DOI: 10.1162/106365601750190398, 2001.

[32] N. Hansen, S. D. Müller, P. Koumoutsakos, "Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)," Evolutionary Computation, vol. 11 no. 1,DOI: 10.1162/106365603321828970, 2003.

[33] N. Hansen, S. Kern, "Evaluating the CMA evolution strategy on multimodal test functions," Parallel Problem Solving from Nature—PPSN VIII, vol. 3242, pp. 282-291, DOI: 10.1007/978-3-540-30217-9_29, 2004.

[34] S. F. Duffy, J. L. Palko, 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, DOI: 10.1115/1.2906663, 1993.

[35] B. Faucher, W. R. Tyson, "On the determination of Weibull parameters," Journal of Materials Science Letters, vol. 7 no. 11, pp. 1199-1203, DOI: 10.1007/BF00722337, 1988.

[36] B. Bergman, "Estimation of Weibull parameters using a weight function," Journal of Materials Science Letters, vol. 5 no. 6, pp. 611-614, DOI: 10.1007/BF01731525, 1986.

[37] J. E. Gentle, Random Number Generation and Monte Carlo Methods, 2006.

[38] S. Acitas, C. H. Aladag, 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.