## 5 Simulation and Data Analysis
### Simulation Study  
Without losing generality, in the simulation, we set \(s=2\); that is, there are two latent risks. To begin with, we generate competing risks sample data from the Burr-XII distributions randomly by the inverse transform sampling method. In this simulation, we set \(a_{1}=0.5\), \(a_{2}=1\), and \(b=1\). The sample size \(n\) takes values of \(n=30\), \(n=50\), and \(n=100\), respectively. In addition, we generate the left endpoints and interval lengths from exponential distributions with mean \(\frac{1}{a}\) and \(\frac{1}{\beta}\), respectively. We consider the following combinations for the parameter pairs \((\alpha,\beta)\): \((0.5,0.5)\), \((0.5,1)\), and \((1,0.5)\). By comparing actual lifetimes with the corresponding interval, we can obtain the middle-censored data. Then, we can calculate the MLEs and 95% ACIs for \(a_{1}\), \(a_{2}\), and \(b\) by employing the _optim_ function or the EM algorithm. For Bayesian estimation, regarding the informative prior, to guarantee that the expected value of the prior distribution of \(b\) coincides with the true value, we set that \(b\) has an exponential prior distribution with \(\lambda=1\). In addition, we select the hyper-parameters as \(\theta_{1}=0.5\) and \(\theta_{2}=1\). This method to select hyper-parameters has been used by some scholars such as [8; 21]. In practical scenarios, researchers often lack prior knowledge, necessitating the use of non-informative priors in Bayesian methodologies. For parameters \(a_{j},j=1,2,\ldots,s\), where Gamma distributions are adopted as the prior, it is common practice to set the shape parameter \(\theta_{j}\) to a decimal value close to 0 in the calculation process. This approach aligns with established methodologies in existing literature, such as references [8; 20]. For parameter \(b\), we introduce an improper uniform prior density, defined over the positive real numbers (\(b>0\)). Then, the generalized non-informative prior density function of \(b\) and the PDF of the joint posterior distribution function can be expressed as  
\[ \begin{array}{c}\pi_{0}^{0}(b)=1,b>0,\\ \\ \pi^{0}(\mathbf{a},b|data)\propto b^{\sigma+n_{1}}e^{-ba^{*}}(\frac{1}{a^{*}}) ^{n_{2}}\prod_{j=1}^{s}a_{j}^{m_{j}+\theta_{j}-1}\\ \\ \times\prod_{i=1}^{n_{1}}[t_{i}^{b-1}(1+t_{i}^{b})^{-(a^{*}+1)}]\prod_{i=n_{1} +1}^{n_{1}+n_{2}}[(1+t_{i}^{b})^{-a^{*}}-(1+t_{i}^{b})^{-a^{*}}].\end{array}  \tag{37}\]  
Through the derivation combined (32), (33), (35), and (37), in the process of Gibbs sampling, we just need to set \(\lambda=0\) in the procedure. It is worth noting that when using a non-informative prior, the method for generating initial values of \(a_{1}\), \(a_{2}\), and \(\lambda\) needs to be adjusted in the first step of Gibbs sampling to avoid generating 0. We then derive the Bayesian estimates and the associated HPD credible intervals by Gibbs sampling. The sampling process consists of \(N=\) 11,000 iterations, during which the first \(d=1000\) iterations are seen as the burn-in period and are discarded. After repeating the process above for \(M=\) 10,000 times, we can compare the performance of all point estimators numerically via bias and mean squared error (MSE). In addition, the outcomes of interval estimators can be assessed by means of average width (AW) and coverage percentages (CP). In this study, the code was executed in R langurange (version 4.3.2) on a Lenovo XiaoXinPro 16 ACH laptop. It is equipped with an 8-core AMD Ryzen 7 5800H with Radeon Graphics processor running at a base frequency, 16 GB of RAM with a speed of 3200 MT/s, and a 512 GB solid-state drive. For EM algorithm, taking into account both the time cost and accuracy, we need to select a proper value for \(\epsilon\). We run one of the schemes and compare the estimation accuracy (according to the MSEs) under different \(\epsilon\) values to identify an appropriate value. Then, we execute the estimation for the other schemes using this determined \(\epsilon\) value.