Bias-Corrected Maximum Likelihood Estimation and Bayesian Inference for the Process Performance Index Using Inverse Gaussian Distribution

: In this study, the estimation methods of bias-corrected maximum likelihood (BCML), bootstrap BCML (B-BCML) and Bayesian using Jeffrey’s prior distribution were proposed for the inverse Gaussian distribution with small sample cases to obtain the ML and Bayes estimators of the model parameters and the process performance index based on the lower speciﬁcation process performance index. Moreover, an approximate conﬁdence interval and the highest posterior density interval of the process performance index were established via the delta and Bayesian inference methods, respectively. To overcome the computational difﬁculty of sampling from the posterior distribution in Bayesian inference, the Markov chain Monte Carlo approach was used to implement the proposed Bayesian inference procedures. Monte Carlo simulations were conducted to evaluate the performance of the proposed BCML, B-BCML and Bayesian estimation methods. An example of the active repair times for an airborne communication transceiver is used for illustration.


Introduction
Process capability analysis has been widely used to identify how well the outputs from an in-control process meet the requirements, specifications and expectations of customers. In practice, process capability analysis methods aim to continuously monitor the process quality via utilizing capability indices to assure if the quality of products is capable of meeting the specifications and supply information based on the product design and process quality improvement. The results of process capability analysis can be the basis of cost reduction, which is attributable to the decrease in product failures; see Pan and Wu [1]. Kane [2] presented the relations of the process capability indices, C p , C pu , C pl and C pk , which are, respectively, defined as Stats 2022, 5 where µ is the population mean, σ is the population standard deviation, L is the lower specification limit and U is the upper specification limit. Kane [2] also indicated that these indices can be a complementary system of measures for evaluating the process performance. Kocherlakota [3] and Kotz and Johnson [4] studied the statistical theories for various process capability indices. Comprehensive discussions about using process capability analysis methods for quality control can be found in the works of Rodriguez [5], Palmer and Tsui [6] and Montgomery [7]. The process performance index is one of the widely used process capability indices for evaluating the quality of lifetime data. Let X denote the product lifetime random variable and L denote a specified threshold of a lower bound. The conforming rate can be defined as δ = P(X > L). The process performance index can have a close connection with the conforming rate of δ. The statistical properties regarding using the process capability indices of C p , C pu , C pl and C pk under the normality assumption have been thoroughly studied in the literature. Hong et al. [8] presented analytical procedures for the ML estimation and hypothesis testing to evaluate the process performance index when the lifetimes of products follow the Pareto distribution. Ahmadi et al. [9] used generalized order statistics to conduct inferential procedures for the process performance index under an exponential distribution. Lee et al. [10] proposed optimal inferential procedures to evaluate the process performance index based on progressively type-II censored samples that were taken from the Burr type XII distribution. Lee et al. [11] proposed an inferential procedure by using uniformly minimum variance unbiased estimation and a hypothesis-testing method to evaluate the process performance index based on type-II censored samples that were taken from the two-parameter exponential distribution. Lee et al. [12] proposed Bayesian inference procedures to assess the process performance index based on progressively type-II censored samples under the Rayleigh distribution. Ahmadi et al. [13] proposed ML estimation methods to estimate the process performance index based on progressively first-failure censored samples when the lifetime data follow a Weibull distribution. Wu and Chiu [14] obtained fourteen different estimates of the process performance index for the two-parameter exponential distribution under a multiple type-II censoring scheme. After a simulation study of a performance comparison, three estimates among them are recommended by Wu and Chiu [14] to develop hypothesis-testing procedures for the process performance index. Wu and Lin [15] proposed a statistical inference for the process performance index based on type-II exponentially distributed samples. Montgomery [7] recommended using the process performance index to evaluate the quality of products. Zhu et al. [16] proposed an inferential procedure by using the ML estimation method to evaluate the process performance index based on power-normal distribution samples. They also discussed the drawbacks of using the exact Fisher information matrix with a delta method under the power-normal distribution to obtain an approximate confidence interval (ACI) of the model parameters. All aforementioned works on process performance index are summarized in Table 1 for easy reference. Table 1. Summary of estimation methods for process performance index.

Method Sampling Distribution
Hong et al. [8] MLE Type-II Pareto Ahmadi et al. [9] Statistical inference Generalized order statistics exponential Lee et al. [10] Optimal inferences Progressively type-II Burr type XII Lee et al. [11] UMVUE Type-II Two-parameter exponential Lee et al. [12] Bayesian inference Progressively type-II Rayleigh Ahmadi et al. [13] MLE Progressively first-failure Weibull Wu and Chiu [14] Fourteen estimates A multiple type-II Two-parameter exponential Wu and Lin [15] Statistical inference Type-II Exponential Zhu et al. [16] MLE Random sample Power-normal Among widely used lifetime distributions, the inverse Gaussian (IG) distribution, also known as the Wald distribution, has been renowned as a versatile lifetime model with sound physical interpretation. The ML and Bayesian estimation methods were commonly used to estimate the IG distribution parameters. Banerjee and Bhattacharyya [17] investigated the Bayesian inferential approach for estimating the IG distribution parameters for the application of equipment failure data. Amry [18] studied the Bayesian inference of the IG parameters using Jeffrey's prior under a quadratic loss function. More information about using the IG distribution for engineering applications can also be found in the books by Chhikara and Folks [19] and Johnson et al. [20]. Sun and Ye [21] discussed the frequentist validity of posterior quantiles for a two-parameter exponential family that includes the IG distribution as a member. Rostamian and Nematollahi [22] studied the stress-strength reliability using the ML estimation method via using an expectationmaximization algorithm and the Bayesian estimation method based on progressively type-II censored IG distribution samples. A survival analysis of the IG distribution based on using Bayesian and Fiducial approaches has been studied by Jayalath and Chhikara [23]. Bera and Jana [24] developed a bootstrap interval of stress-strength reliability assuming that the stress and strength variables are IG distributed.
Sundaraiyer [25] proposed inferential procedures to obtain the ML estimator and bootstrap ACI of the process capability index proposed by Clements [26] when the quality variable follows the IG distribution. We use the term of MLE to denote the ML estimator and ML estimate here and after. Investigating how to reduce the estimation bias for the process performance index based on IG distributed samples and obtain a reliable ACI for the process performance index is helpful for the applications of quality control. Balay [27] used BCMLEs of the generalized inverse Lindley distribution parameters to compute the generalized process capability index C pyk . The C pyk was firstly proposed and studied by Maiti et al. [28]. In addition, Balay [27] obtained a bootstrap ACI of C pyk . Considering the merit of IG distribution being a versatile lifetime model with sound physical interpretation, the IG distribution can be an alternative to the generalized inverse Lindley distribution for reliability analysis applications. The purposes of this article can be twofold. Firstly, we would like to propose analytical procedures to obtain the BCMLEs, whose bias is O(n −2 ) based on the bias correction method proposed by Cordeiro and Klein [29], B-BCMLEs and Bayes estimators via using the Jeffery's prior distribution for the IG distribution parameters and the investigated process performance index. Secondly, we would like to establish the procedures taken to obtain an ACI and the highest posterior density interval (HPDI) for the one-sided version of C pyk . The HPDI is the interval with the shortest length on the posterior density at the given confidence level. We use the term BE to denote the Bayes estimator and Bayes estimate here and after. Because the posterior distribution in the Bayesian estimation procedure is complicated, the Markov chain Monte Carlo (MCMC) approach is used to overcome the computational difficulty when generating random samples from the posterior distribution. To our knowledge, these two aforementioned purposes have not yet been studied in the literature.
The rest of this paper is organized as follows. We address the process capability indices and define the process performance index based on the one-sided version of C pyk in Section 2. In Section 3, we derive the inferential procedures to obtain the CK-BCMLE and B-BCMLE of the model parameters through using the bias correction method proposed by Cordeiro and Klein [29] and the bootstrap method, respectively. A bootstrap algorithm is suggested to obtain the B-BCMLEs of the model parameters and the bootstrap ACI of the process performance index. Moreover, the Bayesian estimation procedure is developed via using the Jeffery's prior distribution, and an MCMC hybrid algorithm of Gibbs sampling and the Metropolis-Hastings algorithm is provided to obtain the BEs of the IG distribution parameters and process performance index. Monte Carlo simulations are conducted in Section 4 to evaluate the performance of the proposed estimation methods. A real data set with 46 active repair times (in hours) for an airborne communication transceiver is given in Section 5 for illustration. Concluding remarks are given in Section 6.

The Generalized Process Performance Index
Based on C p and C pk , two generalized process capability indices can be defined as follows: where ξ = (µ − T)/σ and T is the process target. When the quality measurements of products follow a normal distribution, C p , C pk , C pm and C pmk have been the four most widely used process capability indices in practical applications. However, non-normally distributed quality measurements, which have a skewed distribution, can be found in many works by Clements [26], Gunter [30], Constable and Hobbs [31], Mukherjee and Singh [32], Tang et al. [33] and Chen et al. [34]. Among the aforementioned studies, Clements [26] proposed two generalized versions of process capability indices that are defined by where X γ is the γth quantile of the quality characteristic measure, X. When X follows a normal distribution, C p(q) and C pk(q) reduce to C p and C pk , respectively. C p(q) and C pk(q) are the two most popular process capability indices used to determine the quality of products if the distribution of the quality variable is not normal. Maiti et al. [28] proposed a new generalized version of the process capability index, where F(x) = P(X ≤ x) is the cumulative distribution function (CDF). α 1 and α 2 are the specified lower and upper tailed probabilities of F, respectively. For lifetime products, we are more concerned with if the lifetime is higher than the lower specification limit. Hence, the one-sided version of C pyk was used in this study for the process performance index. We denote the one-sided version of C pyk by C L in this study for simplification. The C L is defined by

The Inverse Gaussian Distribution and Estimation Methods
In this section, the IG distribution is addressed in Section 3.1. Moreover, we propose the analytical procedure in Section 3.2 to obtain the BCMLEs of the model parameters based on the bias correction method proposed by Cordeiro and Klein [29]. This bias correction method is named the CK-BCML method here and after. The bootstrap procedure used to obtain B-BCMLEs of the IG distribution parameters is proposed in Section 3.3. Moreover, the delta method based on the Fisher information matrix is used to obtain an ACI of C L in Section 3.4. In Section 3.5, we use Jeffrey's prior distribution and the proposed MCMC algorithm to develop the Bayesian estimation procedures for obtaining the BEs of the IG distribution parameters and C L . Moreover, the HDPI of C L is also obtained through the corresponding MCMC chain of C L .

Maximum Likelihood Estimation
The probability density function (PDF) and CDF of the IG distribution are defined by and respectively, where Φ(·) is the CDF of the standard normal distribution, Θ = (θ 1 , θ 2 ) = (µ, λ), µ is the mean of the IG distribution and λ is a reciprocal measure of dispersion. Symbolically, denote the IG distribution by X ∼ IG(Θ). The moment generating function, variance, skewness and kurtosis can be obtained, respectively, by Moreover, we can obtain the mean and variance for the reciprocal of X by and The IG distribution is a positively skewed distribution due to the skewness being always positive. Johnson et al. [20] indicated that the IG distribution is a member of the exponential family and is unimodal. The mode of the IG distribution is located at Based on the CDF of the IG distribution, the C L can be defined by where α 1 = 0.0027 can be a reference number for practical applications.
Let X = (X 1 , X 2 , · · · , X n ) denote a random sample taken from the IG distribution. The log-likelihood function based on the realization x = (x 1 , x 2 , · · · , x n ) of X can be expressed by The first derivatives of with respect to µ and λ can be obtained, respectively, by Let 1 = 0 and 2 = 0; we can express the MLEs of µ and λ bỹ The asymptotic distribution ofΘ M = (μ M ,λ M ) can be expressed bỹ where A ≡ A(Θ) = I −1 (Θ) is the inverse of the Fisher information matrix given in Appendix B. From Appendix B, and b and From Table A3 in Appendix C, the matrix B can be expressed by Define V A = vec(A) as a vectorization operator that produces a single column vector by stacking all right columns of A below the adjacent left one. Hence, Under regular conditions of partial derivatives of the log-likelihood function, η ij , η ijk and ij , Cordeiro and Klein [29] showed that the bias ofΘ is O(n −2 ) and can be expressed as Denote the obtained BCMLE of Θ byΘ C = (μ C ,λ C ) T ; then,Θ C can be approximated bỹ . We can show that and the CK-BCMLEΘ CK can be obtained bỹ

Bootstrap Bias-Correction Method
The bootstrap methods have been widely used to obtain an ACI for the model parameter. Readers can refer to the book of Efron and Tibshirani [35] to receive more comprehensive information about using bootstrap methods for statistical inference. The bootstrap biascorrection maximum likelihood (B-BCML) method can be implemented through using the following steps:

Initial
Step: Obtain the MLEΘ M of Θ based on the working sample x = (x 1 , x 2 , · · · , x n ).
Step 1: Generate a bootstrap sample x * = (x * 1 , x * 2 , · · · , x * n ) from the IG(Θ M ) and obtain the MLEΘ * M of Θ based on the bootstrap sample x * . Step 2: Repeat Step 1 B times and label the obtained MLEs by {Θ * M,j , j = 1, 2, · · · , B}. The bias ofΘ can be obtained bỹ Step 3: The B-BCML estimate, B-BCMLE, can be obtained bỹ The above bootstrap method is also called the parametric bootstrap method.

The Approximate Confidence Interval of C L
Replacing the unknown Θ in C L byΘ M ,Θ CK andΘ B obtained above, the corresponding MLEs of C L can be presented bỹ respectively. Utilizing the delta method, the asymptotic distribution of the MLE,C L , can be shown as the normal distribution; that is, where ∇C L is the gradient of C L and defined by with δ T = (δ 1 , δ 2 ), and  1,B ,δ 2,B ). The plug-in asymptotic variance estimates ofC L,h can bẽ Therefore, the approximate (1 − γ) × 100% CIs of C L can be and based on the typical ML estimation method, the bias correction method proposed by Cordeiro and Klein [29] and the bootstrap bias-correction method, respectively, where

Bayesian Estimation Method
The Fisher information matrix that was given in Section 3.1 has Hence, the Jeffrey's prior distribution of Θ can be obtained by The Jeffery's prior distribution is a parameter-free non-informative prior distribution. Hence, we can use the Jeffery's prior distribution to obtain the BEs of µ, λ and C L with less subjective assumptions. From the posterior distribution of Θ, given x, can be expressed by Hence, the conditional distribution of µ, given λ and x, can be obtained by and the conditional distribution of λ, given µ and x, can be derived as . It can be noted that the conditional distributions of π(µ|λ, x) and π(λ|µ, x) do not have complete analytic forms. In this study, the MCMC hybrid algorithm of Gibbs sampling and the Metropolis-Hastings algorithm can be used to obtain the BEs of µ and λ.
Step 2-Update λ: For h ≥ 1, generate u ∼ U(0, 1) and Step 3-Update C L : Step 4-Obtain BEs: Repeat Step 1 to Step 3 N times, where N is a large integer. Perform s burn-in operation by removing the leading N 1 (≤ N) Markov chains for each parameter. Considering the square loss function for implementing the Bayesian estimation, the BE of each parameter can be the sample mean of the remaining (N − N 1 ) Markov chains. Denote the BEs of µ, λ and C L byμ BE ,λ BE andĈ L,BE , respectively.

Monte Carlo Simulations
In this section, Monte Carlo simulations were conducted to evaluate the performance of the estimation methods of the typical ML, CK-BCML, B-BCML and Bayesian. Random samples with a size of n = 30 and 50 were generated from IG(µ = 8, λ = 5) and IG(µ = 10, λ = 8), respectively, to obtain the MLEs, CK-BCMLEs, B-BCMLEs and BEs of µ and λ. The MLE, CK-BCMLE, B-BCMLE and BE of the process performance index were also obtained based on the plug-in method for L = 0.5, 0.6, 0.8, and 1. When µ = 8 and λ = 5, the true value of the process performance index can be obtained by C L = 1.0043, 0.9957, 0.9644 and 0.9173 for L = 0.5, 0.6, 0.8, and 1, respectively, and the corresponding parts per million are 2876, 7130, 22,625 and 45,939. When µ = 10 and λ = 8, the true value of the process performance index can be obtained by C L = 1.0098, 1.0089, 1.0033 and 0.9898 for L = 0.5, 0.6, 0.8 and 1, respectively, and the corresponding parts per million are 138, 568, 3389 and 10,068.
Bootstrap repetition, B = 500, was used to implement the parametric bootstrap biascorrection method. For Bayesian estimation, firstly, we generated N = 51,000 Markov iterations to implement the MCMC method, and the leading N 1 = 1000 Markov iterations were removed for the burn-in operation. Secondly, a spacing operation by selecting one of every ten iterations was used to reduce the autocorrelation in each Markov chain. Finally, 5000 Markov chains were used to obtain the BEs of the parameters µ, λ and PPI, respectively.
The measures of relative bias (rbias) and relative root mean squared error (rRMSE) were evaluated using 10,000 iteration runs. Assume that the target parameter is θ and the obtained estimates areθ j , j = 1, 2, · · · , 10,000; the rbias and rsMSE are defined by and rRMSE = 1 θ × ∑ 10,000 The rbias and rRMSE are scale-free measures. All simulation results for estimating model parameters are reported in Tables 2 and 3.  From Tables 2 and 3, we can find that the estimation quality of the MLE ofμ is good. The derivation processes in Section 3.2 show that the CK-BCMLE of µ is the same as the typical MLE; that is,μ M =μ CK . Hence, the values of rbias and rRMSE ofμ M are also close to that of the B-BCMLEμ B . However, we can find that the BCMLEsλ CK andλ B outperform the MLEλ M , with a smaller rbias and rRMSE when the sample size is small. Because the proposed Bayes estimation method is developed with the parameter-free non-informative prior distribution of Jeffery, the performance of BE can be compared with its competitor of MLE. In Tables 2 and 3, we also find that the rbias ofμ BE is small but larger than the the bias ofμ M . The rRMSE ofλ BE is larger but close to the rRMSE ofλ M . Based on the findings in this study, using gradient methods for optimization makes the ML estimation method less reliable for estimating the process performance index. The obtained ACI of the process performance index via using the MLEs with the delta method and exact Fisher information matrix is conservative. We will study the CPs of the ACI and HPDI of the process performance index to show the good performance of the proposed Bayesian estimation method.  To verify the performance of the delta and Bayesian estimation methods for the interval inference of PPI under small sample size cases, we established the 95% ACI and HPDI for the PPI at the lower specification limit L = 0.5, 0.6, 0.8 and 1. Moreover, the mean lower bounds and mean upper bounds of the 95% ACI and HPDI of the PPI and their CPs were evaluated based on 10,000 iterations. We use the terms of mLB and mUB to denote the mean lower bounds and mean upper bounds, respectively. All simulation results are summarized in Tables 4 and 5.
In view of Tables 4 and 5, we find that the condition of a small sample has an impact on the quality of the interval inference for the PPI. When n ≤ 50, the CP of the 95% ACI of C L is below its nominal value. In particular, for the case of small L and n = 30, the CP of the 95% ACI of C L seriously underestimates the nominal value 0.95 for the delta method with the ML, CK-BCML and B-BCML methods. The delta method with the ML method performs worst among all competitors. When the sample size increases, the estimation performance of the delta method with the ML, CK-BCML and B-BCML methods can be significantly improved. The proposed Bayesian estimation method outperforms all of the competitors based on the delta method in terms of the CP when obtaining a HDPI for the process performance index; even the sample size and L are small. In summary, we recommend using the proposed Bayesian estimation method to obtain the HPDI of the process performance index for IG distribution.

An Example
A maintenance data set concerning the active repair times for an airborne communication transceiver is used to illustrate the applications of the proposed estimation methods. This data set contains 46 repair times in hours as follows: 0.  4.5, 4.7, 5.0, 5.4, 5.4, 7.0, 7.5, 8.8, 9.0, 10.3, 22.0, 24.5. This data set of active repair times was initially investigated by Chhikara and Folks [36]. They obtained the MLEs µ M = 3.607 andλ M = 1.659 and used the Kolmogorov-Smirnor test to show that the IG distribution can characterize this data set well. The data set was also analyzed based on Bayesian estimation methods by many studies after Chhikara and Folks [36]; for example, Sinha [37], Betrò and Rotondi [38] and Jayalath and Chhikara [23].
Using the proposed methods in Section 3 for the data set of active repair times with L = 0.2, we can obtain the MLEs, CK-BCMLEs, B-BCMLEs and BEs of µ, λ and C L . The proposed Bayesian estimation method was implemented with N = 50,000 and the first N 1 = 1000 generated estimates were removed for teh burn-in operation. The spacing operation, which selects one of every ten generated estimates, was used for cutting the Markov chain to reduce the first-order autocorrelation in every Markov chain. To check the quality of the MCMC method, the Markov chain plots based on 5000 generated estimates ofμ BE andλ BE are given in Figures A1 and A2, respectively. The first-order autocorrelation coefficients based on the 5000 generated estimates ofμ BE andλ BE are 0.046 and 0.06, respectively. We note that two first-order autocorrelation coefficients are close to 0. These findings indicate that the generated Markov chains are almost independent chains. All of the obtained estimates are reported in Table 6. From Table 6, we can find that the B-BCMLE and BE of µ are slightly larger than the MLE and CK-BCMLE of µ. The MLE and BE of λ are larger than two BCMLEs of λ. The BE of C L is the smallest one among the four obtained estimates. The parts per million based on the MLE, CK-BCMLE, B-BCMLE and BE of C L are 6232, 8160, 7813 and 6071, respectively. Table 7 reports the 95% ACI and HPDI of C L via using the delta method and Bayesian estimation method, respectively. Two proposed BCML methods generate close ACIs. We also find that the lower limit of the ACI based on the ML method is slightly larger than the lower limits of the ACIs based on the two proposed BCML methods and the lower limits of HPDI based on the proposed Bayesian method.

Concluding Remarks
Considering the restriction of the sample resource when evaluating the quality of lifetime products, we used the bias correction estimation method proposed by Cordeiro and Klein [29] and the bootstrap bias correction method to improve the estimation quality of the typical ML estimation method for estimating the process performance index under the IG distribution. We derived the exact forms of CK-BCMLEs and provided an algorithm used to obtain the B-BCMLEs for the IG model parameters and process performance index. The delta method was used to obtain an ACI for the process performance index. Moreover, a Bayesian estimation procedure was proposed to obtain the BEs of the IG model parameters and the process performance index. The HPDI of the process performance index was obtained via using the proposed MCMC hybrid algorithm of the Gibbs sampling and Metropolis-Hastings algorithm.
An intensive simulation study was conducted to compare the performance of the proposed CK-BCML, B-BCML and Bayesian estimation methods with the typical ML estimation method. Simulation results show that the ACIs based on the delta method with the MLEs, CK-BCMLEs and B-BCMLEs performed less satisfactory, with a seriously underestimated CP when the sample size and lower specification limit were small. The proposed Bayesian estimation method can be an alternative method, other than the delta method, used to obtain a reliable HPDI for the process performance index when the sample size is small. Because the Jeffrey's prior distribution is a parameter-free prior, the proposed Bayesian estimation method is less subjective and easy to be implemented. A data set composed of 46 active repair times for an airborne communication transceiver was used to illustrate the applications of the proposed methods.
For saving the testing time and cost, censoring schemes are often adopted in engineering applications to collect censored lifetime data. Extending two proposed BCML estimation methods and the Bayesian estimation procedures for the IG distribution with different censoring schemes is interesting and can be studied in the future.

Data Availability Statement:
The data presented in Section 4 were obtained from Chhikara and Folks [36].
Acknowledgments: This research was funded by the Ministry of Science and Technology, Taiwan, grant number MOST 110-2221-E-032-034-MY2; National Natural Science Foundation of China, grant number 52174060. We thank you for the excellent suggestions from reviewers and editors to improve the quality of paper.

Conflicts of Interest:
The authors declare no conflict of interest. Figure A1. The Markov chain with 5000 iterations ofμ BE after burn-in. Figure A2. The Markov chain with 5000 iterations ofλ BE after burn-in.

Appendix B. Fisher Information Matrix
The second derivatives of with respect to µ and λ can be obtained, respectively, by 11 Let η ij = E( ij ), i, j = 1, 2. It is trivial to show that η 12 = η 21 = 0. Moreover, we can obtain the following results: and The Fisher information matrix can be presented by

Appendix C. For MLE Bias Correction
Using the equations of ij and η ij , i, j = 1, 2, the analytic forms of η ijk and η (k) ij , i, j, k = 1, 2 are obtained and reported in Tables A1 and A2. The entries of the sub-matrices B (1) and B (2) are evaluated and reported in Table A3.