Unbiased weighted least squares estimation of Weibull modulus with small samples

Abstract

Keywords: Ceramics; Composite Materials; Fracture Strength; Reliability; Weibull Distribution; Weibull Modulus; Weighted Linear Least Squares; Unbiased Estimation

Introduction

The Weibull distribution is one of the most widely used probability distributions in life testing and reliability studies. The distribution was developed by Weibull [1] and has found use in many different areas of application such as wind-speed data analysis [2–5], unemployment duration analysis, survival analysis as well as reliability analysis [6]. One major application area is the modelling statistical variation in the fracture strength of many materials such as advanced ceramics, metallic matrix composites and ceramic matrix composites [7–9], and strength of micro and nano structures [10]; it is also used to describe the

fracture toughness behavior of steels in ductile-brittle transition region [11], and fatigue behavior of metals [12,13]. The two-parameter Weibull distribution function is given by (1) where m is the shape parameter or Weibull modulus, s0 is the scale parameter or characteristic strength of the distribution, and F is the fracture probability of the material at or below uniaxial tensile stress s. Weibull modulus, m, reflects the inherent variability in strength of a material: the higher

*Corresponding author. *E-mail address: b.birgoren@etu.edu.tr This paper was recommended for publication in revised form by Editor-in-Chief Ahmet Selim Dalkilic Published by Yıldız Technical University Press, İstanbul, Turkey © Author. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).

Sigma J Eng Nat Sci, Vol. 43, No. 5, pp. 1683−1692, October, 2025

the m the lower the dispersion of fracture stress. It is generally independent of the size of the material, assuming the material composition and microstructure remain consistent. For example, if test specimens with different cross-sectional areas used in uniaxial compression and extension tests, their Weibull modulus remains the same. The scale parameter, on the other hand, is closely related to the mean fracture stress, and it changes with the cross-sectional area. The Weibull parameters m and s0 are estimated from a sample of strength measurements. High experimental costs limit the number of samples to be tested. Therefore, methods achieving high estimation precision with small sample sizes are highly desirable [14]. In practice, the Ordinary Least Squares (OLS) method and the Maximum Likelihood Estimation (MLE) method have been the most commonly used methods for estimating these parameters in general reliability applications [15]. In materials science, the OLS method remains the most popular due to its simplicity and ease of use [16], while the MLE method is considered a standard method in important references such as [17] due to its well-known distributional optimality properties in large samples [15]. Both of the methods produce biased estimators, therefore there are many studies for bias correction. One approach is using empirical correction factors, [18] is an example for MLE, and [19–21] are some recent studies for OLS; and the other approach is to satisfy unbiasedness through developing a probability index to be used in OLS. Davies [16,22] provides a detailed literature review and a thorough discussion of how such indices can be developed. Despite this focus on the two-parameter model, recent research has explored the potential of the three-parameter Weibull distribution in materials science, although its practical application remains limited [23–27]. One disadvantage of OLS is that it treats each measurement equally under the assumption that the variance of the error term is a constant. However, this assumption is violated for the OLS method used for Weibull parameter estimation. Several authors have proposed the Weighted Least Squares (WLS) method as a remedy, proposing different weight factors [28–33]. Their studies showed that the WLS methods show significantly superior performance as compared to the OLS methods. On the other hand, the WLS methods also result in biased estimators, and there has been no study in the literature for unbiased estimation of Weibull parameters by WLS to the best of our knowledge. Therefore, the main purpose of this study is to develop one such method for the unbiased estimation of the Weibull modulus. As in the OLS case, empirical correction factors and developing a new probability index are two viable approaches. This study prefers the latter because it allows decreasing estimation variance both through weight factors and fine-tuning coefficients used in the probability index. Finally, all the Monte Carlo simulations are coded and run in the R programming language which uses MersenneTwister random number generator as the default generator

whose cycle period is 219937-1 [34]. The following sections are organized as follows: After a concise discussion of the OLS and WLS methods, a systematic analysis is introduced to develop a new probability index for unbiased estimation of Weibull modulus using Monte Carlo simulations as the basic tool. The probability index, proposed as a result of this analysis, is then compared with the MLE, the OLS and the WLS methods using different criteria. Finally, implementation on a practical example is discussed. Ordinary Least Squares (OLS) Method Eq. (1) becomes a straight line by a double logarithmic transformation: (2) In order to apply the OLS method, the measurements are ranked from the smallest to the largest. F-values are assigned according to the rank i of a measurement, si denoting the ith smallest: s1 ≤ s2 ≤ ... ≤ sn. Then, in Eq. (2), F(s) can be replaced by F(si), and ln s can be replaced by ln si. The most commonly used probability estimators of F(si), which are simply denoted as Fi, are median rank, (3a) mean rank, (3b) and hazen rank, (3c) There are many other estimators proposed for this purpose [16]. Considering the familiar form of a regression equation, Y = aX + b, the left side of Eq. (2) corresponds to Y, ln s corresponds to X, m corresponds to a, and -m ln s0 corresponds to b. Using si and Fi pairs in Eq. (2), a and b are obtained by the OLS procedure. Then the Weibull param. eter estimates are calculated as  = a and Weighted Least Squares (WLS) Method Weighted least squares (WLS) extends ordinary least squares regression by incorporating weights into the sum of squared residuals. These weights, denoted by wi for the ith observation, downplay the influence of data points si with higher variance or lower reliability. The regression model in Eq. (2) has non-constant error variance. As a result use of the weighted least squares method gives estimates with better statistical properties: Bergman [30], Faucher and Tyson [29] and Hung [31] proposed using the following weight factors in the WLS method, respectively: (4a)

Sigma J Eng Nat Sci, Vol. 43, No. 5, pp. 1683−1692, October, 2025

where hazen rank is used in Eq. (4a), and median rank is used in Eq. (4b) and (4c). Lu et al. [28] showed that Z = -ln (1 - F(s)) in Eq. (2) follows a standard exponential distribution. Then Eq. (2) takes the form of

and wi are computed by Eq. (4b). As a side note, wi are set to 1 for OLS computations. In both OLS and WLS /m is a pivotal statistic. In Monte Carlo simulations, therefore, choosing m=1 will suffice for all practical purposes. Similarly, the scale parameter will be fixed as s0 = 1.

(5) where Zi is the ith order statistic of the standard exponential variable. Using the expectation and variance of this variable, Lu et al. [28] proposed the following approximate formula for the weights: (6)

A New Probability Index Formula for Unbiased WLS This section and following sections develop a new probability index producing an unbiased estimator of m, that is, E() = m. Numerous probability indices have been proposed for F in addition to Eq. (3a)-(3c). Davies[16] has recently provided a comprehensive list of them. These probability indices are all developed for the ordinary least squares method and conform to the general form of (8)

Note that this formula does not involve a probability index. Also, they went on to compare these four WLS methods: They pointed out that the denominator of Hung’s weight factor in Eq. (4c) is constant, thus Hung’s weight factor is the same as Bergman’s in Eq (4a). They concluded that Bergman’s weight factor (as well as that of Hung) generates a larger Mean Squared Error (MSE) than the others. Faucher and Tyson’s weight factor in Eq. (4b) and their method show similar levels of performance. As a result, we decided to use Faucher and Tyson’s weight factors (Eq. 4b) in this study due to its ease of use. These weight factors were derived independently from median ranks [29] and should approximately work for any probability rank value between 0 and 1. WLS computations are performed by various statistics packages such as MINITAB® and SPSS®. Microsoft Excel® spreadsheets can also be used for this purpose using the following closed-form formulae: (7a) and (7b)

where 0 ≤ a ≤ 1 and b > -a (because 1 - Fi in Eq. (2) has to be between 0 and 1 for all i = 1, ..., n to avoid negative terms in the logarithms). For example, a = 0.3 and b = 0.4 for median rank in Eq. (3a). The coefficients a and b are usually constant [16]; however, some studies aiming at minimizing the coefficient of variation of  [35,36] or making  unbiased [16,37,38] proposed a and b values as functions of n. One problem with this approach is that the values change inconsistently with changing n which leads to the lack of expressions for F which are valid for all materials engineering relevant values of n [16]. Therefore, this study aims to compute a as functions of n, and b as a constant such that WLS produces  values that are unbiased and have minimum variation. The basic simulation procedure employed in this study involves generating a sample of n values from a Weibull distribution with parameters m=1 and s0 = 1; estimating  using WLS with certain a and b values in Eq. (8) (or any index function with two constants such as in Eq. (10) as will be discussed below) and repeating this R times to compute the sample mean and the standard deviation of . In this procedure the sample mean of , ave, can be considered as a stochastic function g(a, b, n, R). For a given n, finding the values of a and b producing an unbiased estimator (ave = 1) is equivalent to finding the root of g(a, b, n, R) -1. The root-finding procedure employed to compute the a-b combination in Eq. (8) is (or in Eq. (10) as will be discussed shortly) as follows: While keeping b at a certain value b', a is increased incrementally by δ, and at any step going from a' to a' + δ, if ave, estimated by the basic simulation procedure, goes from a value less than 1 to a value

Sigma J Eng Nat Sci, Vol. 43, No. 5, pp. 1683−1692, October, 2025

greater than 1, or the vice versa, this indicates that there exists an unbiased estimator for b = b' and the values of a' + δ / 2 and b' can be considered as the approximate values producing an unbiased estimator. The precision of the procedure depends on the choice of δ. This is repeated for all b values in a given search interval in increments of δ: ....b' - δ, b', b' + δ, ... As an initial attempt to produce approximate a and b values in Eq. (8) a small group sample sizes n = 6, 15, 30, 75, 100 and 150, are selected. In Eq. (8) b is usually between 0 and 1 [16], hence 0 ≤ a ≤ 1 and 0 ≤ b ≤ 1. The root-finding procedure was repeated for each sample size with δ =

Comparison With Other Methods

It is tricky to select methods for comparison purposes, because there is no WLS method developed for unbiased estimation of the Weibull modulus m. We selected Davies’ [16] probability index formula as the first method of comparison; it is designed for unbiased estimation of the mean of m (he proposed others for median and mode) to be used with OLS. While OLS can be seen as a special case of WLS, where all the weight are one, OLS and WLS are treated as

two separate methods in this study. Davies [16] showed that the formula results in unbiased estimates with low variation. As the second and third methods of comparison, the MLE method and Faucher and Tyson’s WLS method [29] are selected. While their estimates are biased, they may still produce estimates close to the true value, which can be measured by the Mean Squared Error (MSE): MSE = E( - m)2. MSE is also related with bias: MSE = Var() + (E() - m)2, where the second part of the summation is the square of the bias of . Note that this part is zero for the proposed method and Davies’ method because of their unbiasedness. The basic simulation procedure described previously is used with the selected methods for comparison and selected sample sizes; that is, parameters are set as m=1 and s0 = 1 for Weibull random number generation; the simulation run number is set as R=106. The MSE values as well as the sample means of  are shown in Table 1. First, let’s look at the precision of the unbiasedness: the first column of the proposed method in Table 1 should be ideally 1, but there are slight deviations which are due to the simulation error and the errors incurred during the root finding procedure. Davies’ method is designed for n ≥ 10, this explains the relatively high deviations for n = 6 and 8. Before comparing the methods let’s note that for two estimators, 1 and 2, the ratio of their MSEs are called relative efficiency, defined as RE[1, 2] = MSE(1)/ MSE(2). A RE value of 0.80 implies that the necessary sample size for the method using 1 is 80% of that needed for 2 to achieve approximately equal MSE value [15]. This may translate as the approximate percent saving in testing costs when the method of 1 is preferred to that of 2. When the MSE’s of the two unbiased methods are compared, Table 1 shows that Davies’ method has significantly larger values than the proposed method for all

Sigma J Eng Nat Sci, Vol. 43, No. 5, pp. 1683−1692, October, 2025

sample sizes. It is to the degree that the levels achieved by the proposed method can be achieved with an additional five test specimens (e.g. n = 15 vs n = 20) for n ≥ 20, and with an additional 10 specimen for n ≥ 35; this result can be seen by comparing MSE values for different n values, and also agrees with the RE column (note that RE produces an approximate result). Note that, since these methods are unbiased, MSE of  is equivalent to the sample variance of , (strictly speaking, they are approximately equal due to small deviations from unbiasedness). While the MLE method has been a popular alternative to OLS and WLS, it suffers from severe overestimation bias in small samples. For instance, in Table 1, the MLE estimate for n = 10 (1.1698) overestimates the unknown Weibull modulus by 16.98%. Table 2 decomposes MSE values in Table 1 into bias and variance components, revealing that the bias becomes negligible (below 2%) for sample sizes exceeding n = 70; note that MSE = (Bias)2+Variance (here, for n = 10, 0.1517 = 0.16982 + 0.1229). Table 2 highlights the substantial bias in MLE estimates for small samples significantly impacts the overall MSE.

The MSE values of the MLE method are at least 10% larger than those of the proposed method for n ≤ 30. The difference is so severe that the levels achieved by the proposed method for n ≤ 12 can be achieved by the MLE method with an addition of at least three test specimens. However, as the sample size increases, the MLE methods starts to show better performance as a result of its optimality properties in large samples. For n ≥ 60, the MLE method starts to perform better, and the MSE for the proposed method becomes 14% larger than that of the MLE method at n = 150. When compared with the proposed method, Faucher and Tyson’s method has smaller MSE values for all sample sizes, and the bias is not severe as can be seen in Table

2. The proposed method’s MSE is about 6% above that of

Faucher and Tyson’s method for n ≤ 15, but as the sample size increases, this percentage goes below 2% after n ≥ 35 and below 1% after n ≥ 60. These results do not indicate that Faucher and Tyson’s method is better than the proposed method, because it has significant bias. It underestimates  and approximately 5% underestimation should be expected for n = 6-20. The percentage of underestimation

Table 1. Sample mean and MSE values of  for the compared methods Proposed (P)

Sigma J Eng Nat Sci, Vol. 43, No. 5, pp. 1683−1692, October, 2025

becomes less severe as the sample size increases and goes below 2% for n ≥ 60. Practical Example An application of the proposed method was carried out on the results of an experimental study [40]: 19 identical composite specimens were prepared from quasi-isotropic carbon-epoxy sheets with (0°)3 configuration, 0.89 mm thickness, and 295 gr/m2 weight and the tension experiments were carried out using an Instron 8516+ universal testing machine. The measured fracture strength values are presented in Table 3.

To enhance the application study in terms of sample size, in addition to the original data set given in Table 3, two random subsets of sizes 6, 10, and 15 were drawn from the same data set. These subsets consist of {522, 532.7, 476.5, 521.6, 439, 507.3} for n = 6, {513.6, 552,519, 521.6, 439,450.9, 463.5, 497.5, 476.5, 477} for n =10, and {532.7, 502.5, 442, 519, 502.7, 477, 552, 522, 439, 513.6, 521.6, 450.9, 476.5, 507.3, 463.5} for n = 15. To implement the proposed method, Eq. (7) is used where Fi is computed by Eq. (11) and Eq. (12). The following Table 4 gives the estimated Weibull modulus values. As n increases, the MSE decreases as can be seen from Table 1,

Table 2. Decomposition of MSE values of  for biased estimation methods MLE (M)

Table 3. Fracture strength of carbon-epoxy composite material specimens (megapascals) Test No.

Table 4. Weibull modulus estimates () for the fracture strength data n

Sigma J Eng Nat Sci, Vol. 43, No. 5, pp. 1683−1692, October, 2025

hence estimates with larger sample sizes have higher precision. Therefore, the estimates for n = 19 should be closer to the true, unknown Weibull modulus value. The estimates for n = 19 are close to each other due to their small MSE values (check the values for n = 20 in Table 1). However, the smallest MSE’s belong to the proposed WLS method and Faucher and Tyson’s WLS method, thus their estimates, 17.6756 and 16.916 can be argued to be closer to the true m value. For smaller n values the estimates are farther away from each other, due to increased MSE and increased bias. Consider n = 6: MSE is a measure of expected distance of an estimate from the true value, and they are 5-10 times as large as those for n = 20. Also, due to the substantial overestimation problem in the MLE method, which can be seen in Table 2, MLE estimates are larger than the others in Table 4.

Conclusion

This research introduced a novel probability index that, for the first time, offers unbiased WLS estimators of the Weibull modulus with a satisfactory degree of precision. To this end, a systematic study was carried out to compute parameters of a probability index to be used with the WLS method for the unbiased estimation of the Weibull modulus. First Faucher and Tyson’s weights were selected among several weighting schemes. Then two parameters of a generic probability index were computed by using a step-by-step approach. The result is a simple index formula presented in Eq. (10) wherein the single parameter a, produced from simulations, is fitted as a function of the sample size in Eq. (11). Use of two parameters in this study allowed variance reduction of the estimator as well as satisfying unbiasedness. The proposed method requires only the use of the new probability index along with Faucher and Tyson’s weights. It is the only WLS method with unbiased Weibull modulus estimator in the literature. When compared with a recent OLS method proposing a new probability index producing unbiased Weibull modulus estimators, it exhibited significantly smaller variance. MSE is used for comparison among methods producing biased estimators. Results showed that the proposed method performs better than the MLE method for sample sizes less than 60, and its performance is close to that of Faucher and Tyson’s method, but with slightly larger MSE values. The advantage of the proposed method is most prominent in the range of sample sizes smaller than 125 where the percent bias in the aforementioned biased methods is greater than 1% and they should be used with caution. In materials science, high experimental costs limit the number of samples to be tested; and in general, smaller the sample size are preferred. In general, sample sizes are less than 60 in such experiments, which makes the proposed method advantageous in terms of its smaller MSE and unbiasedness property. While this paper focuses on

materials science applications, it should be pointed out that the proposed method can be used in other application areas such as wind speed analysis and reliability analysis. In conclusion, the implementation of the proposed WLS method is almost as simple as that of the OLS method and we hope that this helps materials scientists and engineers who generally prefer OLS and MLE in practice.

Data Availability Statement

The authors confirm that the data that supports the findings of this study are available within the article. Raw data that support the finding of this study are available from the corresponding author, upon reasonable request.

Conflict Of Interest

The author declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Ethics

There are no ethical issues with the publication of this manuscript.

Statement On The Use Of Artificial Intelligence

Artificial intelligence was not used in the preparation of the article.

Share and Cite