Estimation of the Shape Parameter of a Wear-Out Failure Period for a Three-Parameter Weibull Distribution in a Small Sample

We propose a method to estimate a shape parameter for a three-parameter Weibull distribution. The proposed method first derives an unbiased estimator for the shape parameter independent of the location and scale parameters and then estimates the shape parameter using a minimum-variance linear unbiased estimator. Since the proposed method is expressed using a hyperparameter, its optimal hyperparameter is searched using Monte Carlo simulations. The recommended hyperparameter used for estimating the shape parameter depends on the sample size, and this causes no problems since the sample size is known when data is obtained. The proposed method is evaluated using a bias and a root mean squared error, and the results are very promising when the population shape parameter is 2 or more in the Weibull distribution representing the wear-out failure period. A numerical dataset is analyzed to demonstrate the practical use of the proposed method.


Introduction
Weibull distribution is often used in lifetime and reliability studies. The probability density function and cumulative distribution function of the three-parameter Weibull distribution are expressed as follows: where η > 0, m > 0, and γ < x are the scale, shape, and location parameters, respectively. Although maximum likelihood estimation (MLE) is often used for parameter estimation, it may not be possible to be used for a three-parameter Weibull distribution. Therefore, the parameter estimation method that extends MLE have proposed for a three-parameter Weibull distribution. Hall & Wang (2005) proposed the Bayesian likelihood (BL) method, which multiplies MLE by an empirical prior. Cousineau (2009a) proposed the weighted maximum likelihood estimation (w-MLE) method that added the three weights to MLE. Nagatsuka et al. (2013) proposed the location and scale parameters free maximum likelihood estimators (LSPF-MLE) method, which estimates the shape parameter using the independent statistics of the location and scale parameters. The parameter estimation method of the three-parameter Weibull distribution has also been studied by other researchers (Smith, 1985;Lawless, 2002;Murthy et al., 2004;Nagatsuka & Kamakura, 2005;Cousineau, 2009b;Ng et al., 2012;Sugiyama et al., 2019).
We propose a method to estimate the shape parameter, independent of the location and scale parameters. Therefore, the shape parameter can be estimated first for the three-parameter Weibull distribution. Estimating the shape parameter with a high accuracy is expected to improve the estimation accuracy of the location and scale parameters.
In Section 2, we derive an unbiased estimator of the shape parameter that does not depend on the location and scale parameters. In Section 3, we search for the recommended hyperparameters for each sample size, and compare the effectiveness of the proposed method with the existing methods. Section 4 describes an attempt at estimating the shape parameter using a numerical example. Our conclusions are presented in Section 5.

Unbiased Estimator of the Shape Parameter Independent of the Location and Scale Parameters
Let X 1 , . . . , X n be random samples from the three-parameter Weibull distribution, and let X (i) denote the i-th order statistic, γ < X (1) < · · · < X (n) (X (i) X ( j) for some i j), i = 1, . . . , n (n ≥ 2). Using the inverse transform method, the following formula is established: where U (i) is the i-th order statistics of the uniform distribution, 0 ≤ U (1) < · · · < U (n) ≤ 1 (U (i) U ( j) , for some i j).
The following results are obtained by double logarithmic conversion of both sides: Subtracting the formula in which i = 1 is substituted into (4) from the formula in which i = 2, . . . , n is substituted into (4), it leads to: The right-hand side of (5) is independent of η. The mean of both sides in (5) is given by: where The expected value of (6) is the function of only i and n. The unbiased estimator of 1/m is expressed as follows: Furthermore, let be the empirical prior of γ, where U is the standard uniform distribution, and c > 0 is a hyperparameter. Both sides of (9) have the same range (0, ∞). The hyperparameter c is chosen to obtain an optimal estimator. Substituting (9) into (8), we obtain: Using the minimum-variance linear unbiased estimator (Ahsanullah et al., 2013), m is estimated by: where 1 = (1 · · · 1) is (n − 1)-dimensional vector, var log( Since η and γ are not included in (11), m can be estimated first among the three parameters. Note that S appears to be expressed with γ, but substituting (13) and (14) into (12) gives an expression without γ. Although (13) and (14) seem to depend on m, they are canceled by substituting into (11), andm can then be estimated without depending on m. Since a proof by the expansion of (11) is complicated calculations, we show thatm has the same value regardless of the setting of m on the right-hand side of (13) and (14) in Appendix B. The optimal hyperparameter c in (10) depends on the sample size; thus, its optimal value is searched using Monte Carlo simulations in Section 3.

Search for Optimal c
We set five cases (γ = 0, η = 1, m = 1, 2, 3, 4, 5) of sample size n = 5, 10, 15, 20 to search for the optimal value of the hyperparameter c. Table 1 shows the Monte Carlo simulation study results of the bias and root mean squared error (RMSE = variance + bias 2 ) based on 10 000 Monte Carlo simulations for each hyperparameter c. We used the software R version 3.6.2 (R Core Team, 2019) for the Monte Carlo simulations. It can be seen that the optimal hyperparameter c values vary greatly with respect to n, and also change depending on m, however, the differences are not large. Since n is known and m is unknown when data is obtained, considering the bias and RMSE from the Monte Carlo simulation study results, we use the values of c = 1.5, 2.3, 3.4, and 4.3 for n = 5, 10, 15, and 20, respectively.

Existing Shape Parameter Estimation Methods
We compare the shape parameter estimated by the proposed method to those by three existing methods. The first method is the BL method proposed by Hall & Wang (2005), which multiplies MLE by an empirical prior. The location and shape parameters are estimated by maximizing The second method is the w-MLE method proposed by Cousineau (2009a), which extends MLE technique by incorporating three weights into the solutions to the maximum likelihood equations. The author states that the weight of W 3 is not a constant but changes during the search as new values of the estimating shape parameter are explored. However, the paper used the linear interpolation, because of the lack of a fast computational algorithm. The w-MLE method used for comparison in this paper is based on Cousineau (2009a), but the weight of W 3 is calculated by the following five steps based on Nagatsuka et al. (2013). (i) A temporary shape parameter,m 0 , is estimated by MLE for the two-parameter Weibull distribution in (X (i) − X (1) ), i = 2, . . . , n. (ii) Uniform random samples, (u 1 , . . . , u n ), are generated from the standard uniform distribution. (iii) Usingm 0 and (u 1 , . . . , u n ), we calculate (iv) Steps (i)-(iii) are repeated independently 30 000 times. (v) Let W 3 be the median of Y 3 . In the Monte Carlo simulations, Steps (i)-(v) are performed every time random samples from the three-parameter Weibull distribution. The location and shape parameters were estimated by minimizing International Journal of Statistics and Probability Vol. 9, No. 6; The third method is the LSPF-MLE method (Nagatsuka et al., 2013), in which the shape parameter is estimated by maximizing where Z (i) = X (i) −X (1) X (n) −X (1) . As a feature of each estimation method, the BL and w-MLE methods first estimate the location and shape parameters simultaneously without the scale parameter. The LSPF-MLE method first estimates the shape parameter without the location and scale parameters, which is also the feature of the proposed method.

Effectiveness of the Proposed Method
To demonstrate the validity of the estimations of the shape parameter, five cases (γ = 0, η = 1, m = 1, 2, 3, 4, 5) of sample size n = 5, 10, 15, 20 were investigated. We compared the shape parameter estimated by the proposed method with the shape parameters estimated by the BL, w-MLE, and LSPF-MLE methods in the Monte Carlo simulations. When the shape parameters estimated by the BL, w-MLE, and LSPF-MLE methods were 10 or more, we treated them as estimated at 10 in the Monte Carlo simulations. Since the Monte Carlo simulations were set with m ≤ 5, the bias and RMSE became too large using the estimation result of 10 or more. However, even if the shape parameter estimated by the proposed method was 10 or more, the estimation result was used to calculate the bias and RMSE. The bias and RMSE were calculated from 10 000 Monte Carlo simulations for each set of configurations, and the results are summarized in Table 2. The RMSE of the shape parameter estimated by the proposed method was the smallest of the four methods at m = 2, 3, 4, 5 in the Weibull distribution representing the wear-out failure period. Table 2 also shows the frequencies of the shape parameter estimation of 9.999 or more. In the maximization or minimization of the function using the software R, the frequencies of 9.999 or more were shown because a value slightly less than 10 may be output, when searching in the range where 10 is the upper limit. The shape parameters estimated by the BL and LSPF-MLE methods often exhibited high frequencies of being estimated at 9.999 or more. Figure 1 shows the histogram of the Monte Carlo simulation study results for m = 5 and n = 10, 20. The biases of the shape parameters estimated by the BL and LSPF-MLE methods seem to be small when m = 5, but it can be seen from the histogram that the bias became small by setting the maximum value to 10. The proposed method gave estimates that were closest to the population shape parameter at m ≥ 2 in the Weibull distribution representing the wear-out failure period.  Vol. 9, No. 6; The proposed method used the values of c = 2.3, 3.4, and 4.3 for n = 10, 15, and 20, respectively. Freq: Frequencies with which the shape parameter was estimated at 9.999 or more.

Estimation of the Location and Scale Parameters
The location and scale parameters were estimated by the BL, w-MLE, and LSPF-MLE methods using the shape parameter estimated by the proposed method. To distinguish the location and scale parameters estimated by those methods from the original methods, in this paper they are referred to as the modified BL, modified w-MLE, and modified LSPF-MLE methods, respectively. The location and scale parameters estimated by the BL, modified BL, w-MLE, modified w-MLE, LSPF-MLE, and modified LSPF-MLE were compared in the Monte Carlo simulations with the same conditions as in Section 3.3. The bias and RMSE of the location and scale parameters were calculated from 10 000 Monte Carlo simulations for each set of configurations, and the results are summarized in Tables 3 and 4, respectively. Using the shape parameter estimated by the proposed method, the RMSE of the location and scale parameters are often smaller than the RMSE of the original methods at m = 2, 3, 4, 5. By estimating the shape parameter with a high accuracy, the estimations of the location and scale parameters became highly accurate, too.
We used the software R for estimating the shape parameters by the four methods (the proposed, BL, w-MLE, and LSPF-MLE methods), and a sample code is shown in Appendix A. From the Monte Carlo simulation study results at n = 20 in Section 3.1, we adopted c = 3.4 as the hyperparameter of the proposed method. The shape parameter was estimated to be 4.048 by the proposed method and 4.772 by the w-MLE method, respectively. Note that since the w-MLE method calculates the median by the Monte Carlo simulations, the estimated value will change slightly each time. The shape parameters were not estimated by the BL and LSPF-MLE methods, because BL(m, γ) and LSPF(m) increased monotonically with m. Since the upper limit in the sample code was set to 10, they were estimated at 10 or less, but it can be seen that the shape parameter was estimated to be infinite, when the initial value and the upper limit were changed. A graph can be drawn to make sure they are monotonically increasing. Since the shape parameters estimated by the proposed and w-MLE methods had a negative bias from the Monte Carlo simulation study results at m ≥ 3 and n = 20, the population shape parameter was considered to be m ≥ 4. As m ≥ 2 was assumed from the numerical example, we can recommend the use of the shape parameter estimated by the proposed method.

Conclusions
We proposed an unbiased estimator of the shape parameter that did not depend on the location and scale parameters. The proposed method was expressed including the hyperparameter, and the optimal hyperparameter depended on the sample size. The procedure of the proposed method first looks for the optimal hyperparameter from the data obtained, and then estimates the shape parameter. The Monte Carlo simulation study results showed that the RMSE of the shape parameter estimated by the proposed method was smallest of the four methods when the population shape parameter was 2 or more in the Weibull distribution representing the wear-out failure period. In particular, the difference between the RMSE of the shape parameter estimated by the proposed method and the RMSE of the shape parameters estimated by the other three methods increased with m. The proposed method is recommended, when the population shape parameter of 2 or more is assumed, as in the numerical example.