Bayesian Analysis Using Joint Progressive Type-II Censoring Scheme

The Burr-XII distribution, originally introduced by Burr [1], has found extensive applications in various domains, including lifetime modeling for reliability analysis, addressing life-testing challenges, and devising acceptance sampling plans, as exemplified in the works of Abbasi et al. [2] and other researchers. It has also been effectively utilized in the analysis of observational data across diverse fields such as meteorology, finance, and hydrology, as showcased in studies conducted by Chen et al. [3], Ali and Jaheen [4], Burr [1] and Lio et al. [5]. Moreover, Shao et al. [6] delved into the modeling of extreme events using the three-parameter Burr-XII distribution (TPBXIID), notably in the context of flood-frequency analysis.

Our decision to employ the TPBXIID stems from its remarkable adaptability, spanning a wide spectrum of shapes, from highly skewed to nearly symmetric. This versatility renders it a valuable model for datasets that do not adhere to standard shapes. Notably, the distribution’s three parameters, denoted as $\alpha$, $\gamma$, and $\theta$, offer straightforward interpretations, simplifying the analysis of statistical outcomes and enabling comparisons across different datasets. As a distribution limited to non-negative values, the Burr-XII distribution is frequently harnessed for modeling data related to lifetimes, sizes, or quantities. It consistently demonstrates an excellent fit to empirical datasets and, in specific data types, is known to outperform other commonly employed distributions, including the Weibull distribution.

Furthermore, the Burr-XII distribution serves as a generalized form encompassing several other distributions, including the Lomax (Pareto II), Burr-XII, and log-logistic distributions. In summary, the three-parameter Burr-XII distribution emerges as a versatile tool for modeling non-negative data, which are characterized by a broad spectrum of shapes. Its flexibility and interpretability render it a popular choice in statistical modeling.

Cook and Johnson [7] applied the Burr model to attain superior fits for a uranium survey dataset, while Zimmer et al. [8] delved into the statistical and probabilistic properties of the Burr-XII distribution and its relationships with other distributions commonly employed in reliability analyses. Additionally, Tadikamalla [9] extended the two-parameter Burr-XII distribution by introducing an additional scale parameter, resulting in the TPBXIID. This extension has sparked increased interest in the applications of the Burr-XII distribution.

Tadikamalla also established mathematical connections among Burr-related distributions, revealing that the Lomax distribution constitutes a special instance of the Burr-XII distribution, and the compound Weibull distribution represents a generalization of the Burr distribution. Furthermore, he demonstrated that several widely used distributions, including Weibull, logistic, log-logistic, normal, and lognormal distributions, can be viewed as specific cases of the Burr-XII distribution by appropriately configuring the distribution parameters.

In essence, the TPBXIID proves to be highly adaptable, encompassing two shape parameters and one scale parameter in the distribution function, thereby allowing it to represent a diverse range of distribution shapes. The TPBXIID can be characterized through its cumulative distribution function (CDF) and probability density function (PDF), as expressed, respectively, in Equations (1) and (2). where $\gamma$ and $\theta$ represent the shape parameters, and $\alpha$ serves as the scale parameter. Notably, when $\theta >1$, the density function exhibits an upside-down bathtub shape (unimodal) with the mode located at $x=\alpha {\left[(\theta -1)/(\theta \gamma +1)\right]}^{\frac{1}{\theta }}$, while it assumes an L-shaped form when $\theta =1$.

Recently, Mead and Afify [10] ventured into defining and examining the properties and applications of the five-parameter Burr-XII distribution, which is referred to as the Kumaraswamy exponentiated Burr-XII. Moreover, Shafqat et al. [11] explored the utilization of moving average control charts in the context of Burr X and inverse Gaussian, and Aslam et al. [12] studied a new generalized Burr-XII distribution with real-life applications.

The joint censoring approach proves to be a valuable and practical method for comparing life tests of products originating from various units within the same facility. Consider a scenario where two production lines operate within the same facility, generating products. In this setup, two independent samples of sizes m and n can be selected from each production line and subjected to simultaneous life-testing experiments. To optimize resource utilization, reduce costs, and save time, researchers often employ a joint progressive Type-II censoring scheme (JP-II-CS). This approach is instrumental in terminating life testing when a predetermined number of failures are observed.

Numerous studies in the literature have explored JP-II-CS and inference methods associated with it. For example, Rasouli and Balakrishnan [13] introduced likelihood inference techniques for two exponential distributions based on JP-II-CS. Doostparast et al. [14] delved into Bayes estimation under the linear exponential loss function using JP-II-CS data. Balakrishnan et al. [15] provided likelihood inference procedures for k exponential distributions under JP-II-CS, while Mondal and Kundu [16] focused on the point and interval estimation of Weibull parameters within the context of JP-II-CS.

Goel and Krishna [17] explored likelihood and Bayesian inference for k Lindley populations under a joint Type-II censoring scheme. Krishna and Goel [18] conducted a study on Lindley populations utilizing JP-II-CS. Additionally, Goel and Krishna [19] discussed statistical inference for two Lindley populations under a balanced JP-II-CS. Bayoud and Raqab [20] investigated classical and Bayesian inferences for two Topp–Leone models under JP-II-CS, while Chen and Gui [21] addressed the statistical inference of the generalized inverted exponential distribution in the context of JP-II-CS.

Pandey and Srivastava [22] focused on Bayesian inference for two log-logistic populations under JP-II-CS, and Qiao and Gui [23] tackled the statistical inference of the Weighted Exponential Distribution under similar censoring conditions.

Recently, Hassan et al. [24] delved into the statistical inference of the Burr Type III distribution under joint progressively Type II censoring. Kumar and Kumari [25] explored Bayesian and likelihood estimation techniques for two inverse Pareto populations under joint progressive censoring conditions.

According to Rasouli and Balakrishnan [13], JP-II-CS is described as follows. Let $X_1\dots ,X_m$, be lifetimes of m units for product A, and they are supposed to be independent and identically distributed (iid) random variables from TPBXIID with a CDF given by and the PDF is In a similar manner, consider a set of n lifetimes denoted as $Y_1,Y_2,\dots ,Y_n$ for the product B. These lifetimes correspond to n units and are treated as iid random variables following the TPBXIID. The CDF for this distribution is given by: and PDF is Where ${\theta }_1$, ${\gamma }_1$, ${\theta }_2$ and ${\gamma }_2$ are shape parameters and ${\alpha }_1$ and ${\alpha }_2$ are scale parameters. In this scenario, let $K=m+n$ denote the total sample size and ${\lambda }_1\le \dots \le {\lambda }_K$ indicate the order statistics of the K random variables $\{X_1,X_2,\dots ,X_m;Y_1,Y_2,\dots ,Y_n\}$. The JP-II-CS method is applied as follows: when the first failure occurs, $R_1$ units are randomly removed from the remaining $K-1$ surviving units. The same process is repeated for the second failure, where $R_2$ units are randomly withdrawn from the remaining $K-R_1-2$ surviving units, and so on. At the $rth$ failure, all remaining $R_r=K-r-{\sum }_{i=1}^{r-1}{R}_i$ surviving units are withdrawn from the experiment. The JP-II-CS is represented by $R=(R_1,R_2,\dots ,R_r)$, and the total number of failures r is predetermined before conducting the experiment. Suppose that $R_i=S_i+T_i$, $i=1,\dots ,r,$ where $S_i$ and $T_i$ represent the number of units withdrawn at the time of $ith$ failure related to X and Y samples, respectively. These values are unknown and random variables. The data observed in this form will consist of $(H,\lambda ,S)$, where $H=(H_1,H_2,\dots ,H_r)$, where $H_i=1$ or 0 if ${\lambda }_i$ comes from X or Y failure, respectively, $\lambda =({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_r)$ with $r<K$ and $S=(S_1,S_2,\dots ,S_r)$.

In this study, we employ a JP-II-CS strategy to formulate statistical inferences and assess two independent samples from the TPBXIID. We derive point and interval estimators using Bayesian and maximum likelihood estimation (MLE) techniques. Subsequently, we calculate asymptotic confidence intervals (ACI) based on the observed information matrix. These confidence intervals (CIs) are computed through various bootstrap techniques, including Bootstrap-P (Boot-P), Bootstrap-T (Boot-T), Bias-Corrected Bootstrap (Boot-BC), and Bias-Corrected Accelerated Bootstrap (Boot-BCa) methods. We assume a gamma prior distribution for both the shape and scale parameters. Employing the Metropolis–Hastings (M-H) method, we obtain Bayes estimates and credible intervals (CRIs) for the informative prior under both the squared error (SE) and linear exponential (LINEX) loss functions. To evaluate the effectiveness of these diverse approaches, we conduct Monte Carlo simulations and analyze real-world data.

The paper is structured in the following way: Section 2 outlines the derivation of the MLEs for the unknown parameters of TPBXIID. In Section 3, there is a presentation of ACIs that depend on the MLEs. Section 4 discusses different bootstrap CIs. The Bayesian analysis is performed in Section 5. To illustrate the estimation methods developed in this paper, we analyze real datasets in Section 6. The results of simulation are presented in Section 7. Finally, a brief conclusion can be found in Section 8.

Based on the work of Rasouli and Balakrishnan [13], the likelihood function for $(H,\lambda ,S)$ can be expressed as follows: where $L(.)L({\alpha }_1,{\alpha }_2,{\theta }_1,{\theta }_2,{\gamma }_1,{\gamma }_2;H,\lambda ,S)$, $w_1\le w_2\le \dots \le w_r$,   $\overline{F}=1-F$,   $\overline{Q}=1-Q$,   ${\sum }_{i=1}^r{s}_i=m-m_r$,   ${\sum }_{i=1}^r{t}_i=n-n_r$,   ${\sum }_{i=1}^r{s}_i+{\sum }_{i=1}^r{t}_i={\sum }_{i=1}^r{R}_i$ and $C=D_1{D}_2$ with The likelihood function in (7) becomes The log-likelihood function for the TPBXIID, as related to Equation (8), is as follows: Deriving the first derivative of Equation (9) concerning ${\alpha }_1,{\alpha }_2,{\theta }_1,{\theta }_2,{\gamma }_1,\mathrm{and},{\gamma }_2$ and setting each one as follows: From (14) and (15), we obtain the MLE $\widehat{\gamma }$ as To acquire the MLEs for ${\alpha }_1,{\alpha }_2,{\theta }_1,\mathrm{and},{\theta }_2$, the resolution of Equations (10)–(13) is necessary. Nevertheless, attaining analytical solutions for these equations proves highly demanding, making it arduous to obtain closed-form expressions for each parameter. Consequently, resorting to a numerical approach, such as the Newton–Raphson iteration method, becomes imperative to estimate the parameter values.

The second derivative of log-likelihood function with respect to ${\alpha }_1,{\alpha }_2,{\theta }_1,\mathrm{and},{\theta }_2$, gives and where

The asymptotic variances–covariances of the MLEs for parameters ${\alpha }_1,{\alpha }_2,{\theta }_1,{\theta }_2,{\gamma }_1, \mathrm{and} {\gamma }_2$ are given by elements of the inverse of the Fisher information matrix (FIM) where the FIM is obtained by taking the expectation of minus Equations (18)–(29) and defined as where $i,j=1,2,3,4,5,6 \mathrm{and} ({\delta }_1,{\delta }_2,{\delta }_3,{\delta }_4,{\delta }_5,{\delta }_6)=({\alpha }_1,{\alpha }_2,{\theta }_1,{\theta }_2,{\gamma }_1\mathrm{and}{\gamma }_2)$.

The mathematical expressions for the expectations referred to in Equation (31) lack precision and can be located in Cohen’s work [26], where he computed the asymptotic variance–covariance matrices for various sample types.

Therefore, we rely on the estimated asymptotic variance–covariance matrix for the MLEs. hence, From the likelihood function in (9), we have $I_{i,j}\left({\delta }_i\right)=0$ if $i\ne j$. Consequently, we have The asymptotic normality of the MLEs can be used to compute the ACIs for parameter ${\delta }_i$, so that we can express the approximate $(1-\eta )100\%$ ACIs for ${\delta }_i$ as where $Z_{\eta /2}$ is the percentile of the standard normal distribution with right-tail probability $\eta /2$.

Bootstrap confidence intervals are a powerful statistical tool used to estimate the uncertainty associated with a parameter or statistic of interest in data analysis. They are particularly valuable when the underlying probability distribution of the data is complex or unknown. The bootstrap technique involves repeatedly resampling the observed data, with replacement, to create a large number of “bootstrap samples.” From these samples, new estimates of the parameter of interest are calculated. By analyzing the distribution of these estimates, bootstrap confidence intervals provide a range within which the true parameter value is likely to fall. This approach is widely used in various fields, such as finance, biology, and machine learning, to gain insight into the robustness and stability of statistical estimates, making it an indispensable tool for researchers and analysts seeking to make informed decisions in the presence of uncertainty.

In this context, we have the Boot-P, introduced by Efron [27], and the Boot-T proposed by Hall [28], along with two variations, Boot-BC and Boot-BCa, which are rooted in the concept presented by DiCiccio and Efron [29]. In 2017, Ghazal and Hasaballah [30] explored bootstrap and MCMC methods based on unified hybrid censored data. Below, we will introduce these four methods.

In this section, we employ the MCMC approach within a Gibbs sampler framework that incorporates a nested M-H algorithm. We utilize this approach to generate parametric samples representing the unknown parameters ${\alpha }_1$, ${\alpha }_2$, ${\theta }_1$, ${\theta }_2$, ${\gamma }_1$ and ${\gamma }_2$ from their respective marginal posteriors, allowing us to derive Bayes estimates for these parameters. The Gibbs sampler is a recursive sampling technique employed to simulate samples from the full conditional posterior distributions, while the M-H algorithm is used to generate samples from arbitrary distributions (Hastings [32]; Metropolis et al. [33]). In this context, we generate N samples using the MCMC technique, with the initial M values discarded during the burn-in period. The remaining N–M sample values are subsequently utilized for further Bayesian analysis.

We assume that the parameters ${\alpha }_1$, ${\alpha }_2$, ${\theta }_1$, ${\theta }_2$, ${\gamma }_1$ and ${\gamma }_2$ have independent gamma prior distributions as where $i=1,2,3,4,5,6 \mathrm{and} ({\delta }_1,{\delta }_2,{\delta }_3,{\delta }_4,{\delta }_5,{\delta }_6)=({\alpha }_1,{\alpha }_2,{\theta }_1,{\theta }_2,{\gamma }_1,{\gamma }_2).$

Where the hyperparameters $a_i$ and $b_i$, $i=1,2,3,4,5,6$ are supposed to be known and non-negative. By using the prior distribution for ${\alpha }_1$, ${\alpha }_2$, ${\theta }_1$, ${\theta }_2$, ${\gamma }_1$ and ${\gamma }_2$, we obtain the joint prior distribution as follows: Based on (8) and (42), the joint posterior density function of ${\alpha }_1$, ${\alpha }_2$, ${\theta }_1$, ${\theta }_2$, ${\gamma }_1$ and ${\gamma }_2$ given $(H,\lambda ,S)$ is written as We observed that (43) is not amenable to analytical solutions due to the formidable challenge of deriving closed-form expressions for the marginal posterior distributions of individual parameters. Consequently, we recommend the utilization of the MCMC method to approximate (43). Several studies have extensively explored the MCMC technique, including works by Chen and Shao [34] and Ghazal and Hasaballah [30,35,36]. From (43), the conditional posterior density function of ${\alpha }_1$, ${\alpha }_2$, ${\theta }_1$, ${\theta }_2$, ${\gamma }_1$ and ${\gamma }_2$ can be obtained as the following proportionality: to simplify, we used ${\pi }_1^{*}\left({\alpha }_1\right)$, ${\pi }_2^{*}\left({\alpha }_2\right)$, ${\pi }_3^{*}\left({\theta }_1\right)$, ${\pi }_4^{*}\left({\theta }_2\right)$, ${\pi }_5^{*}\left({\gamma }_1\right)$ and ${\pi }_6^{*}\left({\gamma }_2\right)$ instead of ${\pi }_1^{*}\left({\alpha }_1\right|{\alpha }_2,{\theta }_1,{\theta }_2,{\gamma }_1,{\gamma }_2,H,\lambda ,S)$, ${\pi }_2^{*}\left({\alpha }_2\right|{\alpha }_1,{\theta }_1,{\theta }_2,{\gamma }_1,{\gamma }_2,H,\lambda ,S)$, ${\pi }_3^{*}\left({\theta }_1\right|{\alpha }_1,{\alpha }_2,{\theta }_2,{\gamma }_1,{\gamma }_2,H,\lambda ,S)$, ${\pi }_4^{*}\left({\theta }_2\right|{\alpha }_1,{\alpha }_2,{\theta }_1,{\gamma }_1,{\gamma }_2,H,\lambda ,S)$, ${\pi }_5^{*}\left({\gamma }_1\right|{\alpha }_1,{\alpha }_2,{\theta }_1,{\theta }_2,{\gamma }_2,H,\lambda ,S)$ and ${\pi }_6^{*}\left({\gamma }_2\right|{\alpha }_1,{\alpha }_2,{\theta }_1,{\theta }_2,{\gamma }_1,H,\lambda ,S),$ respectively: It is evident that the full conditional posterior density function of ${\gamma }_1$, as provided in Equation (48), takes the form of a gamma density with a shape parameter of $(a_5+m_r)$ and a scale parameter of $\left(a_5+m_r,b_5+{\sum }_{i=1}^r(h_i+s_i)ln\left[1+{\left(\frac{{\lambda }_i}{{\alpha }_1}\right)}^{{\theta }_1}\right]\right)$. Similarly, the full conditional posterior density function of ${\gamma }_2$, as shown in Equation (49), follows a gamma distribution with a shape parameter of $(a_6+n_r)$ and a scale parameter of $\left(a_6+n_r,b_6+{\sum }_{i=1}^r(t_i-h_i+1)ln\left[1+{\left(\frac{{\lambda }_i}{{\alpha }_2}\right)}^{{\theta }_2}\right]\right)$. Consequently, samples of ${\gamma }_1$ and ${\gamma }_2$ can be readily generated using any gamma distribution generation method.