Robust Reliability Estimation for Lindley Distribution—A Probability Integral Transform Statistical Approach

: In the modeling and analysis of reliability data via the Lindley distribution, the maximum likelihood estimator is the most commonly used for parameter estimation. However, the maximum likelihood estimator is highly sensitive to the presence of outliers. In this paper, based on the probability integral transform statistic, a robust and e ﬃ cient estimator of the parameter of the Lindley distribution is proposed. We investigate the relative e ﬃ ciency of the new estimator compared to that of the maximum likelihood estimator, as well as its robustness based on the breakdown point and inﬂuence function. It is found that this new estimator provides reasonable protection against outliers while also being simple to compute. Using a Monte Carlo simulation, we compare the performance of the new estimator and several well-known methods, including the maximum likelihood, ordinary least-squares and weighted least-squares methods in the absence and presence of outliers. The results reveal that the new estimator is highly competitive with the maximum likelihood estimator in the absence of outliers and outperforms the other methods in the presence of outliers. Finally, we conduct a statistical analysis of four reliability data sets, the results of which support the simulation results.


Introduction
Reliability is defined as the ability of a system or component to perform its required functions under stated conditions for a specified period of time [1]. In reliability theory, various aspects of reliability, probability, statistics and stochastic modeling are studied in combination with engineering principles in the design and scientific understanding of failure mechanisms [2]. Reliability analysis has been utilized to analyze data from various fields, including engineering, medicine, biology, ecology, economics, sociology and the social sciences [3][4][5]. In some areas, reliability data analysis is also referred to as lifetime, failure-time, survival or event-time data analysis [3]. In the analysis of reliability data, reliability properties are often defined using the mean time to failure, reliability function and failure rate function.
Parametric statistical distributions are often used to model and analyze reliability data. The advantages of applying a parametric model in the analysis of reliability data are as follows-a parametric distribution can be described concisely based on only a few parameters rather than having to report an entire curve and also provides smooth estimates of failure time distributions [3]. Several useful parametric models that are often considered for reliability analysis are the exponential, Weibull, gamma and lognormal distributions [3][4][5]. Due to its mathematical simplicity, the exponential distribution is recognized as the most popular and widely used model for reliability data analysis [3,4].
Let X 1 , X 2 , . . . , X n be a random sample of size n from the Lindley distribution with the PDF as shown in Equation (1). As reported by Ghitany et al. [6], the ML estimator for the parameter θ is given by the following:θ where X is the sample mean. Note that, as explained by Ghitany et al. [6], the method of moments (MOM) estimator for the parameter θ is similar to that obtained by the ML estimator.
By numerically solving the non-linear equation shown in Equation (9), the OLS estimatesθ OLS is obtained.
WLS estimates for parameter can be obtained by first minimizing the following function with respect to θ: (n + 1) 2 (n + 2) i(n − i + 1) F(x (i) ; θ) − i n + 1 2 . (11) Apart from that, the WLS estimatesθ WLS can also be obtained by solving the following non-linear equation: Here, the variable ∆(x, θ) is given by Equation (10).

M-Estimators
M-estimators are generalized ML estimators that provide tools for measuring the robustness of ML estimator. As stated by Huber [33], an estimator δ n is defined either by is called an M-estimator. ρ is a measurable function on X × Θ and ψ(x, δ) = (∂/∂δ)ρ(x, δ) is the derivative of the function ρ with respect to δ (when it exists). Note that if ρ(x, δ) = − log f (x; δ), thenδ n is the ordinary ML estimator.

Efficiency Measure: Asymptotic Relative Efficiency
The ML estimator is well known to be efficient. For this reason, the ML estimator is useful in providing a quantitative benchmark for the measure of efficiency. Note that the ML estimator for the parameter of the Lindley distribution θ, which is given in Equation (7), is asymptotically normal with mean θ and variance θ 2 (θ + 1) 2 /(n(θ 2 + 4θ + 2)), that is, For any competing estimator of parameter θ, say,θ 0 the asymptotic relative efficiency (ARE) is defined as the ratio of the asymptotic variance of the ML estimator to the asymptotic variance of the competing estimator, which can be written as follows: In other words, the ARE measures the relative efficiency of the competing estimatorθ 0 as compared withθ ML . When the ARE value gets closer to 1, the efficiency of the estimatorθ 0 becomes closer to the efficiency of the ML estimator.

Robustness Measures: Influence Function and Breakdown Point
The breakdown point (BP), which is useful for assessing the robustness of a statistical approach, measures the degree of sensitivity of an estimator to data contamination. The BP is defined as the largest proportion of contamination that can be tolerated by an estimator before breaking down [20]. An estimator with a higher BP is more robust against data contamination. Note that there are two types of BP, a lower breakdown point (LBP) and an upper breakdown point (UBP). In the present context of θ estimation, the LBP indicates the largest proportion of lower contamination that can be tolerated by an estimator before forcingθ → ∞ and the UBP is the largest proportion of upper contamination that can be tolerated by an estimator before forcingθ → 0 . Note, however, that since contamination of the upper end of the distribution is of greater interest in most typical applications, we emphasize only the UBP in this study. As mentioned by Huber [33], an estimator that has an unbounded ψ function has a BP equal to 0. Note that the ψ function for the ML estimator is ψ(x, For a finite sample, if a single observation x i → ∞ , then n i=1 x i → ∞ and consequently,θ ML → 0 . This result shows that the function ψ(x,θ ML ) is unbounded in x and consequently suggests thatθ ML has a UBP equal to 0. It can be observed that even an extreme value of a single contaminated data in the upper tail of the observations would contribute to the unreliable performance ofθ ML .
Another approach for measuring robustness is the use of an influence function (IF). An estimator is considered to have desirable robustness if it has a bounded IF [33][34][35]. According to Hampel et al. [34] (p. 101), the IF of an estimator δ n that satisfies Equation (14) is defined by the following: where δ(F) denotes the solution δ n of Equation (14) with samples generated from the CDF F. For the ML estimator, since the function ψ(x,θ ML ) is unbounded in x, its IF is also unbounded. Therefore, it is clear thatθ ML is not a robust estimator of parameter θ. For the case when ρ is not differentiable, the IF can be obtained by using the following expression [34] (p. 84): where ∆ x denotes the probability measure that puts mass 1 at point x.

New Robust M-Estimator for the Parameter of the Lindley Distribution
In this section, a new robust M-estimator is proposed based on the probability integral transform statistic. We also discuss the ARE and robustness of this new estimator.

Probability Integral Transform Statistic Estimator
Let X 1 , X 2 , . . . , X n be a random sample from a Lindley distribution. Since the CDF of the Lindley distribution in Equation (2) is continuous and strictly increasing, it can be seen that the random variables F(X 1 ), F(X 2 ), . . . , F(X n ) follow a standard uniform distribution, that is, F(X)~U(0, 1). The new robust M-estimator for parameter θ, which we refer to as the probability integral transform statistic (PITS) estimator, is defined by: where τ > 0 denotes a tuning parameter that will be used later to adjust the balance between efficiency and robustness. Notice that when τ = 1, [1 + (θX i /(1 + θ))]e −θX i = 1 − F(X i ), is a random variable which follows the standard uniform distribution. Assume that u 1 , u 2 , . . . , u n is a random sample from a standard uniform distribution. Based on the strong law of large numbers, it can be shown that n −1 n i=1 u τ i converges to E[u τ ] = 1/(τ + 1) as n → ∞ with probability 1. Therefore, the PITS estimator of parameterθ is defined as the solution with respect to θ to the following equation: Note that Equation (19) can be solved numerically using a method such as the Newton-Raphson, bisection or secant method. It is clear that the PITS estimator for the parameter of the Lindley distribution given in Equation (19) is a class of M-estimator with: Lemma 1. For any fixed τ > 0, Equation (19), that is, H n,t (θ) = 1/(τ + 1), has exactly one solution.

ARE of the PITS Estimator
By applying Equation (15), the ARE of the PITS estimator for parameter θ is given by: To compute the ARE of the PITS estimator, that is,θ PITS , we must find the asymptotic distribution ofθ PITS . One way to obtain the asymptotic distribution ofθ PITS is by following Corollary 2.5 in Chapter 3 of Huber [33]. However, in the present context of the PITS estimator, this corollary cannot be applied since the function is rather complicated to obtain. Thus, the variance ofθ PITS cannot be determined. As an alternative we apply the Monte Carlo simulation method to estimate the variance ofθ PITS , that is, Var(θ PITS ).
As mentioned above, the balance between efficiency and robustness of the PITS estimator can be adjusted by changing the value of the tuning parameter τ. Note that as the value of τ increases, the ARE decreases. Simply put, when the value of τ increases, the PITS estimator attains robustness but loses its relative efficiency. By taking a τ value close to 0, the ARE of the PITS estimator can be made arbitrarily close to 1. Table 1 presents the AREs of the PITS estimator and the corresponding tuning parameter values obtained from 10,000 simulation runs.

BP of PITS Estimator
To obtain the UBP and LBP of the PITS estimator, we apply the argument presented in Finkelstein et al. [29]. The UBP and LBP of the PITS estimator are presented in Theorem 1.
Proof of Theorem 1. For any integer 1 ≤ k ≤ n, the estimatorθ PITS is defined as: For simplicity, let It can be shown that Therefore, h(x) is strictly decreasing with respect to x, for x > 0. It is also worth noting that h(x) > 0 and h(x) < 1 for any x > 0.
We proceed to prove for UBP first. Assume that x 1 , x 2 , . . . , x k takes on values that approach ∞. Let x min = min{x 1 , . . . , x k }. For any ε > 0, suppose Then, it can be shown that The second inequality above is due to the fact that x > log(x) for any x > 0. It follows that since h(x) is strictly decreasing with respect to x and is always positive for x > 0. Therefore, by definition, From Equation (23), we can write since h(x) < 1 for any x > 0. This is valid if and only if k < nε + nτ/(τ + 1). Thus, the finite sample UBP is nτ/(τ + 1) /n, where · is the ceiling function. By letting n → ∞ in the finite sample UBP, we find that UBP equals τ/(τ + 1). For LBP, suppose that x 1 , x 2 , . . . , x k takes on values that approach 0. Let x max = max{x 1 , . . . , x k }. For any 0 < ε < k/n, suppose Then, it follows that If ε ≥ k/n, then the above inequality is true for all From Equation (23), we can write since h(x) > 0 for any x > 0. This is valid if and only if k < nε + n/(τ + 1). Therefore, the finite sample LBP is n/(τ + 1) /n and by taking n → ∞ , the LBP is equal to 1/(τ + 1).
Based on Theorem 1, Table 2 lists the UBPs and LBPs of the PITS estimator for different ARE levels. We can see that as the ARE level decreases, the UBP increases, which suggests that the robustness of the PITS estimator is increasing against upper contamination. On the other hand, the LBP decreases as the ARE level decreases, which means that the robustness of the PITS estimator decreases against lower contamination.

IF of PITS Estimator
According to Marona [35], the IF of an estimator is an asymptotic version of its sensitivity curve (SC). In general, the SC measures the sensitivity of an estimator to the location of the outlier x 0 for a particular random sample. Since using Equation (16) to obtain the IF of the PITS estimator is rather complicated, we apply its SC as an approximation of IF. The SC of the PITS estimatorθ PITS for the random sample x 1 , x 2 , . . . , x n is defined as a function of the location of the outlier x 0 , which can be written as follows: whereθ PITS (x 1 , x 2 , . . . , x n ) is the PITS estimator using the samples x 1 , x 2 , . . . , x n . To obtain the SC of the PITS estimator, we generated a random sample of size n = 50 from the Lindley distribution with parameter θ = 0.5, 1, 2, 3. Then, a single outlier x 0 is added to the sample data. The value of x 0 changed from 0 to 50 in increments of 1. Note that the 99th percentile are often used as upper boundaries for outlier identification, where any data point that exceed the 99th percentile can be considered as a potential outlier. For each parameter θ = 0.5, 1, 2 and 3, the 99th percentile of Lindley distribution are 12.4940, 5.9902, 2.8330 and 1.8222, respectively. Thus, the value of x 0 that exceed these 99th percentiles of Lindley distribution is large enough to be considered as a potential outlier. We obtained the SC(x 0 ) of the PITS estimator with several different AREs for each value of x 0 , as shown in Figure 1.

Simulation Study
In this section, we investigate the performance of the ML, OLS, WLS and PITS (98%, 90%, 80%, 70% and 60% AREs) estimators in the absence and presence of outliers via a simulation study. In Sections 6.1 and 6.2, we present the framework and results of the simulation study, respectively. Then, in Section 6.3, we provide some guidelines for selecting the appropriate ARE of the PITS estimator for practical application purposes.

Simulation Framework
The simulation procedure goes as follows: Step 1: Generate a random variable from the Lindley distribution for two sample size settings, that is, small (n = 30, 50, 70) and large (n = 100, 300, 500), with parameter θ = 0.5, 1, 2, 3. Based on Figure 1, we can see that as x 0 increases, the SC of the PITS estimators with 98%, 90%, 80%, 70% and 60% AREs converge to a certain limit, which indicates that the curve is bounded. This suggests that the PITS estimator has a bounded IF and is robust against the location of a single outlier. As the level of ARE decreases, the limit of each curve moves closer to 0, which suggests that the PITS estimator becomes more robust.

Simulation Study
In this section, we investigate the performance of the ML, OLS, WLS and PITS (98%, 90%, 80%, 70% and 60% AREs) estimators in the absence and presence of outliers via a simulation study. In Sections 6.1 and 6.2, we present the framework and results of the simulation study, respectively. Then, in Section 6.3, we provide some guidelines for selecting the appropriate ARE of the PITS estimator for practical application purposes.
Step 2: Randomly select some observations and replace them with outliers generated from the Lindley distribution with parameter 0.05θ. Note that by multiplying the true value of parameter θ with 0.05, the Lindley distribution will have heavier upper tail and will produce larger values of random variables, which are interpreted as outliers. For the small sample sizes, generate outliers for several fixed numbers, m = 0, 1, 3, 5. For the large sample sizes, simulate outliers for several fixed proportions, ε = 0%, 1%, 5%, 10%.
Step 5: Calculate the performance of each estimator using the percentage relative root mean square error (RRMSE). For a given true value of parameter θ, the RRMSE is given by the following: whereθ i is the estimated parameter for the i-th (i = 1, 2, . . . , N) sample and N is the number of simulated samples. A smaller RRMSE value indicates that the estimator is more accurate and precise. Thus, any estimation method that minimizes the RRMSE provides the best estimation of the parameter θ.

Simulation Results
Tables 3-8 list the results of the simulation study based on the obtained RRMSEs. We summarize these results as follows: 1.
In the cases of both small and large sample sizes (Tables 3-8), we found the following: • When there are no outliers (m = 0 and ε = 0%) in the data, both the ML and PITS (98% ARE) estimators perform similarly and slightly outperform the OLS and WLS estimators.

•
In the presence of outliers, the performance of the ML estimator is much worse than that of the other estimators. As the degree of contamination increases, the performance of the ML estimator deteriorates significantly.

•
The OLS and WLS estimators are quite robust and offer some protection against outliers.

•
As the sample size increases, for all the cases considered, the performance of all the estimators improve, with the RRMSE values becoming smaller.

2.
In the case of a small sample size (Tables 3-5), we found the following: • When the number of outliers is small (m = 1), the PITS (90% ARE) estimator performs best, that is, slightly better than the OLS and WLS estimators.
• When the number of outliers is moderate (m = 3), the OLS, WLS and PITS (70% and/or 60% AREs) estimators perform almost equally well and are considered to be the best methods for this particular case.

•
When the number of outliers is large (m = 5), the PITS (60% ARE) estimator performs best, that is, slightly better than the OLS and WLS estimators.

3.
In the case of a large sample size (Tables 6-8), we found the following: • When the proportion of outliers is small (ε = 1%), for n = 100, the PITS (90% ARE) estimator performs best, that is, slightly better than the OLS and WLS estimators. For n = 300, 500, the OLS, WLS and PITS (90% and 80% AREs) estimators perform almost equally well and are considered to be the best methods for this particular case.

•
When the proportion of outliers is large (ε = 10%), the performance of the PITS (60% ARE) estimator also surpasses that of other methods. The best method for each case is written in bold.

17.05
The best method for each case is written in bold.

14.21
The best method for each case is written in bold. The best method for each case is written in bold. The best method for each case is written in bold. The best method for each case is written in bold.
Overall, in the presence of outliers, the use of the ML estimator should be avoided since it is not robust and provides no protection against outliers. It is also interesting to note that we found the OLS and WLS estimators to be quite robust and able to offer some protection against outliers. However, the proposed PITS estimator provides the most flexible approach as it can be applied in both the absence and presence of outliers.

Some Guidelines for Selecting the Appropriate ARE of PITS Estimator
Based on the results of a comprehensive simulation study, here we provide some guidelines for selecting the appropriate ARE of PITS estimator in practical application:

1.
When there are no outliers in the data, the PITS (98% ARE) estimator should be applied for estimating parameter θ.

Applications and Discussion
In this section, we report four applications of the Lindley distribution as a reliability model using real data sets and compare the performance of the ML, OLS, WLS and PITS (with several AREs) estimators. Note that the AREs of PITS estimator are determined based on the proportion of outliers found in the data. The first data set (Data Set 1) consists of the time to failure of 18 electronic devices reported by Wang [36]. The second data set (Data Set 2) represents the survival times of 44 patients suffering from head and neck cancers (treated using radiotherapy and chemotherapy), which were initially reported by Efron [37] (see also Reference [38]). The third data set (Data Set 3) consists of an uncensored data set of the remission times of a random sample of 128 bladder cancer patients, which was obtained from Lee and Wang [5]. Finally, the fourth data set (Data Set 4) represents the length of stay of 300 patients suffering from breast cancer, which can be found in Reference [39]. All four data sets are given in Appendix A.
For all the data sets considered, to identify the presence of outliers, we applied the generalized boxplot method [40], which is suitable for skewed and/or heavy-tailed distributions. In practical applications, we suggest to apply this method for the purpose of outlier detection for the data that follows the Lindley distribution. Table 9 provides descriptive statistics for all the data sets and Figure 2 shows generalized boxplots of these data sets. To compare the performance of the considered methods in estimating the parameter of the Lindley distribution, we applied the Kolmogorov-Smirnov (K-S) test, Akaike information criterion (AIC) and Bayesian information criterion (BIC) to assess goodness of fit. The best method is determined by that yielding the highest p-value of the K-S test and the smallest value of K-S statistic,  To compare the performance of the considered methods in estimating the parameter of the Lindley distribution, we applied the Kolmogorov-Smirnov (K-S) test, Akaike information criterion (AIC) and Bayesian information criterion (BIC) to assess goodness of fit. The best method is determined by that yielding the highest p-value of the K-S test and the smallest value of K-S statistic, AIC and BIC. Table 10 lists the estimated parameters and goodness of fits obtained for the Lindley distribution for all the data sets. As presented in Table 10, we can observe that the PITS estimator provides a better estimation of the Lindley parameter than the ML, OLS and WLS estimators based on its smallest K-S statistic and highest p-value of the K-S test. On the other hand, the ML estimator is found to be the best method for estimating the parameter θ based on its smallest value of AIC and BIC. This is because the ML estimator maximized the likelihood function of the Lindley model and as a result the AIC and BIC will always favor the ML estimator. Thus, in our case here, we could say that the AIC and BIC are biased measures for goodness of fit assessment. To further support our claim regarding the AIC and BIC, we provide another example where the data is simulated from Lindley distribution with θ = 2 for sample size n = 100. Then, we random select 5% of the observations and replace them with outliers. Based on this data, we compare the performances of all methods in estimating the parameter of the Lindley distribution. The K-S test, AIC and BIC are utilized to assess goodness of fit. The result of the comparative study is presented in Table 11. The best method is written in bold. The best method is written in bold.
Based on the result from Table 11, it can be seen that the estimated parameter found based on ML estimator (θ = 0.81958) is much deviated from the true value of θ = 2. However, if the AIC and BIC are used for determining the best model, it can be observed that the best fitted Lindley distribution is found when the parameter θ is estimated using ML estimator. This further support our claim that AIC and BIC will always favor the ML estimator although the outliers are present in the data. On the other hand, K-S statistic and p-value of K-S test are able to determine the best method for estimating the parameter θ which is PITS (85% ARE). Since the AIC and BIC would provide a biased measures for goodness of fit assessment, thus, in this study, the best Lindley model is only determined based on the smallest value of K-S statistic and the highest p-value of K-S test.
As we mentioned before, based on the smallest K-S statistic and highest p-value of the K-S test, the PITS estimator is found to provide a better estimation of the Lindley parameter than the ML, OLS and WLS estimators. This result is supported by the fitted PDF shown in Figure 3. In addition, due to the substantial increase in the p-value of the K-S test, we note that the fittings of the Lindley distribution for data Sets 2, 3 and 4 are significantly improved when the PITS estimator is used to estimate parameter θ, as compared to the ML estimator. This is due to the presence of outliers far beyond the rest of the data, as shown in Figure 3b-d. However, the OLS and WLS estimators are also found to be quite reliable for estimating parameter θ despite the presence of outliers in the data.
Since the reliability measures based on the Lindley distribution depend on parameter θ, it is important to employ the most suitable method for estimating parameter θ. Based on its application to real data sets, we have shown that the proposed PITS estimator is a viable alternative for estimating the parameter of the Lindley distribution especially when outliers are present in the data.
beyond the rest of the data, as shown in Figure 3b-d. However, the OLS and WLS estimators are also found to be quite reliable for estimating parameter θ despite the presence of outliers in the data.
Since the reliability measures based on the Lindley distribution depend on parameter θ, it is important to employ the most suitable method for estimating parameter θ. Based on its application to real data sets, we have shown that the proposed PITS estimator is a viable alternative for estimating the parameter of the Lindley distribution especially when outliers are present in the data.

Conclusions
In this paper, we proposed a new robust and efficient estimator of the parameter of the Lindley distribution based on the PITS. The advantage of the PITS estimator is that it is conceptually simple and easy to compute. An assessment of the robustness of the PITS estimator based on the BP and IF revealed that the PITS estimator has a high BP and bounded IF, which means that this estimator offers reasonable protection against outliers. In a simulation study, we compared the performance of the PITS estimator with those of several well-known estimators, namely ML, OLS and WLS. The results of the simulation indicated that the performance of the PITS estimator was similar to that of the ML estimator in the absence of outliers and outperforms all the other methods in the presence of outliers. We also note that the OLS and WLS estimators are quite robust and outperformed the ML estimator in the presence of outliers. Four real data sets were applied for which the parameter of the Lindley distribution was estimated. The results demonstrated that the PITS estimator provides a better fitting of this model than the other methods in terms of smallest K-S statistic and highest p-value of the K-S test. Finally, all the abbreviations are listed in Abbreviations part and the R commands for the PITS estimator are available in Appendix B.