Bayesian Estimation Based on Sequential Order Statistics for Heterogeneous Baseline Gompertz Distributions

: A composite dynamic system (CDS) is composed of multiple components. Each component failure can equally induce higher loading on the surviving components and, hence, enhances the hazard rate of each surviving component. The applications of CDS and the reliability evaluation of CDS has earned more attention in the recent two decades. Because the lifetime quality of components could be inconsistent, the lifetimes of components in the CDS is considered to follow heterogeneous baseline Gompertz distributions in this study. A power-trend hazard rate function is used in order to characterize the hazard rate of the CDS. In order to overcome the difﬁculty of obtaining reliable estimates of the parameters in the CDS model, the Bayesian estimation method utilizing a hybrid Gibbs sampling and Metropolis-Hasting algorithm to implement the Markov chain Monte Carlo approach is proposed for obtaining the Bayes estimators of the CDS parameters. An intensive simulation study is carried out to evaluate the performance of the proposed estimation method. The simulation results show that the proposed estimation method is reliable in providing reliability evaluation information for the CDS. An example regarding the service system of small electric carts is used for illustration.


Literature Review
The composite dynamic system (CDS) is composed of n(≥ 1) identical components. The components fail one-by-one over time until r failures are observed (r ≤ n), and then the CDS is defined as malfunctioned. Each component failure in the CDS is assumed to equally distribute higher loading on the surviving components and enhance the hazard rate of the surviving components. Today, the CDS is one widely used electronic device and it is also named r-out-of-n failure system. The CDS reduces to a series system as r = 1 and reduces to a parallel system as r = n. The corresponding ordered failure times from the r failed components are also named sequential order statistics (SOSs) from the r-out-of-n failure system.
In the past two decades, researchers have paid attention on studying the reliability of a CDS. Kamps [1] proposed the concept of generalized order statistics as a unified approach for different models of ordered random variables. Cramer and Kamps [2] studied the reliability of (n − r + 1)-out-of-n failure system for the exponential distribution. Cramer and Kamps [3] investigated the maximum likelihood estimation method, uniformly minimal variance unbiased estimation method, and best linear unbiased estimation method based on SOS samples. Moreover, they used these three estimation methods in order to obtain the estimates of parameters in the exponential distribution. With no any restrictions being imposed on the model parameters, Cramer and Kamps [4] derived the marginal distribution. Balakrishnan, Jiang, Tsai, Lio, and Chen [24] proposed inferential methods for the CDSs, in which components in a CDS follow a baseline Burr type XII distribution with a power-trend hazard rate function. Hashempour and Doostparast [23] evaluated the reliability of comparing heterogeneous exponential distributions that are based on independent multiply SOS samples with using the conditionally proportional hazard rate assumption. Bedbur, Burkschat, and Kamps [25] released the common proportional assumption and allowed for the hazard rates of the baseline lifetime distributions in situations; the failures can have an impact on the entire shape of the hazard rate of surviving components. Hashempour and Doostparast [26] studied Bayesian inference methods on multiple SOS samples for heterogeneous exponential distributions. They also studied using the generalized likelihood ratio test in order to test the homogeneity property. Baratnia and Doostparast [27] proposed an extension of SOS to model the system lifetimes with independent, but heterogeneous, components for the distribution family that was studied by Burkschat and Navarro [16].

Motivation and Organization
Numerous reliability inference studies regarding CDS have been investigated in literature that is based on SOS samples. Among them, the parametric hazard rate model that was proposed by Balakrishnan, Jiang, Tsai, Lio, and Chen [24] is flexible and simple for obtaining the CI of the baseline survival function of CDS components. Balakrishnan, Jiang, Tsai, Lio, and Chen [24] have presented the analytical statistical properties for the maximum likelihood estimation method. However, two problems are found when implementing the maximum likelihood estimation method that was proposed by Balakrishnan, Jiang, Tsai, Lio, and Chen [24]. First, the MLEs of the model parameters could be difficult to obtain if the sample size is small. Second, the lifetime components in many instances could be manufactured with heterogeneity conditions due to ill-trained operators, while using heterogeneous materials in production or other reasons.
Assume that the target baseline lifetime distribution of the components in a CDS is Gompertz with λ = 2 (scale parameter ) and γ = 3 (shape parameter), and the producer of CDS buys components from one specific supplier, say Supplier A. If the quality of components from Supplier A is consistent, the histogram of the homogeneous data set can be displayed in Figure 1. It is often that the quality of components from the supplier is inconsistent. In order to catch the impact of the inconsistent quality of components on the CDS, we generated 1000 Gompertz observations with scale parameter λ and shape parameter γ, but both of the parameters are random. λ follows the Gamma distribution with shape parameter 2 and scale parameter 1 and γ follows the Gamma distribution with shape parameter 3 and scale parameter 1; that is, the mean value of the Gamma distribution is 2 for λ and 3 for γ. These two mean values are the true parameters in the homogeneous case. Figure 1 shows two more histograms that are based on two heterogeneous data sets that are generated from the aforementioned heterogeneous case. That means that the quantile of the heterogeneous case could be different from the quantile that is based on the homogeneous case. Subsequently, we could obtain different reliability assessment results if we ignore the heterogeneity condition and undertake reliability inference that is based on the homogeneity assumption. It is an important issue to study the impact of the heterogeneity condition on the reliability of a CDS based on SOS samples.
To our best knowledge, the inference for the CDS parameter when the components follow heterogeneous Gompertz distributions is still an open question. It is difficult to evaluate the CDS's reliability under a heterogeneity condition. One method is to assume the CDS's lifetime can be characterized by a mixture model. However, two difficulties are found when using a mixture model for the CDS's lifetime. The first difficulty is that the number of candidate distributions for mixture is unknown. The second is that the quality of parameter estimation becomes very difficult to control. In this study, we consider that the model parameters are random and they follow a joint prior distribution to handle the heterogeneity condition. Subsequently, a Bayesian estimation procedure using a hybrid Gibbs sampling and Metropolis-Hasting (hGSMH) algorithm is proposed in order to overcome the computation intractability to implement the Markov chain Monte Carlo (MCMC) approach for obtaining the BEs of the CDS parameters. Denote the proposed estimation procedure by the hGSMH-MCMC method here and after. The rest of this study is organized, as follows: the statistical model of heterogeneous baseline Gompertz distributions with SOS samples is investigated in Section 2. The powertrend hazard rate function is considered for the CDS to obtain the likelihood function. Moreover, the posterior distribution of the statistical model that is based on the SOS sample from the heterogeneous baseline Gompertz distributions is analytically obtained. The proposed hGSMH-MCMC method is also proposed to obtain the BEs of parameters in the CDS. An intensive simulation study is conducted in Section 3 in order to evaluate the quality of the proposed hGSMH-MCMC method in terms of the measures of bias and mean square error (MSE) of estimates. Moreover, an example regarding the service system of small electric carts is used for illustrating the applications of the proposed method. The simulation results show that the proposed estimation method is reliable for obtaining the BEs of parameters in the CDS. Section 4 provides some concluding remarks.

The Statistical Model
Let the lifetimes of components in a CDS follow a baseline Gompertz distribution, whose probability density function (PDF), cumulative density function (CDF), and hazard rate function are defined, respectively, by and where λ > 0 is the scale parameter and γ > 0 is the shape parameter of the baseline distribution. The Gompertz distribution can reduce to an exponential distribution as γ → 0. Many authors, including Read [28], Makany [29], Rao and Damaraju [30], Franses [31], and Wu and Lee [32], have studied the statistical methodologies and characterizations of the Gompertz distribution since it was introduced by Gompertz [33] to characterize the human mortality and establish actuarial tables. Garg, Rao and Redmond [34] investigated the properties of this distribution. Moreover, he proposed the maximum likelihood estimation method to estimate the distribution parameters of Gompertz. Gordon [35] established the maximum likelihood procedure for the mixture of two Gompertz distributions. Chen [36] proposed an exact confidence interval for γ and an exact joint confidence region for both parameters by using the properties of the conventional exponential distribution. Wu et al. [37] proposed un-weighted and weighted least-square estimates for λ and γ based on first failure-censored samples. Willekens [38] presented the connections among the Gompertz, Weibull, and other type I extreme value distributions. Assume that, when the jth component fails, surviving components in the CDS can equally share the loading. It can be shown that the lifetimes of surviving components still follow a Gompertz distribution, but they have different hazard rates. Following the equalsharing condition for the surviving components in the CDS, the hazard rate of surviving components is increased. The PDF, CDF, and hazard rate function can be defined by the The proportional hazard model assumption that is given in Equation (4) is equiv- For simplifying the notations, let y j = x (j) , j = 1, 2, . . . , r, y = (y 1 , y 2 , . . . , y r ) denote the realizations of the SOS and Θ = (λ, γ, α 1 , . . . , α r ). Following the inference procedure of Beutner [11], the likelihood function that is based on the sample y can be presented by L(Θ|y) = a n,r where a n,r = n!/(n − r)! and m j = (n − j + 1)α j − (n − j)α j+1 − 1 for j = 1, 2, . . . , r − 1. Based on Equations (1), (2), and (5), Balakrishnan, Jiang, Tsai, Lio, and Chen [24] presented the log-likelihood function by L(Θ|y) = a n,r Use the power-trend hazard rate model of α j = θ j , j = 1, 2, . . . , r in Equation (6). Subsequently, Θ can be reduced to Θ = (λ, γ, θ), Equation (6) can be represented as Equation (7) and the log-likelihood function of can be presented as Equation (8): and (Θ) = log(L(Θ|y)) ∝ r(r + 1) 2 log(θ) + r log(λ) + γ r ∑ j=1 y j + log(ω 1,y (Θ)) + log(ω 2,y (Θ)), (8) where ω 1,y (Θ) = e − λ γ ∑ r−1 j=1 (m j +1)(e γy j −1) and ω 2,y (Θ) = e − λ γ θ r (n−r+1)(e γyr −1) . The MLEs of γ and θ that are based on SOS samples do not have explicit form and only the MLE of λ, given γ and θ, can be obtained with an explicit form as the following: It could be difficult to set up the heterogeneity condition in the CSD model to conduct the maximum likelihood estimation procedure. One feasible method is to characterize the heterogeneity among components in the CDS by assuming that the models parameters are random, and then to obtain the BEs of the model parameters for reliability evaluation. Denote the BEs of λ, γ and θ byλ B ,γ B , andθ B , respectively. Let the prior distribution of Θ be defined by where and where η 1 , η 2 , δ 1 , and δ 2 are hyper-parameters. The Gamma distribution has been widely used as the prior distributions of the parameters of the Gompertz distribution, see, for example, Soliman et al. [39], Dey et al. [40], and Chacko and Mohan [41]. Some of the authors considered using non-informative prior distribution as the prior distributions of the parameters of the Gompertz distribution in order to implement Bayesian estimation, see, for example, Ismail [42] and Feroze and Aslam [43]. When the Gamma prior distribution has a large variance, the Gamma prior distribution approaches a non-informative distribution.
In this study, we consider using Gamma distribution as the prior distributions of the scale and shape parameters of the Gompertz distribution, respectively. Moreover, a MCMC approach is proposed in order to overcome the computation complexity that is caused by using the Gamma prior distribution. Using the derivation in Appendix A, the posterior distribution of Θ, given the sample y, can be obtained by where After algebraic computation, we can obtain the conditional posterior distributions of θ, λ and γ by and respectively. It can be found that π 1,y (θ|λ, γ), π 2,y (λ|γ, θ) and π 3,y (λ|θ, γ) are not conjugate distributions. Hence, it is difficult to update all of the parameters via using the marginal posterior distributions of π 1,y (θ|λ, γ), π 2,y (λ|γ, θ) and π 3,y (λ|θ, γ) in order to implement the Gibbs sampling algorithm. Assume that λ (i) , γ (i) and θ (i) are the values of λ, γ and θ at the previous state, respectively. Let θ * be generated from a proposal, say q 1 (θ * |θ (i) ). Given λ (i) and γ (i) , it is trivial to show that the ratio of is simple for computation, due to κ 1 only depending on the right hand side of Equation (18). The Metropolis-Hastings algorithm uses the ratio of κ 1 in order to update θ. Likewise, assume that λ * and γ * are generated from the proposals of q 2 (λ * |λ (i) ) and q 3 (γ * |γ (i+1) ), respectively, and the updated value of θ is θ (i+1) . Given γ (i) and θ (i+1) , we can update λ by the ratio of Let the updated value of λ be λ (i+1) . Given λ (i+1) and θ (i+1) , we can update γ by the ratio of respectively. The steps to generate the Markov chains of θ, λ and γ via using the proposed hGSMH-MCMC method are presented in Procedure 1. Please note that the proposed hGSMH-MCMC method is also a Metropolis-Hastings-within Gibbs approach. The proposed hGSMH-MCMC method combines the merits of the Metropolis-Hastings and Gibbs algorithms. Hence, the proposed hGSMH-MCMC method is simple to use as using the typical Metropolis-Hastings algorithm, but it has a higher update rate than the typical Metropolis-Hastings algorithm.
Step 1: if i = N go to Step 4; and go to Step 2 otherwise, where N is a big positive number.
Step 4: stop. Figure 2 provides the flowchart of Algorithm 1. Based on the primacy of mathematical convenience, we consider the two most commonly used loss functions, the squared loss and absolute loss functions, to obtain the BEs of the CSD parameters in this study. The BEs ofθ B ,λ B andγ B based on the squared loss function can be the sample mean of the Markov . . , N}, respectively. The performance of the proposed hGSMH-MCMC method will be evaluated through using the Monte Carlo simulations shown in Section 3.
The confidence interval inference of the baseline survival function, can provide important information for assessing the reliability of the CDS regarding the probability of a CDS can survive longer than a specific time of x 0 . The Markov chain of S 0 (x 0 ) based on the Markov chains of λ and γ can be obtain and denoted by {S Subsequently, the empirical distribution of parameter can be established based on using { (i) , i = N 1 + 1, N 1 + 2, · · · , N}. In this study, can be λ, γ, θ, or S 0 (x 0 ). The credible interval, as denoted by (L, U), which admits an interpretation of (1 − 2α) × 100% posterior probability of covering the true , can be obtained by the following Procedure 2: find an (1 − 2α) × 100% credible interval of .

Initial
Step: let i = 1. If i < B, then go to Step 1; and go to Step 5 otherwise, where B is a large positive integer.
Step 2: find the αth and (1-α)th quantiles, (d L ) and (d U ) , from the sorted Markov chains N 1 ) and x denotes the largest integer smaller than or equal x.

Markov Carlo Simulations
An intensive Monte Carlo simulation study is conducted in order to evaluate the performance of the proposed hGSMH-MCMC method. Size n random samples of lifetimes are generated from the baseline Gompertz distribution with parameters λ = 2 and γ = 3. The censoring rate is r = c × n, where c = 0.3, 0.5, 0.6 and 0.7. In this study, we consider n = 30, 50 and 100 for simulations. Traditionally, it is difficult to obtain the MLEs of the CDS parameters for the homogeneous cases when the sample size is small. However, it could be difficult to collect a large sample in order to evaluate the quality or reliability of the CDS in practical applications.
In order to overcome this difficulty, we can consider using Gamma prior distributions with a big variance to implement the proposed hGSMH-MCMC method. Subsequently, the obtained BEs can be closed to the MLEs of the CDS parameters due the prior distributions in the proposed Bayesian estimation method being close to non-informative prior distributions. The aforementioned consideration and the proposed hGSMH-MCMC method can be used to obtain the BEs of the CDS parameters under the heterogeneity condition. We consider the following five scenarios of prior parameter combinations shown in Table 1. Scenarios S1 and S2 indicate that the Gamma prior distributions have a mean to match the values of λ = 2 and γ = 3, respectively; that is, we have information to set up the prior distribution in order to release the impact of heterogeneity condition. Scenario S3 indicates that the prior distributions are close to non-informative prior distributions due to two prior distributions of λ and γ have a big variance. Hence, the obtained BEs that are based on Scenario S3 are closed to the MLEs of the CDS parameters. Because the parameters have different scales, the relative bias (rBias) and relative square root of MSE (rsqMSE), which are defined as follows, are considered as the quality measures to evaluate the performance of the BEs that are obtained via using the proposed hGSMH-MCMC method: and The proposals of normal distributions, N λ (i−1) , 1 and N γ (i−1) , 1 are used to generate new values of λ (i) and γ (i) . Let λ * ∼ N λ (i−1) , 1 . If λ * < 0, then λ (i) = λ (i−1) and λ (i) = λ * otherwise. Likewise, let γ * ∼ N γ (i−1) , 1 . If γ * < 0, then γ (i) = γ (i−1) and γ (i) = γ * otherwise. The estimation of power-trend hazard rate parameter, θ, is really a problem in CDS reliability inference. Because we use the power-trend hazard rate model that was proposed by Balakrishnan, Jiang, Tsai, Lio, and Chen (2015) in this study, the true θ should be closed to 1. This property makes that the domain of θ can be easily determined.
In this study, we generate the new θ to update θ (i) from θ * ∼ N θ (i−1) , 0.0001 . If θ * > 1.01, θ (i) = θ (i−1) and update θ (i) by θ * otherwise. The assumption of the true θ is close to 1 is reasonable when the power-trend hazard rate model is used. N =10,000 Markov chains are generated in the proposed hGSMH-MCMC method and the first N 1 = 1000 Markov chains are removed for burn-in. Moreover, 1000 BEs are used in order to obtain the values of the rBias and rsqMSE. That is, the proposed hGSMH-MCMC method with N =10,000 and N 1 = 1000 are repeated 1000 times to obtain the values of the rBias and rsqMSE for each BE. In order to catch the heterogeneity condition assumption, we use 1000 generated values of λ and γ from the Gamma distributions with rate parameter 1 and shape parameter λ and γ, respectively, as the true parameters of λ and γ in order to evaluate the values of rBias and rsqMSE ofλ B andγ B . All of the simulation results are reported in Tables 2-7, in which Si-S and Si-A means the BE is obtained based on the square or absolute loss functions for Scenario Si, respectively, for i = 1, 2, 3.  Tables 2-4 report the values of rBias and rsqMSE that are based onθ B ,λ B , andγ B for n = 30, 50 and 100, respectively. From Tables 2-4, we first find the following results:

1.
When compared with the BEs in S1 and S2, the BEs in S3 are less reliable. Because the BEs in S3 are closed to the MLEs of θ, λ, and γ and the maximum likelihood estimation ignores the heterogeneity condition. Hence, most of the values of rBias and rsqMSE of the BEs shown in Scenario S3 of Table 2 to Table 4 are larger than that in Scenarios S1 and S2. We can claim that the maximum likelihood estimation cannot work well in obtaining reliable estimates of the CDS parameters if the heterogeneity condition is true.

2.
The BE that is based on the square loss function is more competitive than the BE based on the absolute loss function when the heterogeneity condition is true. Most of the values of rBias and rsqMSE of the BE based on the squared loss function are smaller than those based on the absolute loss function.

3.
A sample size of 50 or larger with at least 15 failure times can be used in order to obtain reliable BEs of the CDS parameters via using the proposed hGSMH-MCMC method. The BE of the baseline survival function is also evaluated via using the proposed hGSMH-MCMC method in order to obtain more information for evaluating the quality of the proposed hGSMH-MCMC method on the evaluation for the baseline survival function. Moreover, we use the Markov chains ofγ B andθ B and Equation (9) in order to obtain the estimate of λ, as denoted byλ. Tables 5-7 report all of the simulation results. In the scan of Tables 5-7, we also find the following results: 1.
The BEs of the survival function in Scenario S3 are less reliable when compared with that in Scenario S1 and S2. Moreover, the heterogeneity condition significantly affects the quality of the estimation method and makes the BE of the scale parameter less stable than the BE of the shape parameter. Because the baseline survival function is a function of the scale and shape parameters, the heterogeneity condition significantly affects the quality of the BE of the baseline survival function and makes the values of the rBias and rsqMSE of the BE of the baseline survival function significantly larger than the BEs of the scale or shape parameters.

2.
It is very difficult to use the gradient method to maximize the log-likelihood function in order to obtain the MLEs of the CDS parameters even under the homogeneity condition because of the divergence problem in numerical computation. We can find that most of the values of the rBias and rsqMSE ofλ are larger than that ofλ B in Tables 2-4. This is additional evidence to indicate that the proposed hGSMH-MCMC method outperforms the maximum likelihood estimation to obtain reliable estimates of the CDS parameters.
In summary of the simulation results, we recommend using the proposed hGSMH-MCMC method to replace the maximum likelihood estimation in order to obtain reliable estimates of the CDS parameters and implement reliability inference for the CDS.

An Example
In the real world, there are many CDS examples. For example, a light emitting diode (LED) panel is composed of an array of LEDs. Each LED failure will enhance the risk of the LED panel failure due to the surviving LEDs need to share the current. Another example regards the small electronic service system. A data set of the first failure times in month of small electric carts, which are used for internal transportation and delivery in a large manufacturing facility, is used in order to illustrate the proposed hGSMH-MCMC method. Zimmer, Keats, and Wang [44] have studied this example and used Burr XII distribution to model the first failure times of small electric carts. Because the shortest first failure time of small electric cart is smaller than 1, we cut the shortest first failure times and take logarithm transformation for the first failure times in the small electric carts data set. Let X denote the first failure times of small electric carts and t i = log(x i ), i = 1, 2, . . . , 19. Table 8 reports the 19 logarithm-transformed first failure times of small electric carts.
The MLEs of λ and γ in the Gompertz distribution can be obtained by λ = 0.0594 and γ = 1.0674, respectively, via using the R package maxLik. Denote the first failure time of small electric carts in Table 8 by t j and let u j = j−0.5 n , j = 1, 2, . . . , 19. Let Gomp.t j = F −1 0 u j |λ,γ denote the quantile of u j by replacing λ =λ and γ =γ. Figure 4 displays the quantile-quantile (QQ) plot of the logarithm-transformed first failure times of small electric carts in Table 8 with the quantiles of Gomp.t j , j = 1, 2, . . . , 19 and indicates that the Gompertz distribution with λ = 0.0594 and γ = 1.0674 can characterize this data set well. The Kolmogorov-Smirnov (K-S) statistic is 0.0526 with a p-value close to 1. Hence, the K-S test supports the Gompertz distribution with λ = 0.0594 and γ = 1.0674 can characterize this data set well.  Assume that n = 30 small electric carts are used for internal transportation and delivery in a large manufacturing facility. Each small electric cart failure can equally induce higher loading on surviving small electric carts and, hence, enhance the hazard rate of each surviving small electric carts. The service system is shut down if half or r = 15 small electric carts in the service system are malfunctioned. While using θ = 1.002 a SOS sample with heterogeneity condition is generated and reported in Table 9, in which λ is generated from the Gamma distribution with shape parameter 0.0594 and rate parameter 1; and, γ is generated from the Gamma distribution with a shape parameter 1.0674 and rate parameter 1. The proposed hGSMH-MCMC method is used in order to obtain the BEs of θ, λ and γ. We would like to establish Markov chains for each parameter with burn-in to remove the leading Markov chains. During the numerical computation procedure, we find that the autocorrelation among the generated Markov chins disappear slowly. We generated 500,000 Markov chains and removed the first 1000 Markov chains for burn-in in order to reduce the autocorrelation among the Markov chains in this example. Subsequently, take 1 for every 100 Markov chains to reduce the autocorrelation among the Markov chains. Figures 5-7 show the Markov chains of θ, λ and γ, respectively. The acceptance rates are 1, 0.8857 and 1 for θ, λ and γ respectively. All the acceptance rates are high. Based on the square loss function, the obtained BEs areθ B = 1.0003 ,λ B = 0.0898 andγ B = 0.9520, respectively. Because the heterogeneity condition is considered in the example, the sample means of 100 generated data sets of λ and γ are 0.0680 and 1.0845, respectively. The Markov chains of the baseline survival function for x 0 = 0.8 can be obtained via using Equation (24) and are reported in Figure 8. The BE ofŜ B (x 0 = 0.8) = 0.9056 and the 95% credible interval of S 0 (x 0 = 0.8) is (0.8148, 0.9666).

Conclusions
In this study, we proposed a hGSMH-MCMC method in order to obtain the BEs of the CDS parameters when the lifetimes of components follow a Gompertz distribution with a heterogeneity condition. Moreover, the hazard rate of the CDS is characterized by the power-trend hazard rate function. A Bayesian estimation method using the MCMC approach is proposed in order to overcome the difficulty during obtaining reliable estimates of the parameters in the CDS model.
A hybrid Gibbs sampling and Metropolis-Hasting algorithm are used to implement the MCMC approach. The performance of the proposed hGSMH-MCMC method is evaluated via using an intensive simulation study. The simulation results show that the proposed estimation method is reliable in providing reliability evaluation information for the CDS when the quality of the components has a heterogeneity condition. The simulation results show that the proposed method can provide reliable BEs for the CDS parameters.
The quality of the BE for the power-trend hazard rate parameter that is based on using the proposed hGSMH-MCMC method can be further improved to reduce the autocorrelation. When compared with the impact of the autocorrelation on the Markov chains of the scale and shape parameters of the Gompartz distribution, the impact of autocorrelation on the Markov chains of the power-trend hazard rate parameter is more difficult to reduce, even through a slim operation. How to obtain a reliable estimate of the power-trend hazard rate parameter is still an open question in the parameter estimation study of the CDS. The posterior predictive checking method that was proposed by Gelman, Carlin, Stern, Dunson, Vehtari and Rubin [45] is a good method for implementing predictive checks. However, the observations in a SOS sample are dependent and non-identically distributed. It cannot use the predictive checking method that was proposed by Gelman, Carlin, Stern, Dunson, Vehtari, and Rubin [45] to deal with the SOS samples. How to develop an improved predictive checking method to deal with SOS samples is also an open question. These two topics will be studied in the near future.