Parameter Estimation for Composite Dynamical Systems Based on Sequential Order Statistics from Burr Type XII Mixture Distribution

: Considering the impact of the heterogeneous conditions of the mixture baseline distribution on the parameter estimation of a composite dynamical system (CDS), we propose an approach to infer the model parameters and baseline survival function of CDS using the maximum likelihood estimation and Bayesian estimation methods. The power-trend hazard rate function and Burr type XII mixture distribution as the baseline distribution are used to characterize the changes of the residual lifetime distribution of surviving components. The Markov chain Monte Carlo approach via using a new Metropolis–Hastings within the Gibbs sampling algorithm is proposed to overcome the computation complexity when obtaining the Bayes estimates of model parameters. A numerical example is generated from the proposed CDS to analyze the proposed procedure. Monte Carlo simulations are conducted to investigate the performance of the proposed methods, and results show that the proposed Bayesian estimation method outperforms the maximum likelihood estimation method to obtain reliable estimates of the model parameters and baseline survival function in terms of the bias and mean square error.


Introduction
The components in a n-component composite system often fail sequentially when the system is on duty continuously until the system is declared as a failure system. An ncomponent composite system is an n-component composite dynamical system (CDS) if the failure of a component induces higher work loading on the surviving components of system. Nowadays, the CDS is one of the most widely used electronic devices. The lifespan of a CDS is usually defined as the r th (r ≤ n)-ordered component failure time due to the system structure designed or efficiency concerned. The CDS is also named an r-out-of-n failure system, which includes a series (r = 1) and parallel (r = n) systems of n components, respectively. Because each component failure can induce higher loading on the surviving components, the risk of the system is increased along with the increase of the number of failure components. One working assumption for characterizing the increased risk of the CDS is to assume that all surviving components in the system can equally share the loading of stress. Therefore, the observed ordered component failure times from the CDS system of equal load-share are called sequential order statistics (SOSs). Thereafter, the CDS will be simply called CDS. All acronyms in this study are given in Table 1. Many researchers paid attention to the study of the reliability of CDSs in the past twenty years. Kamps [1] proposed a unified approach based on the concept of generalized order statistics (GOSs) for different models. When the components follow an exponential distribution, Cramer and Kamps [2] investigated the reliability of the r-out-of-n failure system. Cramer and Kamps [3] used SOS samples to study the methods of maximum likelihood estimation, best linear unbiased estimation, and uniformly minimal variance unbiased estimation. They used the three estimation methods to obtain the estimates of the Exponential distribution parameters. Cramer and Kamps [4] provided the marginal distribution function of SOS and GOSs with no restrictions on the parameters. Their findings can be connected to the revelation transform for the beta distribution and Meijer's G-functions. Revathy and Chandrasekar [5] proposed statistical methods to obtain minimum risk equivalent estimates of the parameters in the exponential distribution, Gamma distribution, and location-scale families. Zhuang and Hu [6] investigated the stochastic properties of SOS and conducted multivariate stochastic comparisons based on SOS samples. Balakrishnan et al. [7] studied the reliability of CDS under the a conditionally proportional hazard rate model using semiparametric methods. They also studied the impact of using different link functions on the maximum likelihood estimates (MLEs) based on SOS samples.
Burkschat [8] extended the SOS concept to study exchangeable random variables under a weaker condition than the independent and identically distributed assumption. Burkschat's idea is able to apply to general coherent systems, such as the coherent systems of which components have dependent and non-identically distributed failure time distributions. Beutner and Kamps [9] used SOS samples to obtain the estimates of scale parameters for different CDSs when the lifetimes of components have a general location-scale family. Deshpande et al. [10] developed a general semiparametric multivariate family of distributions. Their family can explicitly characterize CDSs using proportional conditional hazards, and they proposed a nonparametric test to test whether a failure happens earlier than the warrant or happens independently. Beutner [11] used a profile score process and multivariate intensity ratio to analyze the test statistics and investigate their asymptotical distribution. Bedbur [12] studied the uniformly most powerful unbiased tests with the assumption of conditionally proportional hazard rate for different hypotheses, in which the SOS samples are taken from a regular exponential family.
Schenk et al. [13] obtained the Bayes estimates of the Exponential distribution parameters based on multiple type-II censored SOS samples. They used inverse Gamma distributions as prior distributions to implement Bayesian estimation. Based on SOS sam-ples, Navarro and Burkschat [14] presented a new lifetime representation for coherent systems. They assumed that the component lifetimes in a coherent system have a specific dependence-type relationship. Balakrishnan et al. [15] investigated the maximum likelihood estimation method when the SOS samples were obtained from a step stress test. Burkschat and Navarro [16] examined mixture representations of the residual lifetime and the inactivity time of system with failure-dependent components. Burkschat and Torrado [17] examined the reversed hazard rate property of SOS and presented the conditions of reversed hazard rate ordering and the decreasing reversed hazard rate property of SOS. Park [18] proposed an expectation-maximization algorithm to obtain the MLE of the survival function for an equally load-sharing system. Burkschat et al. [19] studied statistical inference methods based on type-I censored SOS samples.
Sutar and Naik-Nimbalkar [20] examined the load-sharing phenomenon in a CDS under accelerated life testing. They also proposed a method for testing the dependence of component lifetimes. Esmailian and Doostparast [21] studied the maximum likelihood estimation method based on Weibull progressive SOS samples with the conditionally proportional hazard rate. Shafay et al. [22] used Bayesian estimation methods to predict future SOSs for one-and two-parameter Exponential distributions. Hashempour and Doostparast [23] used a Fisher information matrix to obtain the approximate confidence intervals based on observed multiply system lifetimes. Bedbur et al. [24] proposed a multisample model for the SOS samples of general step-stress exponents. Their proposed methods can be applied for any absolutely continuous lifetime distribution.
Balakrishnan et al. [25] studied the reliability of a CDS when the baseline distribution of the component lifetimes in the CDS is Burr type-XII distribution (BXIID). They used the power-trend hazard rate (PTHR) function to characterize the increased hazard rate due to each component failure. Hashempour and Doostparast [26] studied the reliability of comparing heterogeneous exponential distributions based on multiply SOS samples and the assumption of a conditionally proportional hazard rate. Bedbur et al. [27] released the commonly used proportional assumption and allowed the hazard rates of the baseline lifetime distributions to have situations such that each component failure can influence the hazard rate of surviving components. When the component lifetimes have heterogeneous Exponential distributions, Hashempour and Doostparast [28] studied Bayesian estimation methods for multiple SOS samples. They used a generalized likelihood ratio test to test the homogeneity property. Baratnia and Doostparast [29] considered an extension of SOS to characterize the system lifetimes with independent but heterogeneous components for the distribution family studied by Burkschat and Navarro [30]. Barakat et al. [31] developed two pivotal quantities and established prediction intervals for future lifetimes from a two-parameter Exponential distribution based on GOSs with a general scheme. Katzur and Kamps [32] used a Bayesian two-class discrimination approach with known prior probabilities for classification based on the recorded data with a SOS scheme.

Motivation and Organization
One CDS application is in light-emitting diode (LED) panels. The LED panel is composed of an array of LEDs. Because the current on the LED panel is unchanged, each LED failure will enhance the stress loading of each surviving LED until the LED panel malfunctions. Tsai et al. [33] proposed a Bayesian estimation method to evaluate the reliability of the CDS with heterogeneity condition, and Tsai et al. [34] proposed a Bayesian method and a profile maximum likelihood estimation method to infer the reliability of a video graphics array (VGA) when the adapters come from dual suppliers based on interval-censored samples. When dual suppliers support components for a CDS company simultaneously, using a mixture distribution is better than using a single distribution for characterizing the baseline component lifetimes for the study of CDS reliability.
Among the existing relevant studies about the reliability inference of CDS based on SOS samples, we find that the parametric PTHR model proposed by Balakrishnan et al. [25] is flexible and simple to use. Balakrishnan et al. [25] presented the analytical statistical properties for their proposed maximum likelihood estimation method. It is important to consider the impact of the heterogeneous quality of component lifetimes from dual suppliers on the reliability inference of CDS with a PTHR function. Because the likelihood function based on the SOS sample drawn from a mixture-BXIID is very complicated, searching the MLEs of the CDS model parameters using gradient computation algorithms could be intractable. Therefore, the Bayesian estimation method is considered in this study to obtain the Bayes estimates of the CDS model parameters based on mixture-BXIID SOS samples. For Bayesian estimation, the algorithms of Gibbs sampling (GS), important sampling (IS), and Metropolis-Hastings (M-H) have been widely used to implement the Markov chain Monte Carlo (MCMC) approach. When the analytic form of the conditional posterior distribution is available, the GS MCMC method is efficient to obtain Bayes estimates. However, it often fails to obtain the analytic form of the conditional posterior distribution. In such situations, we can consider using either the IS MCMC or M-H MCMC method to replace the GS MCMC method to obtain Bayes estimates. However, we find that both GS MCMC and IS MCMC methods are difficult to directly apply in this study.
On the basis of the aforementioned discussions, the impact of the heterogeneous condition of the mixture baseline distribution on the parameter estimation of CDS is considered in this study. Moreover, maximum likelihood estimation and Bayesian estimation methods are proposed to infer the model parameters and the baseline survival function of CDS. In Section 2, we will show that the likelihood functions for the maximum likelihood estimation method in this study are very complicated. This fact could result in divergence results when using gradient iterative methods to find the MLEs of the model parameters. Another difficulty is getting the analytic form of the posterior distribution to find reliable Bayes estimates of the parameters in a CDS based on SOS samples. In order to overcome the computation complexity when obtaining the Bayes estimates of model parameters, the MCMC approach is used for Bayesian estimation, and we propose a new M-H within the GS algorithm to implement the MCMC approach. To the best of our knowledge, no existing works study the baseline distribution function of components in CDS by using SOS samples from a mixture of BXIIDs. We use the term "mixture-BXIID" to denote the mixture distribution of the two-component BXIIDs and keep the term BXIID to denote the typical Burr type-XII distribution hereafter. The major contribution of this study is to propose a maximum likelihood estimation and Bayesian estimation methods for estimating the model parameters and the baseline survival function of mixture-BXIID in the CDS. Moreover, the highest credible interval (HCI) and equal-tailed credible interval (ECI) of the baseline survival function are also obtained for the mixture-BXIID.
The rest of this paper is organized as follows. In Section 2, the BXIID and mixture-BXIID with the PTHR function are investigated based on SOS samples. The maximum likelihood estimation method and Bayesian estimation method are established. Moreover, we present the MH-GS MCMC approach to implement the Bayesian estimation method. In Section 3, the proposed MH-GS MCMC approach is illustrated with a numerical example. The parameter determination of the proposed MH-GS MCMC approach and the process via using Markov chains to establish the HCI and ECI of the baseline survival function of mixture-BXIID are also discussed in Section 3. In Section 4, Monte Carlo simulations are conducted to study the estimation performance of the proposed maximum likelihood estimation method and Bayesian estimation method in terms of the quality measures of bias and mean square error (MSE). Some concluding remarks are given in Section 5.

The Lifetime Distribution and Statistical Methods
Let the CDS components from a supplier have lifetimes that follow a baseline twoparameter BXIID, whose probability density function (PDF), cumulative density function (CDF), and hazard rate function are, respectively, defined by and where δ > 0 is scale parameter and β > 0 is the outer-shape parameter. This is a special case of three-parameter BXIID with inner-shape 1. Therefore, the survival function of each basic component is defined by Because the three-parameter BXIID with two shape parameters could result in a model identification problem for parameter estimation, we only keep the outer-shape and scale parameters for the mixture-BXIID model in this study when dual suppliers occur. Denote the BXIID defined in Equations (1)-(4) by BXIID(β, δ) hereafter. The format of the mixture-BXIID model for dual suppliers will be addressed later.
When all components come from a single supplier, each component failure means that the CDS equally redistributes the stress loading to all surviving components and changes their lifetime distributions. Under this situation, the lifetimes of the surviving components in the CDS follow BXIID with different hazard rates at each failure. At the jth SOS, the PDF, CDF, and hazard rate function of the lifetimes of surviving components . The PTHR model proposed by Balakrishnanet et al. [25], is considered to model the increased hazard rate due to component failure, where α j (≥ 1) denotes the multiple of baseline hazard rate for the hazard rate after the jth failure. The CDF corresponding to h j (x) in Equation (5) is When dual suppliers support the components simultaneously, the failure times of all baseline components should be modeled by the mixture-BXIID of BXIID(β, δ 1 ) and BXIID(β, δ 2 ). Then, the PDF of the baseline mixture-BXIID can be addressed by where We can represent Equation (6) by Let F i0 ≡ F i0 (x | β, δ i ) for i = 1, 2 be the corresponding CDFs. The CDF of the mixture-BXIID can be presented by Denote the mixture-BXIID defined by Equations (7) and (8) by mBXIID(β, δ 1 , δ 2 , p) hereafter.
Let y = x (1) , x (2) , · · · , x (r) = (y 1 , y 2 , · · · , y r ) denote the realizations of the SOS sample of the first r component failure times from the mBXIID(β, δ 1 , δ 2 , p). The likelihood function based on the sample y can be defined by where a n, . The likelihood function of Equation (9) can be represented by and the log-likelihood function can be presented by + log p f 10 y j + q f 20 y j } + m r log(1 − pF 10 (y r ) − qF 20 (y r )) + log(p f 10 (y r ) + q f 20 (y r )).
After straightforward algebraic computation, it can be shown that the first derivatives of (θ) with respect to β, δ 1 , δ 2 , θ and p are given as follows: where and ∂m r ∂θ = (n − r + 1)rθ r−1 .
Because Equations (12)- (16) are complicated, it is very difficult to obtain the second derivative of the log-likelihood function with respect to β, δ 1 , δ 2 , θ, and p, respectively, so as to obtain the Fisher information matrix. Therefore, it is intractable to find the asymptotical confidence intervals of the model parameters using the Fisher information matrix. In this paper, we propose a MH-GS MCMC method to obtain the Bayes estimates of the model parameters and the baseline survival function. Denote the Bayes estimators of β, λ, δ, θ, and p byβ B ,λ B ,δ B ,θ B , andp B , respectively. Consider the following joint prior distribution: where and π 5 (p) = 1, 0 < p < 1.
That is, we consider that the shape parameter is β ∼ Gamma(η 1 , κ 1 ) and the scale parameters are δ 1 ∼ Gamma(η 2 , κ 2 ) and δ 2 ∼ Gamma(η 3 , κ 3 ); θ follows a Uniform distribution over (b 1 , b 2 ), denoted by θ ∼ U(b 1 , b 2 ), and p ∼ U(0, 1). Because the PTHR function is used in this study, it can be shown that b 1 = 1, and b 2 is a constant larger than but close to 1. The determination of b 2 will be studied in Section 3.
Method I: The MH-GS MCMC method.
Step 1: Implementing Step 2 N times, i.e., 1 ≤ i ≤ N, where N is a big positive number.
Step 2: For i = i + 1, use Step 2.1 to Step 2.5 to update β (i) , δ 2 , θ (i) and p (i) . The acceptance criterion is based on the M-H algorithm, see in [35].
Step 3: Go to Step 2 until i = N.
Step 4: Denote the obtained Markov chains by β (i) , i = 1, 2, · · · , N , δ (i) 2 , i = 1, 2, · · · , N , θ (i) , i = 1, 2, · · · , N , and p (i) , i = 1, 2, · · · , N , respectively. Based on the square error loss function, the Bayes estimator can be the sample mean of the obtained Markov chain after burn-in, N 0 , and is denoted bŷ We suggest the proposed transition probability function, ω j (θ * j | θ (i) j ), as a symmetric function about θ (i) j for j = 1, 2, 3, 4, 5 to reduce computational loading. For example, Normal and Uniform distributions with means of θ (i) j are two widely used transition probability functions. Based on the Markov chains of β (i) , i = 1, 2, · · · , N , δ (i) 2 , i = 1, 2, · · · , N , θ (i) , i = 1, 2, · · · , N , and p (i) , i = 1, 2, · · · , N , we can establish the Markov chain of the baseline survival distribution If we use a non-informative prior distribution to obtain the Bayes estimates, the major contribution of the posterior distribution comes from the likelihood function. Subsequently, the obtained Bayes estimates of the parameters are closed to the MLEs. In this study, we can select the hyperparameters such that the variances of the Gamma prior distributions of β, δ 1 , and δ 2 are large enough to make the prior distribution non-informative. Then, the obtained Bayes estimates of parameters are closed to the MLEs. This method can overcome the possible divergence problem by using gradient methods of solving the likelihood equations of = 0, and ∂ (Θ) ∂p = 0 to obtain the MLEs of the model parameters.

Discussions
In this section, a numerical SOS sample of n = 30 and r = 15 was generated from the baseline of mBXIID(β = 5, δ 1 = 5, δ 2 = 10, p = 0.7) to demonstrate the implementation of the proposed MH-GS MCMC method. The PTHR model with θ = 1.02 is considered in this example. The scenario of this numerical example for CDS is to assume that the CDS is claimed as failure if 15 or half of the components fail; that is, the CDS is claimed as failed if half of the components in the CDS fail or the hazard rate increases to approximately 34.6%, h 15 (x) = θ 15 h 0 (x) = (1.346)h 0 (x), from the baseline hazard rate. In order to generate a SOS sample from the mixture baseline distribution of mBXIID(β = 5, δ 1 = 5, δ 2 = 10, p = 0.7), we need to generate one more random sample {u 1 , u 2 , ..., u r } from U(0, 1). If u i < p, then the ith SOS is generated from the BXIID (β = 5, δ 1 = 5); otherwise, the ith SOS is generated from the BXIID (β = 5, δ 2 = 10). The 15 generated SOSs from the mBXIID (β = 5, δ 1 = 5, δ 2 = 10, p = 0.7) are reported in Table 2. To determine the parameters of N and N 1 in the proposed MH-GS MCMC method for practical use, we first generate N = 11,000 Markov chains for the MH-GS MCMC method, in which the first N 1 = 1000 chains are removed for burn-in. The posterior distribution and proposals for the MH-GS MCMC method are set up as follows: the prior distribution of π(Θ) defined in Equation (17) is considered with the hyperparameters η 1 = η 2 = η 3 = 1, k 1 = β, k 2 = δ 1 , k 3 = δ 2 , b 1 = 1 and constant b 2 , where b 2 is the value when the hazard rate of CDS at the time of rth failure increases to 2 times of the baseline hazard rate. Let the one-step-ahead parameters of β * , δ * 1 , δ * 2 , p * , and θ * be β (i) , δ 2 , p (i) , and θ (i) , respectively. The normal distributions of β * ∼ N(β (i) , 1 2 ), δ * 1 ∼ N(δ (i) as proposals to generated the values of β * , δ * 1 , δ * 2 , p * , and θ * in Method I. If β ( * ) < 0, then i for i = 1, 2; and if p * < 0 or p * > 1, then p * = p (i) . We found that the Markov chains of β are highly correlated as their first-order autocorrelation coefficient is about 0.80, and the first-order autocorrelation coefficients of the other Markov chains of δ 1 , δ 2 , θ and p are smaller than 0.15. Therefore, we need to trim the Markov chains by spacing to reduce the autocorrelation among the Markov chains of β.
To study the adequate interval length for spacing by simulations, we regenerate Markov chains of length N = 31,000 and remove the first N 1 = 1000 chains for burn-in based on the SOS sample in Table 2 by using the proposed MH-GS MCMC method and check the first-order autocorrelation coefficient of the Markov chains of β by spacing with the interval lengths of 1, 3, 5, 10, 15, and 20, respectively. We find that the spacing of interval length 10 can significantly reduce the first-order autocorrelation coefficient among the Markov chains of β below 0.2. Based on the simulation findings, all Markov chains are slimmed by taking one of every 10 Markov chains to reduce the autocorrelation of Markov chains. Figures 1-5 present the 3000 slimmed Markov chains of β, δ 1 , δ 2 , θ and p. In the scan of Figures 1-5, we find that the proposed MH-GS MCMC method performs well to obtain the Bayes estimates of the CDS model parameters. The Bayes estimates areβ B = 4.471, δ 1 = 4.826,δ 2 = 9.706,θ B = 1.023, andp B = 0.60. All Markov chains well cover their true parameter. All the acceptance rates are high and reported by 0.715, 0.823, 0.95, 0.708, and 0.582 for the Markov chains of β, δ 1 , δ 2 , θ, and p, respectively. The first-order autocorrelation coefficients are 0.192, −0.001, 0.009, −0.004, and 0.058 for the slimmed Markov chains of β, δ 1 , δ 2 , θ, and p, respectively. Let x = 0.5, it can be shown that the baseline survival function is S g0 (x) = 1 − G 0 (x) = 0.669 for x = 0.5. The Markov chains of S g0 (x = 0.5) can be established through using the slimmed Markov chains of β, δ 1 , δ 2 , θ, p and Equation (8). The Bayes estimate of S g0 (x = 0.5) isŜ g0,B (x) = 0.709.
The Markov chains of S g0 (x) can be used to establish the empirical distribution of S g0,B (x), and then the HCI and ECI of S g0 (x) can be obtained. Based on the SOS sample in Table 2, the 95% HCI of S g0 (x = 0.5) is (0.588, 0.820) and the 95% ECI of S g0 (x = 0.5) is (0.581, 0.815). The coverage probabilities of the 95% HCI and ECI of S g0 (x = 0.5) are evaluated using the proposed MH-GS MCMC method with 1000 iterations. We obtain coverage probabilities of 0.961 and 0.974 for the HCI and ECI of S g0 (x = 0.5), respectively. The simulation results indicate that the coverage probability of HCI is closer to the nominal coverage probability than the coverage probability of ECI. Therefore, the HCI outperforms the ECI for the proposed MH-GS MCMC method in terms of the coverage probability.

Monte Carlo Simulations
To explore the estimation quality of the proposed Bayesian estimation method using an MH-GS MCMC approach, an intensive simulation study was conducted, and the proposed Bayesian estimation method is used to obtain the Bayes estimates of the CDS model parameters based on the SOS samples. As the parameter setting in Section 3, we consider the baseline of mBXIID (β = 5, δ 1 = 5, δ 2 = 10, p = 0.7) for simulations with the parameters of (n, r) = (30, 15), (50, 25), (80, 30), (100, 30), θ = 1.02, and p = 0.2, 0.3, 0.5. The CDS is defined as failure if half or over 30 of the components fail. We are also interested to evaluate the estimation quality of the proposed Bayesian estimation method for different mixture probabilities. We use N = 31,000, N 1 = 1000, and interval length 10 for spacing to implement the MH-GS MCMC method. Moreover, the posterior distribution of π(Θ) in Equation (18) with the hyperparameters η 1 = η 2 = η 3 = 1, k 1 = β, k 2 = δ 1 , k 3 = δ 2 , b 1 = 1 and constant b 2 is used to obtain the Bayes estimates of the model parameters, where b 2 is the value when the hazard rate of CDS at the time of rth failure increases to 2 times of the baseline hazard rate. The proposals of , and θ * ∼ U(1, b 2 ) are used to generated the values of β * , δ * 1 , δ * 2 , θ * and p * in Method I. If β * < 0, then i for i = 1, 2; and if p * < 0 or p * > 1, then p * = p (i) . We repeat the proposed Bayesian estimation method 10,000 times to obtain 10,000 Bayes estimates of β, δ 1 , δ 2 , θ and p, then the 10,000 Bayes estimates are used to evaluate the bias and MSE for each Bayes estimate. Because the scales of the parameters are different, we use relative bias and relative square root of MSE, which are free of the scale, as the performance measures to study the estimation quality of the proposed Bayesian estimation method. The relative bias stands for the bias over the true parameter, and the relative square root of MSE stands for the square root of MSE over the true parameter. We use rBias and rsqMSE to denote the relative bias and relative square root of MSE for simplicity hereafter. All simulation results based on the proposed Bayesian estimation method are reported in Tables 3-5. In order to compare the estimation performance of the proposed Bayesian estimation method with the maximum likelihood estimation method, the MLEs of β, δ 1 , δ 2 , θ, and p are also obtained and denoted byβ,δ 1 ,δ 2 ,θ, andp, respectively. Moreover, the rBias and rsqMSE of each MLE are evaluated based on 10,000 obtained MLEs. Because there is a high probability that the gradient methods cannot successfully solve log-likelihood equations to obtain the MLEs of model parameters, we use the proposed Bayesian estimation method with a non-informative prior distribution to obtain approximate MLEs. Let κ 1 = κ 2 = κ 3 = 100 and η 1 = η 2 = η 3 = 1 in Equation (17) to make the Gamma distributions of β, δ 1 , and δ 2 have a big variance, and then the prior distribution of Equation (17) becomes non-informative. All simulation results based on the maximum likelihood estimation method are reported in Tables 6-8. In view of Tables 3-8, we find the following results. Tables 3-5 are more reliable and outperform the MLEs in  Tables 5-7 in terms of the measures of bias and MSE. The rBias and rsqMSE of the Bayes estimate are smaller than the rBias and rsqMSE of MLE. Tables 3-5 indicate that the estimation quality of the Bayes estimate for the mixture proportion is worse than the estimation quality of the Bayes estimates for other parameters. Moreover, the censoring rate significantly influences the estimation quality of Bayes estimates. 3. The maximum likelihood estimation method is less reliable for estimating the mixture proportion and shape parameter. Moreover, the acceptance rate is low for the Markov chains of the parameters of θ and p if non-informative prior is used to implement the Bayesian estimation methods. The rBias and rsqMSE of the MLEs of β and p in Tables 6-8 are large even when the sample size increases to 100 or the censoring rate is high. The findings mean that it could fail to use the maximum likelihood estimation method to obtain reliable estimates of the model parameters if the component lifetimes in CDS follow a baseline mixture distribution.

Conclusions
In this study, we proposed a maximum likelihood estimation method and a Bayesian estimation procedure using the MH-GS MCMC approach to obtain reliable Bayes estimators of the CDS model parameters and baseline survival function based on SOS samples when the component lifetimes follow the baseline of mixture-BXIID. The PTHR function is used to characterize the increased system hazard caused by component failures. For each component failure, the stress loading can be equally redistributed to the surviving components in the CDS. The proposed methods include the Bayesian estimation for a single baseline of BXIID for CDS as special case when p = 1.
The implementation of the proposed Bayesian estimation method with the MH-GS MCMC approach is presented systematically in Method I. The proposed Bayesian estimation method is studied by utilizing a numerical example. An intensive simulation study is conducted to evaluate the performance of the proposed Bayesian estimation method and the maximum likelihood estimation method. Based on the simulation results in Section 4, we find that the proposed Bayesian estimation method outperforms the maximum likelihood estimation method to obtain reliable estimates of the CDS model parameters and its baseline survival function in terms of the relative bias and MSE.
The maximum likelihood estimation method does not perform well in obtaining reliable MLEs of the shape parameter and mixture proportion, and it is also intractable to obtain the Fisher information matrix to construct an approximate confidence interval for the baseline survival function of mixture-BXIID. Therefore, one major contribution of this study is to propose a simple Bayesian estimation method to overcome the divergence problem of the typical maximum likelihood estimation method to obtain approximate MLEs using a non-informative prior distribution. Moreover, the proposed HCI and ECI methods based on using the Markov chains of the baseline survival function can help practitioners to carry out interval inference for the baseline survival function.
The simulation results in Section 4 also show that the estimation quality of the Bayes estimate for the mixture proportion is worse than the estimation quality of the Bayes estimates for other parameters. Moreover, the censoring rate significantly influences the estimation quality of Bayes estimates. When the censoring rate is low, the rBias and rsqMSE of the Bayes estimate could be less stable and fail to decrease along the sample size. There is also room to improve the estimation quality of using the proposed Bayesian estimation method to estimate the shape parameter of the mixture-BXIID. It is also important to extend the proposed Bayesian estimation method to other mixture distributions. These topics can be studied in the future.