Bayesian and Non-Bayesian Inference to Bivariate Alpha Power Burr-XII Distribution with Engineering Application

Abstract

In this research, we present a new distribution, which is the bivariate alpha power Burr-XII distribution, based on the alpha power Burr-XII distribution. We thoroughly examine the key features of our newly developed bivariate model. We introduce a new class of bivariate models, which are built with the copula function. The statistical properties of the proposed distribution, such as conditional distributions, conditional expectations, marginal distributions, moment-generating functions, and product moments were studied. This was accomplished with two datasets of real data that came from two distinct devices. We employed Bayesian, maximum likelihood estimation, and least squares estimation strategies to obtain estimated points and intervals. Additionally, we generated bootstrap confidence intervals and conducted numerical analyses using the Markov chain Monte Carlo method. Lastly, we compared this novel bivariate distribution’s performance to earlier bivariate models, to determine how well it fit the real data.

1. Introduction

To achieve a reliable system, one method is to analyze the failures that occur in the production systems that are being observed. In bivariate modeling, the relationship between two variables is modeled. It enables one to discern between the interrelationship structure and the marginal distributions of the variable. Many authors have put forth bivariate models, analyzed them, and shown how broadly applicable they are across a variety of scientific fields. Over the past ten years, many articles have explored the idea of developing novel bivariate models based on a multiform copula function. A copula function is a helpful tool for describing the dependency structure of variables, regardless of the variable’s marginals. They have many different dependence models, such as the Gaussian, t, Clayton, and Gumbel copulas, which can be used to characterize different reliance patterns. Copula definitions are given by the unit square. Using a copula function, random numbers with a certain dependence structure can be generated. See Nelsen [1] and Flores [2] for examinations of several bivariate Weibull models derived from different copula functions. Verrill et al. [3] proposed a novel Gaussian–Weibull bivariate model. Many authors have since focused their statistical analysis on bivariate models. For example, El-Sherpieny and Almetwally [4] created a generalized Rayleigh distribution in two dimensions using the Clayton copula function as a basis and utilizing the Ali-Mikhail-Haq and Farlie-Gumbel-Morgenstern copulas, Qura et al. [5] showed a bivariate power Lomax distribution, and Almetwally and Muhammed [6] created a new Fréchet distribution in two dimensions. Almetwally et al. [7] proposed the use of the $FGM$ copula function to create a bivariate Weibull distribution. Ahmed et al. [8] introduced a new type of bivariate model using the copula function. By using various copulas, Samanthi and Sepanski [9] examined a few families of bivariate Kumaraswamy distributions. FGM copulas have been studied in different fields, such as finance (Cossette et al. [10]), economics (Patton [11]), and reliability engineering (Navarro et al. [12]; Navarro and Durante [13]).

The ability to study bivariate distributions with copulas, to capture asymmetric and non-linear connections between variables, makes them widely employed.Copulas provide a general class of dependence models that can be used to model a large range of dependent configurations. Copulas also simplify the computation of certain probability measures, such as the likelihood of a joint occurrence, the likelihood of two occurrences being dependent on one another, and the likelihood of crossing particular thresholds. Because copulas are flexible in simulating the dependence structure where non-linear and asymmetric dependencies between variables are often observed, they are an effective tool for financial and insurance applications. For a two-dimensional copula function, Sklar [14] defines $pdf$ and $cdf$. For a bivariate copula, if X and Y are two random variables for the distribution functions $F\left(x\right)$ and $F\left(y\right)$, then

$$F(x,y)=C(F(x;{\mathsf{\Lambda }}_1),F(y;{\mathsf{\Lambda }}_2))$$

(1)

and

$$f(x,y)=f(x;{\mathsf{\Lambda }}_1)f(y;{\mathsf{\Lambda }}_2)c(F(x;{\mathsf{\Lambda }}_1),F(y;{\mathsf{\Lambda }}_2)).$$

(2)

Other copula functions, such as the Farlie–Gumbel–Morgenstern copula, have been defined based on the aforementioned equations. This copula is a well-known parametric family copula, initially demonstrated by Gumbel [15]. The combination $cdf$ and $pdf$ of the Farlie–Gumbel–Morgenstern copula are given by

$$C(\tau ,\nu )=\tau \nu \{1+\delta [(1-\tau )(1-\nu )\left]\right\}$$

(3)

and

$$c(\tau ,\nu )=1+\delta (1-2\tau )(1-2\nu ),$$

(4)

where the variable’s parameters are represented by the vectors $\tau$ = $f(x;{\mathsf{\Lambda }}_1)$ and $\nu$ = $f(y;{\mathsf{\Lambda }}_2)$, and ${\mathsf{\Lambda }}_1$ and ${\mathsf{\Lambda }}_2$ denote the copula parameters that assume values in $[0,1]$.

The Burr-XII distribution introduced by Burr [16] has attracted attention for its remarkable applications in a variety of areas, including reliability, failure time modeling, and acceptance sampling schemes. For example, Wang and Keats [17] employed the maximum likelihood method to estimate the parameters of the Burr-XII distribution using interval estimators. The Burr-XII distribution was applied by Abdel-Ghaly et al. [18] in software reliability growth modeling, microelectronics, and reliability. A statistical analysis of the Burr-XII distribution and its relation to other distributions were explored by Zimmer et al. [19]. Moors [20] established a confidence range for the form parameter under the Bur-XII distribution’s failure censored plan. Maximum likelihood was the method employed by Wu and Yu [21] to calculate the point estimator of the parameters of the Burr-XII distribution. Bayesian predictive density of order statistics based on finite mixture models was studied by AL-Hussaini [22]. AL-Hussaini and Ahmad [23] developed a Bayesian interval prediction of future records and, according to Kumar et al. [24], derived the Marshall–Olkin extended Burr-XII distribution ratio and inverse moments. One of these generalizations is the alpha power Burr-XII (APB-XII) distribution [25], which has provided better results with real datasets in comparison to many versions of the exponential distribution. A random variable X is said to have an alpha power Burr-XII (APB-XII) distribution with parameters $\lambda$ $>0$ and $\beta$ $>0$ as a shape parameter, and $\alpha$ $>0$ as a scale parameter, if its cdf and pdf are given by

$$F\left(x\right)={\displaystyle \frac{{\alpha }^{1-{(1+x^{\beta })}^{-\lambda }}-1}{(\alpha -1)}},$$

(5)

and

$$f\left(x\right)={\displaystyle \frac{log\left[\alpha \right]}{(\alpha -1)}}\lambda \beta x^{\beta -1}{(1+x^{\beta })}^{-(1+\lambda )}{\alpha }^{1-{(1+x^{\beta })}^{-\lambda }}.$$

(6)

This paper introduces the BAPB-XII distribution, an innovative model derived from the alpha power Burr-XII distribution, enhancing the available options in bivariate statistical analysis. It thoroughly explores the statistical properties of this distribution, including conditional and marginal distributions, moment-generating functions, and product moments, which are crucial for understanding its behavior across different applications. To achieve robust parameter estimation, the study employs multiple estimation techniques, such as Bayesian methods, maximum likelihood estimation, and least squares estimation, adding flexibility and reliability. The practical applicability of the BAPB-XII model was demonstrated through real-world data analysis, specifically focusing on Burr measurements from different machines and sheet failure rates, showcasing the model’s effectiveness on life testing data. Additionally, this paper conducted a comparative performance evaluation with existing bivariate models, underscoring the BAPB-XII distribution’s advantages and suitability for fitting real data, thereby providing valuable insights for researchers and practitioners.

The BAPB-XII distribution finds numerous applications across various fields. It is particularly suitable for life testing data analysis, playing a critical role in reliability engineering and quality control by effectively modeling failure rates. This distribution proves practical in assessing failure rates, as illustrated through its application in the evaluation of sheet failures, which is essential for understanding real-world failure mechanisms. Additionally, the BAPB-XII distribution is designed to handle bivariate data, making it valuable in contexts where two interdependent variables are analyzed, such as in engineering and manufacturing. Its statistical characteristics allow for in-depth exploration of properties, enabling its application in fields requiring advanced statistical modeling. Furthermore, comparative analysis with traditional bivariate distributions highlights the potential of the BAPB-XII distribution for broader applications in statistical analysis and modeling.

This paper is structured as follows. The distribution of BAPB-XII is defined in Section 2. The statistical characteristics of the BAPB-XII distribution are given in Section 3. The maximum likelihood and Bayesian estimations are performed in Section 4. Section 5 introduces the simulation study. In Section 6, the application on real data is reported. The manuscript ends with the conclusions Section.

2. The Bivariate Alpha Power Burr-XII Model

The bivariate alpha power Burr-XII (BAPB-XII) model is stated in Equations (1) and (2), and the copula is specified in Equations (3) and (4), with x and y being two random variables after the APB-XII function as in Equations (5) and (6), used to obtain the bivariate distribution. The bivariate APB-XII distribution with the Farlie–Gumbel–Morgenstern copula function’s joint $pdf$ and $cdf$ are provided as follows

$$F(x,y)=u\left(x\right)v\left(y\right)[1+\delta [1-u\left(x\right)\left]\right[1-v\left(y\right)\left]\right]$$

(7)

and

$$\begin{array}{cc}\hfill f(x,y)& ={\lambda }_1{\lambda }_2{\beta }_1{\beta }_2{\displaystyle \frac{log\left[{\alpha }_1\right]log\left[{\alpha }_2\right]}{({\alpha }_1-1)({\alpha }_2-1)}}{x}^{{\beta }_1-1}{y}^{{\beta }_2-1}{(1+x^{{\beta }_1})}^{-(1+{\lambda }_1)}{(1+y^{{\beta }_2})}^{-(1+{\lambda }_2)}\hfill \\ \hfill & {\alpha }_1^{1-{(1+x^{{\beta }_1})}^{-{\lambda }_1}}{\alpha }_2^{1-{(1+y^{{\beta }_2})}^{-{\lambda }_2}}\left[1+\delta \left(1-2u\left(x\right)\right)\left(1-2v\left(y\right)\right)\right],\hfill \end{array}$$

(8)

where ${\alpha }_1,{\beta }_1,{\lambda }_1,{\alpha }_2,{\lambda }_2,{\beta }_2>0$, $\delta \in (-1,1)$, $x,y>0$

$$u\left(x\right)={\displaystyle \frac{{\alpha }_1^{1-{(1+x^{{\beta }_1})}^{-{\lambda }_1}}-1}{{\alpha }_1-1}}\mathrm{and}v\left(y\right)={\displaystyle \frac{{\alpha }_2^{1-{(1+y^{{\beta }_2})}^{-{\lambda }_2}}-1}{{\alpha }_2-1}}.$$

Figure 1 illustrates various shapes of the pdf of BAPB-XII in Equation (7) for different parameter values, demonstrating the flexibility of the joint pdf in modeling bivariate skewed data. One can observe that the joint pdf of the BAPB-XII distribution has different shapes depending on the values of its parameters, which reflects the ability of the joint pdf to model bivariate skewed data.

3. Properties of the New Model

This section explores and presents several significant statistical properties of the BAPB-XII distribution.

3.1. Marginal and Conditional Distributions

The marginal density functions of X and Y are

$$f_1(x;{\alpha }_1,{\beta }_1,{\lambda }_1)={\lambda }_1{\beta }_1{\displaystyle \frac{log\left[{\alpha }_1\right]}{{\alpha }_1-1}}{x}^{{\beta }_1-1}{(1+x^{{\beta }_1})}^{-({\lambda }_1+1)}{\alpha }_1^{1-{(1+x^{{\beta }_1})}^{-{\lambda }_1}},x>0,{\alpha }_1,{\beta }_1,{\lambda }_1>0$$

(9)

and

$$f_2(y;{\alpha }_2,{\beta }_2,{\lambda }_2)={\lambda }_2{\beta }_2{\displaystyle \frac{log\left[{\alpha }_2\right]}{{\alpha }_2-1}}{y}^{{\beta }_2-1}{(1+y^{{\beta }_2})}^{-({\lambda }_2+1)}{\alpha }_2^{1-{(1+y^{{\beta }_2})}^{-{\lambda }_2}},y>0,{\alpha }_2,{\beta }_2,{\lambda }_2>0,$$

(10)

respectively.

The distribution of the conditional probability of X given Y is

$$f(x\mid y)={\lambda }_1{\beta }_1{\displaystyle \frac{log\left[{\alpha }_1\right]}{{\alpha }_1-1}}{x}^{{\beta }_1-1}{(1+x^{{\beta }_1})}^{-({\lambda }_1+1)}{\alpha }_1^{1-{(1+x^{{\beta }_1})}^{-{\lambda }_1}}[1+\delta [1-2u\left(x\right)][1-2v\left(y\right)]],$$

(11)

So, the conditional probability of Y given x is defined as

$$f(y\mid x)={\lambda }_2{\beta }_2{\displaystyle \frac{log\left[{\alpha }_2\right]}{{\alpha }_2-1}}{y}^{{\beta }_2-1}{(1+y^{{\beta }_2})}^{-({\lambda }_2+1)}{\alpha }_2^{1-{(1+y^{{\beta }_2})}^{-{\lambda }_2}}[1+\delta [1-2u\left(x\right)][1-2v\left(y\right)]].$$

(12)

The conditional $cdf$ of X given Y, and $Y$ given $X$, are given, respectively, by

$$F(x\mid y)=u\left(x\right)[1+\delta [1-u\left(x\right)\left]\right[1-v\left(y\right)\left]\right]$$

(13)

and

$$F(y\mid x)=v\left(y\right)[1+\delta [1-u\left(x\right)\left]\right[1-v\left(y\right)\left]\right].$$

(14)

The conditional method which is outlined below can be used to create a bivariate sample from the APB-XII distribution.

-

Generate independent variables Z and W from uniform ($0,1$) distribution.

-

Let $X=Q_{APB-XII}\left(Z\right)={\left({\left(1-\left[{\displaystyle \frac{log\left(Z\right(\alpha -1)+1)}{log\left(\alpha \right)}}\right]\right)}^{{\displaystyle \frac{-1}{\lambda }}}-1\right)}^{{\displaystyle \frac{1}{\beta }}}.$

-

Let $F\left(y\right|x)=W$ to find y by a numerical method.

-

To obtain $(x_i,y_i)$, $i=1,\dots ,n$, repeat the above steps n times.

3.2. Product Moments

The mth and lth product moments of a BAPB-XII distribution for the random variables X and Y can be calculated as follows:

$${\mu }_{m,l}=E\left(x^m{y}^l\right)=\underset{0}{\overset{\infty }{\int }}\underset{0}{\overset{\infty }{\int }}{x}^m{y}^{l}f(x,y)dxdy,$$

(15)

then,

$$\begin{array}{cc}\hfill {\stackrel{`}{\mu }}_{m,l}& =\varphi \sum _{k=0}^1\sum _{p=0}^k\sum _{q=0}^k\sum _{c=0}^p\sum _{G=0}^q\sum _{s=0}^{\infty }\sum _{j=0}^{\infty }\sum _{Q=0}^s\sum _{r=0}^j\xi \Omega {\delta }^k{\left({\displaystyle \frac{1}{{\alpha }_1-1}}\right)}^p{\left({\displaystyle \frac{1}{{\alpha }_2-1}}\right)}^q\hfill \\ \hfill & \left[{\displaystyle \frac{{(log\left({\alpha }_1\right))}^s}{s!}}\right]\left[{\displaystyle \frac{{(log\left({\alpha }_2\right))}^z}{j!}}\right]B\left({\displaystyle \frac{m}{{\beta }_1}}+1,-{\lambda }_1(Q+1)\right)B\left({\displaystyle \frac{l}{{\beta }_2}}+1,-{\lambda }_2(r+1)\right),\hfill \end{array}$$

(16)

where$\Omega =$$\left({\displaystyle \genfrac{}{}{0pt}{}{1}{k}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{p}{c}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{q}{G}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{Q}{s}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{r}{Z}}\right){(c+1)}^s{(G+1)}^j$, $\varphi ={\lambda }_1{\lambda }_2{\beta }_1{\beta }_2{\displaystyle \frac{log\left[{\alpha }_1\right]log\left[{\alpha }_2\right]}{({\alpha }_1-1)({\alpha }_2-1)}}$ and $\xi ={(-2)}^{p+q}{(-1)}^{{\displaystyle \frac{m+{\beta }_1}{{\beta }_1}}}{(-1)}^{p-G+p-c+Q+r}{(-1)}^{{\displaystyle \frac{l+{\beta }_2}{{\beta }_2}}}.$

3.3. Moment Generating Function

The random variables X and Y have a $pdf$ that be defined in Equation (8), then the moment generating function ($MGF$) for X and Y is obtained as

$$M_{x,y}(t_1,t_2)=\sum _{m=0}^{\infty }\sum _{l=0}^{\infty }{\displaystyle \frac{t_1^m{t}_2^l}{m!t!}}\underset{0}{\overset{\infty }{\int }}\underset{0}{\overset{\infty }{\int }}{x}^m{y}^{l}f(x,y)dxdy,$$

(17)

then,

$$M_{x,y}(t_1,t_2)=\sum _{m=0}^{\infty }\sum _{l=0}^{\infty }{\displaystyle \frac{t_1^m{t}_2^l}{m!t!}}{\stackrel{`}{\mu }}_{m,l}.$$

(18)

Using the product moment defined in the previous section in Equation (18), the $MGF$of $X,Y$ is written as

$$\begin{array}{cc}\hfill M_{x,y}(t_1,t_2)& =\sum _{m=0}^{\infty }\sum _{l=0}^{\infty }\sum _{k=0}^1\sum _{p=0}^k\sum _{q=0}^k\sum _{c=0}^p\sum _{G=0}^q\sum _{s=0}^{\infty }\sum _{z=0}^{\infty }\sum _{Q=0}^s\sum _{r=0}^z\varphi \Omega {\delta }^k{\displaystyle \frac{t_1^m{t}_2^l}{m!t!}}\hfill \\ \hfill & \times {\left({\displaystyle \frac{1}{{\alpha }_1-1}}\right)}^p{\left({\displaystyle \frac{1}{{\alpha }_2-1}}\right)}^q\left[{\displaystyle \frac{{(log\left({\alpha }_1\right))}^s}{s!}}\right]\left[{\displaystyle \frac{{(log\left({\alpha }_2\right))}^z}{z!}}\right]\hfill \\ \hfill & \times B\left({\displaystyle \frac{m}{{\beta }_1}}+1,-{\lambda }_1(Q+1)\right)B\left({\displaystyle \frac{l}{{\beta }_2}}+1,-{\lambda }_2(r+1)\right),\hfill \end{array}$$

(19)

where $\Omega =$ $\left({\displaystyle \genfrac{}{}{0pt}{}{1}{k}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{p}{c}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{q}{G}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{Q}{s}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{r}{Z}}\right){(G+1)}^z{(c+1)}^s,$ $\varphi ={\lambda }_1{\lambda }_2{\beta }_1{\beta }_2{\displaystyle \frac{log\left[{\alpha }_1\right]log\left[{\alpha }_2\right]}{({\alpha }_1-1)({\alpha }_2-1)}}$ and $\xi ={(-2)}^{p+q}$${(-1)}^{{\displaystyle \frac{m+{\beta }_1}{{\beta }_1}}}{(-1)}^{p-G+p-c+Q+r}{(-1)}^{{\displaystyle \frac{l+{\beta }_2}{{\beta }_2}}}$.

3.4. Survival and Hazard Rate Functions

The joint survival function was proved by Osmetti and Chiodini [26], who applied X and Y as random variables with survival functions ($1-F\left(x\right))$ and ($1-F\left(y\right))$, which is more practical.

The survival functions of the marginal distributions are defined as

$$S\left(x\right)=1-F\left(x\right)={\displaystyle \frac{{\alpha }_1-{\alpha }_1^{1-{(1+x^{{\beta }_1})}^{-{\lambda }_1}}}{{\alpha }_1-1}}$$

(20)

and

$$S\left(y\right)=1-F\left(y\right)={\displaystyle \frac{{\alpha }_2-{\alpha }_2^{1-{(1+x^{{\beta }_2})}^{-{\lambda }_2}}}{{\alpha }_2-1}}.$$

(21)

The joint survival function for a copula is expressed as

$$S(x,y)=C\left(S\right(x),S(y\left)\right).$$

(22)

Hence, the survival function of the BAPB-XII distribution is

$$\begin{array}{cc}\hfill S(x,y)& =\left({\displaystyle \frac{{\alpha }_1-{\alpha }_1^{1-{(1+x^{{\beta }_1})}^{-{\lambda }_1}}}{{\alpha }_1-1}}\right)\left({\displaystyle \frac{{\alpha }_2-{\alpha }_2^{1-{(1+y^{{\beta }_2})}^{-{\lambda }_2}}}{{\alpha }_2-1}}\right)\hfill \\ \hfill & \left[1+\delta \left(1-{\displaystyle \frac{{\alpha }_1-{\alpha }_1^{1-{(1+x^{{\beta }_1})}^{-{\lambda }_1}}}{{\alpha }_1-1}}\right)\left(1-{\displaystyle \frac{{\alpha }_2-{\alpha }_2^{1-{(1+y^{{\beta }_2})}^{-{\lambda }_2}}}{{\alpha }_2-1}}\right)\right].\hfill \end{array}$$

(23)

The bivariate hazard rate function is given by

$$h(x,y)={\displaystyle \frac{f(x,y)}{S(x,y)}},$$

(24)

then, it is obtained as

$$\begin{array}{cc}\hfill h(x,y)=& {\lambda }_1{\lambda }_2{\beta }_1{\beta }_{2}log\left[{\alpha }_1\right]log\left[{\alpha }_2\right]x^{{\beta }_1-1}{y}^{{\beta }_2-1}{\alpha }_1^{1-{(1+x^{{\beta }_1})}^{-{\lambda }_1}}{\alpha }_2^{1-{(1+y^{{\beta }_2})}^{-{\lambda }_2}}{(1+x^{{\beta }_1})}^{-(1+{\lambda }_1)}\hfill \\ \hfill \times & {(1+y^{{\beta }_2})}^{-(1+{\lambda }_2)}\left[1+\delta \left(1-2u\left(x\right)\right)\left(1-2v\left(y\right)\right)\right]\hfill \\ \hfill \times & {\left[1+\delta \left(1-{\displaystyle \frac{{\alpha }_1-{\alpha }_1^{1-{(1+x^{{\beta }_1})}^{-{\lambda }_1}}}{{\alpha }_1-1}}\right)\left(1-{\displaystyle \frac{{\alpha }_2-{\alpha }_2^{1-{(1+y^{{\beta }_2})}^{-{\lambda }_2}}}{{\alpha }_2-1}}\right)\right]}^{-1}\hfill \\ \hfill \times & {\left[\left({\alpha }_1-{\alpha }_1^{1-{(1+x^{{\beta }_1})}^{-{\lambda }_1}}\right)\left({\alpha }_2-{\alpha }_2^{1-{(1+y^{{\beta }_2})}^{-{\lambda }_2}}\right)\right]}^{-1}.\hfill \end{array}$$

(25)

4. Estimation Methods

4.1. Maximum Likelihood Estimation

This subsection derives the maximum likelihood estimate ($MLE$) for the parameter vector $\Theta =$(${\alpha }_1,{\beta }_1,{\lambda }_1,{\alpha }_2,{\beta }_2,{\lambda }_2,\delta ).$ Assume that the n observed values from the BAPB-XII distribution are $(x_{1:n},y_{1:n})$, $(x_{2:n}$, $y_{2:n}),\dots ,(x_{n:n},y_{n:n})$. Following the approach of Kim et al. [27], who discussed the likelihood function for a bivariate model, the likelihood function for a complete sample is given as

$$\begin{array}{cc}\hfill L(\Theta )& =\prod _{i=1}^n{f}_{x,y}(x_i,y_i)=\prod _{i=1}^n{f}_x\left(x_i\right)f_y\left(y_i\right)c(F_x\left(x_i\right),F_y\left(y_i\right)),\hfill \end{array}$$

(26)

then,

$$\begin{array}{cc}\hfill L(\Theta )& =\prod _{i=1}^n{\lambda }_1{\lambda }_2{\beta }_1{\beta }_2{\displaystyle \frac{log\left[{\alpha }_1\right]log\left[{\alpha }_2\right]}{({\alpha }_1-1)({\alpha }_2-1)}}{x}_i^{{\beta }_1-1}{y}_i^{{\beta }_2-1}{(1+x_i^{{\beta }_1})}^{-(1+{\lambda }_1)}{(1+y_i^{{\beta }_2})}^{-(1+{\lambda }_2)}\hfill \\ \hfill & {\alpha }_1^{1-{(1+x_i^{{\beta }_1})}^{-{\lambda }_1}}{\alpha }_2^{1-{(1+y_i^{{\beta }_2})}^{-{\lambda }_2}}\left[1+\delta \left(1-2u\left(x_i\right)\right)\left(1-2v\left(y_i\right)\right)\right],\hfill \end{array}$$

(27)

where $u\left(x_i\right)={\displaystyle \frac{{\alpha }_1^{1-{(1+x_i^{{\beta }_1})}^{-{\lambda }_1}}-1}{{\alpha }_1-1}}$, $v\left(y_i\right)={\displaystyle \frac{{\alpha }_2^{1-{(1+y_i^{{\beta }_2})}^{-{\lambda }_2}}-1}{{\alpha }_2-1}}.$

The log-likelihood function $l(\Theta )$can be written as follows:

$$\begin{array}{cc}\hfill l(\Theta )& =nlog\left({\lambda }_1\right)+nlog\left({\lambda }_2\right)+nlog\left({\beta }_1\right)+nlog\left({\beta }_2\right)+nlog(log\left[{\alpha }_1\right])+nlog(log\left[{\alpha }_2\right])\hfill \\ \hfill & -nlog({\alpha }_1-1)-nlog({\alpha }_2-1)+({\beta }_1-1)\sum _{i=1}^{n}log\left(x_i\right)+({\beta }_2-1)\sum _{i=1}^{n}log\left(y_i\right)\hfill \\ \hfill & -(1+{\lambda }_1)\sum _{i=1}^{n}log(1+x_i^{{\beta }_1})-(1+{\lambda }_2)\sum _{i=1}^{n}log(1+y_i^{{\beta }_2})+log\left({\alpha }_1\right)\sum _{i=1}^n\left(1-{(1+x_i^{{\beta }_1})}^{-{\lambda }_1}\right)\hfill \\ \hfill & +log\left({\alpha }_2\right)\sum _{i=1}^n\left(1-{(1+y_i^{{\beta }_2})}^{-{\lambda }_2}\right)+\sum _{i=1}^{n}log\left(1+\delta [1-2u\left(x_i\right)][1-2v\left(y_i\right)]\right).\hfill \end{array}$$

(28)

The MLEs, $\widehat{\Theta }=(\widehat{{\alpha }_1},\widehat{{\beta }_1},\widehat{{\lambda }_1},\widehat{{\alpha }_2},\widehat{{\beta }_2},\widehat{{\lambda }_2},\widehat{\delta })$, of $\Theta$ are obtained by simultaneously solving a set of nonlinear equations. These equations are derived by partially differentiating Equation (28) with respect to each parameter and setting each resulting expression to zero, as follows:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial l}{\partial {\alpha }_1}}& ={\displaystyle \frac{n}{{\alpha }_{1}log\left({\alpha }_1\right)}}+{\alpha }_1^{-1}\sum _{i=1}^n\left(1-{(1+x^{{\beta }_1})}^{-{\lambda }_1}\right)-{\displaystyle \frac{n}{{\alpha }_1-1}}\hfill \\ \hfill & -\sum _{i=1}^n\left[{\alpha }_1^{-{(1+x_i^{{\beta }_1})}^{-{\lambda }_1}}\left(1-{(1+x_i^{{\beta }_1})}^{-{\lambda }_1}\right)-u\left(x_i\right)\right]\left[{\displaystyle \frac{2\delta {({\alpha }_1-1)}^{-1}\left(1-2v\left(y_i\right)\right)}{1+\delta [1-2u\left(x_i\right)][1-2v\left(y_i\right)]}}\right]=0,\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial l}{\partial {\alpha }_2}}& ={\displaystyle \frac{n}{{\alpha }_{2}log\left({\alpha }_2\right)}}+{\alpha }_2^{-1}\sum _{i=1}^n\left(1-{(1+y_i^{{\beta }_2})}^{-\lambda 2}\right)-{\displaystyle \frac{n}{{\alpha }_2-1}}\hfill \\ \hfill & -\sum _{i=1}^n\left[{\alpha }_2^{-{(1+y_i^{{\beta }_2})}^{-{\lambda }_2}}\left(1-{(1+y_i^{{\beta }_2})}^{-{\lambda }_2}\right)-v\left(y_i\right)\right]\left[{\displaystyle \frac{2\delta {({\alpha }_2-1)}^{-1}\left(1-2u\left(x_i\right)\right)}{1+\delta [1-2u\left(x_i\right)][1-2v\left(y_i\right)]}}\right]=0,\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial l}{\partial {\beta }_1}}& ={\displaystyle \frac{n}{{\beta }_1}}+\sum _{i=1}^{n}log\left(x_i\right)+log\left({\alpha }_1\right)\sum _{i=1}^n{\lambda }_1{x}_i^{{\beta }_1}log\left(x_i\right){\left(x_i^{{\beta }_1}+1\right)}^{-{\lambda }_1-1}-\left({\lambda }_1+1\right)\sum _{i=1}^n{\displaystyle \frac{x_i^{{\beta }_1}log\left(x_i\right)}{x_i^{{\beta }_1}+1}}\hfill \\ \hfill & -\sum _{i=1}^n{x}_i^{{\beta }_1}log\left(x_i\right){\left(x_i^{{\beta }_1}+1\right)}^{-{\lambda }_1-1}{\alpha }_1^{1-{\left(x_i^{{\beta }_1}+1\right)}^{-{\lambda }_1}}\left[{\displaystyle \frac{2\delta {\lambda }_{1}log\left({\alpha }_1\right){({\alpha }_1-1)}^{-1}\left(1-2v\left(y_i\right)\right)}{1+\delta [1-2u\left(x_i\right)][1-2v\left(y_i\right)]}}\right]=0,\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial l}{\partial {\beta }_2}}& ={\displaystyle \frac{n}{{\beta }_2}}+\sum _{i=1}^{n}log\left(y_i\right)+log\left({\alpha }_2\right)\sum _{i=1}^n{\lambda }_2{y}_i^{{\beta }_2}log\left(y_i\right){\left(y_i^{{\beta }_2}+1\right)}^{-{\lambda }_2-1}-\left({\lambda }_2+1\right)\sum _{i=1}^n{\displaystyle \frac{y_i^{{\beta }_2}log\left(y_i\right)}{y_i^{{\beta }_2}+1}}\hfill \\ \hfill & -\sum _{i=1}^n{y}_i^{{\beta }_2}log\left(y_i\right){\left(y_i^{{\beta }_2}+1\right)}^{-{\lambda }_2-1}{\alpha }_2^{1-{\left(y_i^{{\beta }_2}+1\right)}^{-{\lambda }_2}}\left[{\displaystyle \frac{2\delta {\lambda }_{2}log\left({\alpha }_2\right){({\alpha }_2-1)}^{-1}\left(1-2u\left(x_i\right)\right)}{1+\delta [1-2u\left(x_i\right)][1-2v\left(y_i\right)]}}\right]=0,\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial l}{\partial {\lambda }_1}}& ={\displaystyle \frac{n}{{\lambda }_1}}-\sum _{i=1}^{n}log\left(x_i^{{\beta }_1}+1\right)+log\left({\alpha }_1\right)\sum _{i=1}^n{\left(x_i^{{\beta }_1}+1\right)}^{-{\lambda }_1}log\left(x_i^{{\beta }_1}+1\right)\hfill \\ \hfill & -\sum _{i=1}^{n}log\left({\alpha }_1\right){\left(x_i^{{\beta }_1}+1\right)}^{-{\lambda }_1}log\left(x_i^{{\beta }_1}+1\right){\alpha }_1^{1-{\left(x_i^{{\beta }_1}+1\right)}^{-{\lambda }_1}}\left[{\displaystyle \frac{2\delta {({\alpha }_1-1)}^{-1}\left(1-2v\left(y_i\right)\right)}{1+\delta [1-2u\left(x_i\right)][1-2v\left(y_i\right)]}}\right]=0,\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial l}{\partial {\lambda }_2}}& ={\displaystyle \frac{n}{{\lambda }_2}}-\sum _{i=1}^{n}log\left(y_i^{{\beta }_2}+1\right)+log\left({\alpha }_2\right)\sum _{i=1}^n{\left(y_i^{{\beta }_2}+1\right)}^{-{\lambda }_2}log\left(y_i^{{\beta }_2}+1\right)\hfill \\ \hfill & -\sum _{i=1}^{n}log\left({\alpha }_2\right){\left(y_i^{{\beta }_2}+1\right)}^{-{\lambda }_2}log\left(y_i^{{\beta }_2}+1\right){\alpha }_2^{1-{\left(y_i^{{\beta }_2}+1\right)}^{-{\lambda }_2}}\left[{\displaystyle \frac{2\delta {({\alpha }_2-1)}^{-1}\left(1-2u\left(x_i\right)\right)}{1+\delta [1-2u\left(x_i\right)][1-2v\left(y_i\right)]}}\right]=0,\hfill \end{array}$$

(34)

and

$${\displaystyle \frac{\partial l}{\partial \delta }}=\sum _{i=1}^n{\displaystyle \frac{[1-2u\left(x_i\right)][1-2v\left(y_i\right)]}{1+\delta [1-2u\left(x_i\right)][1-2v\left(y_i\right)]}}.$$

(35)

It is necessary to use numerical methods to solve Equations (29)–(35), such as the Newton–Raphson methodology, implemented in mathematical software packages, to find a solution. Alternatively, direct approaches for maximizing the log-likelihood functions in Equation (28) are available in some statistical and mathematical packages such as Mathematica software (NMaximize and FindMaximum functions).

Confidence Interval

The confidence intervals are based on the asymptotic normality of the MLEs for the parameters. To achieve this, we take the negative second derivatives of the log-likelihood function with respect to the parameter vector $\Theta$, which allows us to construct a Fisher information matrix $I(\Theta )$. The asymptotic variance–covariance matrix of the MLE estimates for the parameters is then derived from this Fisher information matrix.

This method assumes that the asymptotic variance–covariance matrix for the parameter vector $\Theta$ is supplied by

$$I\left(\widehat{\Theta }\right)=\left[\begin{array}{ccccccc}{I}_{{\widehat{\alpha }}_1{\widehat{\alpha }}_1}& & & & & & \\ I_{{\widehat{\beta }}_1{\widehat{\alpha }}_1}& I_{{\widehat{\beta }}_1{\widehat{\beta }}_1}& & & & & \\ I_{{\widehat{\lambda }}_1{\widehat{\alpha }}_1}& I_{{\widehat{\lambda }}_1{\widehat{\beta }}_1}& I_{{\widehat{\lambda }}_1{\widehat{\lambda }}_1}& & & & \\ I_{{\widehat{\alpha }}_2{\widehat{\alpha }}_1}& I_{{\widehat{\alpha }}_2{\widehat{\beta }}_1}& I_{{\widehat{\alpha }}_2{\widehat{\lambda }}_1}& I_{{\widehat{\alpha }}_2{\widehat{\alpha }}_2}& & & \\ I_{{\widehat{\beta }}_2{\widehat{\alpha }}_1}& I_{{\widehat{\beta }}_2{\widehat{\beta }}_1}& I_{{\widehat{\beta }}_2{\widehat{\lambda }}_1}& I_{{\widehat{\beta }}_2{\widehat{\alpha }}_2}& I_{{\widehat{\beta }}_2{\widehat{\beta }}_2}& & \\ I_{{\widehat{\lambda }}_2{\widehat{\alpha }}_1}& I_{{\widehat{\lambda }}_2{\widehat{\beta }}_1}& I_{{\widehat{\lambda }}_2{\widehat{\lambda }}_1}& I_{{\widehat{\lambda }}_2{\widehat{\alpha }}_2}& I_{{\widehat{\lambda }}_2{\widehat{\beta }}_2}& I_{{\widehat{\lambda }}_2{\widehat{\lambda }}_2}& \\ I_{\widehat{\delta }{\widehat{\alpha }}_1}& I_{\widehat{\delta }{\widehat{\beta }}_1}& I_{\widehat{\delta }{\widehat{\lambda }}_1}& I_{\widehat{\delta }{\widehat{\alpha }}_2}& I_{\widehat{\delta }{\widehat{\beta }}_2}& I_{\widehat{\delta }{\widehat{\lambda }}_2}& I_{\widehat{\delta }\widehat{\delta }}\end{array}\right],$$

(36)

where the variance–covariance matrix is given by $V\left(\widehat{\Theta }\right)=I^{-1}\left(\widehat{\Theta }\right)$. With the asymptotic normality of the MLE, we can construct a $100(1-\gamma )\%$ asymptotic confidence interval for the parameter vector $\Theta$, as follows:

$${\widehat{\alpha }}_j\pm Z_{\gamma }\sqrt{var\left({\widehat{\alpha }}_j\right)},{\widehat{\beta }}_j\pm Z_{\gamma }\sqrt{var\left({\widehat{\beta }}_j\right)},{\widehat{\lambda }}_j\pm Z_{\gamma }\sqrt{var\left({\widehat{\lambda }}_j\right)},\widehat{\delta }\pm Z_{\gamma }\sqrt{var\left(\widehat{\delta }\right)},j=1,2.$$

(37)

In this expression, $Z_{\gamma }$ represents the percentile of the standard normal distribution corresponding to the right tail with a probability of ${\displaystyle \frac{\gamma }{2}}$. Note that $var\left(\widehat{{\Theta }_k}\right)$ is the square root of the k$th$ diagonal element of the matrix $I^{-1}\left(\widehat{\Theta }\right)$.

4.2. Bayesian Estimation

In this subsection, the Bayesian estimation method is applied to estimate the parameters of the BAPB-XII distribution. For parameter estimation based on a complete sample, we use the symmetric square error loss function, assuming the parameters of the BAPB-XII distribution are independent. We consider gamma-independent priors for the parameters $({\alpha }_j,{\beta }_j,{\lambda }_j,\delta )$for $j=1,2,$ defined as follows:

$$\begin{array}{cc}\hfill \pi \left({\alpha }_j\right)& ={\alpha }_j^{q_j-1}{e}^{-w_j{\alpha }_j}{\alpha }_j,q_j,w_j>0,\hfill \\ \hfill \pi \left({\beta }_j\right)& ={\beta }_j^{r_j-1}{e}^{-u_j{\beta }_j}{\beta }_j,r_j,u_j>0\hfill \\ \hfill & \mathrm{and}\hfill \\ \hfill \pi \left({\lambda }_j\right)& ={\lambda }_j^{o_j-1}{e}^{-p_j{\lambda }_j}{\lambda }_j,o_j,p_j>0,\hfill \end{array}$$

(38)

where the hyperparameters $q_j,w_j,r_j,u_j,o_j$and $p_j$are chosen based on prior assumptions. For the prior distribution of the copula parameter $\delta$, we assign a uniform prior over the interval $(-1,1)$.

The following is the joint prior:

$$\pi (\Theta )\propto {\alpha }_1^{q_1-1}{\beta }_1^{r_1-1}{\lambda }_l^{o_l-1}{e}^{-w_1{\alpha }_l-u_1{\beta }_1-p_l{\lambda }_l}{\alpha }_2^{q_2-1}{\beta }_2^{r_2-1}{\lambda }_2^{o_2-1}{e}^{-w_2{\alpha }_2-u_2{\beta }_2-p_2{\lambda }_2}.$$

(39)

To find out how to elicit the hyperparameters of the independent joint prior, utilize a variance–covariance matrix and estimate them using the likelihood technique. The resulting hyperparameters can be represented by the mean and variance of the gamma priors, see Mohamed et al. [28], Abulebda et al. [29] and Hassan and Chesneau [30].

The joint posterior distribution can be described as follows:

$$\begin{array}{cc}\hfill {\pi }^{*}(\Theta |x,y)& \propto {\alpha }_1^{q_1-1}{\beta }_1^{n+r_1-1}{\lambda }_1^{n+o_1-1}{\alpha }_2^{q_2-1}{\beta }_2^{n+r_2-1}{\lambda }_2^{n+o_2-1}{e}^{-w_1{\alpha }_1-u_1{\beta }_1-p_3{\lambda }_1-n}{e}^{-w_2{\alpha }_2-u_2{\beta }_2-p_2{\lambda }_2}\hfill \\ \hfill & log{\left[{\alpha }_1\right]}^n{({\alpha }_1-1)}^{-n}log{\left[{\alpha }_2\right]}^n{({\alpha }_2-1)}^{-n}\prod _{i=1}^n{\alpha }_2^{1-{(1+y_i^{{\beta }_2})}^{-{\lambda }_2}}{\alpha }_1^{1-{(1+x_i^{{\beta }_1})}^{-{\lambda }_1}}\hfill \\ \hfill & \prod _{i=1}^n{x}_i^{{\beta }_1-1}{y}_i^{{\beta }_2-1}{(1+y_i^{{\beta }_2})}^{-({\lambda }_2+1)}{(1+x_i^{{\beta }_1})}^{-({\lambda }_1+1)}\hfill \\ \hfill & \prod _{i=1}^n\left[1+\delta \left(1-2u\left(x_i\right)\right)\left(1-2v\left(y_i\right)\right)\right].\hfill \end{array}$$

(40)

The posterior conditional density functions of $({\alpha }_j,{\beta }_j,{\lambda }_j,\delta ),$where $j=1,2,$can be generated as follows:

$$\begin{array}{cc}\hfill {\pi }^{*}\left({\alpha }_1\right|{\beta }_1,{\lambda }_1,{\alpha }_2,{\beta }_2,{\lambda }_2,\delta )& ={\alpha }_1^{q_1-1}{e}^{-w_1{\alpha }_1}log{\left[{\alpha }_1\right]}^n{({\alpha }_1-1)}^{-n}\prod _{i=1}^n{\alpha }_1^{1-{(1+x_i^{{\beta }_1})}^{-{\lambda }_1}}\hfill \\ \hfill & \prod _{i=1}^n\left[1+\delta \left(1-2u\left(x_i\right)\right)\left(1-2v\left(y_i\right)\right)\right],\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill {\pi }^{*}\left({\beta }_1\right|{\alpha }_1,{\lambda }_1,{\alpha }_2,{\beta }_2,{\lambda }_2,\delta )& ={\beta }_1^{n+r_1-1}{e}^{-u_1{\beta }_1}\prod _{i=1}^n{\alpha }_1^{1-{(1+x_i^{{\beta }_1})}^{-{\lambda }_1}}{x}_i^{{\beta }_1-1}{(1+x_i^{{\beta }_1})}^{-({\lambda }_1+1)}\hfill \\ \hfill & \prod _{i=1}^n\left[1+\delta \left(1-2u\left(x_i\right)\right)\left(1-2v\left(y_i\right)\right)\right],\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill {\pi }^{*}\left({\lambda }_1\right|{\alpha }_1,{\beta }_1,{\alpha }_2,{\beta }_2,{\lambda }_2,\delta )& ={\lambda }_1^{n+o_1-1}{e}^{-p_1{\lambda }_1}\prod _{i=1}^n{\alpha }_1^{1-{(1+x_i^{{\beta }_1})}^{-{\lambda }_1}}{(1+x_i^{{\beta }_1})}^{-({\lambda }_1+1)}\hfill \\ \hfill & \prod _{i=1}^n\left[1+\delta \left(1-2u\left(x_i\right)\right)\left(1-2v\left(y_i\right)\right)\right],\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill {\pi }^{*}\left({\alpha }_2\right|{\alpha }_1,{\beta }_1,{\lambda }_1,{\beta }_2,{\lambda }_2,\delta )& ={\alpha }_2^{q_2-1}{e}^{-w_2{\alpha }_2}log{\left[{\alpha }_2\right]}^n{({\alpha }_2-1)}^{-n}\prod _{i=1}^n{\alpha }_2^{1-{(1+y_i^{{\beta }_2})}^{-{\lambda }_2}}\hfill \\ \hfill & \prod _{i=1}^n\left[1+\delta \left(1-2u\left(x_i\right)\right)\left(1-2v\left(y_i\right)\right)\right],\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill {\pi }^{*}\left({\beta }_2\right|{\alpha }_1,{\beta }_1,{\lambda }_1,{\alpha }_2,{\lambda }_2,\delta )& ={\beta }_2^{n+r_2-1}{e}^{-u_2{\beta }_2}\prod _{i=1}^n{\alpha }_2^{1-{(1+y_i^{{\beta }_2})}^{-{\lambda }_2}}{y}_i^{{\beta }_2-1}{(1+y_i^{{\beta }_2})}^{-({\lambda }_2+1)}\hfill \\ \hfill & \prod _{i=1}^n\left[1+\delta \left(1-2u\left(x_i\right)\right)\left(1-2v\left(y_i\right)\right)\right],\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill {\pi }^{*}\left({\lambda }_2\right|{\alpha }_1,{\beta }_1,{\lambda }_1,{\alpha }_2,{\beta }_2,\delta )& ={\lambda }_2^{n+o_2-1}{e}^{-p_2{\lambda }_2}\prod _{i=1}^n{\alpha }_2^{1-{(1+y_i^{{\beta }_2})}^{-\lambda 2}}{(1+y_i^{{\beta }_2})}^{-({\lambda }_2+1)}\hfill \\ \hfill & \prod _{i=1}^n\left[1+\delta \left(1-2u\left(x_i\right)\right)\left(1-2v\left(y_i\right)\right)\right]\hfill \end{array}$$

(46)

and

$${\pi }^{*}\left(\delta \right|{\alpha }_1,{\beta }_1,{\lambda }_1,{\alpha }_2,{\beta }_2,{\lambda }_2)=\prod _{i=1}^n\left[1+\delta \left(1-2u\left(x_i\right)\right)\left(1-2v\left(y_i\right)\right)\right].$$

(47)

Under the squared error loss function, which is a symmetrical loss function that assigns equal losses to overestimation and underestimation, the Bayesian estimator of $\Theta$ is the posterior mean $\widehat{\Theta }=E(\Theta |\mathrm{data})$. Given the analytical challenges in evaluating expectations for such loss functions, we utilized the Markov chain Monte Carlo (MCMC) method. This includes key tools like Gibbs sampling and Metropolis-within-Gibbs, which are well-suited for complex models. Among these, the Metropolis–Hastings (MH) algorithm is particularly useful, generating candidate values from a proposal distribution for each iteration of the sampling process. This approach was thoroughly discussed by [31,32]. The MH algorithm provides a transition mechanism for constructing the chain, as follows:

• Sample ${\Theta }^{*}$ from a proposal distribution $q\left({\Theta }^{*}\right|{\Theta }_i)$, with ${\Theta }_i$ as the current state.
• Decide to either accept ${\Theta }^{*}$ as the next sample or retain the current sample ${\Theta }_i$ based on the acceptance probability:

  $${\tau }_{{\Theta }^{*}|{\Theta }_i}=min\left(1,{\displaystyle \frac{q\left({\Theta }_i\right|{\Theta }^{*})p\left({\Theta }^{*}\right)}{q\left({\Theta }^{*}\right|{\Theta }_i)p\left({\Theta }_i\right)}}\right).$$

This acceptance probability ensures that the target density $p(\Theta |\mathrm{data})$ remains invariant, thereby allowing the chain to converge to its stationary distribution from any initial state. Thus, the MH approach guarantees that the chain approximates the desired posterior distribution of $\Theta$.

The approach of Chen et al. [33] is frequently used to construct highest posterior density ($HPD$) intervals for unknown distribution parameters in Bayesian estimation. For instance, two endpoints from the $MCMC$ sampling outputs, the lower $5\%$ and upper $95\%$ percentiles, can be used to calculate a $90\%$ HPD interval. The following is how reliable Bayesian intervals of the parameters are obtained:

• Arrange ${\alpha }_j^{\left[1\right]}<{\alpha }_j^{\left[2\right]}<\dots <{\alpha }_j^{\left[B\right]},{\beta }_j^{\left[1\right]}<{\beta }_j^{\left[2\right]}<\dots <{\beta }_j^{\left[B\right]},{\lambda }_j^{\left[1\right]}<{\lambda }_j^{\left[2\right]}<\dots <{\lambda }_j^{\left[B\right]}$and ${\delta }^{\left[1\right]}<{\delta }^{\left[2\right]}<\dots <{\delta }^{\left[B\right]}$ where $j=1,2$ and B is the length of $MCMC$generated.
• The $100(1-\gamma )\%$ HPD intervals become (${\alpha }_j^{\left[B{\displaystyle \frac{\gamma }{2}}\right]},{\alpha }_j^{\left[B(1-{\displaystyle \frac{\gamma }{2}})\right]}$), (${\beta }_j^{\left[B{\displaystyle \frac{\gamma }{2}}\right]},{\beta }_j^{\left[B(1-{\displaystyle \frac{\gamma }{2}})\right]}$), (${\lambda }_j^{\left[B{\displaystyle \frac{\gamma }{2}}\right]},{\lambda }_j^{\left[B(1-{\displaystyle \frac{\gamma }{2}})\right]}$) and (${\delta }^{\left[B{\displaystyle \frac{\gamma }{2}}\right]},{\delta }^{\left[B(1-{\displaystyle \frac{\gamma }{2}})\right]}$).

4.3. Least Squares Estimation

The least squares method ($LES$) was proposed for estimating the parameters ${\alpha }_1,{\beta }_1,{\lambda }_1,{\alpha }_2,$${\beta }_2,{\lambda }_2$ and$\delta$ by minimizing the function $W({\alpha }_1,{\beta }_1,{\lambda }_1,{\alpha }_2,{\beta }_2,{\lambda }_2,\delta )$ defined as

$$W({\alpha }_1,{\beta }_1,{\lambda }_1,{\alpha }_2,{\beta }_2,{\lambda }_2,\delta )=\sum _{i=1}^n{\left(F(x_i,y_i)-F_n(x_i,y_i)\right)}^2,$$

(48)

where $F(x_i,y_i)$ is the $F_n(x_i,y_i)$, i = $1,2,\dots ,n$ represents the empirical distribution function, evaluated at the observed value of $(x_i,y_i)$, given by

$$F_n(x_i,y_i)={\displaystyle \frac{{\sum }_{i=1}^n\left((X_i\le x_i),(Y_i\le y_i)\right)}{n}}.$$

(49)

We propose to use the percentile bootstrap (Boot-p) method to estimate the confidence intervals for the parameters based on least squares estimates ($LSEs$). To estimate the bootstrap interval, we proceed as follows:

(i)

Obtain the $LSEs$ for the parameters ${\Theta }_i;i=1,\dots ,7$.

(ii)

Generate bootstrap sample from BAPB-XII, using ${\widehat{\Theta }}_i;i=1,\dots ,7$ and then obtain the bootstrap estimate of ${\Theta }_i$ say ${\widehat{\Theta }}_i;i=1,\dots ,7$ respectively.

(iii)

Repeat step $\left(ii\right)$ a large number N of times.

(iv)

From ${\widehat{\Theta }}_{\left(1\right)}^{(*)}\le {\widehat{\Theta }}_{\left(2\right)}^{(*)}\le \cdots \le {\widehat{\Theta }}_{\left(N\right)}^{(*)},$ the bootstrap confidence interval $100(1-\gamma )\%$, $(0<\gamma <1)$ for the parameter ${\Theta }_i;i=1,\dots ,7$ is given by $\left({\widehat{\Theta }}_{\left(v_1\right)}^{(*)},{\widehat{\Theta }}_{\left(v_2\right)}^{(*)}\right)$, where $v_1=N\left({\displaystyle \frac{\alpha }{2}}\right),v_2=N-v_1$. This is performed for the vector parameters $\Theta$.

5. Simulation Study

This section presents simulations studies conducted to study the performance of the maximum likelihood, Bayesian, and least squares estimation approaches. For the parameters of the BAPB-XII distribution, the point estimates were assessed using mean square error (MSE) and average bias (AB). Other criteria such as the average width (AW) and coverage probabilities (CP) were used to asses the performance of $95\%$ ACI, HPD, and Boot-p intervals. The simulation was conducted for 1000 iterations, each with sample size $n=25, 50, 100$, and $150.$ The following algorithm was used to generate random samples from the BAPB-XII distribution and conduct the simulation study:

• For a given sample size n, generate independent random values Z and W from a uniform (0,1) distribution.
• For the parameters ${\alpha }_1,{\beta }_1$, and ${\lambda }_1,$ use the inverse transformation method to generate X from the APB-XII distribution utilizing the equation

  $$X={\left[\left(1-{\left({\displaystyle \frac{log(Z({\alpha }_1-1)+1)}{log\left({\alpha }_1\right)}}\right)}^{{\displaystyle \frac{-1}{{\lambda }_1}}}\right)-1\right]}^{{\displaystyle \frac{1}{{\beta }_1}}}.$$
• For ${\alpha }_2,{\beta }_2,{\lambda }_2$, and $\delta$, the random variable $Y$ is obtained numerically, setting $F\left(Y\right|X)=W$ in Equation (14).
• Obtain the $MLEs$ for the parameters, say ${\widehat{\alpha }}_j,{\widehat{\beta }}_j,{\widehat{\lambda }}_j$, and $\widehat{\delta },j=1,2$.
• Run the $MCMC$ algorithm 10,000 times, with normal proposals, discarding the first 1000 as a burn-in period and take every 3rd value in the thinning process to reduce the dependence between observations.
• Construct the $95\%$ ACI of MLEs and HPDs of Bayesian estimates.
• Obtain the least square estimates and their respective Boot-p intervals for the parameters.
• Repeat these steps 1000 times.

The above algorithm was considered for the following scenarios of parameters: (${\alpha }_1$, ${\beta }_1$, ${\lambda }_1$,${\alpha }_2$, ${\beta }_2$, ${\lambda }_2$) = (0.4, 1.5, 0.5, 0.8, 1.4, 1.5), (1.25, 0.25, 1.5, 0.25, 0.5, 0.5). A specific copula parameter value was assigned to each set of parameters, namely $\delta =0.6$ and $-0.6,$ with the results shown in

Table 1. Simulation results: ${\alpha }_1=0.4, {\beta }_1=1.5, {\lambda }_1=0.5, {\alpha }_2=0.8, {\beta }_2=1.4, {\lambda }_2=1.5$ and $\delta =0.6$.

Param.	n	MLEs				BSEs				LSEs
		MSEs	ABs	AWs	CPs	MSEs	ABs	AWs	CPs	MSEs	ABs	AWs	CPs
${\alpha }_1$	25	2.3433	0.5725	10.4989	0.93	0.1902	0.07	2.1533	1.00	4.1839	0.4059	5.7878	1.00
	50	2.0475	0.5609	8.9687	0.97	0.1062	−0.0353	1.9919	1.00	3.4491	0.2353	5.3441	1.00
	75	1.9159	0.4636	5.3525	0.97	0.157	0.0311	2.0932	1.00	2.2283	0.0358	5.2647	1.00
	100	1.8765	0.4541	3.9192	0.93	0.1941	0.0601	1.9435	0.98	1.9697	−0.0226	4.8453	1.00
	150	1.486	0.2975	2.5405	0.93	0.1037	−0.0183	1.8238	0.98	1.2659	−0.0673	4.5969	1.00
${\beta }_1$	25	0.0023	0.0113	0.4038	1.00	0.0026	0.0133	0.1828	0.94	0.0038	−0.011	0.2573	0.88
	50	0.0012	−0.0034	0.2571	1.00	0.0012	0.015	0.1325	0.96	0.0025	−0.0005	0.1741	0.96
	75	0.0012	−0.0061	0.1945	0.96	0.0007	0.0014	0.1068	0.90	0.0014	0.0023	0.1674	1.00
	100	0.0007	0.0012	0.1838	1.00	0.0009	0.0136	0.0982	0.92	0.0019	−0.0031	0.1487	0.96
	150	0.0006	0.0003	0.1542	0.99	0.0004	0.0036	0.0807	0.96	0.001	−0.0011	0.1358	1.00
${\lambda }_1$	25	0.1651	0.0315	5.1286	0.99	0.0786	0.0322	1.2186	0.98	0.4644	0.0046	2.6705	1.00
	50	0.1464	0.0877	3.5781	0.97	0.0571	−0.0361	1.0259	0.96	0.3865	−0.1445	1.9991	0.96
	75	0.1297	0.0766	2.6555	0.99	0.0382	0.0081	0.9374	0.96	0.3730	−0.3126	2.0522	1.00
	100	0.118	0.0064	2.4095	0.98	0.0424	−0.0606	0.8096	0.94	0.365	−0.2423	1.9311	1.00
	150	0.0962	0.0457	2.0727	0.99	0.0212	−0.0278	0.7965	0.98	0.2993	−0.2908	1.8164	1.00
${\alpha }_2$	25	2.994	1.2804	4.2643	0.98	2.3759	1.3422	2.4483	0.98	1.7454	−0.05	2.1791	1.00
	50	1.3296	1.4686	4.5498	0.96	1.7978	1.0715	2.6455	0.96	1.2565	0.4719	3.2442	1.00
	75	0.7596	1.6672	3.9925	0.95	1.2417	0.8061	1.9475	0.87	0.8388	0.2604	3.8811	1.00
	100	0.8317	1.9868	3.6489	0.99	0.4372	0.3848	1.7647	0.86	0.5425	0.3787	3.6142	1.00
	150	0.3195	2.2979	2.1629	0.93	0.1952	0.1881	1.4078	0.89	0.0396	0.4049	4.2645	1.00
${\beta }_2$	25	0.0074	−0.0124	0.5453	1.00	0.0154	0.0522	0.3884	0.96	0.0115	−0.0431	0.4636	0.96
	50	0.0055	−0.0095	0.4098	0.99	0.0027	0.0098	0.2533	1.00	0.0121	−0.019	0.3887	0.96
	75	0.0052	−0.0345	0.3388	0.99	0.0026	0.0024	0.2237	0.98	0.0063	−0.0294	0.2973	0.92
	100	0.0036	−0.0353	0.3229	0.99	0.0018	−0.0108	0.1979	0.96	0.0044	−0.0283	0.2552	0.88
	150	0.0053	−0.0546	0.2677	0.96	0.0020	−0.0250	0.1587	0.90	0.0058	−0.0192	0.2446	0.84
${\lambda }_2$	25	0.8254	0.658	3.0996	1.00	0.4896	0.6645	1.1496	0.96	0.29	0.0465	1.6117	1.00
	50	0.7217	0.7008	2.1727	0.94	0.3387	0.5237	1.0821	0.92	0.1913	0.2309	1.5838	1.00
	75	0.6975	0.7695	2.0123	0.91	0.2922	0.4565	0.9699	0.88	0.1752	0.1668	1.5679	0.88
	100	0.5907	0.8050	1.8827	0.92	0.1384	0.2822	0.8851	0.86	0.1616	0.2554	1.5066	0.92
	150	0.5597	0.8784	1.5509	0.89	0.0827	0.1553	0.7288	0.86	0.0215	0.214	1.4724	0.89
$\delta$	25	0.159	−0.2242	2.2894	1.00	1.4237	−0.1191	2.3985	0.86	0.6643	−0.2385	2.0012	1.00
	50	0.1454	−0.2692	1.5731	1.00	0.256	−0.2914	1.4189	0.90	0.5225	−0.257	1.8458	0.92
	75	0.172	−0.3054	1.2767	0.95	0.2329	−0.3269	1.2163	0.90	0.3954	−0.1194	1.8312	1.00
	100	0.1256	−0.2518	1.1229	0.89	0.1676	−0.3057	1.0648	0.92	0.1676	−0.477	1.7268	0.96
	150	0.1289	−0.284	0.9076	0.88	0.1557	−0.3164	0.8866	0.86	0.1264	−0.241	1.5231	0.90

,

Table 2. Simulation results: ${\alpha }_1=0.4, {\beta }_1=1.5, {\lambda }_1=0.5, {\alpha }_2=0.8, {\beta }_2=1.4, {\lambda }_2=1.5$ and $\delta =-0.6$.

Param.	n	MLEs				BSEs				LSEs
		MSEs	ABs	AWs	CPs	MSEs	ABs	AWs	CPs	MSEs	ABs	AWs	CPs
${\alpha }_1$	25	0.4782	0.0725	9.0164	0.97	0.1615	0.1194	2.036	1.00	3.2026	0.3257	1.9536	1.00
	50	0.4376	0.0287	7.8716	0.96	0.1569	−0.0002	1.9493	1.00	2.324	0.1755	1.948	0.93
	75	0.4205	0.098	7.2338	0.97	0.1363	−0.0206	1.9346	0.98	2.0636	−0.0344	1.6644	1.00
	100	0.4044	−0.0101	6.476	0.96	0.1157	−0.0031	1.9933	0.98	1.3973	0.0052	1.3039	1.00
	150	0.4019	0.0408	6.5828	0.96	0.0994	−0.0658	1.9727	1.00	0.4595	0.0570	1.2130	1.00
${\beta }_1$	25	0.0015	0.0109	0.2806	0.99	0.003	0.0142	0.1818	0.92	0.0035	−0.0142	0.2049	0.9
	50	0.0010	0.0056	0.2300	0.96	0.0010	0.0032	0.1236	0.92	0.0017	−0.0134	0.1450	0.93
	75	0.0006	0.0031	0.1924	0.98	0.0007	0.0015	0.1046	0.92	0.0015	0.001	0.1061	0.93
	100	0.0006	0.0049	0.1812	0.97	0.0006	0.0001	0.0928	0.96	0.0023	0.0064	0.1096	0.92
	150	0.0004	0.0014	0.1689	0.96	0.0004	0.0003	0.0791	0.96	0.001	0.0047	0.0963	0.93
${\lambda }_1$	25	0.1098	−0.0305	3.2322	1.00	0.0748	−0.049	1.2308	0.98	0.4812	0.1919	1.5217	0.99
	50	0.0838	−0.0272	2.8700	0.99	0.0545	−0.0377	0.9772	0.96	0.4143	−0.1318	1.1668	0.87
	75	0.0591	−0.0056	2.4637	0.99	0.0321	−0.0210	0.9439	0.96	0.3323	−0.1239	1.0057	0.93
	100	0.0675	−0.0247	2.3589	0.98	0.0315	−0.0148	0.9014	1.00	0.4267	−0.2448	0.977	0.87
	150	0.0637	−0.0047	2.2985	0.98	0.0191	0.015	0.8055	0.98	0.2468	−0.2131	0.9581	0.87
${\alpha }_2$	25	1.4015	0.5064	3.617	1.00	1.3098	0.5799	3.0891	0.94	1.5816	0.5873	1.8105	1.00
	50	1.2999	0.5803	2.7236	1.00	1.0206	0.4541	3.2777	0.96	1.5351	0.4857	1.1999	1.00
	75	0.8496	0.627	2.5848	0.95	0.5114	0.5658	2.9715	0.92	1.4983	0.4336	1.1599	1.00
	100	0.6288	0.7831	2.7722	0.93	0.4936	0.2912	2.6104	0.94	0.6728	0.4096	1.1234	1.00
	150	0.5800	0.7836	2.1168	0.96	0.1488	0.2676	3.2153	0.96	0.0979	0.0392	0.897	1.00
${\beta }_2$	25	0.0086	0.054	0.5014	1.00	0.0308	0.0986	0.4851	0.90	0.013	−0.0252	0.3703	0.93
	50	0.0057	0.0542	0.3538	0.99	0.0094	0.0584	0.3158	0.94	0.0079	−0.0379	0.3700	0.93
	75	0.0051	0.0523	0.2892	0.99	0.0065	0.0517	0.2639	0.90	0.0051	−0.0011	0.2624	0.93
	100	0.0048	0.0528	0.2565	0.97	0.0061	0.0547	0.2328	0.86	0.005	0.0036	0.2211	0.93
	150	0.0045	0.0531	0.2042	0.95	0.0021	0.0269	0.194	0.96	0.0033	0.0085	0.2079	0.94
${\lambda }_2$	25	0.1457	0.3238	1.5465	1.00	0.2434	0.4865	0.9309	0.98	0.1006	−0.0044	0.8923	1.00
	50	0.1274	0.3159	1.0044	0.94	0.1779	0.4064	0.8117	0.97	0.0968	0.0812	0.6192	1.00
	75	0.1232	0.2901	0.8818	0.92	0.1555	0.3738	0.6999	0.94	0.0939	−0.0473	0.4731	0.93
	100	0.1142	0.2822	0.7688	0.87	0.1085	0.2759	0.5809	0.92	0.051	0.0230	0.6142	0.93
	150	0.1131	0.2540	0.5946	0.90	0.0462	0.0673	0.5292	0.91	0.0222	0.0189	0.4461	0.93
$\delta$	25	0.2462	0.4232	2.6156	1.00	0.6002	0.3712	2.6413	0.92	0.9064	0.485	1.4711	0.90
	50	0.2410	0.3830	1.6752	1.00	0.3168	0.3542	1.5215	0.88	0.2937	−0.0785	0.8445	0.87
	75	0.2377	0.3609	1.3525	0.98	0.2558	0.3405	1.2028	0.84	0.2778	0.2921	0.4531	0.90
	100	0.2350	0.3184	1.1617	0.87	0.1815	0.3088	1.0855	0.86	0.1373	0.2664	0.6865	0.86
	150	0.2274	0.2195	0.9405	0.88	0.181	0.0403	0.9256	0.98	0.0572	0.1874	0.2514	0.83

,

Table 3. Simulation results: ${\alpha }_1=1.25, {\beta }_1=0.25, {\lambda }_1=1.5, {\alpha }_2=0.25, {\beta }_2=0.5, {\lambda }_2=0.5$ and $\delta =0.6$.

Param.	n	MLEs				BSEs				LSEs
		MSEs	ABs	AWs	CPs	MSEs	ABs	AWs	CPs	MSEs	ABs	AWs	CPs
${\alpha }_1$	25	0.9371	0.4388	3.0563	0.99	0.2427	0.13	1.5769	0.96	1.1263	0.6765	4.2227	1.00
	50	0.9443	0.3934	2.4268	0.98	0.1887	−0.003	1.4498	0.96	1.0586	0.4259	4.0605	1.00
	75	0.4819	0.3377	2.1758	0.98	0.1482	0.0477	1.2657	0.96	0.7659	0.325	4.0567	1.00
	100	0.6099	0.3023	1.5881	0.97	0.1056	0.1265	1.2320	0.98	0.7516	0.2606	3.6988	1.00
	150	0.2627	0.1876	0.9971	0.98	0.1051	0.1167	1.2025	0.96	0.1587	0.1067	2.5077	1.00
${\beta }_1$	25	0.0911	0.0141	1.6759	0.99	0.0889	0.0542	1.1324	0.96	0.1142	−0.0205	1.6611	1.00
	50	0.0691	−0.0306	1.1500	0.94	0.0277	0.0461	0.8289	0.94	0.0825	−0.0946	1.1247	0.96
	75	0.0393	−0.0278	0.8777	0.94	0.0428	0.0021	0.7296	0.94	0.0431	−0.1445	1.0117	1.00
	100	0.0367	−0.0423	0.7587	0.95	0.0274	0.0114	0.6401	0.94	0.0788	−0.1408	0.9014	1.00
	150	0.0198	−0.0307	0.6074	0.97	0.0169	−0.0186	0.5336	0.94	0.0369	−0.0281	0.8266	0.96
${\lambda }_1$	25	0.0659	0.1106	1.6525	1.00	0.0297	−0.0465	0.712	0.94	0.1116	0.1394	1.0527	1.00
	50	0.0523	0.1028	1.2329	1.00	0.0262	−0.0587	0.6561	0.98	0.0843	0.0998	1.1502	1.00
	75	0.0319	0.0704	0.9023	0.98	0.0212	−0.0498	0.5613	0.94	0.0719	0.0901	1.0805	1.00
	100	0.0333	0.0796	0.8457	0.99	0.0161	−0.0438	0.5465	0.94	0.0664	0.0450	1.1277	1.00
	150	0.0219	0.0462	0.6598	0.99	0.0139	−0.0176	0.4776	0.96	0.0594	0.0344	1.0288	1.00
${\alpha }_2$	25	1.2701	0.4109	5.795	0.98	0.1148	−0.1173	1.6113	1.00	2.7087	0.9012	4.3806	1.00
	50	1.3546	0.2996	3.3668	0.93	0.1085	−0.117	1.6075	0.96	1.0732	0.7377	5.0823	1.00
	75	1.2346	0.3079	3.2333	0.93	0.0902	−0.0823	1.5380	0.98	1.0221	0.4272	5.3337	1.00
	100	1.2031	0.2988	2.7805	0.96	0.0722	−0.1082	1.5307	1.00	0.9165	0.3647	5.2027	1.00
	150	0.6466	0.0854	2.9428	0.95	0.0518	−0.1044	1.4713	0.93	0.2587	−0.0317	5.2039	1.00
${\beta }_2$	25	0.0522	−0.0339	1.6577	0.98	0.0468	0.0489	0.9069	0.96	0.1482	−0.2525	1.318	0.94
	50	0.0415	−0.0648	1.1445	0.94	0.0152	−0.0041	0.6434	1.00	0.0626	−0.1347	1.0850	0.92
	75	0.0271	−0.0403	1.0543	0.98	0.0157	−0.0052	0.5532	0.98	0.0552	−0.1031	0.9994	0.94
	100	0.0265	−0.0568	0.9498	0.98	0.0115	0.0037	0.4852	0.98	0.0464	−0.0625	0.8792	1.00
	150	0.0151	−0.0362	0.7625	0.98	0.0108	−0.0231	0.4151	0.96	0.0354	−0.0408	0.7367	1.00
${\lambda }_2$	25	0.4228	0.3374	5.2249	1.00	0.1226	0.1363	1.5515	0.96	0.2805	0.2577	3.0066	0.96
	50	0.3183	0.2711	3.9152	0.98	0.1073	0.1525	1.3674	0.96	0.5123	0.4014	2.9126	1.00
	75	0.2887	0.2853	3.5089	0.98	0.0549	0.1033	1.2778	0.98	0.4595	0.1877	2.8987	1.00
	100	0.2323	0.2751	3.2908	0.99	0.0557	0.1113	1.1700	1.00	0.6896	0.3503	2.9705	1.00
	150	0.1939	0.2414	2.6652	0.98	0.0550	0.1342	1.1504	0.98	0.5485	0.1422	2.8217	1.00
$\delta$	25	0.1458	−0.2005	2.3234	1.00	1.7577	−0.1091	2.2023	0.78	1.0274	−0.5315	2.1481	1.00
	50	0.1318	−0.2163	1.6292	1.00	1.0322	−0.2688	1.8607	0.84	0.3821	−0.2501	1.9757	1.00
	75	0.1322	−0.2463	1.3187	1.00	0.2026	−0.3382	1.2383	0.88	0.2351	−0.276	1.9738	1.00
	100	0.1236	−0.2546	1.1420	0.91	0.2219	−0.3067	1.1263	0.83	0.2479	−0.3629	1.8619	0.92
	150	0.1262	−0.2732	0.9338	0.87	0.1953	−0.2821	0.9186	0.87	0.2191	−0.3012	1.7415	0.92

and

Table 4. Simulation results: ${\alpha }_1=1.25, {\beta }_1=0.25, {\lambda }_1=1.5, {\alpha }_2=0.25, {\beta }_2=0.5, {\lambda }_2=0.5$ and $\delta =-0.6$.

Param.	n	MLEs				BSEs				LSEs
		MSEs	ABs	AWs	CPs	MSEs	ABs	AWs	CPs	MSEs	ABs	AWs	CPs
${\alpha }_1$	25	0.7585	0.3143	6.1793	0.97	0.1709	0.1428	1.3477	0.94	0.9297	0.0549	4.0885	1.00
	50	0.1945	0.2590	3.5390	0.96	0.1375	0.1151	1.5254	0.96	0.7266	0.137	4.0047	1.00
	75	0.4224	0.2399	3.1069	0.98	0.1156	0.0978	1.3762	0.92	0.6202	0.3774	3.7453	1.00
	100	0.5427	0.1622	2.9150	0.97	0.0830	0.0311	1.1283	0.90	0.3202	0.0696	3.5493	1.00
	150	0.1612	0.0792	1.7657	0.95	0.0642	0.0104	1.2615	0.98	0.2115	0.3353	2.7123	1.00
${\beta }_1$	25	0.0467	−0.004	2.0008	0.99	0.0908	0.0953	1.1682	0.96	0.1162	−0.1965	1.2765	0.92
	50	0.0277	−0.012	1.2923	0.98	0.0422	0.0502	0.8809	0.94	0.0918	−0.0496	1.2578	0.96
	75	0.0241	−0.0096	0.9112	0.99	0.0283	0.0189	0.7023	0.98	0.0503	−0.1005	0.9411	0.92
	100	0.0243	−0.0270	0.7885	0.98	0.0225	0.0255	0.6265	0.94	0.0455	−0.0612	0.8899	0.92
	150	0.0201	0.0014	0.6165	0.99	0.0248	0.0257	0.5486	0.92	0.0479	−0.0787	0.7865	0.93
${\lambda }_1$	25	0.0551	0.0778	1.9446	1.00	0.0278	−0.0916	0.6901	0.94	0.0763	0.0568	1.1755	1.00
	50	0.0303	0.0535	1.3408	1.00	0.0202	−0.0338	0.6449	0.96	0.0641	−0.0346	1.0362	1.00
	75	0.0275	0.0481	0.9870	0.99	0.0184	−0.0316	0.5961	0.96	0.0582	0.0927	1.0184	1.00
	100	0.0297	0.0516	0.8636	1.00	0.0168	−0.0532	0.5142	0.94	0.0507	0.0624	1.1168	1.00
	150	0.0162	−0.0031	0.6800	0.99	0.0159	−0.0143	0.5189	0.96	0.0394	0.0019	1.1034	1.00
${\alpha }_2$	25	2.2027	0.7450	6.2183	0.99	0.2786	0.2692	1.8723	0.98	3.1676	0.1508	4.6536	1.00
	50	1.7195	0.5471	5.1358	0.98	0.1862	0.2006	1.8898	0.98	3.1672	−0.0865	4.5475	1.00
	75	1.6457	0.5316	5.0529	0.99	0.1811	0.1990	1.9809	0.98	0.4790	−0.1751	4.4509	1.00
	100	0.9434	0.2831	2.5911	0.97	0.1353	0.1296	1.8849	0.98	0.4517	0.6070	4.4115	1.00
	150	0.6238	0.1984	2.5310	0.97	0.1209	0.1111	1.7161	0.94	0.4177	0.8623	3.944	1.00
${\beta }_2$	25	0.0433	0.0952	1.5546	1.00	0.0980	0.1591	0.9479	0.94	0.1061	0.1292	1.4069	1.00
	50	0.0285	0.0573	1.0280	0.99	0.0434	0.1299	0.7449	0.96	0.1005	0.1219	1.1025	0.96
	75	0.0305	0.0517	1.0424	0.99	0.0393	0.1293	0.6492	0.96	0.0597	0.0543	1.0715	0.96
	100	0.0295	0.0513	0.9770	1.00	0.0343	0.1236	0.5532	0.94	0.0442	0.0257	0.9539	1.00
	150	0.0292	0.0447	0.7588	0.94	0.0329	0.0863	0.4814	0.92	0.0290	0.0235	0.8721	1.00
${\lambda }_2$	25	0.1127	−0.0942	3.0058	1.00	0.1096	−0.1666	1.2351	0.90	0.3040	−0.1229	2.0991	1.00
	50	0.1115	−0.1006	2.2132	1.00	0.102	−0.2089	1.0896	0.86	0.2915	−0.2876	1.9295	1.00
	75	0.1106	−0.0671	2.2781	0.97	0.0906	−0.2214	0.9823	0.88	0.2734	−0.4025	1.8676	1.00
	100	0.1046	−0.0447	2.1764	0.96	0.0721	−0.1733	0.9192	0.94	0.2488	−0.2035	1.6521	1.00
	150	0.0789	−0.1122	1.7172	0.96	0.0585	−0.1998	0.8243	0.88	0.2126	−0.1862	1.6512	1.00
$\delta$	25	0.1405	0.2619	2.5285	1.00	1.9997	0.3410	2.3576	0.93	0.8010	0.3737	1.9747	1.00
	50	0.1289	0.2428	1.6350	1.00	1.2348	0.3334	1.6736	0.84	0.2643	0.3565	1.8570	1.00
	75	0.1221	0.2282	1.3451	1.00	0.7017	0.3122	1.4008	0.84	0.3583	0.2090	1.7154	1.00
	100	0.1102	0.2185	1.1471	0.99	0.2266	0.1521	1.0650	0.86	0.2080	0.1557	1.6364	0.96
	150	0.1079	0.1663	0.9359	0.86	0.1759	−0.0611	0.9023	0.82	0.1997	0.0053	1.4516	0.84

.

• All estimation methods, Bayesian and non-Bayesian approaches, performed well and could be utilized to effectively estimate the parameters of the BAPB-XII model.
• Across parameters, the MSEs generally decreased as the sample size n increased, indicating that all estimation methods became more accurate with larger samples. The BSEs showed lower MSEs than the MLEs and LSEs for most parameters, suggesting that Bayesian estimation might yield more precise parameter estimates.
• The biases for most estimators were small, especially as n increased. The BSEs often had slightly lower biases.
• The AWs tended to decrease with sample size for all estimates, but the BSEs had narrower confidence intervals. However, for the copula parameters $\delta$, the MLEs and LSEs often had larger AWs.
• The CPs were close to the nominal level for all estimation approaches.
• The BSEs generally demonstrated a balanced performance with low MSEs, reduced bias, and acceptable coverage probabilities with increasing sample sizes, indicating they were a reliable choice for parameter estimation in this setup. These findings suggest that BSEs may be preferable for accuracy and precision, whereas MLEs provide a more conservative choice with reliable coverage probabilities.
• Negative initial values of the copula parameter $\delta$ tended to introduce greater bias and variability, especially in the BSEs, where the ABs for $\delta$ tended to be larger when the initial value was negative. Positive initial values led to better stability, lower biases, and higher coverage probabilities for $\delta$ across all estimators.

6. Application of Real Data

In this section, we demonstrate the practical usefulness of the proposed bivariate model through its application to a real-world dataset. The data consist of 50 bivariate observations on burr measurements, where X represents burr measurements (in millimeters) for components with a hole diameter of 12 mm and a sheet thickness of 3.15 mm, and Y represents burr measurements for components with a hole diameter of 9 mm and a sheet thickness of 2 mm. The data were obtained from Dasgupta [34], and involved readings taken from jobs with a single, fixed hole, selected based on a predetermined orientation. These two datasets correspond to burr measurements from the two different machines being compared. Additionally, the model was applied to analyze the relationship between X and Y, accounting for all possible relationships, such as $X<Y$, $X>Y$, and $X=Y$.

As detailed in Section 3.1, the marginal distributions of X and Y were assumed to follow the BAPB-XII distribution. To validate this assumption, we analyzed the goodness of fit for each marginal distribution.

Table 5. MLEs, $-\ell$, K-S, and p-value statistic for X, Y.

Variable	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\lambda }$	$-\mathit{\ell }$	K-S	p-Value
X	$2.9951$	$1.9205$	$35.189$	$56.047$	$0.1066$	$0.6210$
Y	$3.3104$	$1.8187$	$34.389$	$57.635$	$0.1488$	$0.2184$

summarizes the results, providing MLEs for the parameters along with the Kolmogorov–Smirnov (K-S) statistics and corresponding p-value. These values indicate that the marginal distributions fit well for both X and Y.

This conclusion is further illustrated in Figure 2 and Figure 3, where we present the estimated pdfs, survival curves, and probability–probability (pp) plots for the marginal distributions of X and Y. These visualizations demonstrate a close alignment between the BAPB-XII distribution and the observed data, confirming the adequacy of this model for the given burr measurement data.

We now proceed to model the bivariate relationship between X and Y. This section evaluates the performance of several bivariate models to identify the best fit for the data. Specifically, we compared the BAPB-XII distribution, the bivariate Fréchet (BF) distribution [35], the bivariate power Lomax (BL) distribution [5], the bivariate Lomax exponential (BLE) distribution [36], and the bivariate generalized half-logistic (BGHL) distribution [30].

Table 6. MLEs and SEs (in parenthesis) for parameters of bivariate distributions.

Model	${\mathit{\alpha }}_1$	${\mathit{\beta }}_1$	${\mathit{\lambda }}_1$	${\mathit{\alpha }}_2$	${\mathit{\beta }}_2$	${\mathit{\lambda }}_2$	$\mathit{\delta }$	AIC	BIC
BAPB-XII	2.999	1.917	35.259	3.267	1.815	34.434	0.545	−215.283	−201.299
	(4.0900)	(0.3761)	(13.8221)	(3.7410)	(0.3111)	(12.833)	(0.368)
BL	1.511	1.342	0.1409	2.687	0.865	0.7941	1.00	−123.92	−110.5
	(0.6700)	(0.1930)	(0.1655)	(1.3282)	(0.2498)	(0.4080)	(8.1052)
BF	0.092	1.248		0.082	1.197		0.894	−144.684	−135.124
	(0.0112)	(0.1180)		(0.0104)	(0.1160)		(0.6274)
BGHL	0.097	0.917		0.0869	0.8808		1.00	−186.297	−176.736
	(0.0501)	(0.6034)		(0.0439)	(0.5704)		(0.6177)
BLE	75.698	6.7600	0.8156	72.559	11.293	1.447	0.588	−148.248	−134.864
	(12.4850)	(2.4371)	(0.1264)	(12.681)	(4.6140)	(0.2305)	(0.7322)

presents the MLEs of the parameters for each model, along with their respective standard errors (SEs), and compares the goodness-of-fit using the Akaike information criterion (AIC) and Bayesian information criterion (BIC). The results indicate that the BAPB-XII distribution provided the best fit to the data, as evidenced by its lower AIC and BIC values.

The findings of our comparative analysis highlight the superior efficacy of the BAPB-XII model in accurately representing the bivariate relationship between X and Y compared to the other well-established models, including the bivariate Fréchet, power Lomax, Lomax exponential, and bivariate generalized half-logistic distributions. This result highlights the potential of BAPB-XII as a significant contribution to the current research. To this end, the BAPB-XII model establishes itself as a valuable tool for bivariate data modeling, paving the way for broader applications in future research.

Finally, we employed the estimation methods discussed earlier, utilizing both Bayesian and non-Bayesian approaches, to obtain estimates of the parameter vector $\Theta$. The parameter estimates derived from the MCMC algorithm, performed over 55,000 iterations, provided significant insights into the Bayesian framework for analyzing burr measurements. After discarding the initial 5000 runs as burn-in, we selected every fifth observation during the thinning process and utilized 10,000 bootstrap samples to establish $95\%$ bootstrapped confidence intervals.

Table 7. MLEs, BSEs, LSEs, and corresponding interval estimates for the datasets.

Parameter	MLEs	ACI	BSEs	HPD Intervals	LSEs	Boot-p Intervals
${\alpha }_1$	2.9998	(−5.0167, 11.016)	3.0128	(1.7592, 4.5466)	3.6809	(0.3376, 6.7696)
${\beta }_1$	1.9176	(1.1804, 2.6548)	1.5386	(1.3057, 1.7771)	1.6487	(1.2035, 1.9207)
${\lambda }_1$	35.260	(8.1684, 62.351)	17.915	(12.947, 24.125)	25.929	(11.181, 30.236)
${\alpha }_2$	3.2670	(−4.0660, 10.600)	3.0366	(1.8097, 4.5239)	3.6193	(0.30594, 6.7584)
${\beta }_2$	1.8156	(1.2059, 2.4254)	1.4797	(1.2530, 1.7040)	1.4272	(0.9883, 1.6044)
${\lambda }_2$	34.435	(9.2803, 59.589)	18.025	(12.730, 24.024)	20.539	(9.2123, 24.274)
$\delta$	0.5456	(−0.1773, 1.2684)	0.5362	(−0.4031, 1.2930)	0.3966	(−1.0, 1.0)

summarizes the parameter estimates acquired through MLEs, BSEs, and LSEs, each accompanied by the interval estimations.

We employed different diagnostic tools to assess the convergence of the MCMC chain. Firstly, the graph of the posterior density distributions for the estimated parameters, displayed in Figure 4, appears approximately normal, indicating well-behaved parameter estimates. Secondly, the trace plots displayed in Figure 5 reinforce our findings by illustrating the convergence of parameter values around specific pivot points, demonstrating the stability of the MCMC chains, which is essential for validating the effectiveness of the sampling process. Lastly, Figure 6 depicts the autocorrelation of the first 100 lags obtained from the MCMC algorithm, showing that the autocorrelation values started at positive one and gradually declined towards zero. This behavior indicates a rapid decrease in correlation with greater lags, supporting the effectiveness of the thinning process.

7. Conclusions

In this study, we introduced a novel bivariate distribution, the BAPB-XII distribution, derived from the alpha power Burr-XII distribution and the Farlie–Gumbel–Morgenstern copula function. This model exhibits well-defined statistical properties, making it particularly suitable for analyzing complex life testing data. We thoroughly examined the marginal distributions of the BAPB-XII model, along with its product moments and moment-generating functions, and developed closed-form expressions for these properties, thereby enhancing the model’s theoretical foundation. This rigorous mathematical framework facilitated seamless integration into real-world data analysis, and the efficacy of these methods was assessed through simulation studies, providing robust tools for parameter estimation, which can be applied in various practical scenarios.

A range of estimation methods were employed to evaluate the parameters of the (BAPB-XII) model, including non-Bayesian techniques such as maximum likelihood and least squares, alongside Bayesian estimates obtained using the MCMC approach. We conducted a real-data analysis focused on sheet failure rates to demonstrate the model’s practicality and effectiveness. The analysis revealed that the BAPB-XII model is well-suited for assessing failure rates in situations where the bivariate dependencies between random variables are significant. This case study exemplified the model’s potential utility in understanding and modeling complex, interdependent failure processes.

While the BAPB-XII model shows great promise, there are challenges in bivariate data modeling that warrant future research. Addressing issues related to estimating model parameters under various censoring schemes is essential. Moreover, we aim to explore the integration of this model in accelerating life testing analyses. Investigating the model’s performance across diverse fields—such as reliability engineering, finance, and risk assessment—represents a valuable direction for future studies, ensuring that the BAPB-XII model evolves into a versatile tool in statistical modeling.