Yi Zhang
School of Mathematical Sciences, Guizhou Normal University, Guiyang, China.
DOI: 10.4236/ojs.2026.164012   PDF    HTML   XML   19 Downloads   118 Views   Citations

Abstract

Bivariate survival data with dependent event times and censored observations are common in medical research. Traditional survival analysis methods cannot effectively capture such dependence structures. In this paper, we propose bivariate Clayton Weibull (BCW) and bivariate Clayton exponential (BCE) models based on the Clayton copula under Type II censoring. Parameter estimation is performed using maximum likelihood estimation and Bayesian MCMC methods. Monte Carlo simulations evaluate point estimation accuracy, interval estimation, and dependence parameter estimation under various censoring levels. Results show that bias and mean squared error decrease as effective sample size increases. Under the same censoring level, Bayesian estimation outperforms maximum likelihood estimation. The models are applied to catheter infection data of kidney disease patients, demonstrating their validity and practicality on real medical data.

Keywords

Clayton Copula, Weibull Distribution, Type-II Censoring, MLE, Bayesian Estimation

Share and Cite:

1. Introduction

In clinical research, disease progression often involves multiple correlated event times. For kidney disease patients, the times of first and second infections after catheter implantation not only reflect disease recurrence patterns but may also be dependent [1]. Similarly, in organ transplantation and tumor recurrence, joint analysis of multiple event times is crucial for understanding disease evolution and guiding treatment [2]. However, such bivariate survival data face two challenges: potential non-independent dependence between event times [3], and the presence of censored observations due to limited follow-up [4].

Traditional methods for bivariate survival data, such as frailty models or multivariate risk regression, often impose strong parametric assumptions on the dependence structure. Copula functions offer greater flexibility by separating marginal distributions from the dependence structure. Nelsen systematically presented the mathematical properties and statistical applications of Copulas in his classic work [5]. Clayton first introduced Copula concepts into bivariate life table analysis [6], with his model capturing lower-tail dependence—making it particularly suitable for analyzing associations between early events.

The success of Copula methods in survival analysis has driven their extension to more complex censored data. Under Type II censoring, the trial terminates upon observing a predetermined r-th event—a design widely used in reliability testing and medical follow-up [7]. To support practical applications, Sun and Ding developed the CopulaCenR package, which enables Copula-based regression for bivariate censored data, supporting Clayton, Gumbel, Frank copulas, and Weibull, log-logistic marginal models [8].

The Weibull distribution is widely used in survival analysis because its shape parameter can flexibly fit monotonically increasing, decreasing, or constant hazard rate functions. The exponential distribution, as a special case of the Weibull distribution (with the shape parameter equal to 1), assumes a constant hazard rate and has a more parsimonious model structure [9]. In recent years, El-Sherpieny et al. [10] investigated parameter estimation for bivariate Weibull distributions based on different Copula functions under progressive Type II censoring. Wang and Yan [11] explored reliability and dependence parameter estimation for a three-variable stress-strength model based on the Clayton Copula under progressive Type II censoring, employing method of moments, maximum likelihood estimation, and Bayesian methods for comparative analysis. Although the properties of these two marginal distributions under univariate censored data have been thoroughly studied, a systematic comparison between the bivariate Weibull model and the bivariate exponential model linked by the Clayton Copula under the Type II censoring framework has not yet been reported in the literature. The differences between the two models in terms of parameter estimation accuracy, interval estimation reliability, and capability to characterize dependence structures remain unclear.

Parametric and semi-parametric models, particularly the Cox proportional hazards model, have long been the cornerstone of survival analysis, offering interpretable estimation of covariate effects without baseline hazard assumptions [12]. Recent advances in machine learning, including random survival forests and DeepSurv, have demonstrated superior predictive performance in large clinical datasets by adaptively capturing complex nonlinearities and interactions [13]. However, many machine learning methods struggle with small samples or high censoring rates [14], and semi-parametric approaches lack a direct framework for modeling dependence structures among variables. By contrast, copula-based parametric models explicitly characterize event-time dependence, jointly estimate marginal and dependence parameters, and provide full uncertainty quantification [15]—advantages that are especially critical in small-sample or heavily censored medical studies. Thus, rather than competing, the proposed method complements existing approaches: machine learning excels in prediction with abundant data, while fully parametric copula models offer a principled framework for dependence quantification and interpretable inference in small-sample settings, with each approach occupying a distinct niche under different data regimes and analytical goals.

This paper is situated within the Type-II censoring scheme, where the experiment is terminated upon the occurrence of a predetermined number r of events. This censoring mechanism is particularly valuable in reliability experiments and certain clinical follow-up studies, as it ensures a fixed number of complete observations, thus providing adequate information for parameter estimation. Although the proposed estimation framework for the two Clayton copula-based models is developed under Type-II censoring, it can be straightforwardly generalized to Type-I or random right-censoring settings by adjusting the censoring terms in the likelihood function.

To address this gap, we construct bivariate Clayton Weibull (BCW) and bivariate Clayton exponential (BCE) models under Type II censoring, using both maximum likelihood estimation and Bayesian MCMC methods. Monte Carlo simulations evaluate point estimation accuracy, interval reliability, and Copula dependence estimation under various censoring levels, comparing Weibull and exponential marginal distributions for Type II censored data. The proposed models are then applied to the kidney catheter infection data of McGilchrist and Aisbett to demonstrate their validity and practicality on real medical data, offering a methodological reference for analyzing dependent survival data.

2. Copula Theory

This section delineates the essential components of bivariate copula theory: 1) the definition of two dependence measures, Spearman’s $\rho$ and Kendall’s $\tau$ ; 2) the theoretical framework of the Clayton copula function; and 3) Sklar’s theorem, which underpins the entire copula theory.

Let $X_1$ and $X_2$ be continuous random variables with marginal distribution functions $H_1\left(x_1\right)$ and $H_2\left(x_2\right)$ , respectively, and let their joint distribution be given by the copula function $C\left(u_1,u_2\right)$ , where $u_i=H_i\left(x_i\right)$ for $i=1,2$ . The population Kendall’s tau, denoted $\tau$ , is defined as the probability of concordance minus the probability of discordance, and can be expressed in terms of the copula as:

$\tau =4{\displaystyle {\int }_0^1{\displaystyle {\int }_0^{1}C\left(u_1,u_2\right)\text{d}C\left(u_1,u_2\right)}}-1=4{\displaystyle {\int }_0^1{\displaystyle {\int }_0^{1}C\left(u_1,u_2\right)c\left(u_1,u_2\right)\text{d}{u}_1\text{d}{u}_2}}-1,$

where $c\left(u_1,u_2\right)={\partial }^{2}C\left(u_1,u_2\right)/\partial u_1\partial u_2$ is the copula density function. The population Spearman’s rho, denoted $\rho$ , is defined as the correlation between the probability-transformed variables and is given by:

$\rho =12{\displaystyle {\int }_0^1{\displaystyle {\int }_0^{1}C\left(u_1,u_2\right)\text{d}{u}_1\text{d}{u}_2}}-3.$

The Archimedean copula family is a class of copulas characterized by an explicit generator, among which the Clayton copula is one of the most widely used. Proposed by Clayton [6] in 1978, the Clayton copula has the generator:

$\psi \left(t\right)=\frac{1}{\eta }\left(t^{-\eta }-1\right)$ . For the association parameter $\eta \in \left(0,\infty \right)$ , the cumulative distribution function (cdf) of the bivariate Clayton copula is given by:

$C\left(u_1,u_2\right)={\left(u_1^{-\eta }+u_2^{-\eta }-1\right)}^{-\frac{1}{\eta }},\eta \in \left(0,\infty \right).$

In the limit $\eta \to 0^{+}$ , $C\left(u_1,u_2\right)\to u_1{u}_2$ , corresponding to independence; as $\eta \to \infty$ , the variables become perfectly positively dependent, and the copula approaches the Fréchet-Hoeffding upper bound. Thus, the Clayton copula is characterized by lower-tail dependence, indicating stronger association in the lower tail of the joint distribution.

The probability density function (pdf) of the Clayton copula is given by:

$c\left(u_1,u_2\right)=\left(1+\eta \right){\left(u_1{u}_2\right)}^{-1-\eta }{\left(u_1^{-\eta }+u_2^{-\eta }-1\right)}^{-2-\frac{1}{\eta }},\eta \in \left(0,\infty \right).$

The relationship between $\tau$ and the copula parameter $\eta$ is given by: $\tau =\frac{\eta }{\eta +2}.$ Thus, when $\eta =0$ , we have $\tau =0$ , corresponding to independence

between the variables; when $\eta \to \infty$ , we have $\tau \to 1$ , corresponding to perfect positive dependence.

Sklar’s theorem, proposed by Sklar in 1959 [16], is the foundational core of copula theory. The theorem states that for any random variables $X_1$ and $X_2$ with marginal cumulative distribution functions $H_1\left(x_1\right)$ and $H_2\left(x_2\right)$ , and marginal probability density functions $h_1\left(x_1\right)$ and $h_2\left(x_2\right)$ , there exists a unique copula $C:{\left[0,1\right]}^2\to \left[0,1\right]$ such that their joint distribution function can be written as:

$H\left(x_1,x_2\right)=C\left(H_1\left(x_1\right),H_2\left(x_2\right)\right),$

where $\eta$ is the dependence parameter that quantifies the dependence structure between the variables. If the marginal distributions are continuous, the joint probability density function is given by:

$h\left(x_1,x_2\right)=h_1\left(x_1\right)h_2\left(x_2\right)c\left(H_1\left(x_1\right),H_2\left(x_2\right)\right),$

where $c(\cdot )$ denotes the copula density function.

3. Bivariate Distributions

This section introduces two bivariate distributions constructed by linking the exponential distribution and Weibull distribution with the Clayton copula function separately.

3.1. Bivariate Clayton Exponential Distribution (BCE)

The cumulative distribution function (cdf) and probability density function (pdf) of the exponential distribution are, respectively:

$H\left(x;\lambda \right)=1-{\text{e}}^{-\lambda x},\lambda >0,x>0,$

$h\left(x;\lambda \right)=\lambda {\text{e}}^{-\lambda x},\lambda >0,x>0.$

Its reliability function and hazard rate function are, respectively:

$S\left(x;\lambda \right)={\text{e}}^{-\lambda x},\lambda >0,x>0,$

$r\left(x;\lambda \right)=\lambda ,\lambda >0,x>0.$

The joint pdf and joint cdf of the bivariate Clayton exponential (BCE) model are given, respectively, by:

$\begin{array}{c}{h}_{BCE}\left(x_1,x_2\right)={\lambda }_1{\text{e}}^{-{\lambda }_1{x}_1}{\lambda }_2{\text{e}}^{-{\lambda }_2{x}_2}\left(1+\eta \right){\left(u_1\left(x_1\right)u_2\left(x_2\right)\right)}^{-\eta -1}\\ \times {\left(u_1{\left(x_1\right)}^{-\eta }+u_2{\left(x_2\right)}^{-\eta }-1\right)}^{-2-\frac{1}{\eta }},\end{array}$

$H_{BCE}\left(x_1,x_2\right)={\left(u_1{\left(x_1\right)}^{-\eta }+u_2{\left(x_2\right)}^{-\eta }-1\right)}^{-\frac{1}{\eta }},$

where $u_1\left(x_1\right)=1-{\text{e}}^{-{\lambda }_1{x}_1}$ , $u_2\left(x_2\right)=1-{\text{e}}^{-{\lambda }_2{x}_2}$ , with ${\lambda }_1,{\lambda }_2,\eta \in \left(0,\infty \right)$ .

The respective reliability functions of the marginal distributions are:

$S_E\left(x_1;{\lambda }_1\right)={\text{e}}^{-{\lambda }_1{x}_1},{\lambda }_1>0,x_1>0,$

$S_E\left(x_2;{\lambda }_2\right)={\text{e}}^{-{\lambda }_2{x}_2},{\lambda }_2>0,x_2>0.$

The joint reliability function of the bivariate copula is: $S\left(x_1,x_2\right)=C\left(S\left(x_1\right),S\left(x_2\right)\right)$ . Therefore, the reliability function of the BCE distribution is given by:

$S_{BCE}\left(x_1,x_2\right)={\left[{\left(1-u_1\left(x_1\right)\right)}^{-\eta }+{\left(1-u_2\left(x_2\right)\right)}^{-\eta }-1\right]}^{-\frac{1}{\eta }},{\lambda }_1,{\lambda }_2>0,\eta \in \left(0,\infty \right).$

The expression of the bivariate failure rate function is $r\left(x_1,x_2\right)=\frac{h\left(x_1,x_2\right)}{S\left(x_1,x_2\right)}$ . [17] Consequently, the hazard function of the BCE distribution is given by:

$\begin{array}{l}{r}_{BCE}\left(x_1,x_2\right)\\ =\frac{{\lambda }_1{\text{e}}^{-{\lambda }_1{x}_1}{\lambda }_2{\text{e}}^{-{\lambda }_2{x}_2}\left(1+\eta \right){\left(u_1\left(x_1\right)u_2\left(x_2\right)\right)}^{-1-\eta }{\left(u_1{\left(x_1\right)}^{-\eta }+u_2{\left(x_2\right)}^{-\eta }-1\right)}^{-2-\frac{1}{\eta }}}{{\left[{\left(1-u_1\left(x_1\right)\right)}^{-\eta }+{\left(1-u_2\left(x_2\right)\right)}^{-\eta }-1\right]}^{-\frac{1}{\eta }}},\end{array}$

with ${\lambda }_1,{\lambda }_2>0$ and $\eta \in \left(0,\infty \right)$ .

Figure 1 and Figure 2 show the joint density and joint hazard surfaces of the BCE model under different parameter combinations. As the Copula parameter $\eta$ increases, the lower-tail dependence strengthens, and the peaks of both surfaces concentrate near the origin, exhibiting a sharp peaking pattern. This reflects the amplifying effect of lower-tail dependence on early system failure risk.

Figure 1. Joint density plots of BCE under some parameter values.

Figure 2. Joint hazard rate plots of BCE under some parameter values.

3.2. Bivariate Clayton Weibull Distribution (BCW)

The Weibull distribution is a continuous distribution defined by Fréchet [18] in 1927, and the cdf and the pdf of it are respectively presented as:

$H\left(x;\lambda ,\gamma \right)=1-{\text{e}}^{-{\left(\lambda x\right)}^{\gamma }},\lambda ,\gamma >0,x\ge 0,$

$h\left(x;\lambda ,\gamma \right)=\lambda \gamma {\left(\lambda x\right)}^{\gamma -1}{\text{e}}^{-{\left(\lambda x\right)}^{\gamma }},\lambda ,\gamma >0,x\ge 0,$

where $\gamma$ is the shape parameter and $1/\lambda$ is the scale parameter. When $\gamma =1$ , the distribution becomes an exponential distribution.

Individually, the reliability function and the hazard rate function of it are:

$S\left(x;\lambda ,\gamma \right)={\text{e}}^{-{\left(\lambda x\right)}^{\gamma }},\lambda ,\gamma >0,x\ge 0,$

$r\left(x;\lambda ,\gamma \right)=\lambda \gamma {\left(\lambda x\right)}^{\gamma -1},\lambda ,\gamma >0,x\ge 0.$

Here are the joint pdf and joint cdf of the BCW distribution:

$\begin{array}{c}{h}_{BCW}\left(x_1,x_2\right)={\lambda }_1{\gamma }_1{\left({\lambda }_1{x}_1\right)}^{{\gamma }_1-1}{\text{e}}^{-{\left({\lambda }_1{x}_1\right)}^{{\gamma }_1}}{\lambda }_2{\gamma }_2{\left({\lambda }_2{x}_2\right)}^{{\gamma }_2-1}{\text{e}}^{-{\left({\lambda }_2{x}_2\right)}^{{\gamma }_2}}\\ \times \left(1+\eta \right){\left(u_1\left(x_1\right)u_2\left(x_2\right)\right)}^{-1-\eta }{\left(u_1{\left(x_1\right)}^{-\eta }+u_2{\left(x_2\right)}^{-\eta }-1\right)}^{-2-\frac{1}{\eta }},\end{array}$

$H_{BCW}\left(x_1,x_2\right)={\left[u_1{\left(x_1\right)}^{-\eta }+u_2{\left(x_2\right)}^{-\eta }-1\right]}^{-\frac{1}{\eta }},{\lambda }_1,{\gamma }_1,{\lambda }_2,{\gamma }_2>0,\eta \in \left(0,\infty \right),$

where $u_1\left(x_1\right)=1-{\text{e}}^{-{\left({\lambda }_1{x}_1\right)}^{{\gamma }_1}}$ , $u_2\left(x_2\right)=1-{\text{e}}^{-{\left({\lambda }_2{x}_2\right)}^{{\gamma }_2}}$ .

Separately, the reliability functions belonging to the marginal distributions are:

$S_W\left(x_1;{\lambda }_1,{\gamma }_1\right)={\text{e}}^{-{\left({\lambda }_1{x}_1\right)}^{{\gamma }_1}},{\lambda }_1,{\gamma }_1>0,x_1\ge 0,$

$S_W\left(x_2;{\lambda }_2,{\gamma }_2\right)={\text{e}}^{-{\left({\lambda }_2{x}_2\right)}^{{\gamma }_2}},{\lambda }_2,{\gamma }_2>0,x_2\ge 0.$

Then the reliability function of the BCW distribution is given by the following formula:

$S_{BCW}\left(x_1,x_2\right)={\left[{\left(1-u_1\left(x_1\right)\right)}^{-\eta }+{\left(1-u_2\left(x_2\right)\right)}^{-\eta }-1\right]}^{-\frac{1}{\eta }},{\lambda }_1,{\gamma }_1,{\lambda }_2,{\gamma }_2>0,\eta \in \left(0,\infty \right).$

The expression of the hazard function of the BCW distribution is as follows:

$\begin{array}{c}{r}_{BCW}\left(x_1,x_2\right)={\lambda }_1{\gamma }_1{\left({\lambda }_1{x}_1\right)}^{{\gamma }_1-1}{\text{e}}^{-{\left({\lambda }_1{x}_1\right)}^{{\gamma }_1}}{\lambda }_2{\gamma }_2{\left({\lambda }_2{x}_2\right)}^{{\gamma }_2-1}{\text{e}}^{-{\left({\lambda }_2{x}_2\right)}^{{\gamma }_2}}\\ \times \frac{\left(1+\eta \right){\left(u_1\left(x_1\right)u_2\left(x_2\right)\right)}^{-1-\eta }{\left(u_1{\left(x_1\right)}^{-\eta }+u_2{\left(x_2\right)}^{-\eta }-1\right)}^{-2-\frac{1}{\eta }}}{{\left[{\left(1-u_1\left(x_1\right)\right)}^{-\eta }+{\left(1-u_2\left(x_2\right)\right)}^{-\eta }-1\right]}^{-\frac{1}{\eta }}},\end{array}$

where ${\lambda }_1,{\gamma }_1,{\lambda }_2,{\gamma }_2>0$ , $\eta \in \left(0,\infty \right)$ .

Figure 3 and Figure 4 illustrate the joint probability density and joint hazard rate surfaces of the BCW model under various parameter settings. It can be observed that as $\eta$ increases, the lower-tail dependence between variables intensifies, leading to a more pronounced concentration of the joint density and hazard rate near the origin, with the surface shape evolving from smooth to sharp-peaked.

Figure 3. Joint density plots of BCW under some parameter values.

Figure 4. Joint hazard rate plots of BCW under some parameter values.

4. Inferential Analysis for Uncensored and Censored Samples

This section presents the maximum likelihood estimation (MLE) and Bayesian parameter estimation for the bivariate models. To facilitate clarity, the BCW model is used as an example for the mathematical derivations; the estimation for the BCE model follows similarly.

4.1. MLE with Full Samples

MLE is a classical approach to parameter estimation, which seeks the parameter values that maximize the probability of the observed data. Under complete sampling, let $\left(x_{1i},x_{2i}\right)$ , for $i=1,\cdots ,n$ , be independent and identically distributed observations from a bivariate distribution $H_{X_1,X_2}\left(x_1,x_2\right)$ , with $\theta =\left({\lambda }_1,{\gamma }_1,{\lambda }_2,{\gamma }_2,\eta \right)$ denoting the parameter vector. The likelihood function for $\theta$ based on the complete sample is given by:

$L\left(\theta \right)={\displaystyle \prod _{i=1}^n{h}_{X_1,X_2}\left(x_{1i},x_{2i}\right)}={\displaystyle \prod _{i=1}^n{h}_{X_1}}\left(x_{1i}\right)h_{X_2}\left(x_{2i}\right)c\left(H\left(x_{1i}\right),H\left(x_{2i}\right)\right),$

where $c\left(H\left(x_{1i}\right),H\left(x_{2i}\right)\right)$ is the density function of the Clayton copula that connects the two marginal distributions.

Given observed data pairs $\left(x_{11},x_{21}\right),\cdots ,\left(x_{1n},x_{2n}\right)$ from the Weibull distribution, the likelihood function for parameter $\theta$ is:

$L\left(\theta \right)={\displaystyle \prod _{i=1}^n{\lambda }_1{\gamma }_1}{\left({\lambda }_1{x}_{1i}\right)}^{{\gamma }_1-1}{\text{e}}^{-{\left({\lambda }_1{x}_{1i}\right)}^{{\gamma }_1}}{\lambda }_2{\gamma }_2{\left({\lambda }_2{x}_{2i}\right)}^{{\gamma }_2-1}{\text{e}}^{-{\left({\lambda }_2{x}_{2i}\right)}^{{\gamma }_2}}{c}_{BCW}\left(H\left(x_{1i}\right),H\left(x_{2i}\right)\right),$

where

$c_{BCW}\left(H\left(x_{1i}\right),H\left(x_{2i}\right)\right)=\left(1+\eta \right){\left(u_1\left(x_{1i}\right)u_2\left(x_{2i}\right)\right)}^{-1-\eta }{\left(u_1{\left(x_{1i}\right)}^{-\eta }+u_2{\left(x_{2i}\right)}^{-\eta }-1\right)}^{-2-\frac{1}{\eta }}$ ,

in which $u_1\left(x_{1i}\right)=1-{\text{e}}^{-{\left({\lambda }_1{x}_{1i}\right)}^{{\gamma }_1}}$ , $u_2\left(x_{2i}\right)=1-{\text{e}}^{-{\left({\lambda }_2{x}_{2i}\right)}^{{\gamma }_2}}$ .

Let $l_{BCW}\left(\theta \right)$ denote the log-likelihood function of the BCW distribution, which is given by:

$\begin{array}{c}{l}_{BCW}\left(\theta \right)={\mu }_{BCW}\left(x_{1i}\right)+{\nu }_{BCW}\left(x_{2i}\right)+n\ln \left(1+\eta \right)\\ -\left(1+\eta \right){\displaystyle \sum }_{i=1}^n\left[\ln u_1\left(x_{1i}\right)+\ln u_2\left(x_{2i}\right)\right]\\ -\left(2+\frac{1}{\eta }\right){\displaystyle \sum }_{i=1}^n\ln \left[u_1{\left(x_{1i}\right)}^{-\eta }+u_2{\left(x_{2i}\right)}^{-\eta }-1\right],\end{array}$

where

$\begin{array}{c}{\mu }_{BCW}\left(x_{1i}\right)={\displaystyle \sum }_{i=1}^n\ln h_{X_1}\left(x_{1i}\right)\\ =n\ln {\lambda }_1+n\ln {\gamma }_1+\left({\gamma }_1-1\right){\displaystyle \sum }_{i=1}^n\ln \left({\lambda }_1{x}_{1i}\right)-{\displaystyle \sum }_{i=1}^n{\left({\lambda }_1{x}_{1i}\right)}^{{\gamma }_1},\end{array}$

$\begin{array}{c}{\nu }_{BCW}\left(x_{2i}\right)={\displaystyle \sum }_{i=1}^n\ln h_{X_2}\left(x_{2i}\right)\\ =n\ln {\lambda }_2+n\ln {\gamma }_2+\left({\gamma }_2-1\right){\displaystyle \sum }_{i=1}^n\ln \left({\lambda }_2{x}_{2i}\right)-{\displaystyle \sum }_{i=1}^n{\left({\lambda }_2{x}_{2i}\right)}^{{\gamma }_2}.\end{array}$

By calculating the partial derivatives of $l_{BCW}\left(\theta \right)$ , we obtain:

$\begin{array}{l}\frac{\partial l_{BCW}\left(\theta \right)}{\partial {\lambda }_k}=\frac{n}{{\lambda }_k}+\frac{n\left({\gamma }_k-1\right)}{{\lambda }_k}-{\displaystyle \sum }_{i=1}^n{\gamma }_k{x}_{ki}{\left({\lambda }_k{x}_{ki}\right)}^{{\gamma }_k-1}-\left(1+\eta \right){\displaystyle \sum }_{i=1}^n\frac{{\Lambda }_k\left(x_{ki}\right)}{u_k\left(x_{ki}\right)}\\ +\left(2+\frac{1}{\eta }\right){\displaystyle \sum }_{i=1}^n\frac{\eta {\Lambda }_k\left(x_{ki}\right)\left[u_k{\left(x_{ki}\right)}^{-\eta -1}\right]}{u_1{\left(x_{1i}\right)}^{-\eta }+u_2{\left(x_{2i}\right)}^{-\eta }-1},k=1,2,\end{array}$

$\begin{array}{l}\frac{\partial l_{BCW}\left(\theta \right)}{\partial {\gamma }_k}=\frac{n}{{\gamma }_k}+{\displaystyle \sum }_{i=1}^n\ln \left({\lambda }_k{x}_{ki}\right)-{\displaystyle \sum }_{i=1}^n{\left({\lambda }_k{x}_{ki}\right)}^{{\gamma }_k}\ln \left({\lambda }_k{x}_{ki}\right)-\left(1+\eta \right){\displaystyle \sum }_{i=1}^n\frac{{\Sigma }_k\left(x_{ki}\right)}{u_k\left(x_{ki}\right)}\\ +\left(2+\frac{1}{\eta }\right){\displaystyle \sum }_{i=1}^n\frac{\eta {\Sigma }_k\left(x_{ki}\right)\left[u_k{\left(x_{ki}\right)}^{-\eta -1}\right]}{u_1{\left(x_{1i}\right)}^{-\eta }+u_2{\left(x_{2i}\right)}^{-\eta }-1},k=1,2,\end{array}$

$\begin{array}{c}\frac{\partial l_{BCW}\left(\theta \right)}{\partial \eta }=\frac{n}{1+\eta }-{\displaystyle \sum }_{i=1}^n\left[\ln u_1\left(x_{1i}\right)+\ln u_2\left(x_{2i}\right)\right]\\ +\frac{1}{{\eta }^2}{\displaystyle \sum }_{i=1}^n\ln \left[u_1{\left(x_{1i}\right)}^{-\eta }+u_2{\left(x_{2i}\right)}^{-\eta }-1\right]\\ +\left(2+\frac{1}{\eta }\right){\displaystyle \sum }_{i=1}^n\frac{u_1{\left(x_{1i}\right)}^{-\eta }\ln u_1\left(x_{1i}\right)+u_2{\left(x_{2i}\right)}^{-\eta }\ln u_2\left(x_{2i}\right)}{u_1{\left(x_{1i}\right)}^{-\eta }+u_2{\left(x_{2i}\right)}^{-\eta }-1},\end{array}$

where

${\Lambda }_k\left(x_{ki}\right)=\frac{\partial u_k\left(x_{ki}\right)}{\partial {\lambda }_k}={\gamma }_k{x}_{ki}{\left({\lambda }_k{x}_{ki}\right)}^{{\gamma }_k-1}{\text{e}}^{-{\left({\lambda }_k{x}_{ki}\right)}^{{\gamma }_k}},k=1,2,$

${\Sigma }_k\left(x_{ki}\right)=\frac{\partial u_k\left(x_{ki}\right)}{\partial {\gamma }_k}={\left({\lambda }_k{x}_{ki}\right)}^{{\gamma }_k}\ln \left({\lambda }_k{x}_{ki}\right){\text{e}}^{-{\left({\lambda }_k{x}_{ki}\right)}^{{\gamma }_k}},k=1,2.$

The maximum likelihood estimates are obtained by numerically solving the nonlinear system resulting from setting the partial derivatives to zero.

4.2. MLE under Type-II Censoring

In practice, complete lifetime data are often unavailable due to limited observation periods or cost constraints, leading to censored data. Type-II censoring is a typical mechanism: in an experiment with $n$ samples, the test terminates when the $r$ -th failure ($r\le n$ ) is observed, yielding $r$ complete failures and $n-r$ censored observations.

Consider bivariate random samples $\left(x_{1\left(i:n\right)},x_{2\left(i:n\right)}\right)$ for $i=1,\cdots ,n$ from a bivariate Weibull distribution, with joint cdf $H\left(x_1,x_2\right)$ and pdf $h\left(x_1,x_2\right)$ . Let $x_{1\left(1:n\right)}<\cdots <x_{1\left(n:n\right)}$ be the order statistics of $X_1$ , and let $x_{2\left[i:n\right]}$ be the concomitant of the $i$ -th order statistic [19]. Under Type-II censoring, only the first $r$ ($r<n$ ) pairs $\left(x_{1\left(i:n\right)},x_{2\left[i:n\right]}\right)$ for $i=1,\cdots ,r$ are observed. The likelihood function is given by Balakrishnan and Kim [20] and Kim et al. [21]:

$L^{\text{*}}\left(\theta \right)=\frac{n!}{\left(n-r\right)!}{\left[1-H_{X_1}\left(x_{1r:n}\right)\right]}^{n-r}{\displaystyle \prod _{i=1}^r{h}_{X_1}}\left(x_{1i:n}\right)h_{X_2}\left(x_{2\left[i:n\right]}\right)c\left(H_{X_1}\left(x_{1i:n}\right),H_{X_2}\left(x_{2\left[i:n\right]}\right)\right),$

where $c\left(H_{X_1}\left(x_{1i:n}\right),H_{X_2}\left(x_{2\left[i:n\right]}\right)\right)$ denotes the Clayton copula density evaluated at the marginal cumulative distribution functions $H_{X_1}\left(x_{1i:n}\right)$ and $H_{X_2}\left(x_{2\left[i:n\right]}\right)$ . The log-likelihood function of the BCW distribution, denoted by ${\ell }_{BCW}^{\text{*}}\left(\theta \right)$ , is then given by:

$\begin{array}{c}{l}_{BCW}^{*}\left(\theta \right)\propto -\left(n-r\right){\left({\lambda }_1{x}_{1\left(r:n\right)}\right)}^{{\gamma }_1}+{\mu }_{BCW}^{*}\left(x_{1\left(i:n\right)}\right)+{\nu }_{BCW}^{*}\left(x_{2\left[i:n\right]}\right)\\ +r\ln \left(1+\eta \right)-\left(1+\eta \right){\displaystyle \sum }_{i=1}^r\left[\ln u_1\left(x_{1\left(i:n\right)}\right)+\ln u_2\left(x_{2\left[i:n\right]}\right)\right]\\ -\left(2+\frac{1}{\eta }\right){\displaystyle \sum }_{i=1}^r\ln \left[u_1{\left(x_{1\left(i:n\right)}\right)}^{-\eta }+u_2{\left(x_{2\left[i:n\right]}\right)}^{-\eta }-1\right],\end{array}$

where

$\begin{array}{c}{\mu }_{BCW}^{*}\left(x_{1\left(i:n\right)}\right)={\displaystyle \sum }_{i=1}^r\ln h_{X_1}\left(x_{1\left(i:n\right)}\right)\\ =r\ln {\lambda }_1+r\ln {\gamma }_1+\left({\gamma }_1-1\right){\displaystyle \sum }_{i=1}^r\ln \left({\lambda }_1{x}_{1\left(i:n\right)}\right)-{\displaystyle \sum }_{i=1}^r{\left({\lambda }_1{x}_{1\left(i:n\right)}\right)}^{{\gamma }_1},\end{array}$

$\begin{array}{c}{\nu }_{BCW}^{*}\left(x_{2\left[i:n\right]}\right)={\displaystyle \sum }_{i=1}^r\ln h_{X_2}\left(x_{2\left[i:n\right]}\right)\\ =r\ln {\lambda }_2+r\ln {\gamma }_2+\left({\gamma }_2-1\right){\displaystyle \sum }_{i=1}^r\ln \left({\lambda }_2{x}_{2\left[i:n\right]}\right)-{\displaystyle \sum }_{i=1}^r{\left({\lambda }_2{x}_{2\left[i:n\right]}\right)}^{{\gamma }_2}.\end{array}$

Differentiating $l_{BCW}^{\text{*}}\left(\theta \right)$ partially with respect to $\theta$ leads to:

$\begin{array}{c}\frac{\partial l_{BCW}^{*}\left(\theta \right)}{\partial {\lambda }_1}=-\left(n-r\right){\gamma }_1{x}_{1\left(r:n\right)}{\left({\lambda }_1{x}_{1\left(r:n\right)}\right)}^{{\gamma }_1-1}+\frac{r}{{\lambda }_1}+\frac{r\left({\gamma }_1-1\right)}{{\lambda }_1}\\ -{\displaystyle \sum }_{i=1}^r{\gamma }_1{x}_{1\left(i:n\right)}{\left({\lambda }_1{x}_{1\left(i:n\right)}\right)}^{{\gamma }_1-1}-\left(1+\eta \right){\displaystyle \sum }_{i=1}^r\frac{{\Lambda }_1\left(x_{1\left(i:n\right)}\right)}{u_1\left(x_{1\left(i:n\right)}\right)}\\ +\left(2+\frac{1}{\eta }\right){\displaystyle \sum }_{i=1}^r\frac{\eta {\Lambda }_1\left(x_{1\left(i:n\right)}\right)\left[u_1{\left(x_{1\left(i:n\right)}\right)}^{-\eta -1}\right]}{u_1{\left(x_{1\left(i:n\right)}\right)}^{-\eta }+u_2{\left(x_{2\left[i:n\right]}\right)}^{-\eta }-1},\end{array}$

$\begin{array}{c}\frac{\partial l_{BCW}^{*}\left(\theta \right)}{\partial {\lambda }_2}=\frac{r}{{\lambda }_2}+\frac{r\left({\gamma }_2-1\right)}{{\lambda }_2}-{\displaystyle \sum }_{i=1}^r{\gamma }_2{x}_{2\left[i:n\right]}{\left({\lambda }_2{x}_{2\left[i:n\right]}\right)}^{{\gamma }_2-1}-\left(1+\eta \right){\displaystyle \sum }_{i=1}^r\frac{{\Lambda }_2\left(x_{2\left[i:n\right]}\right)}{u_2\left(x_{2\left[i:n\right]}\right)}\\ +\left(2+\frac{1}{\eta }\right){\displaystyle \sum }_{i=1}^r\frac{\eta {\Lambda }_2\left(x_{2\left[i:n\right]}\right)\left[u_2{\left(x_{2\left[i:n\right]}\right)}^{-\eta -1}\right]}{u_1{\left(x_{1\left(i:n\right)}\right)}^{-\eta }+u_2{\left(x_{2\left[i:n\right]}\right)}^{-\eta }-1},\end{array}$

$\begin{array}{c}\frac{\partial l_{BCW}^{*}\left(\theta \right)}{\partial {\gamma }_1}=-\left(n-r\right){\left({\lambda }_1{x}_{1\left(r:n\right)}\right)}^{{\gamma }_1}\ln \left({\lambda }_1{x}_{1\left(r:n\right)}\right)+\frac{r}{{\gamma }_1}+{\displaystyle \sum }_{i=1}^r\ln \left({\lambda }_1{x}_{1\left(i:n\right)}\right)\\ -{\displaystyle \sum }_{i=1}^r{\left({\lambda }_1{x}_{1\left(i:n\right)}\right)}^{{\gamma }_1}\ln \left({\lambda }_1{x}_{1\left(i:n\right)}\right)-\left(1+\eta \right){\displaystyle \sum }_{i=1}^r\frac{{\Sigma }_1\left(x_{1\left(i:n\right)}\right)}{u_1\left(x_{1\left(i:n\right)}\right)}\\ +\left(2+\frac{1}{\eta }\right){\displaystyle \sum }_{i=1}^r\frac{\eta {\Sigma }_1\left(x_{1\left(i:n\right)}\right)\left[u_1{\left(x_{1\left(i:n\right)}\right)}^{-\eta -1}\right]}{u_1{\left(x_{1\left(i:n\right)}\right)}^{-\eta }+u_2{\left(x_{2\left[i:n\right]}\right)}^{-\eta }-1},\end{array}$

$\begin{array}{c}\frac{\partial l_{BCW}^{*}\left(\theta \right)}{\partial {\gamma }_2}=\frac{r}{{\gamma }_2}+{\displaystyle \sum }_{i=1}^r\ln \left({\lambda }_2{x}_{2\left[i:n\right]}\right)-{\displaystyle \sum }_{i=1}^r{\left({\lambda }_2{x}_{2\left[i:n\right]}\right)}^{{\gamma }_2}\ln \left({\lambda }_2{x}_{2\left[i:n\right]}\right)\\ -\left(1+\eta \right){\displaystyle \sum }_{i=1}^r\frac{{\Sigma }_2\left(x_{2\left[i:n\right]}\right)}{u_2\left(x_{2\left[i:n\right]}\right)}+\left(2+\frac{1}{\eta }\right){\displaystyle \sum }_{i=1}^r\frac{\eta {\Sigma }_2\left(x_{2\left[i:n\right]}\right)\left[u_2{\left(x_{2\left[i:n\right]}\right)}^{-\eta -1}\right]}{u_1{\left(x_{1\left(i:n\right)}\right)}^{-\eta }+u_2{\left(x_{2\left[i:n\right]}\right)}^{-\eta }-1},\end{array}$

$\begin{array}{c}\frac{\partial l_{BCW}^{*}\left(\theta \right)}{\partial \eta }=\frac{r}{1+\eta }-{\displaystyle \sum }_{i=1}^r\left[\ln u_1\left(x_{1\left(i:n\right)}\right)+\ln u_2\left(x_{2\left[i:n\right]}\right)\right]\\ +\frac{1}{{\eta }^2}{\displaystyle \sum }_{i=1}^r\ln \left[u_1{\left(x_{1\left(i:n\right)}\right)}^{-\eta }+u_2{\left(x_{2\left[i:n\right]}\right)}^{-\eta }-1\right]\\ +\left(2+\frac{1}{\eta }\right){\displaystyle \sum }_{i=1}^r\frac{u_1{\left(x_{1\left(i:n\right)}\right)}^{-\eta }\ln u_1\left(x_{1\left(i:n\right)}\right)+u_2{\left(x_{2\left[i:n\right]}\right)}^{-\eta }\ln u_2\left(x_{2\left[i:n\right]}\right)}{u_1{\left(x_{1\left(i:n\right)}\right)}^{-\eta }+u_2{\left(x_{2\left[i:n\right]}\right)}^{-\eta }-1},\end{array}$

where

${\Lambda }_1\left(x_{1\left(i:n\right)}\right)=\frac{\partial u_1\left(x_{1\left(i:n\right)}\right)}{\partial {\lambda }_1}={\gamma }_1{x}_{1\left(i:n\right)}{\left({\lambda }_1{x}_{1\left(i:n\right)}\right)}^{{\gamma }_1-1}{\text{e}}^{-{\left({\lambda }_1{x}_{1\left(i:n\right)}\right)}^{{\gamma }_1}}，$

${\Lambda }_2\left(x_{2\left[i:n\right]}\right)=\frac{\partial u_2\left(x_{2\left[i:n\right]}\right)}{\partial {\lambda }_2}={\gamma }_2{x}_{2\left[i:n\right]}{\left({\lambda }_2{x}_{2\left[i:n\right]}\right)}^{{\gamma }_2-1}{\text{e}}^{-{\left({\lambda }_2{x}_{2\left[i:n\right]}\right)}^{{\gamma }_2}},$

${\Sigma }_1\left(x_{1\left(i:n\right)}\right)=\frac{\partial u_1\left(x_{1\left(i:n\right)}\right)}{\partial {\gamma }_1}={\left({\lambda }_1{x}_{1\left(i:n\right)}\right)}^{{\gamma }_1}\ln \left({\lambda }_1{x}_{1\left(i:n\right)}\right){\text{e}}^{-{\left({\lambda }_1{x}_{1\left(i:n\right)}\right)}^{{\gamma }_1}},$

${\Sigma }_2\left(x_{2\left[i:n\right]}\right)=\frac{\partial u_2\left(x_{2\left[i:n\right]}\right)}{\partial {\gamma }_2}={\left({\lambda }_2{x}_{2\left[i:n\right]}\right)}^{{\gamma }_2}\ln \left({\lambda }_2{x}_{2\left[i:n\right]}\right){\text{e}}^{-{\left({\lambda }_2{x}_{2\left[i:n\right]}\right)}^{{\gamma }_2}}.$

Equating the partial derivatives to zero gives a nonlinear system in $\theta$ . Since $\widehat{\theta }$ has no closed-form solution, it is obtained numerically using nonlinear optimization.

4.3. Bayesian Estimation

In contrast to frequentist methods that treat parameters as fixed, Bayesian methods consider them as random variables described by prior distributions. Priors, based on historical or domain knowledge, are updated with sample data to produce posterior distributions, thereby combining prior and data information. For reliability analysis with scarce data, this framework enhances estimation accuracy and stability, making prior selection critical. We adopt independent Gamma priors for ${\lambda }_1,{\gamma }_1,{\lambda }_2,{\gamma }_2,\eta$ , chosen for their positive support $\left(0,\infty \right)$ and distributional flexibility [22]. Under the mean squared error loss function, the Bayesian posterior estimates are given by the following priors:

${\pi }_j\left({\theta }_j\right)\propto {\theta }_j^{{\alpha }_j-1}{\text{e}}^{-{\beta }_j{\theta }_j},{\alpha }_j,{\beta }_j>0,j=1,\cdots ,5.$

Then the joint prior distribution is:

$\begin{array}{l}\pi \left(\theta \right)\propto {\lambda }_1^{{\alpha }_1-1}{\gamma }_1^{{\alpha }_2-1}{\lambda }_2^{{\alpha }_3-1}{\gamma }_2^{{\alpha }_4-1}{\eta }^{{\alpha }_5-1}\\ \times {\text{e}}^{-\left({\beta }_1{\lambda }_1+{\beta }_2{\gamma }_1+{\beta }_3{\lambda }_2+{\beta }_4{\gamma }_2+{\beta }_5\eta \right)},{\alpha }_j,{\beta }_j>0,j=1,\cdots ,5.\end{array}$

The hyperparameters $\left({\alpha }_j,{\beta }_j\right)$ are calibrated using the moment matching approach. Suppose that for each parameter ${\delta }_j$ ($j=1,\cdots ,5$ ), a total of $k$ MLE estimates ${\left\{{\widehat{\delta }}_j^{\left(i\right)}\right\}}_{i=1}^k$ are obtained from $k$ independent simulated datasets. The sample mean and sample variance of these estimates are then calculated as:

${\overline{\delta }}_j=\frac{1}{k}{\displaystyle \sum }_{i=1}^k{\widehat{\delta }}_j^{\left(i\right)},s_j^2=\frac{1}{k-1}{\displaystyle \sum }_{i=1}^k{\left({\widehat{\delta }}_j^{\left(i\right)}-{\overline{\delta }}_j\right)}^2.$

By matching the theoretical moments of the Gamma prior distribution—specifically, its mean ${\alpha }_j/{\beta }_j$ and variance ${\alpha }_j/{\beta }_j^2$ —to their sample counterparts ${\overline{\delta }}_j$ and $s_j^2$ , we obtain the following closed-form solutions for the hyperparameters:

${\alpha }_j=\frac{{\overline{\delta }}_j^2}{s_j^2},{\beta }_j=\frac{{\overline{\delta }}_j}{s_j^2},j=1,\cdots ,5.$

The correspondence between the parameters and the ${\delta }_j$ notation is established as: ${\delta }_1={\lambda }_1$ , ${\delta }_2={\gamma }_1$ , ${\delta }_3={\lambda }_2$ , ${\delta }_4={\gamma }_2$ , ${\delta }_5=\eta$ .

Let ${\pi }^{\text{*}}\left({\lambda }_1,{\gamma }_1,{\lambda }_2,{\gamma }_2,\eta |\text{data}\right)$ denote the joint posterior distribution. By Bayes’ theorem, the posterior is proportional to the prior times the likelihood. For the BCW model, the joint posterior is given by:

${\pi }_{BCW}^{\text{*}}\left({\lambda }_1,{\gamma }_1,{\lambda }_2,{\gamma }_2,\eta |\text{data}\right)\propto \text{prior}\times \text{likelihood}.$

Under the squared error loss function, the Bayes estimator is the posterior mean. However, due to the lack of a closed-form solution, we resort to MCMC methods for posterior sampling. MCMC constructs a Markov chain with the target posterior as its stationary distribution. After convergence, the sampled draws are used to compute the Bayesian estimates. The full conditional posterior distributions for each parameter are derived as follows:

$\begin{array}{c}{\pi }_{BCW}^{\text{*}}\left({\lambda }_1|{\gamma }_1,{\lambda }_2,{\gamma }_2,\eta ,\text{data}\right)\propto {\lambda }_1^{{\alpha }_1+r{\gamma }_1-1}{\text{e}}^{-{\beta }_1{\lambda }_1}{\displaystyle \prod _{i=1}^r{\text{e}}^{-{\left({\lambda }_1{x}_{1\left(i:n\right)}\right)}^{{\gamma }_1}}}\\ \times c_{BCW}\left(F_{X_1}\left(x_{1\left(i:n\right)}\right),F_{X_2}\left(x_{2\left[i:n\right]}\right)\right),\end{array}$

$\begin{array}{c}{\pi }_{BCW}^{\text{*}}\left({\gamma }_1|{\lambda }_1,{\lambda }_2,{\gamma }_2,\eta ,\text{data}\right)\propto {\lambda }_1^{{\alpha }_1+r{\gamma }_1-1}{\gamma }_1^{{\alpha }_2+r-1}{\text{e}}^{-{\beta }_2{\gamma }_1}{\displaystyle \prod _{i=1}^r{x}_{1\left(i:n\right)}^{{\gamma }_1-1}{\text{e}}^{-{\left({\lambda }_1{x}_{1\left(i:n\right)}\right)}^{{\gamma }_1}}}\\ \times c_{BCW}\left(F_{X_1}\left(x_{1\left(i:n\right)}\right),F_{X_2}\left(x_{2\left[i:n\right]}\right)\right),\end{array}$

$\begin{array}{c}{\pi }_{BCW}^{\text{*}}\left({\lambda }_2|{\lambda }_1,{\gamma }_1,{\gamma }_2,\eta ,\text{data}\right)\propto {\lambda }_2^{{\alpha }_3+r{\gamma }_2-1}{\text{e}}^{-{\beta }_3{\lambda }_2}{\displaystyle \prod _{i=1}^r{\text{e}}^{-{\left({\lambda }_2{x}_{2\left[i:n\right]}\right)}^{{\gamma }_2}}}\\ \times c_{BCW}\left(F_{X_1}\left(x_{1\left(i:n\right)}\right),F_{X_2}\left(x_{2\left[i:n\right]}\right)\right),\end{array}$

$\begin{array}{c}{\pi }_{BCW}^{\text{*}}\left({\gamma }_2|{\lambda }_1,{\gamma }_1,{\lambda }_2,\eta ,\text{data}\right)\propto {\gamma }_2^{{\alpha }_4+r-1}{\text{e}}^{-{\beta }_4{\gamma }_2}{\displaystyle \prod _{i=1}^r{\left({\lambda }_2{x}_{2\left[i:n\right]}\right)}^{{\gamma }_2-1}{\text{e}}^{-{\left({\lambda }_2{x}_{2\left[i:n\right]}\right)}^{{\gamma }_2}}}\\ \times c_{BCW}\left(F_{X_1}\left(x_{1\left(i:n\right)}\right),F_{X_2}\left(x_{2\left[i:n\right]}\right)\right),\end{array}$

${\pi }_{BCW}^{\text{*}}\left(\eta |{\lambda }_1,{\lambda }_2,{\gamma }_2,\eta ,\text{data}\right)\propto {\eta }^{{\alpha }_5-1}{\text{e}}^{-{\beta }_5\eta }{\displaystyle \prod _{i=1}^r{c}_{BCW}}\left(H_{x_1}\left(x_{1\left(i:n\right)}\right),H_{x_2}\left(x_{2\left[i:n\right]}\right)\right).$

Given the complexity of the full conditional distributions, direct sampling is infeasible. We thus employ the Metropolis-Hastings (MH) algorithm, using a normal proposal distribution, to approximate the posterior distribution through an acceptance-rejection sampling scheme. The MH algorithm with a normal proposal is a well-established approach in Bayesian inference for obtaining parameter estimates and credible intervals [23].

5. Interval Estimation

This section presents two types of interval estimates for the BCW model parameters ${\lambda }_1,{\gamma }_1,{\lambda }_2,{\gamma }_2,\eta$ : confidence intervals and credible intervals. The construction procedures are detailed below.

5.1. Asymptotic Confidence Intervals

Based on the asymptotic normality property of maximum likelihood estimation, confidence intervals for the unknown parameters can be constructed using the Fisher information matrix. Let $\widehat{\theta }={\left({\widehat{\lambda }}_1,{\widehat{\gamma }}_1,{\widehat{\lambda }}_2,{\widehat{\gamma }}_2,\widehat{\eta }\right)}^{\top }$ denote the maximum likelihood estimator. The Fisher information matrix $I\left(\theta \right)$ , evaluated at $\widehat{\theta }$ , quantifies the information contained in the observed data. The asymptotic variance-covariance matrix of the estimator is obtained by inverting $I\left(\widehat{\theta }\right)$ , i.e., $V\left(\widehat{\theta }\right)=I^{-1}\left(\widehat{\theta }\right)$ . The explicit form of $I\left(\widehat{\theta }\right)$ is given by:

$I\left(\widehat{\theta }\right)=-\left(\begin{array}{ccccc}{I}_{{\widehat{\lambda }}_1{\widehat{\lambda }}_1}& & & & \\ I_{{\widehat{\lambda }}_2{\widehat{\lambda }}_1}& I_{{\widehat{\lambda }}_2{\widehat{\lambda }}_2}& & & \\ I_{{\widehat{\gamma }}_1{\widehat{\lambda }}_1}& I_{{\widehat{\gamma }}_1{\widehat{\lambda }}_2}& I_{{\widehat{\gamma }}_1{\widehat{\gamma }}_1}& & \\ I_{{\widehat{\gamma }}_2{\widehat{\lambda }}_1}& I_{{\widehat{\gamma }}_2{\widehat{\lambda }}_2}& I_{{\widehat{\gamma }}_2{\widehat{\gamma }}_1}& I_{{\widehat{\gamma }}_2{\widehat{\gamma }}_2}& \\ I_{\widehat{\eta }{\widehat{\lambda }}_1}& I_{\widehat{\eta }{\widehat{\lambda }}_2}& I_{\widehat{\eta }{\widehat{\gamma }}_1}& I_{\widehat{\eta }{\widehat{\gamma }}_2}& I_{\widehat{\eta }\widehat{\eta }}\end{array}\right).$

In accordance with the asymptotic normality of MLEs, the $100\left(1-\alpha \right)\text{\%}$ confidence interval for each parameter ${\theta }_k$ is formulated as:

${\widehat{\theta }}_k\pm z_{\alpha /2}\sqrt{V_{kk}},k=1,\cdots ,5,$

where ${\widehat{\theta }}_k$ is the $k$ -th parameter estimate, $z_{\alpha /2}$ is the upper $\alpha /2$ quantile of the standard normal distribution, and $V_{kk}$ is the $k$ -th diagonal element of $V\left(\widehat{\theta }\right)$ , representing the asymptotic variance of ${\widehat{\theta }}_k$ .

5.2. Confidence Intervals

Bayesian credible intervals for the individual parameters are constructed via the following numerical procedure [24]. Let ${\theta }^{\left(j\right)}$ denote the parameter vector sampled at the $j$ -th MCMC iteration. After discarding the burn-in period, extract the marginal samples for each parameter, yielding sequences of length $L$ :

${\lambda }_1^{\left[1\right]},{\lambda }_1^{\left[2\right]},\cdots ,{\lambda }_1^{\left[L\right]},{\gamma }_1^{\left[1\right]},{\gamma }_1^{\left[2\right]},\cdots ,{\gamma }_1^{\left[L\right]},{\lambda }_2^{\left[1\right]},{\lambda }_2^{\left[2\right]},\cdots ,{\lambda }_2^{\left[L\right]},{\gamma }_2^{\left[1\right]},{\gamma }_2^{\left[2\right]},\cdots ,{\gamma }_2^{\left[L\right]},{\eta }_1^{\left[1\right]},{\eta }_1^{\left[2\right]},\cdots ,{\eta }_1^{\left[L\right]},$

where $L$ is the effective number of MCMC samples after burn-in. Sorting the $L$ samples of each parameter in ascending order yields the corresponding order statistics:

${\tilde{\lambda }}_i^{\left(1\right)}\le {\tilde{\lambda }}_i^{\left(2\right)}\le \cdots \le {\tilde{\lambda }}_i^{\left(L\right)},{\tilde{\gamma }}_i^{\left(1\right)}\le {\tilde{\gamma }}_i^{\left(2\right)}\le \cdots \le {\tilde{\gamma }}_i^{\left(L\right)},{\tilde{\eta }}_i^{\left(1\right)}\le {\tilde{\eta }}_i^{\left(2\right)}\le \cdots \le {\tilde{\eta }}_i^{\left(L\right)}\left(i=1,2\right).$

Based on these ordered samples, the percentile method is adopted to construct the $100\left(1-\alpha \right)\text{\%}$ symmetric Bayesian credible intervals. For a given significance level $\alpha$ , the lower and upper bounds are taken as the $L_{\alpha /2}$ -th and $L_{\left(1-\alpha /2\right)}$ -th order statistics, respectively, where $\lfloor \cdot \rfloor$ denotes the floor function. Thus, the credible intervals for the individual parameters are:

${\lambda }_i:\left({\tilde{\lambda }}_i^{\left(L_{\alpha /2}\right)},{\tilde{\lambda }}_i^{\left(L_{\left(1-\alpha /2\right)}\right)}\right),{\gamma }_i:\left({\tilde{\gamma }}_i^{\left(L_{\alpha /2}\right)},{\tilde{\gamma }}_i^{\left(L_{\left(1-\alpha /2\right)}\right)}\right),\eta :\left({\tilde{\eta }}^{\left(L_{\alpha /2}\right)},{\tilde{\eta }}^{\left(L_{\left(1-\alpha /2\right)}\right)}\right)\left(i=1,2\right).$

Setting $\alpha =0.05$ yields the 95% symmetric credible intervals for all unknown parameters.

6. Numerical Simulation Experiments

In this section, Monte Carlo simulations are conducted to evaluate the statistical performance of the BCE and BCW models, with a focus on point estimation accuracy, interval estimation reliability, and the effectiveness of copula dependence parameter estimation. Following the conditional distribution method proposed by Nelsen [5], bivariate random samples are generated under the Clayton copula framework, and the performance of the models is systematically compared under different censoring levels.

Two sets of true parameter values are considered:

1) ${\lambda }_1=0.9,{\gamma }_1=0.7,{\lambda }_2=1.0,{\gamma }_2=0.6,\eta =1.3$ ;

2) ${\lambda }_1=3.3,{\gamma }_1=1.7,{\lambda }_2=1.5,{\gamma }_2=2.5,\eta =2.4$ ;

3) ${\lambda }_1=1.8,{\gamma }_1=1.2,{\lambda }_2=2.2,{\gamma }_2=1.8,\eta =1.8$ .

Under fixed sample sizes $n=45$ and $n=140$ , we evaluate both complete and Type-II censored samples, with censoring levels $r=27,36,45$ for $n=45$ and $r=105,130,140$ for $n=140$ , where $r$ is the number of observed failures. For each parameter-censoring combination, 250 independent replications are carried out to ensure robust simulation results. Parameter estimation is carried out using Markov chain Monte Carlo (MCMC) methods within the Bayesian framework, with MCMC chains generated using the coda package in R, each running for 12,000 iterations and discarding the first 2000 as burn-in to ensure convergence to the stationary distribution. For comparison, point estimates are also obtained via MLE using the Newton-Raphson algorithm implemented in the maxLik package.

The simulation results for the parameters ${\lambda }_1$ , ${\gamma }_1$ , ${\lambda }_2$ , ${\gamma }_2$ , $\eta$ , along with the dependence measures $\tau$ and $\rho$ , are summarized in Tables 1-3 for both MLE and Bayesian estimation approaches. The performance of the estimators is evaluated using the following metrics: MSEs (mean squared errors), Abias (average bias), LACI (average length of asymptotic confidence intervals), and LCCI (average length of credible intervals).

Table 1. The estimation results for parameters ${\lambda }_1=0.9$ , ${\gamma }_1=0.7$ , ${\lambda }_2=1.0$ , ${\gamma }_2=0.6$ , $\eta =1.3$ .

n	r		MLE						Bayes
			BCE			BCW			BCE			BCW
			Abias	MSEs	LACI	Abias	MSEs	LACI	Abias	MSEs	LCCI	Abias	MSEs	LCCI
45	27	${\lambda }_1$	0.0451	0.0361	0.6890	0.0922	0.0970	1.0432	0.0066	0.1262	0.5541	0.0474	0.2308	0.9461
		${\gamma }_1$				0.0433	0.0207	0.4782				0.0339	0.1512	0.4969
		${\lambda }_2$	0.0498	0.0489	0.9239	0.1229	0.1930	1.6507	0.0198	0.1535	0.6997	0.1427	0.3612	1.5309
		${\gamma }_2$				0.0323	0.0105	0.3727				0.0122	0.0936	0.3690
		$\eta$	0.1409	0.2474	1.6691	0.1120	0.3030	2.0198	0.0077	0.1800	1.0512	−0.0534	0.5507	2.0223
		$\rho$	0.0136	0.0107	0.3566	0.0020	0.0141	0.4276	−0.0010	0.0404	0.2414	−0.0408	0.1458	0.4849
		$\tau$	0.0138	0.0069	0.2874	0.0054	0.0089	0.3446	−0.0003	0.0324	0.1920	−0.0277	0.1113	0.3775
	36	${\lambda }_1$	0.0430	0.0243	0.5890	0.0710	0.0586	0.9017	0.0189	0.1149	0.5051	0.0406	0.2067	0.8233
		${\gamma }_1$				0.0169	0.0104	0.3824				0.0175	0.0978	0.3801
		${\lambda }_2$	0.0578	0.0440	0.7170	0.1289	0.1570	1.3154	0.0093	0.1309	0.5869	0.0919	0.2945	1.1936
		${\gamma }_2$				0.0184	0.0066	0.3068				0.0028	0.0682	0.2952
		$\eta$	0.0960	0.1660	1.4820	0.1043	0.2348	1.7866	0.0007	0.1771	1.0038	−0.0005	0.4722	1.8056
		$\rho$	0.0089	0.0075	0.3239	0.0058	0.0102	0.3831	−0.0026	0.0405	0.2318	−0.0207	0.1210	0.4239
		$\tau$	0.0092	0.0049	0.2601	0.0075	0.0066	0.3083	−0.0016	0.0324	0.1842	−0.0130	0.0932	0.3324
	45	${\lambda }_1$	0.0229	0.0167	0.5104	0.0338	0.0417	0.7943	0.0130	0.1087	0.4587	0.0282	0.1886	0.7739
		${\gamma }_1$				0.0284	0.0079	0.3191				0.0048	0.0742	0.3016
		${\lambda }_2$	0.0149	0.0220	0.5593	0.0239	0.0714	1.0310	0.0123	0.1173	0.5046	0.0558	0.2339	1.0080
		${\gamma }_2$				0.0172	0.0051	0.2687				0.0068	0.0664	0.2601
		$\eta$	0.0641	0.1403	1.3475	0.0440	0.1740	1.5657	−0.0047	0.1695	0.9575	0.0471	0.4105	1.6381
		$\rho$	0.0034	0.0068	0.3006	−0.0046	0.0093	0.3490	−0.0036	0.0388	0.2222	−0.0035	0.0951	0.3728
		$\tau$	0.0046	0.0043	0.2406	−0.0012	0.0058	0.2790	−0.0024	0.0310	0.1765	−0.0005	0.0750	0.2951
140	105	${\lambda }_1$	−0.1474	0.0499	0.2419	−0.2382	0.0961	0.3469	0.0060	0.0846	0.3289	0.0199	0.1324	0.4919
		${\gamma }_1$				0.0312	0.0074	0.2070				0.0021	0.0522	0.2191
		${\lambda }_2$	−0.0897	0.0206	0.3502	−0.1554	0.0489	0.5335	0.0146	0.1054	0.4071	0.0186	0.1827	0.7020
		${\gamma }_2$				0.0459	0.0041	0.1790				0.0016	0.0431	0.1737
		$\eta$	−0.0195	0.0474	0.8259	−0.0921	0.0716	0.9019	−0.0214	0.1976	0.8307	0.0184	0.2474	1.0086
		$\rho$	−0.0087	0.0026	0.1728	−0.0278	0.0042	0.2371	−0.0085	0.0472	0.1988	−0.0009	0.0541	0.2296
		$\tau$	−0.0062	0.0016	0.1442	−0.0210	0.0026	0.1801	−0.0062	0.0374	0.1570	0.0001	0.0435	0.1829
	130	${\lambda }_1$	−0.0545	0.0114	0.2745	−0.0583	0.0193	0.4225	0.0095	0.0760	0.2941	−0.0874	0.0451	0.6605
		${\gamma }_1$				0.0097	0.0026	0.1817				0.0167	0.0139	0.4268
		${\lambda }_2$	−0.0333	0.0093	0.3187	−0.0240	0.0212	0.5723	0.0041	0.0858	0.3328	−0.0139	0.0029	0.2082
		${\gamma }_2$				0.0154	0.0015	0.1559				0.0442	0.0256	0.6309
		$\eta$	0.0182	0.0448	0.7606	0.0107	0.0504	0.8881	0.0135	0.1948	0.7660	0.0241	0.1182	1.3526
		$\rho$	0.0003	0.0023	0.1756	−0.0019	0.0025	0.2094	−0.0003	0.0447	0.1771	−0.0012	0.0014	0.1631
		$\tau$	0.0009	0.0015	0.1381	−0.0008	0.0016	0.1612	0.0004	0.0357	0.1409	−0.0002	0.0012	0.1455
	140	${\lambda }_1$	0.0084	0.0050	0.2845	0.0163	0.0135	0.4527	0.0057	0.0708	0.2813	0.0097	0.1128	0.4461
		${\gamma }_1$				0.0046	0.0019	0.1743				0.0001	0.0457	0.1710
		${\lambda }_2$	0.0114	0.0075	0.3164	0.0247	0.0225	0.5875	0.0128	0.0812	0.3146	0.0152	0.1460	0.5787
		${\gamma }_2$				0.0071	0.0013	0.1495				0.0019	0.0384	0.1478
		$\eta$	0.0287	0.0328	0.7431	0.0438	0.0570	0.8815	0.0314	0.1939	0.7530	0.0064	0.2136	0.8715
		$\rho$	0.0038	0.0016	0.1651	0.0054	0.0026	0.2037	0.0040	0.0428	0.1714	−0.0025	0.0482	0.2001
		$\tau$	0.0035	0.0010	0.1320	0.0051	0.0017	0.1579	0.0038	0.0344	0.1368	−0.0013	0.0386	0.1593

Table 2. The estimation results for parameters ${\lambda }_1=3.3$ , ${\gamma }_1=1.7$ , ${\lambda }_2=1.5$ , ${\gamma }_2=2.5$ , $\eta =2.4$ .

n	r		MLE						Bayes
			BCE			BCW			BCE			BCW
			Abias	MSEs	LACI	Abias	MSEs	LACI	Abias	MSEs	LCCI	Abias	MSEs	LCCI
45	27	${\lambda }_1$	0.1622	0.4818	2.5603	0.1020	0.1777	1.4899	−0.0028	0.2652	1.5362	0.0323	0.3847	1.5216
		${\gamma }_1$				0.1024	0.1170	1.1653				0.0594	0.2840	1.1040
		${\lambda }_2$	0.0734	0.1043	1.3295	0.0296	0.0221	0.5701	0.0169	0.1660	0.8694	0.0225	0.1533	0.5953
		${\gamma }_2$				0.1581	0.2753	1.7357				0.0820	0.4141	1.6066
		$\eta$	0.2300	0.4574	2.2547	0.1854	0.7702	3.1927	0.0060	0.2455	1.4281	0.0395	0.7809	3.1474
		$\rho$	0.0125	0.0038	0.2205	−0.0021	0.0080	0.3174	−0.0013	0.0265	0.1570	−0.0158	0.0948	0.3554
		$\tau$	0.0147	0.0036	0.2137	0.0032	0.0070	0.3035	−0.0008	0.0253	0.1482	−0.0100	0.0844	0.3238
	36	${\lambda }_1$	0.1609	0.3293	2.1586	0.0849	0.1093	1.2845	0.0317	0.2587	1.4567	0.0362	0.3361	1.2974
		${\gamma }_1$				0.0395	0.0603	0.9091				0.0294	0.2233	0.8995
		${\lambda }_2$	0.0860	0.0914	1.0640	0.0331	0.0153	0.4306	0.0067	0.1588	0.7625	0.0164	0.1143	0.4429
		${\gamma }_2$				0.0790	0.1322	1.3267				0.0253	0.3032	1.2735
		$\eta$	0.1636	0.3158	2.0606	0.1806	0.5649	2.7681	0.0098	0.2527	1.3820	0.0747	0.6716	2.7623
		$\rho$	0.0086	0.0030	0.2071	0.0035	0.0054	0.2755	−0.0010	0.0271	0.1514	−0.0064	0.0751	0.3008
		$\tau$	0.0103	0.0028	0.1999	0.0070	0.0049	0.2648	−0.0005	0.0259	0.1431	−0.0026	0.0695	0.2784
	45	${\lambda }_1$	0.0795	0.2077	1.8059	0.0303	0.0848	1.1608	0.0225	0.2661	1.3527	0.0234	0.2970	1.2024
		${\gamma }_1$				0.0645	0.0430	0.7597				0.0086	0.1793	0.7202
		${\lambda }_2$	0.0286	0.0474	0.8154	0.0815	0.0084	0.3584	0.0118	0.1423	0.6645	0.0088	0.0894	0.3702
		${\gamma }_2$				0.0767	0.0043	1.1006				0.0184	0.2712	1.0629
		$\eta$	0.0969	0.2708	1.8997	0.0570	0.3785	2.3735	0.0005	0.2306	1.3270	0.1360	0.6215	2.4944
		$\rho$	0.0026	0.0027	0.1978	−0.0058	0.0046	0.2528	−0.0017	0.0250	0.1462	0.0031	0.0636	0.2595
		$\tau$	0.0043	0.0024	0.1894	−0.0026	0.0040	0.2398	−0.0012	0.0238	0.1381	0.0057	0.0605	0.2442
140	105	${\lambda }_1$	−0.5788	0.6435	0.9635	−0.4152	0.3104	0.6225	0.0025	0.3110	1.2039	0.0139	0.1980	0.7457
		${\gamma }_1$				0.0707	0.0405	0.4871				0.0039	0.1232	0.5242
		${\lambda }_2$	−0.1930	0.0843	0.4842	−0.0864	0.0142	0.2144	0.0188	0.1538	0.6013	0.0005	0.0700	0.2614
		${\gamma }_2$				0.1690	0.0820	0.7339				0.0089	0.1853	0.7664
		$\eta$	−0.0309	0.0938	1.1520	−0.1585	0.2233	1.3761	−0.0135	0.2587	1.1388	0.0369	0.3861	1.5484
		$\rho$	−0.0064	0.0011	0.1433	−0.0246	0.0029	0.1836	−0.0037	0.0292	0.1274	−0.0005	0.0391	0.1659
		$\tau$	−0.0054	0.0010	0.1312	−0.0216	0.0025	0.1620	−0.0030	0.0275	0.1200	0.0005	0.0377	0.1571
	130	${\lambda }_1$	−0.1224	0.1046	1.0006	−0.0874	0.0451	0.6605	0.0300	0.2685	1.0544	0.0136	0.1792	0.6791
		${\gamma }_1$				0.0167	0.0139	0.4268				0.0200	0.1287	0.4465
		${\lambda }_2$	−0.0312	0.0160	0.4676	−0.0139	0.0029	0.2082	0.0070	0.1239	0.4884	0.0043	0.0570	0.2159
		${\gamma }_2$				0.0442	0.0256	0.6309				0.0237	0.1689	0.6423
		$\eta$	0.0464	0.0868	1.0829	0.0241	0.1182	1.3526	0.0266	0.2687	1.0783	0.0005	0.3638	1.3745
		$\rho$	0.0023	0.0009	0.1254	−0.0012	0.0014	0.1631	0.0006	0.0290	0.1171	−0.0043	0.0400	0.1506
		$\tau$	0.0029	0.0008	0.1168	−0.0002	0.0012	0.1455	0.0011	0.0277	0.1112	−0.0030	0.0377	0.1420
	140	${\lambda }_1$	0.0245	0.0669	1.0064	0.0226	0.0290	0.6689	0.0182	0.2495	0.9941	0.0045	0.1755	0.6723
		${\gamma }_1$				0.0155	0.0112	0.4124				−0.0013	0.1094	0.4050
		${\lambda }_2$	0.0110	0.0135	0.4570	0.0050	0.0028	0.2072	0.0166	0.1143	0.4534	0.0008	0.0539	0.2070
		${\gamma }_2$				0.0053	0.0261	0.6043				0.0002	0.1587	0.5975
		$\eta$	0.0496	0.0623	1.0554	0.0402	0.0996	1.3319	0.0508	0.2785	1.0634	0.0358	0.3416	1.3335
		$\rho$	0.0034	0.0007	0.1071	0.0012	0.0011	0.1428	0.0031	0.0290	0.1135	0.0003	0.0356	0.1428
		$\tau$	0.0037	0.0006	0.1055	0.0019	0.0010	0.1361	0.0035	0.0279	0.1083	0.0011	0.0341	0.1356

Table 3. The estimation results for parameters ${\lambda }_1=1.8$ , ${\gamma }_1=1.2$ , ${\lambda }_2=2.2$ , ${\gamma }_2=1.8$ , $\eta =1.8$ .

n	r		MLE						Bayes
			BCE			BCW			BCE			BCW
			Abias	MSEs	LACI	Abias	MSEs	LACI	Abias	MSEs	LCCI	Abias	MSEs	LCCI
45	27	${\lambda }_1$	0.0894	0.1435	1.3845	0.0861	0.1115	1.1663	0.0019	0.1886	0.9495	0.0212	0.1926	0.9304
		${\gamma }_1$				0.0726	0.0587	0.8186				0.0252	0.1257	0.6002
		${\lambda }_2$	0.1069	0.2285	1.9942	0.0587	0.0862	1.1552	0.0188	0.2272	1.2573	0.0324	0.2010	0.9517
		${\gamma }_2$				0.1039	0.1141	1.1827				0.0290	0.1705	0.8204
		$\eta$	0.1839	0.3287	1.9293	0.1468	0.4861	2.5498	0.0087	0.2125	1.2267	0.0333	0.1821	1.3299
		$\rho$	0.0141	0.0061	0.2829	−0.0001	0.0106	0.3736	−0.0012	0.0338	0.1975	−0.0014	0.0288	0.2139
		$\tau$	0.0153	0.0048	0.2478	0.0047	0.0079	0.3259	−0.0004	0.0291	0.1698	−0.0007	0.0248	0.1839
	36	${\lambda }_1$	0.0870	0.0970	1.1735	0.0699	0.0686	1.0093	0.0264	0.1811	0.8868	0.0184	0.1838	0.8275
		${\gamma }_1$				0.0280	0.0300	0.6459				0.0180	0.1172	0.5222
		${\lambda }_2$	0.1277	0.2059	1.5701	0.0687	0.0645	0.8842	0.0084	0.2204	1.1111	0.0259	0.1768	0.7728
		${\gamma }_2$				0.0558	0.0634	0.9346				0.0104	0.1427	0.7172
		$\eta$	0.1272	0.2275	1.7405	0.1389	0.3630	2.2234	0.0057	0.2145	1.1831	0.0110	0.1942	1.2836
		$\rho$	0.0090	0.0048	0.2614	0.0046	0.0075	0.3283	−0.0018	0.0343	0.1908	−0.0004	0.0315	0.2054
		$\tau$	0.0101	0.0037	0.2282	0.0076	0.0057	0.2868	−0.0009	0.0294	0.1639	0.0002	0.0270	0.1768
	45	${\lambda }_1$	0.0441	0.0638	1.0008	0.0283	0.0523	0.9063	0.0183	0.1788	0.8210	0.0039	0.1674	0.7665
		${\gamma }_1$				0.0471	0.0222	0.5425				0.0074	0.0963	0.4385
		${\lambda }_2$	0.0378	0.1042	1.2110	0.0030	0.0352	0.7348	0.0167	0.2041	0.9805	0.0090	0.1430	0.6513
		${\gamma }_2$				0.0558	0.0445	0.7998				0.0152	0.1456	0.6440
		$\eta$	0.0794	0.1936	1.5939	0.0501	0.2564	1.9292	−0.0013	0.1999	1.1319	0.0377	0.1997	1.2518
		$\rho$	0.0029	0.0043	0.2466	−0.0058	0.0067	0.3019	−0.0026	0.0322	0.1837	0.0038	0.0314	0.1969
		$\tau$	0.0046	0.0032	0.2136	−0.0021	0.0048	0.2607	−0.0016	0.0276	0.1577	0.0038	0.0271	0.1702
140	105	${\lambda }_1$	0.0087	0.0249	0.6604	0.0050	0.0190	0.5663	0.0081	0.1395	0.5926	0.0127	0.1347	0.5381
		${\gamma }_1$				0.0169	0.0115	0.3761				0.0038	0.0748	0.3360
		${\lambda }_2$	0.0082	0.0579	0.8960	0.0038	0.0203	0.5167	0.0170	0.1791	0.7896	−0.0018	0.1269	0.4976
		${\gamma }_2$				0.0313	0.0235	0.5405				0.0052	0.1078	0.4757
		$\eta$	0.0063	0.0590	0.9646	−0.0029	0.1008	1.2127	−0.0063	0.1663	0.8323	0.0153	0.1891	0.9783
		$\rho$	−0.0024	0.0015	0.1556	−0.0062	0.0027	0.1957	−0.0027	0.0274	0.1357	0.0004	0.0295	0.1558
		$\tau$	−0.0012	0.0011	0.1336	−0.0041	0.0020	0.1681	−0.0019	0.0233	0.1164	0.0009	0.0255	0.1344
	130	${\lambda }_1$	0.0234	0.0232	0.5873	0.0190	0.0197	0.5295	0.0140	0.1288	0.5392	0.0115	0.1215	0.4914
		${\gamma }_1$				0.0097	0.0069	0.3188				0.0126	0.0783	0.2935
		${\lambda }_2$	0.0282	0.0412	0.7366	0.0140	0.0144	0.4423	0.0048	0.1571	0.6670	0.0117	0.1042	0.4210
		${\gamma }_2$				0.0158	0.0153	0.4655				0.0162	0.1008	0.4226
		$\eta$	0.0295	0.0553	0.9018	0.0287	0.0892	1.1167	0.0184	0.1754	0.7956	−0.0057	0.1965	0.9103
		$\rho$	0.0016	0.0014	0.1430	−0.0003	0.0022	0.1767	0.0011	0.0282	0.1274	0.0032	0.0320	0.1471
		$\tau$	0.0021	0.0010	0.1233	0.0008	0.0016	0.1524	0.0014	0.0242	0.1097	−0.0022	0.0273	0.1264
	140	${\lambda }_1$	−0.0051	0.0180	0.5517	−0.0074	0.0155	0.5143	0.0068	0.1222	0.5181	0.0025	0.1177	0.4839
		${\gamma }_1$				0.0087	0.0062	0.2941				−0.0002	0.0682	0.2700
		${\lambda }_2$	0.0001	0.0294	0.6767	−0.0017	0.0113	0.4185	0.0208	0.1496	0.6311	0.0001	0.0986	0.3988
		${\gamma }_2$				0.0195	0.0136	0.4431				0.0032	0.0986	0.4028
		$\eta$	0.0026	0.0487	0.8723	−0.0043	0.0655	1.0637	0.0350	0.1800	0.7851	0.0080	0.1827	0.8883
		$\rho$	−0.0024	0.0013	0.1408	−0.0045	0.0017	0.1721	0.0038	0.0280	0.1241	−0.0007	0.0290	0.1423
		$\tau$	−0.0014	0.0009	0.1208	−0.0030	0.0013	0.1477	0.0037	0.0242	0.1072	−0.0001	0.0250	0.1226

In terms of estimation performance, the Bayesian method generally surpasses MLE across all scenarios, producing smaller bias, lower MSE, and shorter interval estimates, indicating that prior information effectively improves estimation accuracy. Regarding model comparison, the BCW model exhibits greater flexibility and robustness than the BCE model due to its shape parameters, and achieves superior performance across all metrics despite its higher parameter dimensionality. As sample size increases from 45 to 140, estimation precision improves for all parameters, with reduced bias, MSE, and interval lengths. Complete samples yield the best estimates, and reducing the censoring rate from 40% to 20% consistently improves precision, with the Clayton copula parameter η being the most sensitive to censoring.

7. Empirical Study

To evaluate the applicability and effectiveness of the proposed model, a case study is conducted using the kidney disease patient data from McGilchrist and Aisbett [25]. The dataset comprises 30 patients on portable dialysis who required long-term treatment. Catheter site infections are common during dialysis, affecting treatment continuity and potentially causing severe complications. For each patient, two time points are recorded: the time from catheter implantation to the first infection ($X_1$ ) and the time from the second implantation to the second infection ($X_2$ ), both in days. The data are presented in Table 4.

Table 4. Inter-infection times in kidney disease patients.

$X_1$	8	23	22	447	30	24	7	511	53	15	7	141	96	149	536
	152	402	13	39	12	113	132	34	2	130	17	185	292	22	15
$X_2$	16	13	28	318	12	245	9	30	196	154	333	8	38	70	25
	362	24	66	46	40	201	156	30	25	26	4	117	114	159	108

Figure 5 shows a bivariate scatter plot matrix of $\left(X_1,X_2\right)$ . Diagonal panels display histograms with kernel density estimates; the off-diagonal panel shows the joint scatter plot with a locally smoothed fit. Both variables are right-skewed, with observations concentrated at lower values and extended tails, consistent with a Weibull distribution. The joint plot indicates clear dependence between $X_1$ and $X_2$ , violating independence and supporting the use of the Clayton copula.

Figure 5. Scatter plot matrix of the data.

Table 5 shows that the exponential distribution is rejected for $X_1$ at the 5% significance level, whereas the Weibull distribution passes the Kolmogorov-Smirnov (K-S) test for both $X_1$ and $X_2$ . The Akaike information criterion (AIC) and Hannan-Quinn information criterion (HQIC) values further confirm that the Weibull model outperforms the exponential model for $X_1$ . Therefore, the Weibull distribution is chosen as the marginal distribution for both infection interval datasets.

Based on the data characteristics and goodness-of-fit results, the Weibull distribution is chosen as the marginal distribution and the Clayton copula is used to model the dependence structure, i.e., the BCW model is adopted for further analysis.

Table 5. Goodness-of-fit test results for sample data.

			Est	StEr	K-S	p-value	AIC	HQIC
$X_1$	Exponential	$\lambda$	0.0083	0.0015	0.2577	0.0372	349.7309	352.1792
	Weibull	$\lambda$	0.0099	0.0026	0.1456	0.5484	346.9844	351.8809
		$\gamma$	0.7514	0.1055
$X_2$	Exponential	$\lambda$	0.0101	0.0018	0.1721	0.3364	337.7678	340.2160
	Weibull	$\lambda$	0.0104	0.0022	0.1458	0.5462	339.5149	344.4114
		$\gamma$	0.9319	0.1326

To evaluate the BCW model’s fitting performance on the kidney disease patient infection data, it is compared with several models from [26], including BCPL, BGPL, and BFPL. Table 6 summarizes the parameter estimates and goodness-of-fit metrics for each model. In the table, Est is the parameter estimate, and StEr is the standard error, measuring estimation precision. The negative log-likelihood $l$ indicates model fit, with smaller values being better. AIC and HQIC are model selection criteria that balance fit and complexity; smaller values indicate a better model.

Table 6. Parameter estimates and goodness-of-fit for different models.

		${\lambda }_1$	${\gamma }_1$	${\lambda }_2$	${\gamma }_2$	$\eta$	$l$	AIC	HQIC
BCW	Est	0.0100	0.7396	0.0105	0.9139	0.4792	338.4692	686.9384	689.1797
	StEr	0.0026	0.1044	0.0022	0.1325	0.3806
BFPL	Est	0.5538	0.1596	0.6633	0.1020	0.8184	338.9452	687.8905	694.8964
	StEr	0.0657	0.0528	0.0802	0.0398	1.0106
BCPL	Est	0.5489	0.1634	0.6546	0.1068	0.4187	338.5119	687.0239	694.0299
	StEr	0.0657	0.0541	0.0807	0.0419	0.3633
BGPL	Est	0.5525	0.1686	0.4333	0.4501	1.0463	371.7617	753.5235	760.5295
	StEr	0.0661	0.0607	0.0505	0.0984	0.1100
BFXL	Est	0.0162		0.0193		0.8682	361.3897	728.7794	732.9830
	StEr	0.0021		0.0026		0.7539
CBPL	Est	0.0157		0.0181		0.2274	359.1277	724.2553	728.4589
	StEr	0.0020		0.0025		0.1154
BGXL	Est	0.0181		0.0365		1.0219	432.1889	870.3779	874.5815
	StEr	0.0024		0.0031		0.0242
BFGMV	Est	0.7521	100.5887	0.9298	95.8087	0.3984	338.9352	687.8704	694.8764
	StEr	0.1054	25.8636	0.1323	19.9343	0.4979
BFGMF	Est	22.7499	0.7287	27.9035	0.8661	0.6121	340.7127	691.4254	698.4314
	StEr	6.1236	0.0979	6.2978	0.1175	0.5614
BFGMG	Est	0.6761	0.0056	0.9288	0.0094	0.3781	339.4908	688.9817	695.9877
	StEr	0.1478	0.0018	0.2089	0.0028	0.4839
BFGMGE	Est	0.6595	0.0062	0.9224	0.0096	0.4100	339.5462	689.0924	696.0984
	StEr	0.1481	0.0017	0.2174	0.0024	0.4849

Table 6 shows that the BCW model has the lowest negative log-likelihood, AIC, and HQIC among all competing models, indicating the best fit. Moreover, its parameter estimates have relatively small standard errors, suggesting stability. Overall, the BCW model exhibits the best fitting performance on this dataset, confirming the effectiveness of the proposed model.

To evaluate the BCW model’s estimation performance under Type-II censoring, three scenarios are considered: complete samples ($r=30$ ) and censored samples ($r=25,20$ ). Table 7 presents the MLE and Bayesian results. As censoring increases (effective sample size decreases from 30 to 20), the parameter estimates show small fluctuations, but their standard errors increase, indicating greater estimation uncertainty. Under the same censoring level, Bayesian standard errors are generally smaller than those of MLE, demonstrating higher estimation accuracy under small-sample Type-II censoring. Moreover, the copula parameter $\eta$ is positive across all censoring levels, confirming a positive dependence between the two infection intervals.

Table 7. MLE and Bayesian estimation results for the BCW model.

			$r=30$		$r=25$		$r=20$
			Est	StEr	Est	StEr	Est	StEr
BCW	MLE	${\gamma }_1$	0.7431	0.1051	0.7215	0.1143	0.7027	0.1268
		${\lambda }_1$	0.0101	0.0026	0.0102	0.0030	0.0097	0.0034
		${\gamma }_2$	0.9145	0.1329	0.8784	0.1361	0.8381	0.1534
		${\lambda }_2$	0.0106	0.0022	0.0109	0.0026	0.0110	0.0032
		$\eta$	0.4468	0.3814	0.4627	0.4835	0.3341	0.7586
	Bayes	${\gamma }_1$	0.7344	0.0912	0.7044	0.1011	0.6783	0.1044
		${\lambda }_1$	0.0102	0.0022	0.0103	0.0024	0.0099	0.0026
		${\gamma }_2$	0.8798	0.1191	0.8503	0.1254	0.7948	0.1265
		${\lambda }_2$	0.0106	0.0021	0.0107	0.0022	0.0109	0.0026
		$\eta$	0.6738	0.2368	0.7468	0.2772	0.7376	0.3090

Figure 6 shows the MCMC trace plots for each parameter under three Type-II censoring scenarios. The sampling sequences of all parameters display stable random fluctuations in the later iterations, with no obvious trends, clustering, or drift, indicating good convergence of the Markov chains. This confirms that posterior sampling is effective under Type-II censored data, and the Bayesian estimates are statistically reliable.

Figure 6. MCMC trace plots.

8. Conclusion

In this paper, we construct the bivariate Clayton Weibull (BCW) and bivariate Clayton exponential (BCE) models under Type-II censoring, with parameters estimated via maximum likelihood estimation (MLE) and Bayesian Markov chain Monte Carlo (MCMC) methods. Monte Carlo simulations reveal that Bayesian estimation consistently outperforms MLE, that larger sample sizes enhance estimation precision, and that the BCW model is more robust than the BCE model. When applied to kidney infection data, BCW yields the best fit, as evidenced by the lowest AIC and HQIC values among competing models. The estimated shape parameters for both infection times are below 1, indicating a declining risk of infection over time, and the positive Clayton copula parameter confirms a positive dependence between the two event times, thereby supporting patient risk stratification. The simulation results further inform the selection of sample sizes and follow-up cutoffs in clinical trial design. While the Clayton copula is specifically designed to capture lower-tail dependence and medical data may present more heterogeneous dependence patterns, the proposed estimation framework can be straightforwardly extended to other copula families to accommodate diverse data structures. Overall, the proposed methodology offers a practical and robust statistical framework for analyzing dependent survival data under Type-II censoring.