EM Algorithm for Mixture Distributions Model with Type-I Hybrid Censoring Scheme

: An expectation–maximization (EM) likelihood estimation procedure is proposed to obtain the maximum likelihood estimates of the parameters in a mixture distributions model based on type-I hybrid censored samples when the mixture proportions are unknown. Three bootstrap methods are applied to construct the conﬁdence intervals of the model parameters. Monte Carlo simulations are conducted to evaluate the performance of the proposed methods. Simulation results show that the proposed methods can perform well to obtain reliable point and interval estimation results. Three examples are used to illustrate the applications of the proposed methods.


Introduction
The finite mixture probability models have been applied to model heterogeneity in a dataset. For example, Peel and McLachlan [1] proposed a mixture of t distributions with the application of clustering multivariate data that contains atypical observations. For more information, the reader may refer to the book by McLachlan and Peel [2]. Reliability inference is one of the important issues in industrial management to maintain the good quality of productions. In many practical cases, source items are coming from different suppliers, which have been certified, to support the company's production line. Hence, it is imperative to investigate the mixture type probability model for life testing to deal with the heterogeneity from different resources in industrial reliability analysis. For example, Sultan et al. [3] explored the mixture of two inverse Weibull distributions and provided some important properties, and Razali and Salih [4] provided a mixture of two Weibull distributions to analyze the lifetimes of electronic components. Ruhi et al. [5] investigated the mixture probability model using Weibull distributions as components and provided a case study.
Recent technology advancement has successfully enhanced product reliability. Spending a long time collecting complete lifetime samples from life tests to conduct reliability inference becomes unrealistic under the considerations of the scheduled test time and budget. Different censoring schemes have been adopted in the literature for life testing to conduct reliability inferences; see Balakrishnan and Aggarwala [6] and Balakrishnan and Cramer [7] for comprehensive information. Among all censoring schemes developed in industry life testing, the type-I and type-II censoring schemes have earned more attention because of easy implementation. Place n items on a life test at the same initial time and denote all failure times by X 1 , X 2 , · · · , X n , respectively. Without loss of generality, let the initial time be 0. The life test continues up to the scheduled time τ 1 for the type-I censoring scheme; while the life test under the type-II censoring scheme is carried out until a preassigned number, r(≤ n), of failures is observed. Let X (1) < X (2) < · · · < X (n) be Moreover, let τ denote the censored time of the TIHC scheme; it can be shown that: Hence, the TIHC sample with the scheduled r and τ can be denoted by X = {X 1 , X 2 , · · · , X D , rep(τ, n − D)}, where rep(τ, n − D) indicates n − D items withdrawn at time τ. The EM algorithm is simple to use, and hence the EM algorithm has been widely applied to obtain the MLEs of the model parameters for censoring data and missing data. When the mixture proportions and the data are type-I hybrid censored, it is not easy to obtain the MLE of the model parameters. If the mixture proportions are known, the EM algorithm surely has a better performance than the cases for which the mixture proportions are unknown. The rest of this article will be structured as follows. Section 2 addresses the likelihood function using the TIHC sample from the finite mixture distribution. The EM algorithm to find the MLEs of the model parameters and survival function, as well as three bootstrap procedures for confidence intervals, are presented in Section 2. In Section 3, a simulation study will be conducted under the mixture model of two-component distributions with Weibull, generalized exponential (GE) and GR distributions. The simulation results will be discussed in Section 3. Moreover, one numerical TIHC sample is used for illustration. Section 4 addresses applications by using three examples and Section 5 provides concluding remarks.
In practical applications, the lifetime quality of items could not be consistent. For example, if the items are provided by multiple suppliers and the lifetime quality of items from those suppliers are inconsistent. In the aforementioned situation, a mixture distributions model can be a good candidate for characterizing the lifetimes of items.

Likelihood Function
Let the lifetimes be taken from a mixture distributions model with k groups and the PDFs and CDFs of the k groups be g j (x|θ j ) and G j (x|θ j ) for j = 1, 2, · · · , k. The PDF and CDF of the mixture distributions model can be represented by: and respectively, where 0 ≤ δ j ≤ 1 is the proportion parameter to indicate the likelihood of the random variable x from the jth distribution, ∑ k j=1 δ j = 1 and Θ = (θ 1 , θ 2 , · · · , θ k , δ 1 , δ 2 , · · · , δ k ). Replacing f (x|Θ) and F(x|Θ) with Equations (3) and (4), Equations (1) and (2) can be respectively represented by: and The normal equation is obtained as ∇ (Θ) ≡ 0, where ∇ (Θ) is a gradient of (Θ) with respect to Θ, and 0 is the zero vector with the same dimension of ∇ (Θ). The MLEΘ of Θ is the solution of ∇ (Θ) ≡ 0. Because of the complicity of the normal equation, the solution procedure is not tractable. The following EM algorithm is suggested to find the MLEΘ of Θ.

EM Algorithm
For i = 1, 2, · · · , n, let i be the indicator of the lifetime of the ith tested item being censored or not and defined as: If X i is taken from the jth distribution, then the contribution to the joint PDF can be presented by: Define latent variable by y ij = 1 if x i is from the jth distribution; otherwise, y ij = 0 for i = 1, 2, · · · , n and j = 1, 2, · · · , n. Then, the likelihood function of Equation (5) and log-likelihood function of Equation (6) can be respectively represented as follows: and The following two steps of the EM algorithm can be iteratively used to obtain the MLEs of the model parameters. E-step: The conditional expected value of (Θ)|x) is calculated and labeled by: (9) where η ij = E[y ij |x] indicates the expected value of the ith sample observation from the jth distribution under right censored at τ, η = {η ij , i = 1, 2, · · · , n, j = 1, 2, · · · , k} and δ = {δ j , j = 1, 2, · · · , k}. During the E-step, for each i in {1, 2, · · · , n}, η ij can be derived by utilizing the conditional probability concept and the result is given as:

Applications
Three common used lifetime distributions of Weibull, GE and GR are used for illustration. Without loss of generality, it is assumed that the mixture distributions model has two components. For simplification, let H ≡ H(Θ, η, δ) here and after.

The Mixture Weibull Distributions
Consider two Weibull distributions as the members for mixture, the PDF and CDF are respectively given by: where α j is shape parameter, λ j is scale parameter and j = 1, 2. Hence, and The MLEs of α j and λ j , j = 1, 2 are the solutions of the system by letting the Equations (16) and (17) be 0.

The Mixture GE Distributions
Consider two GE distributions for mixture, the PDF and CDF are respectively given by where α j is shape parameter, λ j is scale parameter and j = 1, 2. Hence, and The MLEs of α j and λ j , j = 1, 2 are the solutions of the system by letting the Equations (18) and (19) to 0.

The Mixture GR Distributions
Consider two GR distributions for mixture, the PDF and CDF are respectively given by where α j is shape parameter, λ j is scale parameter and j = 1, 2. Hence, and The MLEs of α j and λ j , j = 1, 2 are the solutions of the system by letting the Equations (20) and (21) be 0.
Analogously, the results of mixture Weibull, mixture GE and mixture GR distributions can be extended to the mixture distributions with two different distributions; that is, the mixture Weibull and GE distributions, the mixture Weibull and GR distribution and the mixture GR and GE distributions. We skip the details here to save pages.

Maximum Likelihood Estimate of Survival Function
From Equation (4), the survival function at given time x 0 > 0 is defined by LetΘ denote the MLE of Θ by using the proposed EM algorithm. Then, the MLE of S(x 0 ) can be obtained by the plug-in method and be labeled by: It is difficult to obtain the confidence intervals of Θ and S(x 0 ) based on the sampling distribution ofΘ andŜ(x 0 ). To overcome this difficulty, the parametric bootstrap percentile (PBP) procedure and two bootstrap correction methods are proposed to obtain approximate confidence intervals of Θ and S(x 0 ), respectively.
(a) PBP procedure Step 1.Under the TIHC scheme with r and τ 1 , a TIHC sample, , rep(τ, n − D)}, is collected from a mixture distribution; Step 2.The MLEsΘ based on TIHC sample from Step 1 is obtained by utilizing the proposed EM algorithm and the MLEŜ(x 0 ) is obtained via the plug-in method; Step 3.A bootstrap TIHC sample with r and τ 1 is drawn from the same mixture distribution with the parameters Θ substituted byΘ. Let the bootstrap TIHC sample be denoted by Label all MLEs byΘ * j andŜ * j (x 0 ) for j = 1, 2, · · · , B; Step 6.Letω j , j = 1, 2, · · · , B be a B bootstrap MLEs obtained in Step 5 for a parameter ω considered. The empirical distribution,Ĝ B ω , of {ω j , j = 1, 2, · · · , B} can be obtained. Given 0 < γ < 1, the γth empirical quantile, (Ĝ B ω ) −1 (γ), is defined as the γB th order statistic,ω * ( γB ) , ofω j , j = 1, 2, · · · , B, where γB is the greatest positive integer less than or equal to γB. Then an 1 − γ bootstrap percentile confidence interval of ω can be obtained by More information regarding PBP procedure, readers may refer to the books by Efron [16], Efron and Tibshirani [17] and Shao and Tu [18]. The two bootstrap correction methods are presented as follows: (b) Hybrid bootstrap percentile (HBP) procedure Shao and Tu [18] proposed the HBP procedure is provided as follows: Step 1.Repeat Step 1 to Step 5 of the PBP procedure. Letω be the MLE utilizing the original TIHC sample. Moreover, letω * j , j = 1, 2, · · · , B be the B bootstrap MLEs of ω; Step 2.LetĤ B be the empirical distribution ofφ j , which is defined byφ Step 3.LetĤ −1 B (γ) be the γth quantile ofĤ B , where 0 < γ < 1. An 1 − γ HBP confidence interval for ω can be obtained as the following interval: The second bootstrap correction method that can be found in the books by Efron and Tibshirani [17], as well as in Shao and Tu [18], is addressed as follows: (c) Bootstrap bias-corrected percentile (BCP) procedure Step 1.Implement the Step 1 of the HBP procedure. Let the empirical distribution, G B ω , be constructed based on the bootstrap MLEsω * j , j = 1, 2, · · · , B; Step 2.An 1 − γ bootstrap BCP confidence interval for parameter ω can be obtained by:

Simulation Study
To evaluate the accuracy of the proposed EM algorithm to obtain the MLEs of model parameters using the censored sample under the TIHC scheme with r and τ 1 , an intensive simulation study will be conducted in this section. The simulation study will utilize three distributions that include Weibull, GE and GR distributions as the baseline component distributions for the two-component mixture distributions model; that is, k = 2. All computation work is obtained via using R codes. R is an open software and can be freely downloaded. Population parameters (α 1 , λ 1 , α 2 , λ 2 ) = (2, 1.5, 4, 2), mixture proportions (δ 1 , δ 2 ) = (0.6, 0.4) and (0.75, 0.25), as well as the pre-required failure numbers r = 0.6n and 0.75n, are considered in this simulation study for the sample sizes of n = 500 and 1000. τ 1 is taken as the maximum of the 85th percentiles of the two distributions for mixture; that is, For each set of simulation inputs, the simulation study was conducted 1000 iterations. The maximal iteration for the convergence of EM algorithm is 50 and B = 10,000 is used for bootstrap sampling. To evaluate the quality of the MLEω of parameter ω, relative bias rBias = Bias ω and relative squared root of mean square error rsMSE = sMSE ω are used, where: The simulation results for the MLEs of the model parameters are displaced in Tables 1-3.    To evaluate the interval estimator, the proposed bootstrap procedures of PBP, HBP and BCP described in Section 2.4 with the B = 10, 000 bootstrap sample are used. For each simulation run, the empirical distributions of the MLEs of the distribution parameters and survival functions are constructed based on 10,000 bootstrap samples. Moreover, the confidence intervals of the distribution parameters and survival function are obtained based on the bootstrap methods of PBP, HBP and BCP. The quality of the PBP, HBP and BCP bootstrap methods is evaluated based on the coverage probability (CP) with 1000 confidence intervals. Part of the simulated CPs for 90% confidence intervals of survival function values based on sample size n = 1000 and r = 0.8n are given as follows: 1. For the mixture Weibull distributions model with λ 1 = 1.5, α 1 = 2, λ 2 = 2, α 2 = 4, δ = 0.6 and S(x 0 ) = 0.5 at x 0 = 1.098, the CP = 0.894, 0.895 and 0.897 by using PBP, HBP and BCP procedures, respectively. Figure 1 shows that the lengths of the obtained PBP confidence interval is more consistent than that of the obtained HBP and BCP confidence intervals. The meam length of 1000 confidence intervals based on the PBP, HBP and BCP methods are 0.0256, 0.0286 and 0.0218, respectively. 2. For the mixture GE distributions model with λ 1 = 1.5, α 1 = 2, λ 2 = 2, α 2 = 4, δ = 0. 6 and S(x 0 ) = 0.5 at x 0 = 0.862, the CP = 0.898, 0.897 and 0.899 by using PBP, HBP and BCP procedures, respectively. Figure 2 shows that the lengths of the obtained PBP confidence interval are more consistent than those of the obtained HBP and BCP confidence intervals. The meam lengths of 1000 confidence intervals based on the PBP, HBP and BCP methods are 0.0425, 0.0449 and 0.0358, respectively. 3. For the mixture GR distributions model with λ 1 = 1.5, α 1 = 2, λ 2 = 2, α 2 = 4, δ = 0. 6 and S(x 0 ) = 0.5 at x 0 = 0.708, the CP = 0.894, 0.894 and 0.896 by using PBP, HBP and BCP procedures, respectively. Figure 3 shows that the lengths of the obtained PBP confidence interval is more consistent than that of the obtained HBP and BCP confidence intervals. The meam lengths of 1000 confidence intervals based on the PBP, HBP and BCP methods are 0.0424, 0.045 and 0.0357, respectively.
Based on the simulation results, the PBP, HBP and BCP methods are competitive, and the PBP method slightly outperforms the other two due to the length of the PBP confidence interval is shorter than that of the other two bootstrap methods.  The simulation study is also conducted to evaluate the PDF and CDF for the three mixture distributions models, which include the mixture model with GR and GE distribution, the mixture model with GR and Weibull distributions, and the mixture model with Weibull and GE distributions. Additional simulation parameter inputs include population parameters (α 1 , λ 1 , α 2 , λ 2 ) = (2, 1, 5, 1) and mixture proportion, (δ 1 , δ 2 ) = (0.2, 0.8), (0.50, 0.50). We skip the outputs to save pages. Generally, all the empirical PDFs and CDFs are close to their true function curves.
A simulated example is used to show the applications of the proposed EM-MLE method. The two datasets in Table 4 present the strength of carbon fibers measured in GPa for single carbon fibers and impregnated 1000-carbon fiber tows; see Murthy et al. [19] Kundu and Raqab [20]. These two datasets were originally reported by Badar and Pries [21] and later studied by Kundu and Raqab [20]. The first and second datasets contain 63 and 69 strength measurements, respectively, which are tested under tension at gauge lengths of 10 mm and 20 mm. The histograms and empirical density plots are given in Figure 4. The Kolmogorov-Smirnov test is used to check the goodness-of-fit of using the Weibull distribution to model these two datasets. The first dataset provides the MLEsλ 1 = 2.651 andα 1 = 5.502, and the second dataset provides the MLEsλ 2 = 3.31 andα 1 = 5.048. Replacing the parameters with their MLEs, we obtain the Kolmogorov-Smirnov distances D = 0.067 with p-value = 0.965 for the first dataset and D = 0.09 with p-value = 0.775 for the second dataset. The testing results indicate that the Weibull distribution can characterize these two datasets well.  Table 4. These two datasets are not TIHC samples. Assuming a future pool of these two datasets could be coming from a mixture distributions model with Weibull(α 1 = 5.502, λ 1 = 2.651) and Weibull(α 2 = 5.048, λ 2 = 3.315) and δ = 0.6. We regenerate a TIHC sample with r = 0.8n and n = 500, τ 1 = 1.95, and x 0 = 3.76. We are interested in estimating the parameters of the mixture model of two Weibull distributions and the survival function at 0.5(or at x 0 = 3.76). Utilizing the proposed EM algorithm for the maximum likelihood estimation method and the regenerated sample. The generated TIHC sample is given in Table 5. The obtained MLEs areλ 1 = 2.672,α 1 = 5.435,λ 2 = 3.306,α 2 = 5.261,δ = 0.593 and the estimated survival function is 0.498. Using the three proposed bootstrap methods with 5000 iterations to construct the bootstrap empirical distribution of the survival, we obtain the confidence intervals of (0.482, 0.521) for the parametric PBP procedure, (0.474, 0.512) for the HBC procedure and (0.473, 0.513) for the bootstrap BCP procedure.

Illustrative Examples
Three practical datasets will be used for the further investigation on the mixture modeling and the illustration of the proposed EM-MLE algorithm. To save space, only the first example dataset is reported in Table 6. Table 6 has 20 electronic failure times that were used by Razali and Salih [4] for the mixture distribution based on two Weibull distributions where one has a positive support with a positive left boundary and the other has the whole positive support. The second dataset regarding 153 aircraft windshield failure times in terms of 1000 h was originally given in Blischke and Murthy [19]. This dataset was also used by Ruhi et al. [5] for the twofold Weibull mixture model. The third dataset contains 600 lifetimes of adapters which are used in the liquid-crystal display (LCD) assembling monitors. Tsai et al. [22] has studied this dataset and their suggested model of two Weibull mixture distributions can characterize this dataset well. Table 6 has 20 electronic failure times that were used by Razali and Salih [4] for the mixture distribution based on two Weibull distributions, where one has a positive support with a positive left boundary and the other has the whole positive support. In this example, the mixture model of Weibull and GE distributions is used to model TIHC samples that can be generated from 20 electronic failure times with r = 0.5n, 0.6n, 0.75n, 0.8n, 0.9n, n with n = 20, and a large value of τ 1 is used to allow all observations are not censored. Given a TIHC sample generated by the TIHC with each given r and τ 1 , the MLEs of the mixture distribution parameters were obtained through the proposed EM algorithm. Then the corresponding six mixture distribution PDFs using six different sets of MLEs as the parameter inputs were obtained and drawn in Figure 5.

Example 1
All six MLEs of mixture PDFs are very close except at x 0 = 2, where the estimated PDFs using smaller r values have higher peak. Moreover, the 95% confidence intervals of S(x 0 ) are reported in Table 7 based on the PBP, HBP and BCP procedures, where the LB and UB denote the lower and upper bounds of a bootstrap confidence interval, respectively. All PBP, HBP and BCP confidence intervals at each line of given r are closed. Moreover, two boundaries from the PBP and BCP procedures are closed for each given value of r among three procedures.

Example 2
In this example, aircraft windshield failure times will be used to check the mixture model through the three proposed distributions. Similar to the procedure used by Ruhi et al. [5], we separated 153 observations into two groups; one group has 88 observations and the other group has 65 observations. Ruhi et al. [5] checked Weibull distribution as good-of-fit for both groups. The Kolmogorov-Smirnov test is used to further check the goodness-of-fit of using GE distribution to model these two groups. The first group provides the MLEsλ 1 = 0.7593 andα 1 = 3.648, and the second group provides the MLEŝ λ 2 = 0.7024 andα 2 = 1.9458. Replacing the parameters by their MLEs, we obtain the Kolmogorov-Smirnov distances D = 0.11287 with p-value = 0.2122 for the first group and D = 0.14015 with p-value = 0.1413 for the second group. The testing results indicate that the GE distribution can characterize these two groups well.
In this study, the mixture model with Weibull and GE distributions is suggested to model TIHC samples that were generated from these 153 observations with r = 0.5 n, 0.6 n, 0.75 n, 0.8 n, 0.9 n, n, n = 153 and τ 1 = 2. Given a TIHC sample generated by the TIHC with each given r and τ 1 , the MLEs of mixture distribution parameters were obtained through the proposed EM algorithm. Then the corresponding six mixture distribution PDFs using the calculated six different sets of MLEs as the parameter inputs were obtained.
The 95% confidence intervals of S(t 0 ) are reported in Table 8 based on the PBP, HBP and BCP procedures, where the LB and UB denote the lower and upper bounds of a bootstrap confidence interval, respectively. Table 8 shows that all estimated values for r ≥ n 2 are very closed to the corresponding estimated results by using all 153 observations. Moreover, the 95% confidence intervals by BP and BC are very closed and have shorter lengths than the 95% confidence intervals by HB.

Example 3
This example was based on a past statistical consultation project for a LCD manufacturer in Asia. The manufacturer's information is not disclosed here due to confidentiality. We refer to this manufacturer as the company. The dataset is regarding the failure time of a type of video graphics array (VGA) adapter used for the liquid-crystal display (LCD) assembling monitors. The work done by Tsai et al. [22] was suggested and the company continued to purchase adapters from two suppliers and started to record the exact failure time in terms of hours and the originality of each adapter. In this example, the dataset that contains a sample of n = 600 failure times in terms of days, with their originality, are considered. Therefore, all 600 failure times can be separated into two groups. The first group has 266 observations and the second group has 334 observations. The Kolmogorov-Smirnov test is used to check the goodness-of-fit of using the Weibull distribution to model these two groups. The first group provides the MLEsλ 1 = 58.18 andα 1 = 2.027, and the second group provides the MLEsλ 2 = 39.62 andα 2 = 4.145. Replacing the parameters with their MLEs, we obtain the Kolmogorov-Smirnov distances D = 0.0000+ with p-value > 0.999 for the first group and D = 0.0000+ with p-value > 0.999 for the second group. The testing results indicate that the Weibull distribution can characterize the dataset of two groups well.
In this study, the mixture model with Weibull and Weibull distributions is suggested to model TIHC samples, which were generated based on these 600 observations. The TIHC scheme is r = 0.8n = 480 and τ 1 = 80. Based on the generated sample, the MLEs of mixture distribution parameters were obtained byλ 1 = 58.479,α 1 = 2.01,λ 2 = 39.396 and α 1 = 3.978. The 95% confidence intervals of S(t 0 ) = 0.5 were also obtained and reported in Table 9 based on the PBP, HBP and BCP procedures. Table 9 shows that all 95% bootstrap confidence intervals cover the true survival function.

Conclusions and Remarks
In the real world, source items could come from different suppliers to support a company's production line. The heterogeneity from different resources may cause inconsis-tencies in the production quality. In these circumstances, it would be more reasonable to utilize a mixture distributions model to characterize the product's lifetime quality. Moreover, it is very difficult to collect the entire lifetime sample from the life test because of the finite test schedule and restricted budget. Hence, the TIHC scheme with the given failure number r and life test schedule τ 1 has been proposed for collecting lifetime data from numerous practitioners. In this work, the EM algorithm has been proposed to find the MLEs of the mixture distribution model parameters and the survival function. Three different bootstrap procedures-the PBP, HBC and BCP procedures-are provided to establish the empirical distribution of the MLE of the mixture distribution model parameters and the survival function based on TIHC samples, and to obtain the confidence intervals. The parameter estimation for the mixture distributions model requires a large sample; even the mixture proportions are known with a complete sample. The TIHC scheme renders the sample censored and incomplete. Hence, the proposed EM-MLE requests a very large sample in order to obtain good quality MLEs for the model parameters. This is a drawback of the proposed EM-MLE method, and how to reduce the sample size can be left for future study.
All existing studies on the mixture distributions model have been focused on the mixture of common distributions but with different distribution parameters. In the current work, we extended the EM-MLE method to the mixture distribution model with different distributions as members. The approximate MLE approach proposed by Mendenhall and Hader [8] is a very interesting procedure for mixture distribution modeling. This will be the focus of a future research project. Meanwhile, the Bayesian estimation approach and some other censoring schemes might provide novel future research potential to add to the mixture distributions model. These two topics will be studied in the future.