Discrete Joint Random Variables in Fréchet-Weibull Distribution: A Comprehensive Mathematical Framework with Simulations, Goodness-of-Fit Analysis, and Informed Decision-Making

Abstract

This paper introduces a novel four-parameter discrete bivariate distribution, termed the bivariate discretized Fréchet–Weibull distribution (BDFWD), with marginals derived from the discretized Fréchet–Weibull distribution. Several statistical and reliability properties are thoroughly examined, including the joint cumulative distribution function, joint probability mass function, joint survival function, bivariate hazard rate function, and bivariate reversed hazard rate function, all presented in straightforward forms. Additionally, properties such as moments and their related concepts, the stress–strength model, total positivity of order 2, positive quadrant dependence, and the median are examined. The BDFWD is capable of modeling asymmetric dispersion data across various forms of hazard rate shapes and kurtosis. Following the introduction of the mathematical and statistical frameworks of the BDFWD, the maximum likelihood estimation approach is employed to estimate the model parameters. A simulation study is also conducted to investigate the behavior of the generated estimators. To demonstrate the capability and flexibility of the BDFWD, three distinct datasets are analyzed from various fields, including football score records, recurrence times to infection for kidney dialysis patients, and student marks from two internal examination statistical papers. The study confirms that the BDFWD outperforms competitive distributions in terms of efficiency across various discrete data applications.

1. Introduction

Bivariate probability distributions play a pivotal role in statistical modeling, allowing for the exploration of dependencies between two random variables across a wide spectrum of fields, including engineering, biology, finance, and social sciences. In bivariate analysis, the goal is to capture the interdependence structure between paired data points, whether they represent continuous or discrete variables. This enables better insights into complex systems where outcomes are not mutually exclusive and are often influenced by shared underlying factors.

For continuous data, bivariate distributions such as the bivariate exponential and the bivariate Weibull are commonly employed. These models offer considerable flexibility in capturing relationships between continuous variables and are particularly useful for modeling duration, reliability, and life-span data. Key methods for continuous bivariate modeling include copulas, which enable the decomposition of joint distributions into marginal distributions and a dependence structure. This decomposition is advantageous as it allows the modeling of each variable individually before addressing their joint behavior, providing both interpretive clarity and flexibility in analysis. So far, several continuous bivariate distributions have been developed by several researchers, including Jose et al. [1], Kundu and Gupta [2], Sarhan et al. [3], El-Sherpieny [4], Wagner and Artur [5], Balakrishnan and Shiji [6], El-Bassiouny et al. [7], El-Gohary et al. [8]. One well-known and commonly used bivariate distribution is the bivariate inverse Weibull distribution, developed by Hiba [9], which has been applied to football score data. Later, Rasool and Akbar [10] introduced the bivariate exponentiated extended-Weibull distribution. The bivariate Gumbel-G family of distributions and bivariate odd Weibull-G family of distributions were discussed in the works of Eliwa and El-Morshedy [11,12], respectively.

In contrast, discrete bivariate distributions are suitable for scenarios where data come in countable units, such as the number of occurrences, scores, or classifications, and are frequently encountered in studies related to count data analysis, zero-inflation, and over-dispersion. Techniques for constructing discrete bivariate distributions vary but often involve methods like the trivariate reduction technique, compounding approaches, and copula-based extensions, which allow for greater modeling flexibility and the ability to capture nuanced dependency structures in discrete data. Some researchers have developed various discrete bivariate distributions that are documented in the literature. Encyclopedic surveys of different discrete bivariate distributions can be found in the contributions of Kocherlakota and Kocherlakota [13]. Additionally, see Basu and Dhar [14], Kumar [15], Kemp [16], Lee and Cha [17], Nekoukhou and Kundu [18], and the references cited therein. For discrete bivariate data, the bivariate discrete inverse Weibull distribution is one of the most commonly used models, with marginals that are the discrete inverse Weibull, as proposed by Eliwa and El-Morshedy [19]. Later, El-Morshedy et al. [20] developed a bivariate discrete exponentiated Weibull distribution, and Ali et al. [21] discussed the bivariate discrete Nadarajah and Haghighi distribution in relation to football data and nasal drainage datasets. In recent years, there have been several advancements in bivariate discrete Weibull distributions. For instance, Kundu and Nekoukhou [22] developed the bivariate discrete Weibull distribution. Furthermore, Shibu and Beegum [23] discussed the bivariate discrete modified Weibull distribution, which has marginals that are the discrete modified Weibull.

The selection of an appropriate methodology for both continuous and discrete bivariate data is influenced by factors such as the nature of dependencies, marginal distribution types, and empirical patterns in the data. Nonetheless, the work of Johnson et al. [24] and Balakrishnan and Lai [25] stands out as highly influential and well-regarded in this area of research. Johnson et al. [24] conducted an extensive and detailed survey of various discrete bivariate distributions, while Balakrishnan and Lai [25] offered an in-depth analysis of numerous continuous bivariate distributions. The proposed paper introduces a versatile discrete probability model, developed as an extension of the Weibull model, termed the bivariate discretized Fréchet–Weibull distribution (BDFWD). To establish a foundation for the proposed model, it is essential to first outline the limitations of classical bivariate distributions and the advantages provided by the new approach.

The Novelty and Motivation Behind the BDFWD

In numerous real-world applications, classical discrete bivariate distributions frequently fail to accurately reflect observed data. This limitation has led to a growing interest in creating more flexible bivariate discrete models. The aim of this study is to propose a new adaptable bivariate discrete model, constructed via the maximization method, based on the discretized Fréchet–Weibull distribution (DFWD). According to the approach of Chakraborty [26], the DFWD was initially introduced by Das and Das [27] through the transformed-transformer (T-X) method. The probability mass function (pmf) of $Y\sim \mathrm{DFWD}(\alpha ,\beta ,m,k)$ was derived using the survival function approach to discretization and is expressed as follows

(1)$P[Y=y;\alpha ,\beta ,m,k]=exp\left\{-{\beta }^{\alpha }{\left({\displaystyle \frac{m}{y+1}}\right)}^{\alpha k}\right\}-exp\left\{-{\beta }^{\alpha }{\left({\displaystyle \frac{m}{y}}\right)}^{\alpha k}\right\},$

where $y\in {\mathbf{Z}}_{+}$ and the parameters include $\alpha >0,\beta >0,m>0,k>0$. The parameters $\alpha$ and k are the shape parameters, whereas $\beta$ and m are the scale parameters. The corresponding cumulative distribution function (cdf) can be formulated as

(2)$F_Y(y;\alpha ,\beta ,m,k)=exp\left\{-{\beta }^{\alpha }{\left({\displaystyle \frac{m}{y+1}}\right)}^{\alpha k}\right\};y\in {\mathbf{Z}}_{+}.$

After re-parameterization, $\theta =exp\left\{-{\beta }^{\alpha }{m}^{\alpha k}\right\}$ and $\lambda =\alpha k$, the cdf in Equation (2) can be written as:

(3)$F_Y(y;\theta ,\lambda )={\theta }^{{(y+1)}^{-\lambda }};y\in {\mathbf{Z}}_{+},$

where $0<\theta <1$ and $\lambda >0,\forall (\alpha ,\beta ,m,k)>0$. The DFWD is a versatile discrete model capable of fitting datasets with increasing, decreasing, and bathtub-shaped hazard rate functions. It can be applied to data exhibiting both positive and negative skewness under various kurtosis shapes. Additionally, the DFWD is suitable for modeling both over-dispersion and under-dispersion. Building on these advantages, this paper proposes a flexible bivariate extension of the DFWD, termed the bivariate discretized Fréchet-Weibull distribution (BDFWD). This extension utilizes the maximization methodology for constructing bivariate discrete models, which will be discussed in detail in Section 2. The motivations for introducing the BDFWD model are as follows:

1. Closed-form Expressions: The joint probability mass function, joint survival function, and joint reversed (hazard) rate function are provided in closed forms, enabling straightforward and efficient application across a wide range of fields.

2. Flexibility in Handling Over-dispersion and Zero-inflation: Many real-world datasets exhibit features such as over-dispersion (where observed variability exceeds that expected under a simple distribution) and zero-inflation (an excess of zero counts). Traditional models often face challenges addressing these characteristics, necessitating the development of the BDFWD model. BDFWD is specifically designed to manage both over-dispersion and zero-inflation effectively, making it highly applicable in various practical contexts.

3. Applications to Positively Ordered, Quadrant-dependent, and Discrete Data: The BDFWD model supports the analysis of data with total positivity of order two (TP2) properties and positive quadrant dependence (PQD). The TP2 is a property of bivariate distributions that reflects a specific positive dependence between two random variables, where an increase in one variable is associated with an increase in the other. This indicates a positive association and synergy between the variables. In contrast, PQD indicates that if one variable is higher, the likelihood of the other variable also being higher increases. Both TP2 and PQD are important in fields like finance and reliability analysis, as they help in understanding relationships between dependent variables for effective modeling and prediction.

4. Adaptable Joint Hazard Rate Function: The joint hazard rate function in BDFWD is adaptable to various shapes, depending on parameter configurations, providing flexibility for modeling diverse hazard rate behaviors. This adaptability allows for tailored modeling based on specific data requirements.

5. The stress-strength model relies on the parameters of the bivariate distribution.

6. It can be applied to a maintenance model or a stress model.

7. Effective Marginals for Hazard Rate Analysis: The marginals in BDFWD support robust analysis of different hazard rate profiles, enabling detailed examination and assessment of hazard rates across diverse scenarios.

8. Suitability for Asymmetric Datasets: BDFWD is particularly effective for modeling asymmetric datasets, especially those that are right-skewed or have varying levels of kurtosis. This feature makes it valuable in applications where data exhibit non-symmetric distributions.

9. Model Performance and Competitiveness: Comparative studies demonstrate the BDFWD’s superior performance over other bivariate discrete models in terms of model fit across various datasets. This highlights BDFWD’s robustness and versatility as a tool for modeling dependent discrete data, making it an excellent choice for practitioners.

The remainder of this paper is structured as follows: Section 2 offers a comprehensive development of the BDFWD, including the derivation of the joint probability mass function and joint survival function, along with their graphical representations. Section 3 presents several significant propositions related to the BDFWD. Section 4 describes various statistical and reliability characteristics of the BDFWD. Section 5 addresses the estimation of the model’s unknown parameters using the maximum likelihood method. Section 6 discusses the insights gained from the simulation results. Section 7 demonstrates the application of the model to three bivariate count datasets. Finally, Section 8 concludes the study with a summary of the key findings.

2. Bivariate Discretized Fréchet-Weibull Distribution

A bivariate random vector $\mathbf{Y}=(Y_1,Y_2)$ adheres to a bivariate model if and only if it can be expressed as $Y_i=max\{V_i,V_3\}$; $i=1,2$ where $V_1$, $V_2$, and $V_3$ are independent random variables. The maximization technique is employed to propose a novel flexible bivariate distribution that has not been explored in the statistical literature previously, specifically tailored for two-dimensional discrete datasets. To derive the joint cdf of BDFWD, we have used the cdf of DFWD as stated in Equation (3). Suppose $V_1\sim DFWD({\theta }_1,\lambda )$, $V_2\sim DFWD({\theta }_2,\lambda )$, and $V_3\sim DFWD({\theta }_3,\lambda )$, and they are independently distributed, and if $Y_1=max\{V_1,V_3\}$ and $Y_2=max\{V_2,V_3\}$, then we can say that the bivariate vector $\mathbf{Y}$ has a BDFWD with the parameter vector $\Psi ={({\theta }_1,{\theta }_2,{\theta }_3,\lambda )}^T$. We denote this bivariate discrete distribution by $BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$. If $\mathbf{Y}\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$, then the joint cdf of $\mathbf{Y}$ for $y_1,y_2\in {\mathbf{Z}}_{+}$ and for $z=min\{y_1,y_2\}$, which is given by:

(4)$\begin{array}{ccc}\hfill F_{Y_1,Y_2}(y_1,y_2;\Psi )& =& {\theta }_1^{{(y_1+1)}^{-\lambda }}{\theta }_2^{{(y_2+1)}^{-\lambda }}{\theta }_3^{{(z+1)}^{-\lambda }}\hfill \\ & =& F_{DFWD}(y_1;{\theta }_1,\lambda )F_{DFWD}(y_2;{\theta }_2,\lambda )F_{DFWD}(z;{\theta }_3,\lambda )\hfill \\ & =& \left\{\begin{array}{cc}{F}_{DFWD}(y_1;{\theta }_1{\theta }_3,\lambda )F_{DFWD}(y_2;{\theta }_2,\lambda )\hfill & ;\mathrm{if}0<y_1<y_2<\infty \hfill \\ F_{DFWD}(y_1;{\theta }_1,\lambda )F_{DFWD}(y_2;{\theta }_2{\theta }_3,\lambda )\hfill & ;\mathrm{if}0<y_2<y_1<\infty \hfill \\ F_{DFWD}(y;{\theta }_1{\theta }_2{\theta }_3,\lambda )\hfill & ;\mathrm{if}0<y_1=y_2=y<\infty .\hfill \end{array}\right.\hfill \end{array}$

Proposition 1.

If $(Y_1,Y_2)\sim$ $\mathrm{BDFWD}({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$, then the marginal distributions of $Y_i(i=1,2)$ are DFWD with marginal cdf

(5) $F_{Y_i}\left(y_i\right)=F_{DFWD}(y_i;{\theta }_i{\theta }_3,\lambda );i=1,2.$

Proof. 

The marginal cdf of $Y_i(i=1,2)$, denoted by $F_{Y_i}\left(y_i\right)$, can be obtained as

$\begin{array}{ccc}\hfill F_{Y_i}\left(y_i\right)& =& P[Y_i\le y_i]\hfill \\ & =& P[max\{V_i,V_3\}\le y_i]\hfill \\ & =& P[V_i\le y_i,V_3\le y_i]\hfill \\ & =& {\theta }_i^{{(y_i+1)}^{-\lambda }}{\theta }_3^{{(y_i+1)}^{-\lambda }}\hfill \\ & =& {\left({\theta }_i{\theta }_3\right)}^{{(y_i+1)}^{-\lambda }}=F_{DFWD}(y_i;{\theta }_i{\theta }_3,\lambda ).\hfill \end{array}$

Hence the result. □

2.1. The Joint Probability Mass Function of BDFWD

The corresponding joint probability mass function (pmf) of $\mathbf{Y}\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$ for $y_1,y_2\in {\mathbf{Z}}_{+}$ is given by

(6)$\begin{array}{c}\hfill f_{Y_1,Y_2}(y_1,y_2;\Psi )=\left\{\begin{array}{cc}{f}_1(y_1,y_2;\Psi )\hfill & ;\mathrm{if}0<y_1<y_2<\infty \hfill \\ f_2(y_1,y_2;\Psi )\hfill & ;\mathrm{if}0<y_2<y_1<\infty \hfill \\ f_3(y;\Psi )\hfill & ;\mathrm{if}0<y_1=y_2=y<\infty ,\hfill \end{array}\right.\end{array}$

where

$\begin{array}{ccc}\hfill f_1(y_1,y_2;\Psi )& =& f_{DFWD}(y_1;{\theta }_1{\theta }_3,\lambda )f_{DFWD}(y_2;{\theta }_2,\lambda ),\hfill \\ \hfill f_2(y_1,y_2;\Psi )& =& f_{DFWD}(y_1;{\theta }_1,\lambda )f_{DFWD}(y_2;{\theta }_2{\theta }_3,\lambda ),\hfill \\ \hfill f_3(y;\Psi )& =& F_{DFWD}(y;{\theta }_2,\lambda )f_{DFWD}(y;{\theta }_1{\theta }_3,\lambda )-F_{DFWD}(y-1;{\theta }_2{\theta }_3,\lambda )f_{DFWD}(y;{\theta }_1,\lambda ).\hfill \end{array}$

The expressions of $f_1(y_1,y_2;\Psi ),f_2(y_1,y_2;\Psi )$ and $f_3(y;\Psi )$ for $y_1,y_2\in {\mathbf{Z}}_{+}$ can be easily obtained by using the following relation:

(7)$f_{Y_1,Y_2}(y_1,y_2;\Psi )=F_{Y_1,Y_2}(y_1,y_2;\Psi )-F_{Y_1,Y_2}(y_1-1,y_2;\Psi )-F_{Y_1,Y_2}(y_1,y_2-1;\Psi )+F_{Y_1,Y_2}(y_1-1,y_2-1;\Psi ).$

In Figure 1, the scatter plots of the joint pmf of BDFWD are presented for different parameter values. From the scatter plots, it is clear that the joint pmf of BDFWD can take different shapes depending upon the choice of the parameter values $({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$. It is also seen that the joint mass has a long right tail as compared to its left tail.

2.2. The Joint Survival Function

The joint survival (reliability) function of $\mathbf{Y}\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$ for $y_1,y_2\in {\mathbf{Z}}_{+}$ is given by

(8)$\begin{array}{ccc}\hfill S_{Y_1,Y_2}(y_1,y_2;\Psi )& =& 1-F_{Y_1}(y_1;{\theta }_1{\theta }_3,\lambda )-F_{Y_2}(y_2;{\theta }_2{\theta }_3,\lambda )+F_{Y_1,Y_2}(y_1,y_2;\Psi )\hfill \\ & =& \left\{\begin{array}{cc}{S}_1(y_1,y_2;\Psi )\hfill & ;\mathrm{if}0<y_1<y_2<\infty \hfill \\ S_2(y_1,y_2;\Psi )\hfill & ;\mathrm{if}0<y_2<y_1<\infty \hfill \\ S_3(y;\Psi )\hfill & ;\mathrm{if}0<y_1=y_2=y<\infty ,\hfill \end{array}\right.\hfill \end{array}$

where

$\begin{array}{ccc}\hfill S_1(y_1,y_2;\Psi )& =& 1-F_{DFWD}(y_1;{\theta }_1{\theta }_3,\lambda )-F_{DFWD}(y_2;{\theta }_2{\theta }_3,\lambda )+\hfill \\ & & F_{DFWD}(y_1;{\theta }_1{\theta }_3,\lambda ).F_{DFWD}(y_2;{\theta }_2,\lambda ),\hfill \\ \hfill S_2(y_1,y_2;\Psi )& =& 1-F_{DFWD}(y_1;{\theta }_1{\theta }_3,\lambda )-F_{DFWD}(y_2;{\theta }_2{\theta }_3,\lambda )+\hfill \\ & & F_{DFWD}(y_1;{\theta }_1,\lambda ).F_{DFWD}(y_2;{\theta }_2{\theta }_3,\lambda ),\hfill \\ \hfill S_3(y;\Psi )& =& 1-F_{DFWD}(y;{\theta }_1{\theta }_3,\lambda )-F_{DFWD}(y;{\theta }_2{\theta }_3,\lambda )+F_{DFWD}(y;{\theta }_1{\theta }_2{\theta }_3,\lambda ).\hfill \end{array}$

In Figure 2, the scatter plots of the joint survival function (sf) of BDFWD are presented for different values of the model parameters. This figure explains the change in the shape of the joint sf for varying values of ${\theta }_1,{\theta }_2,{\theta }_3$ and $\lambda$, depicting the flexibility of BDFWD.

3. Key Propositions Regarding the BDFWD

This section derives and thoroughly discusses several useful propositions related to the BDFWD model.

Proposition 2.

If $(X_1,X_2)\sim BFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$, then for the bivariate vector $(Y_1,Y_2)$, defined as $Y_1=\lfloor X_1\rfloor$ and $Y_2=\lfloor X_2\rfloor$, has BDFWD with parameter vector ${({\theta }_1,{\theta }_2,{\theta }_3,\lambda )}^T$, where $\lfloor .\rfloor$ is the floor function.

See Appendix A for the proof.

Remark 1.

This result indicates that the proposed bivariate discrete distribution BDFWD(${\theta }_1,{\theta }_2,{\theta }_3,\lambda$) is a natural discrete version of the bivariate continuous counterpart BFWD(${\theta }_1,{\theta }_2,{\theta }_3,\lambda$).

Proposition 3.

Suppose $(Y_1,Y_2)\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$ then $max\{Y_1,Y_2\}\sim DFWD(\theta ,\lambda )$; where $\theta ={\theta }_1{\theta }_2{\theta }_3$.

See Appendix A for the proof.

Proposition 4.

Suppose $(Y_{i1},Y_{i2})\sim BDFWD({\theta }_{i1},{\theta }_{i2},{\theta }_{i3},\lambda )$ for $i=1,2,\dots ,n$, and they are independently distributed. If $Z_1=max\{Y_{11},Y_{21},Y_{31},\dots ,Y_{n1}\}$ and $Z_2=max\{Y_{12},Y_{22},Y_{32},\dots ,Y_{n2}\}$, then $(Z_1,Z_2)\sim BDFWD\left({\prod }_{i=1}^n{\theta }_{i1},{\prod }_{i=1}^n{\theta }_{i2},{\prod }_{i=1}^n{\theta }_{i3},\lambda \right)$.

See Appendix A for the proof.

Corollary 1.

Let $(Y_{i1},Y_{i2});(i=1,2,\dots ,n)$ be independently and identically distributed bivariate random variables accordingly as $BDFWD({\theta }_1{\theta }_2{\theta }_3,\lambda )$. If $Z_1=max\{Y_{11},Y_{21},Y_{31},\dots ,Y_{n1}\}$ and $Z_2=max\{Y_{12},Y_{22},Y_{32},\dots ,Y_{n2}\}$, then $(Z_1,Z_2)\sim BDFWD\left({\theta }_1^n,{\theta }_2^n,{\theta }_3^n,\lambda \right).$

Proposition 5.

If $(Y_1,Y_2)\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$, then the stress–strength reliability is given by $P[Y_1<Y_2]={\sum }_{y=0}^{\infty }{\left({\theta }_1{\theta }_3\right)}^{{(y+1)}^{-\lambda }}\left[{\left({\theta }_2{\theta }_3\right)}^{{(y+2)}^{-\lambda }}-{\left({\theta }_2{\theta }_3\right)}^{{(y+1)}^{-\lambda }}\right].$

See Appendix A for the proof.

Proposition 6.

Suppose $(Y_1,Y_2)\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$, then it has the total positively of order two $\left(T{P}_2\right)$ property.

See Appendix A for the proof.

Proposition 7.

If $(Y_1,Y_2)\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$, then $Y_1$ and $Y_2$ are positive quadrant dependent $\left(PQD\right)$.

See Appendix A for the proof.

Proposition 8.

If $(Y_1,Y_2)\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$, then the median of $Y_1$ and $Y_2$ can be expressed as $M_{Y_i}={\displaystyle \frac{1}{\lambda }}\left[{\displaystyle \frac{log\left(2\right)}{log\left({\theta }_i{\theta }_3\right)}}\right]-1;i=1,2.$

See Appendix A for the proof.

Proposition 9.

If the random variables $X_1,X_2,X_3$ are independently distributed as $Exponential$$\left(p_i\right);0<p_i<1$ and $U_i={\lfloor X_i\rfloor }^{-{\displaystyle \frac{1}{\lambda }}}-1;\lambda >0i=1,2,3$, then $\left(max\{U_1,U_3\},max\{U_2,U_3\}\right)$ follows BDFWD with parameters $({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$. Also, $p_i=log\left({\displaystyle \frac{1}{{\theta }_i}}\right)$ and $\lfloor .\rfloor$ is the floor function.

See Appendix A for the proof.

Proposition 10.

If the discrete random variables $X_1,X_2,X_3$ are independently distributed as $Geometric(1-{\theta }_i);0<{\theta }_i<1$ and $U_i=\left(X_i^{-{\displaystyle \frac{1}{\lambda }}}\right)-1;\lambda >0;i=1,2,3$, then the pair $\left(max\{U_1,U_3\},max\{U_2,U_3\}\right)$ follows BDFWD with parameters $({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$.

See Appendix A for the proof.

4. Statistical and Reliability Characteristics

4.1. The Joint Probability Generating Function Along with Associated Measures

If the bivariate vector $\mathbf{Y}\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$, then the joint probability generating function (pgf) of $Y_1$ and $Y_2$ for $|t_1|<1$ and $|t_2|<1$ can be written as

(9)$\begin{array}{cc}\hfill G_{Y_1,Y_2}(t_1,t_2)& =E\left(t_1^{Y_1}{t}_2^{Y_2}\right)\hfill \\ \hfill & =\sum _{j,i=0}^{\infty }P[Y_1=i,Y_2=j]t_1^i{t}_2^j\hfill \\ \hfill & =\sum _{j=0}^{\infty }\sum _{i=0}^{j-1}\left[{{\left({\theta }_1{\theta }_3\right)}^{(i+1)}}^{-\lambda }-{{\left({\theta }_1{\theta }_3\right)}^i}^{-\lambda }\right].\left[{{\theta }_2^{(j+1)}}^{-\lambda }-{{\theta }_2^j}^{-\lambda }\right]t_1^i{t}_2^j\hfill \\ \hfill & +\sum _{j=0}^{\infty }\sum _{i=j+1}^{\infty }\left[{{\theta }_1^{(i+1)}}^{-\lambda }-{{\theta }_1^i}^{-\lambda }\right].\left[{{\left({\theta }_2{\theta }_3\right)}^{(j+1)}}^{-\lambda }-{{\left({\theta }_2{\theta }_3\right)}^j}^{-\lambda }\right]t_1^i{t}_2^j\hfill \\ \hfill & +\sum _{i=0}^{\infty }{{\theta }_2^{(i+1)}}^{-\lambda }\left[{{\left({\theta }_1{\theta }_3\right)}^{(j+1)}}^{-\lambda }-{{\left({\theta }_1{\theta }_3\right)}^j}^{-\lambda }\right]t_1^i{t}_2^i\hfill \\ \hfill & -\sum _{i=0}^{\infty }{{\left({\theta }_2{\theta }_3\right)}^{(i+1)}}^{-\lambda }\left[{{\theta }_1^{(i+1)}}^{-\lambda }-{{\theta }_1^i}^{-\lambda }\right]t_1^i{t}_2^i;|t_1|,|t_2|<1.\hfill \end{array}$

Hence, using the above joint pgf of BDFWD, different moments and product moments of BDFWD can be obtained as infinite series. For a bivariate vector $(Y_1,Y_2)$, taking partial derivative of the joint pgf $G_{Y_1,Y_2}(t_1,t_2)$ for “$r$” times with respect to “$t_i$”, gives the $r^{th}$ factorial moment of $Y_i;(i=1,2)$, and this can be expressed as

(10)${{\mu }_{Y_i}^{\prime }}_{\left(r\right)}={\displaystyle \frac{{\partial }^r}{\partial t_i^r}}{G}_{Y_1,Y_2}(t_1,t_2){|}_{t_1=1,t_2=1};i=1,2,$

where ${{\mu }_{Y_i}^{\prime }}_{\left(r\right)}$ is the $r^{th}$ factorial moment of $Y_i;(i=1,2)$. Taking the first partial derivative of Equation (10) yields the following relations:

$\begin{array}{ccc}\hfill E\left(Y_1\right)& =& {{\mu }_{Y_1}^{\prime }}_{\left(1\right)}={\displaystyle \frac{\partial }{\partial t_1}}{G}_{Y_1,Y_2}(t_1,t_2){|}_{t_1=1,t_2=1},\hfill \\ \hfill E\left(Y_2\right)& =& {{\mu }_{Y_2}^{\prime }}_{\left(1\right)}={\displaystyle \frac{\partial }{\partial t_2}}{G}_{Y_1,Y_2}(t_1,t_2){|}_{t_1=1,t_2=1}.\hfill \end{array}$

The double partial derivative of Equation (10) yields:

${{\mu }_{Y_1}^{\prime }}_{\left(2\right)}={\displaystyle \frac{{\partial }^2}{\partial t_1^2}}{G}_{Y_1,Y_2}(t_1,t_2){|}_{t_1=1,t_2=1}\mathrm{and}{{\mu }_{Y_2}^{\prime }}_{\left(2\right)}={\displaystyle \frac{{\partial }^2}{\partial t_2^2}}{G}_{Y_1,Y_2}(t_1,t_2){|}_{t_1=1,t_2=1}.$

Thus, the variances of $Y_1$ and $Y_2$ can be obtained as

$V\left(Y_1\right)={{\mu }_{Y_1}^{\prime }}_{\left(2\right)}+{{\mu }_{Y_1}^{\prime }}_{\left(1\right)}-{\left\{{{\mu }_{Y_1}^{\prime }}_{\left(1\right)}\right\}}^2\mathrm{and}V\left(Y_2\right)={{\mu }_{Y_2}^{\prime }}_{\left(2\right)}+{{\mu }_{Y_2}^{\prime }}_{\left(1\right)}-{\left\{{{\mu }_{Y_2}^{\prime }}_{\left(1\right)}\right\}}^2.$

Joint skewness, say $JSk(Y_1,Y_2)$ can be computed by examining the third central moments and their joint distribution. On the other hand, joint kurtosis, say $JKu(Y_1,Y_2)$, involves the fourth central moments and their combinations, requiring a more complex formula where

$\begin{array}{ccc}\hfill JSk(Y_1,Y_2)& =& {\displaystyle \frac{E[(Y_1-E\left(Y_1\right)(Y_2-E\left(Y_2\right))(Y_1-E\left(Y_1\right)+Y_2-E\left(Y_2\right))]}{{\left[Var\left(Y_1\right)Var\left(Y_2\right)\right]}^{3/2}}},\hfill \\ \hfill JKu(Y_1,Y_2)& =& {\displaystyle \frac{E\left[\right(Y_1-E{\left(Y_1\right)}^4(Y_2-E{\left(Y_2\right)}^4]}{{\left[Var\left(Y_1\right)Var\left(Y_2\right)\right]}^2}}.\hfill \end{array}$

Table 1 presents numerical values for both $JSk(Y_1,Y_2)$ and $JKu(Y_1,Y_2)$ based on different model parameters.

Based on various numerical examples and as demonstrated in Figure 1, the BDFWD model effectively facilitates the discussion and evaluation of asymmetric data with diverse kurtosis shapes. This suggests that the model is versatile in analyzing data distributions with varying degrees of tail heaviness and asymmetry.

4.2. Computational Assessment Using Various Statistical Metrics

In this section, numerical computations were performed to analyze the patterns of the covariance, Pearson’s correlation coefficient, Sperman’s rho ($\rho$), and Kendall’s tau ($\tau$) of a random pair ($Y_1,Y_2$) randomly drawn from BDFWD(${\theta }_1,{\theta }_2,{\theta }_3,\lambda$). From Equation (10), we can have

$E\left(Y_1{Y}_2\right)={\displaystyle \frac{{\partial }^2}{\partial t_1\partial t_2}}{G}_{Y_1,Y_2}(t_1,t_2){|}_{t_1=1,t_2=1},$

which can be used to obtain the covariance of $(Y_1,Y_2)$, using the relation $Cov(Y_1,Y_2)=E\left(Y_1{Y}_2\right)-E\left(Y_1\right)E\left(Y_2\right)$. In Table 2, the above-mentioned descriptive statistic values were calculated by varying the values of ${\theta }_i$ ($i=1,2$ and 3) and $\lambda$, with the remaining three parameter values fixed.

From the above Table 2, the following observations can be noted:

• The values of covariance, correlation, Sperman’s rho, and Kendall’s tau were greater than zero (which is a direct consequence of PQD).

• For $0\le {\theta }_i\le 1$, as the value of ${\theta }_i$ increases ($i=1,2$ and 3), with a fixed value of the remaining three parameters, and the values of covariance, correlation, Sperman’s rho, and Kendall’s tau remained constant.

• For fixed values of ${\theta }_1,{\theta }_2$ and ${\theta }_3$ with $\lambda \to \infty$, the values of covariance, correlation, Sperman’s rho, and Kendall’s tau became constant.

• The values of correlation, Sperman’s rho, and Kendall’s tau were equal for fixed values of ${\theta }_1,{\theta }_2$ and ${\theta }_3$ with $\lambda \to \infty$.

4.3. Median Correlation Coefficient

The median of a bivariate distribution can be determined using marginal medians, geometric medians, conditional medians, or median contours, depending on the distribution’s specific structure and characteristics. The selection of method is guided by the nature of the bivariate relationship and its intended application. In this section, the median of the BDFWD is determined using marginal medians, following the approach of Domma [28], where the median correlation coefficient $M_{Y_1,Y_2}$ is expressed as $M_{Y_1,Y_2}=4F_{Y_1,Y_2}(M_{Y_1},M_{Y_2})-1$, with $M_{Y_1}$ and $M_{Y_2}$ representing the medians of $Y_1$ and $Y_2$, respectively. If $Y_1\sim DFWD$(${\theta }_1{\theta }_3,\lambda$) and $Y_2\sim DFWD$(${\theta }_2{\theta }_3,\lambda$), then

(11)$M_{Y_1,Y_2}=\left\{\begin{array}{c}4F_{DFWD}\left(M_{Y_2};{\theta }_2{\theta }_3,\lambda \right)\times F_{DFWD}\left(M_{Y_1};{\theta }_1{\theta }_3,\lambda \right)-1\mathrm{if}{y}_1\le y_2\hfill \\ 4F_{DFWD}\left(M_{Y_1};{\theta }_1{\theta }_3,\lambda \right)\times F_{DFWD}\left(M_{Y_2};{\theta }_2{\theta }_3,\lambda \right)-1\mathrm{if}{y}_1>y_2,\hfill \end{array}\right.$

where $M_{Y_i}={\displaystyle \frac{1}{\lambda }}\left[{\displaystyle \frac{log2}{log{\theta }_i{\theta }_3}}\right]-1$. The median of a bivariate model, determined through marginal medians, is valuable across numerous fields by providing insights into the central tendencies of variables.

4.4. The Bivariate Hazard Rate Function

For the two-dimensional random variables ($Y_1,Y_2$), with the joint pmf $f(y_1,y_2)$, and the joint sf $S(y_1,y_2)$, Basu [29] defined the bivariate hazard rate function (BHRF) as:

(12)$h(y_1,y_2)={\displaystyle \frac{f(y_1,y_2)}{S(y_1,y_2)}}.$

The BHRF of $\mathbf{Y}\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$ for $y_1,y_2\in {\mathbf{Z}}_{+}$ can be presented as

(13)$\begin{array}{ccc}\hfill h_{Y_1,Y_2}(y_1,y_2;\Psi )& =& \left\{\begin{array}{cc}{h}_1(y_1,y_2;\Psi )\hfill & ;\mathrm{if}0<y_1<y_2<\infty \hfill \\ h_2(y_1,y_2;\Psi )\hfill & ;\mathrm{if}0<y_2<y_1<\infty \hfill \\ h_3(y;\Psi )\hfill & ;\mathrm{if}0<y_1=y_2=y<\infty .\hfill \end{array}\right.\hfill \end{array}$

where the terms $h_1(y_1,y_2;\Psi )$, $h_1(y_1,y_2;\Psi )$ and $h_1(y_1,y_2;\Psi )$ are defined as follows:

$h_1(y_1,y_2;\Psi )={\displaystyle \frac{f_1(y_1,y_2;\Psi )}{S_1(y_1,y_2;\Psi )}},h_2(y_1,y_2;\Psi )={\displaystyle \frac{f_2(y_1,y_2;\Psi )}{S_2(y_1,y_2;\Psi )}},\mathrm{and}{h}_3(y;\Psi )={\displaystyle \frac{f_3(y;\Psi )}{S_3(y;\Psi )}}$

The terms $f_1(y_1,y_2;\Psi )$, $f_2(y_1,y_2;\Psi )$, and $f_3(y;\Psi )$ are the joint pmf of BDFWD as defined in Equation (6). And the terms (the joint sf) of BDFWD $S_1(y_1,y_2;\Psi )$, $S_2(y_1,y_2;\Psi )$, and $S_3(y;\Psi )$ are as defined in Equation (8). In Figure 3, scatter plots of the BHRF for the BDFWD are displayed across various model parameter values, visually showcasing the distinct shapes that emerge with different parameter choices. This effectively demonstrates the flexibility of the BDFWD model.

Suppose we have a parallel system consisting of two components, then the BHRF can be defined as a vector, as discussed by Cox [30], which is denoted as $h\left(\underline{\mathbf{y}}\right)$ and is given using the following equation

(14)$h\left(\underline{\mathbf{y}}\right)=\left(h\left(y\right),h_{12}\left(y_1\right|y_2),h_{21}\left(y_2\right|y_1)\right),$

where $h\left(\underline{\mathbf{y}}\right)$ is the vector form of BHRF for $(Y_1,Y_2)$. This is useful to measure the total life span of a two-component parallel system. Here, the first element $h\left(y\right)$, represents the hazard rate function of the system using the information that both the components survived beyond y. The second element $h_{12}\left(y_1\right|y_2)$ represents the span of the hazard rate function for the first component given that it has survived to an age of $y_1$ and the second component failed at age $y_2$. The third element $h_{21}\left(y_2\right|y_1)$, can also be similarly defined. Let $(Y_1,Y_2)$ be a random vector from BDFWD(${\theta }_1,{\theta }_2,{\theta }_3,\lambda$), then the elements of $h\left(\underline{\mathbf{y}}\right)$ can be defined as stated in Equation (14) where

$\begin{array}{ccc}\hfill h\left(y\right)& =& {\displaystyle \frac{f_T\left(y\right)}{S_T\left(y\right)}}{|}_{T=min(Y_1,Y_2)}={\displaystyle \frac{-\lambda [A(y;{\theta }_1{\theta }_3)+A(y;{\theta }_2{\theta }_3)-A(y;{\theta }_1{\theta }_2{\theta }_3)]}{S_T\left(y\right)}},\hfill \\ \hfill h_{12}\left(y_1\right|y_2)& =& {\displaystyle \frac{\frac{{\partial }^2}{\partial y_1\partial y_2}S(y_1,y_2)}{\frac{\partial }{\partial y_2}S(y_1,y_2)}}{|}_{y_1>y_2}={\displaystyle \frac{\lambda A(y_1;{\theta }_1)}{1-F_{DFWD}(y_1;{\theta }_1,\lambda )}},\hfill \\ \hfill h_{21}\left(y_2\right|y_1)& =& {\displaystyle \frac{\frac{{\partial }^2}{\partial y_1\partial y_2}S(y_1,y_2)}{\frac{\partial }{\partial y_1}S(y_1,y_2)}}{|}_{y_1<y_2}={\displaystyle \frac{\lambda A(y_2;{\theta }_2)}{1-F_{DFWD}(y_2;{\theta }_2,\lambda )}},\hfill \end{array}$

where $A(y;\theta )=ln\left(\theta \right){\theta }^{{(y+1)}^{-\lambda }}$.

4.5. The Bivariate Reversed Hazard Rate Function

For the two-dimensional random variables ($Y_1,Y_2$), the bivariate reversed hazard rate function (BRHRF) is given as

(15)$h^{*}(y_1,y_2)={\displaystyle \frac{f(y_1,y_2)}{F(y_1,y_2)}}.$

where $f(y_1,y_2)$ is the joint pmf and $F(y_1,y_2)$ is the joint cdf. The BRHRF of $\mathbf{Y}\sim BDFWD({\theta }_1,$${\theta }_2,{\theta }_3,\lambda )$ for $y_1,y_2\in {\mathbf{Z}}_{+}$ can be presented by:

(16)$\begin{array}{ccc}\hfill h_{Y_1,Y_2}^{*}(y_1,y_2;\Psi )& =& \left\{\begin{array}{cc}{h}_1^{*}(y_1,y_2;\Psi )\hfill & ;\mathrm{if}0<y_1<y_2<\infty \hfill \\ h_2^{*}(y_1,y_2;\Psi )\hfill & ;\mathrm{if}0<y_2<y_1<\infty \hfill \\ h_3^{*}(y;\Psi )\hfill & ;\mathrm{if}0<y_1=y_2=y<\infty .\hfill \end{array}\right.\hfill \end{array}$

where the terms $h_1(y_1,y_2;\Psi )$, $h_1(y_1,y_2;\Psi )$, and $h_1(y_1,y_2;\Psi )$ are defined as follows:

$h_1(y_1,y_2;\Psi )={\displaystyle \frac{f_1(y_1,y_2;\Psi )}{F_1(y_1,y_2;\Psi )}},h_2(y_1,y_2;\Psi )={\displaystyle \frac{f_2(y_1,y_2;\Psi )}{F_2(y_1,y_2;\Psi )}},\mathrm{and}{h}_3(y;\Psi )={\displaystyle \frac{f_3(y;\Psi )}{F_3(y;\Psi )}}$

The terms $f_1(y_1,y_2;\Psi )$, $f_2(y_1,y_2;\Psi )$, and $f_3(y;\Psi )$ are the joint pmf of BDFWD as defined in Equation (6). And the terms of the joint cdf of BDFWD $F_1(y_1,y_2;\Psi )$, $F_2(y_1,y_2;\Psi )$, and $F_3(y;\Psi )$ are as defined in Equation (4).

4.6. Joint Failure Analysis

To examine the monotonicity of the hazard components in the BDFWD, it is necessary to derive a mathematical formula based on the conditional distribution (see, Glaser [31]), say $\eta (.)$, where

(17)$\eta (.;{\theta }_1,{\theta }_2,{\theta }_3,\lambda )=\left\{\begin{array}{c}-ln{\displaystyle \frac{\partial }{\partial y_1}}f\left(y_1\right|Y_2>y_2)\mathrm{if}{y}_1<y_2\hfill \\ -ln{\displaystyle \frac{\partial }{\partial y_2}}f\left(y_2\right|Y_1>y_1)\mathrm{if}{y}_1>y_2,\hfill \end{array}\right.$

and

$\begin{array}{ccc}\hfill f\left(y_1\right|Y_2& >& y_2)={\displaystyle \frac{h_1(y_1,y_2)S_1(y_1,y_2)}{1-F_{DFWD}(y_2;{\theta }_2{\theta }_3,\lambda )}},\hfill \\ \hfill f\left(y_2\right|Y_1& >& y_1)={\displaystyle \frac{h_2(y_1,y_2)S_2(y_1,y_2)}{1-F_{DFWD}(y_1;{\theta }_1{\theta }_3,\lambda )}}.\hfill \end{array}$

By examining the joint probability mass function of the BDFWD alongside its corresponding joint survival and joint hazard rate functions, it becomes evident that ${\displaystyle \frac{\partial }{\partial y_i}}\eta (y_i;{\theta }_1,$${\theta }_2,{\theta }_3,\lambda )<0;$$i=1,2$. Hence, $h_1(y_1,y_2)$ is decreasing in $y_1$ and $h_2(y_1,y_2)$ is decreasing in $y_2$.

5. Parameter Estimation: Maximum Likelihood Method

Let us consider a bivariate sample of size n, of the form $\{(y_{11},y_{21}),(y_{12},y_{22}),\dots ,(y_{1n},y_{2n})\}$ from BDFWD. Now, to estimate the values of unknown parameters (${\theta }_1,{\theta }_2,{\theta }_3,\lambda$) of BDFWD, using the maximum likelihood estimation (MLE) method, let us assume $I_1=\{y_{1j}<y_{2j}\}$, $I_2=\{y_{1j}>y_{2j}\}$ and $I_3=\{y_{1j}=y_{2j}=y_j\}$, such that $I=I_1\cup I_2\cup I_3$. Also, let the number of elements in $I_1$, $I_2$, and $I_3$ be $n_1$, $n_2$, and $n_3$ respectively, such that $n=n_1+n_2+n_3$. Using these notations, for the parameter vector, $\Omega ={({\theta }_1,{\theta }_2,{\theta }_3,\lambda )}^T$, the likelihood function of BDFWD is given by:

(18)$L(\Omega )=\prod _{j=1}^{n_1}{f}_1(y_{1j},y_{2j};\Omega ).\prod _{j=1}^{n_2}{f}_2(y_{1j},y_{2j};\Omega ).\prod _{j=1}^{n_3}{f}_3(y_j;\Omega ).$

Then, the log-likelihood function becomes as

(19)$\begin{array}{ccc}\hfill l(\Omega )& =& \sum _{j=1}^{n_1}ln{f}_1(y_{1j},y_{2j};\Omega )+\sum _{j=1}^{n_2}ln{f}_2(y_{1j},y_{2j};\Omega )+\sum _{j=1}^{n_3}ln{f}_3(y_j;\Omega )\hfill \\ & =& \sum _{j=1}^{n_1}ln\Phi (y_{1j};{\theta }_1{\theta }_3,\lambda )+\sum _{j=1}^{n_1}ln\Phi (y_{2j};{\theta }_2,\lambda )\hfill \\ & +& \sum _{j=1}^{n_2}ln\Phi (y_{1j};{\theta }_1{\theta }_3,\lambda )+\sum _{j=1}^{n_1}ln\Phi (y_{2j};{\theta }_2,\lambda )\hfill \\ & +& \sum _{j=1}^{n_3}ln\left\{{\theta }_2^{{(y_j+1)}^{-\lambda }}\Phi (y_j;{\theta }_1{\theta }_3,\lambda )-{\left({\theta }_2{\theta }_3\right)}^{{\left(y_j\right)}^{-\lambda }}\Phi (y_j;{\theta }_1,\lambda )\right\},\hfill \end{array}$

where $\Phi (y;\theta ,\lambda )=f_{DFWD}(y;\theta ,\lambda )={\theta }^{{(y+1)}^{\lambda }}-{\theta }^{y^{\lambda }}$. The maximum likelihood estimators of the parameters ${\theta }_1,{\theta }_2,{\theta }_3,and \lambda$ can be obtained by solving the partial derivative of Equation (19) with respect to ${\theta }_1,{\theta }_2,{\theta }_3$, and $\lambda$ and then equating the results to zero. The partial derivative equations of Equation (19) with respect to ${\theta }_1,{\theta }_2,{\theta }_3$, and $\lambda$ are as given in Equations (20), (21), (22), and (23), respectively.

(20)$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial l}{\partial {\theta }_1}}& =& \sum _{j=1}^{n_1}{\displaystyle \frac{{\Phi }_1(y_{1j}+1;{\theta }_1{\theta }_3,\lambda )-{\Phi }_1(y_{1j};{\theta }_1{\theta }_3,\lambda )}{\Phi (y_{1j};{\theta }_1{\theta }_3,\lambda )}}+\sum _{j=1}^{n_2}{\displaystyle \frac{{\Phi }_1(y_{1j}+1;{\theta }_1,\lambda )-{\Phi }_1(y_{1j};{\theta }_1,\lambda )}{\Phi (y_{1j};{\theta }_1,\lambda )}}\hfill \\ & +& \sum _{j=1}^{n_3}{\displaystyle \frac{{\theta }_2^{{(y_j+1)}^{-\lambda }}\left[{\Phi }_1(y_j+1;{\theta }_1{\theta }_3,\lambda )-{\Phi }_1(y_j;{\theta }_1{\theta }_3,\lambda )\right]}{{\theta }_2^{{(y_j+1)}^{-\lambda }}\Phi (y_j;{\theta }_1{\theta }_3,\lambda )-{\left({\theta }_2{\theta }_3\right)}^{{\left(y_j\right)}^{-\lambda }}\Phi (y_j;{\theta }_1,\lambda )}}\hfill \\ & -& \sum _{j=1}^{n_3}{\displaystyle \frac{{\left({\theta }_2{\theta }_3\right)}^{{\left(y_j\right)}^{-\lambda }}\left[{\Phi }_1(y_j+1;{\theta }_1,\lambda )-{\Phi }_1(y_j;{\theta }_1,\lambda )\right]}{{\theta }_2^{{(y_j+1)}^{-\lambda }}\Phi (y_j;{\theta }_1{\theta }_3,\lambda )-{\left({\theta }_2{\theta }_3\right)}^{{\left(y_j\right)}^{-\lambda }}\Phi (y_j;{\theta }_1,\lambda )}},\hfill \end{array}$

(21)$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial l}{\partial {\theta }_2}}& =& \sum _{j=1}^{n_1}{\displaystyle \frac{{\Phi }_1(y_{2j}+1;{\theta }_2,\lambda )-{\Phi }_1(y_{2j};{\theta }_2,\lambda )}{\Phi (y_{2j};{\theta }_2,\lambda )}}+\sum _{j=1}^{n_2}{\displaystyle \frac{{\Phi }_1(y_{2j}+1;{\theta }_2{\theta }_3,\lambda )-{\Phi }_1(y_{2j};{\theta }_2{\theta }_3,\lambda )}{\Phi (y_{2j};{\theta }_2{\theta }_3,\lambda )}}\hfill \\ & +& \sum _{j=1}^{n_3}{\displaystyle \frac{\Phi (y_j;{\theta }_1{\theta }_3,\lambda ).{\Phi }_1(y_j+1;{\theta }_2,\lambda )-\Phi (y_j;{\theta }_1,\lambda ).{\Phi }_1(y_j;{\theta }_2{\theta }_3,\lambda )}{{\theta }_2^{{(y_j+1)}^{-\lambda }}\Phi (y_j;{\theta }_1{\theta }_3,\lambda )-{\left({\theta }_2{\theta }_3\right)}^{{\left(y_j\right)}^{-\lambda }}\Phi (y_j;{\theta }_1,\lambda )}},\hfill \end{array}$

(22)$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial l}{\partial {\theta }_3}}& =& \sum _{j=1}^{n_1}{\displaystyle \frac{{\Phi }_1(y_{1j}+1;{\theta }_1{\theta }_3,\lambda )-{\Phi }_1(y_{1j};{\theta }_1{\theta }_3,\lambda )}{\Phi (y_{1j};{\theta }_1{\theta }_3,\lambda )}}+\sum _{j=1}^{n_2}{\displaystyle \frac{{\Phi }_1(y_{2j}+1;{\theta }_2{\theta }_3,\lambda )-{\Phi }_1(y_{2j};{\theta }_2{\theta }_3,\lambda )}{\Phi (y_{2j};{\theta }_2{\theta }_3,\lambda )}}\hfill \\ & +& \sum _{j=1}^{n_3}{\displaystyle \frac{{\theta }_2^{{(y_j+1)}^{-\lambda }}\left[{\Phi }_1(y_j+1;{\theta }_1{\theta }_3,\lambda )-{\Phi }_1(y_j;{\theta }_1{\theta }_3,\lambda )\right]}{{\theta }_2^{{(y_j+1)}^{-\lambda }}\Phi (y_j;{\theta }_1{\theta }_3,\lambda )-{\left({\theta }_2{\theta }_3\right)}^{{\left(y_j\right)}^{-\lambda }}\Phi (y_j;{\theta }_1,\lambda )}}\hfill \\ & -& \sum _{j=1}^{n_3}{\displaystyle \frac{\Phi (y_j;{\theta }_1,\lambda ).{\Phi }_1(y_j;{\theta }_2{\theta }_3,\lambda )}{{\theta }_2^{{(y_j+1)}^{-\lambda }}\Phi (y_j;{\theta }_1{\theta }_3,\lambda )-{\left({\theta }_2{\theta }_3\right)}^{{\left(y_j\right)}^{-\lambda }}\Phi (y_j;{\theta }_1,\lambda )}},\hfill \end{array}$

(23)$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial l}{\partial \lambda }}& =& \sum _{j=1}^{n_1}{\displaystyle \frac{{\Phi }_2(y_{1j};{\theta }_1{\theta }_3,\lambda )-{\Phi }_2(y_{1j}+1;{\theta }_1{\theta }_3,\lambda )}{\Phi (y_{1j};{\theta }_1{\theta }_3,\lambda )}}+\sum _{j=1}^{n_1}{\displaystyle \frac{{\Phi }_2(y_{2j};{\theta }_2,\lambda )-{\Phi }_2(y_{2j};{\theta }_2,\lambda )}{\Phi (y_{2j}+1;{\theta }_2,\lambda )}}\hfill \\ & +& \sum _{j=1}^{n_2}{\displaystyle \frac{{\Phi }_2(y_{1j};{\theta }_1,\lambda )-{\Phi }_2(y_{1j};{\theta }_1,\lambda )}{\Phi (y_{1j}+1;{\theta }_1,\lambda )}}+\sum _{j=1}^{n_2}{\displaystyle \frac{{\Phi }_2(y_{2j};{\theta }_2{\theta }_3,\lambda )-{\Phi }_2(y_{2j}+1;{\theta }_2{\theta }_3,\lambda )}{\Phi (y_{2j};{\theta }_2{\theta }_3,\lambda )}}\hfill \\ & +& \sum _{j=1}^{n_2}{\displaystyle \frac{{\theta }_2^{{(y_j+1)}^{-\lambda }}\left[{\Phi }_2(y_j;{\theta }_1{\theta }_3,\lambda )-{\Phi }_2(y_j+1;{\theta }_1{\theta }_3,\lambda )\right]+\Phi (y_j;{\theta }_1{\theta }_3,\lambda ){\Phi }_2(y_j+1;{\theta }_2,\lambda )}{{\theta }_2^{{(y_j+1)}^{-\lambda }}\Phi (y_j;{\theta }_1{\theta }_3,\lambda )-{\left({\theta }_2{\theta }_3\right)}^{{\left(y_j\right)}^{-\lambda }}\Phi (y_j;{\theta }_1,\lambda )}}\hfill \\ & -& \sum _{j=1}^{n_2}{\displaystyle \frac{{\left({\theta }_2{\theta }_3\right)}^{{(y_j+1)}^{-\lambda }}\left[{\Phi }_2(y_j;{\theta }_1,\lambda )-{\Phi }_2(y_j+1;{\theta }_1,\lambda )\right]+\Phi (y_j;{\theta }_1,\lambda ){\Phi }_2(y_j+1;{\theta }_2{\theta }_3,\lambda )}{{\theta }_2^{{(y_j+1)}^{-\lambda }}\Phi (y_j;{\theta }_1{\theta }_3,\lambda )-{\left({\theta }_2{\theta }_3\right)}^{{\left(y_j\right)}^{-\lambda }}\Phi (y_j;{\theta }_1,\lambda )}}\hfill \\ \hfill \end{array}$

where ${\Phi }_1(y;\theta ,\lambda )=y^{-\lambda }{\theta }^{y^{-\lambda }-1}$ and ${\Phi }_2(y;\theta ,\lambda )=ln\left(\theta \right)y{\theta }^{y^{-\lambda }}$. The solution of these four equations is not easy to solve. To obtain the MLEs, we need a numerical technique. Here in this study, the stats4 package in R-programming was used to obtain the MLEs of the parameters ${\theta }_1,{\theta }_2,{\theta }_3$ and $\lambda$ of BDFWD. We can derive the asymptotic confidence interval of BDFWD for the estimated parameter vector $\widehat{\Omega }={(\widehat{{\theta }_1},\widehat{{\theta }_2},\widehat{{\theta }_3},\widehat{\lambda })}^T$ by using the variance–covariance matrix $\mathcal{V}\left(\widehat{\Omega }\right)$, which was approximated using the observed Fisher’s information matrix $\mathcal{I}\left(\widehat{\Omega }\right)$. The elements of $\mathcal{I}\left(\widehat{\Omega }\right)$ were derived using the second partial derivative of the log-likelihood function for each of the parameters. The observed Fisher’s information matrix $\mathcal{I}\left(\widehat{\Omega }\right)$ is given by

$\mathcal{I}\left(\widehat{\Omega }\right)=-\left[\begin{array}{cccc}{I}_{11}& I_{12}& I_{13}& I_{14}\\ I_{21}& I_{22}& I_{23}& I_{24}\\ I_{31}& I_{32}& I_{33}& I_{34}\\ I_{41}& I_{42}& I_{43}& I_{44}\end{array}\right],$

and the $\mathcal{V}\left(\widehat{\Omega }\right)$ is then given by

$\mathcal{V}\left(\widehat{\Omega }\right)=\mathcal{I}{\left(\widehat{\Omega }\right)}^{-1}=-{\left[\begin{array}{cccc}{I}_{11}& I_{12}& I_{13}& I_{14}\\ I_{21}& I_{22}& I_{23}& I_{24}\\ I_{31}& I_{32}& I_{33}& I_{34}\\ I_{41}& I_{42}& I_{43}& I_{44}\end{array}\right]}^{-1}=\left[\begin{array}{cccc}{V}_{11}& V_{12}& V_{13}& V_{14}\\ V_{21}& V_{22}& V_{23}& V_{24}\\ V_{31}& V_{32}& V_{33}& V_{34}\\ V_{41}& V_{42}& V_{43}& V_{44}\end{array}\right],$

where

$\begin{array}{cccc}\hfill I_{11}& ={\displaystyle \frac{{\partial }^{2}l}{\partial {\theta }_1^2}}{|}_{(\widehat{{\theta }_1},\widehat{{\theta }_2},\widehat{{\theta }_3},\widehat{\lambda })},\hfill & \hfill I_{22}& ={\displaystyle \frac{{\partial }^{2}l}{\partial {\theta }_2^2}}{|}_{(\widehat{{\theta }_1},\widehat{{\theta }_2},\widehat{{\theta }_3},\widehat{\lambda })},\hfill \\ \hfill I_{33}& ={\displaystyle \frac{{\partial }^{2}l}{\partial {\theta }_3^2}}{|}_{(\widehat{{\theta }_1},\widehat{{\theta }_2},\widehat{{\theta }_3},\widehat{\lambda })},\hfill & \hfill I_{44}& ={\displaystyle \frac{{\partial }^{2}l}{\partial {\theta }_4^2}}{|}_{(\widehat{{\theta }_1},\widehat{{\theta }_2},\widehat{{\theta }_3},\widehat{\lambda })},\hfill \\ \hfill I_{12}& =I_{21}={\displaystyle \frac{{\partial }^{2}l}{\partial {\theta }_1\partial {\theta }_2}}{|}_{(\widehat{{\theta }_1},\widehat{{\theta }_2},\widehat{{\theta }_3},\widehat{\lambda })},\hfill & \hfill I_{13}& =I_{31}={\displaystyle \frac{{\partial }^{2}l}{\partial {\theta }_1\partial {\theta }_3}}{|}_{(\widehat{{\theta }_1},\widehat{{\theta }_2},\widehat{{\theta }_3},\widehat{\lambda })},\hfill \\ \hfill I_{14}& =I_{41}={\displaystyle \frac{{\partial }^{2}l}{\partial {\theta }_1\partial {\theta }_4}}{|}_{(\widehat{{\theta }_1},\widehat{{\theta }_2},\widehat{{\theta }_3},\widehat{\lambda })},\hfill & \hfill I_{23}& =I_{32}={\displaystyle \frac{{\partial }^{2}l}{\partial {\theta }_2\partial {\theta }_3}}{|}_{(\widehat{{\theta }_1},\widehat{{\theta }_2},\widehat{{\theta }_3},\widehat{\lambda })},\hfill \\ \hfill I_{24}& =I_{42}={\displaystyle \frac{{\partial }^{2}l}{\partial {\theta }_2\partial {\theta }_4}}{|}_{(\widehat{{\theta }_1},\widehat{{\theta }_2},\widehat{{\theta }_3},\widehat{\lambda })}\mathrm{and}\hfill & \hfill I_{34}& =I_{43}={\displaystyle \frac{{\partial }^{2}l}{\partial {\theta }_3\partial {\theta }_4}}{|}_{(\widehat{{\theta }_1},\widehat{{\theta }_2},\widehat{{\theta }_3},\widehat{\lambda })}.\hfill \end{array}$

Using the diagonal elements of the $\mathcal{V}\left(\widehat{\Omega }\right)$, we can derive the $(1-\xi )100\%$ confidence interval of the estimated parameter vector $\widehat{\Omega }={(\widehat{{\theta }_1},\widehat{{\theta }_2},\widehat{{\theta }_3},\widehat{\lambda })}^T$ in the following form:

$\widehat{{\Omega }_k}\pm z_{{\displaystyle \frac{\xi }{2}}}\sqrt{V_{k,k}};k=1,2,3,4,$

where $V_{k,k}$ is the $k^{th}$ diagonal element of $\mathcal{V}\left(\widehat{\Omega }\right)$ and $z_{{\displaystyle \frac{\xi }{2}}}$ is the upper ${{\displaystyle \frac{\xi }{2}}}^{th}$ percentile of the standard normal distribution.

6. Estimator Performance in Simulation Framework

In Proposition 2, it is clearly shown that the proposed $BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$ is a natural discrete version of the continuous $BFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$. The following algorithm can be used to generate a random sample from the bivariate distribution $BDFWD$ as follows

• Step 1:

  Generate $U_1\sim FWD({\theta }_1,\lambda )$, $U_2\sim FWD({\theta }_2,\lambda )$, and $U_3\sim FWD({\theta }_3,\lambda )$, using the inverse transformation method of cdf.

• Step 2:

  Obtain $X_1=max\{U_1,U_3\}$ and $X_2=max\{U_2,U_3\}$.

• Step 3:

  The pair $(Y_1,Y_2)$ is the desired random sample from BDFWD(${\theta }_1,{\theta }_2,{\theta }_3,\lambda$), where $Y_1=\lfloor X_1\rfloor$ and $Y_2=\lfloor X_2\rfloor$.

• Step 4:

  Continue Steps 1, 2, and 3 for “n” times to get a random sample of size n from BDFWD.

Here to analyze the performance of the MLEs of BDFWD, a simulation study was carried out using the stats4 package in R-programming software version (4.3.1). Random samples from BDFWD for different sizes ($n=20,40,60,80$, and 100), with known values of the parameters ${\theta }_1,{\theta }_2,{\theta }_3$, and $\lambda$, were generated with 1000 replications for the following two cases:

• (i) 

  Case I: ${\theta }_1=0.12,{\theta }_2=0.74,{\theta }_3=0.09,\lambda =2.21$,

• (ii) 

  Case II: ${\theta }_1=0.2,{\theta }_2=0.5,{\theta }_3=0.6,\lambda =3.0$.

The accuracy and precision of the estimation method are evaluated and validated based on two criteria: bias and mean square error (MSE), where

• $Bias\left(\delta \right)={\displaystyle \frac{1}{1000}}{\sum }_{i=1}^{1000}(\widehat{{\delta }_i}-\delta )$

• $MSE\left(\delta \right)={\displaystyle \frac{1}{1000}}{\sum }_{i=1}^{1000}{(\widehat{{\delta }_i}-\delta )}^2$

where $\widehat{\delta }$ is the estimated value of the parameter $\delta$ using the MLE method.

Based on the simulation results from Table 3, it is observed that as “n” increases, the bias and MSE decrease. This shows the consistency and unbiasedness of the MLEs. Moreover, the magnitude of bias always decreases to zero as $n\to \infty$. From these observations, it can be concluded that the estimates of the parameters using the MLE method can be used quite effectively for data analysis purposes.

7. Data Analysis: Goodness-of-Fit and Informed Decision-Making

This section analyzes three real-life bivariate datasets to establish the proposed distribution’s (BDFWD) applicability in real-life practice. First, each dataset was analyzed for fitting to the parent distribution DFWD. Later, the datasets were analyzed for BDFWD. Further, this section compares five existing competing discrete bivariate distributions. Before analyzing these datasets using BDFWD, at first, the fitting of the marginal distribution (DFWD) on $Y_1$, $Y_2$, and min($Y_1,Y_2$) is checked separately. Here, the concept of scaled TTT plots was considered to get an idea about the shape of the empirical hazard function of the marginals $Y_1$ and $Y_2$. According to Aarset [32], the scaled TTT transform will be convex for the decreasing hazard rate, concave for the increasing hazard rate, diagonally straight line for the constant hazard rate, and tilted S-shaped for the bathtub-shaped hazard rate. After analyzing the shape of the hazard function for the marginals, we proceed to calculate the maximum likelihood estimates (MLEs) of the parameters, the negative log-likelihood (−log l), Akaike information criterion (AIC), Bayesian information criterion (BIC), corrected Akaike information criterion (CAIC), Hannan–Quinn information criterion (HQIC), the chi-square $\left({\chi }^2\right)$, and its corresponding p-value. Then the estimated frequency plots of the marginals for $Y_1$, $Y_2$, and min($Y_1,Y_2$) are also discussed for each case.

After analyzing the marginals, now the data analysis for BDFWD is performed. To illustrate the efficiency of the proposed distribution, comparisons with five other competitive distributions, viz., bivariate discrete Rayleigh (BDR), bivariate discrete Weibull (BDW), bivariate discrete generalized Weibull (BDGW), bivariate discrete inverse Weibull (BDIW), and bivariate discrete generalized inverse Weibull (BDGIW), were performed.

7.1. Application I

The first dataset represents a football match score in an Italian football match (Serie A) from 1996 to 2011 between ACF Fiorentina ($Y_1$) and Juventus ($Y_2$) and is reported in Lee and Cha [17]. The data are as reported in Table 4. These data basically consist of 26 pairs of observations ($Y_1$,$Y_2$), where $Y_1$ represents the score obtained by the team ACF Fiorentina and $Y_2$ represents the score secured by the team Juventus, when both the teams played against each other during the period of 1996 to 2011. The dataset I is given in Table 4.

The values of Pearson’s correlation coefficient, Spearman’s rho $\left(\rho \right)$, and Kendall’s tau $\left(\tau \right)$ between $Y_1$ and $Y_2$ for this dataset are 0.815, 0.772, and 0.736 respectively. Clearly, it can be concluded that there is a positive dependence between $Y_1$ and $Y_2$. In Figure 4, the TTT plot of the variables $Y_1$ and $Y_2$ is presented, and it is clearly seen that both have bathtub-shaped empirical hazard functions. Thus, this dataset can be considered to be analyzed with BDFWD.

Now we would like to check whether the BDFWD model fits the data. As it is already proved that for the BDFWD, the marginals will be DFWD; hence, at first, we will be testing the marginals only. The values of MLEs, negative log likelihood (−log l), AIC, BIC, CAIC, HQIC, ${\chi }^2$, and its corresponding p-value for $Y_1$, $Y_2$, and min($Y_1,Y_2$) using DFWD are reported in Table 5. Since the p-values for both $Y_1$ and $Y_2$ are significant, hence we have reason to conclude that the marginals give a good fit to DFWD. Thus, it is reasonable to use the BDFWD model to analyze the dataset I. And in Figure 5, the observed and expected frequency plots of the marginals $Y_1$, $Y_2$ and min($Y_1,Y_2$) using the DFWD for the dataset I are presented.

The efficiency of the BDFWD model in fitting the first dataset is compared against five competing models. Table 6 provides the values of MLEs, $-logl$, AIC, BIC, CAIC, and HQIC for the models BDR, BDW, BDGW, BDIW, BDGIW, and BDFWD using dataset I.

Table 6 shows that dataset I fits best with the proposed BDFWD model, as evidenced by the lowest values of $-logl$, AIC, CAIC, BIC, and HQIC for BDFWD. Additionally, Figure 6 presents the observed and estimated joint frequency distribution plot of dataset I for the BDFWD, BDGIW, BDGW, BDIW, BDW, and BDR distributions. The graphical representation further illustrates that dataset I aligns most closely with the BDFWD model.

7.2. Application II

The second dataset is extracted from Barbiero [33] representing the data resulting from a study of thirty-eight kidney dialysis patients. This dataset was initially considered in the study of McGilchrist and Aisbett [34], concerned with the pervasiveness of infection at the catheter insertion point. The dataset consists of observations on two recurrence times to infection for each patient, viz., the first recurrence time to infection ($Y_1$, in days), and the second recurrence time to infection ($Y_2$, in days). Catheter removal is necessary if an infection develops. However, for any reason other than infection that results in the removal of the catheter, it counts as censoring the first recurrence time. Then, after some predetermined period, the catheter is again reinserted. For the second catheter, if infection occurs again, the second recurrence time to infection ($Y_2$, in days) is measured as the time elapsed between the second insertion and the occurrence of the second infection or the second censoring. The removal of the catheter either for any reason other than infection or the follow-up period for the patient ends before the infection occurs. In that case, the second infection recurrence time counts as censored. The dataset II is provided in Table 7.

The values of Pearson’s correlation coefficient, Spearman’s rho $\left(\rho \right)$, and Kendall’s tau $\left(\tau \right)$ for $Y_1$ and $Y_2$ are all equal to 0.1716, indicating a positive dependence between $Y_1$ and $Y_2$ for dataset II. Furthermore, the TTT plot in Figure 7 shows that both variables $Y_1$ and $Y_2$ exhibit bathtub-shaped empirical hazard functions. Therefore, this dataset can be effectively analyzed using the BDFWD model.

The MLEs of the parameters corresponding to DFWD for $Y_1$, $Y_2$, and min($Y_1,Y_2$), along with the corresponding values of -log l, AIC, BIC, CAIC, HQIC, ${\chi }^2$, and its corresponding p-value for $Y_1$, $Y_2$, and min($Y_1,Y_2$) using DFWD are reported in Table 8. Since the p-values for both $Y_1$ and $Y_2$ are significantly high, we have reason to conclude that the marginals give a good fit to DFWD. Thus, it is reasonable to use the BDFWD model to analyze the dataset II. In Figure 8, the observed and expected frequency plots of the marginals $Y_1$, $Y_2$, and min($Y_1,Y_2$) using DFWD for the dataset II are presented.

In Table 9, the values MLEs, −log l, AIC, BIC, CAIC, and HQIC corresponding to BDR, BDW, BDGW, BDIW, BDGIW, and BDFWD using dataset II is presented.

Table 9 indicates that dataset II fits best with the proposed BDFWD distribution, as evidenced by the lowest values of $-logl$, AIC, CAIC, BIC, and HQIC for BDFWD. Additionally, Figure 9 displays the observed and estimated joint frequency distribution plot of dataset II for the BDFWD, BDGIW, BDGW, BDIW, BDW, and BDR distributions. The graphical representation further confirms that dataset II aligns most closely with the BDFWD model.

7.3. Application III

The third dataset represents the marks twenty-four students received in the second semester of their Master’s degree internal test. These students are from the 2022–2023 batch in the Department of Statistics at North Eastern Hill University. The dataset contains marks scored on two papers, namely, Distribution theory and matrix algebra and Statistical analysis of data for survival analysis. The Department of Statistics, North Eastern Hill University, provided the dataset for this collection. The observations for variables $Y_1$ and $Y_2$ represent the marks the students received on the first and second papers, respectively. The tabulated values of $Y_1$ and $Y_2$ are $0,1,2,3$, and 4. These values are classified as follows: 0 represents marks less than $20\%$, 1 represents marks between $20\%$ and $40\%$, 2 represents marks between $40\%$ and $60\%$, 3 represents marks between $60\%$ and $80\%$, and 4 represents marks more than $80\%$. Table 10 displays the dataset III.

The values of Pearson’s correlation coefficient, Spearman’s rho $\left(\rho \right)$, and Kendall’s tau $\left(\tau \right)$ between $Y_1$ and $Y_2$ for the fifth dataset are 0.940, 0.948 and 0.916, respectively, which clearly suggests a positive dependence between $Y_1$ and $Y_2$. In Figure 10, the TTT plots of the variables $Y_1$ and $Y_2$ are presented, and it is clearly seen that both have bathtub-shaped empirical hazard functions. Thus, this dataset can be considered to be analyzed with BDFWD.

The values of MLEs, $-logl$, AIC, BIC, CAIC, HQIC, ${\chi }^2$ and its corresponding p-value for $Y_1$, $Y_2$, and min($Y_1,Y_2$) using DFWD are reported in Table 11. Since the p-values for both $Y_1$ and $Y_2$ are significantly high, we have reason to conclude that the marginals give a good fit to DFWD. Thus, it is reasonable to use the BDFWD model to analyze the dataset III. And in Figure 11, the observed and expected frequency plots of the marginals $Y_1$, $Y_2$, and min($Y_1,Y_2$) using DFWD for the dataset III are presented.

Table 12 presents the values of MLEs, $-logl$, AIC, BIC, CAIC, and HQIC corresponding to the models BDR, BDW, BDGW, BDIW, BDGIW, and BDFWD using dataset III.

Table 12 shows that dataset III fits best with the proposed BDFWD distribution, as indicated by the lowest values of $-logl$, AIC, CAIC, BIC, and HQIC for BDFWD. Furthermore, Figure 12 presents the observed and estimated joint frequency distribution plot of dataset III for the BDFWD, BDGIW, BDGW, BDIW, BDW, and BDR distributions. The graphical representation also demonstrates that dataset III aligns most closely with the BDFWD model.

8. Conclusions

This paper introduced a novel four-parameter discrete bivariate distribution, termed the bivariate discretized Fréchet–Weibull distribution (BDFWD), with its marginals derived from the discretized Fréchet–Weibull distribution (DFWD). Several statistical properties and reliability characteristics of this proposed discrete model were thoroughly derived. Notably, the joint cumulative distribution function, joint probability mass function, joint survival function, bivariate hazard rate function, and bivariate reversed hazard rate function were expressed in straightforward forms, enhancing the practicality of the BDFWD for modeling bivariate discrete data. Moreover, the BDFWD proved effective for modeling asymmetric data across various kurtosis forms and hazard rate shapes, making it suitable for addressing phenomena such as over-dispersion and zero inflation. The properties of TP2 and PQD were validated for the new model, further underscoring its robustness. Additionally, the stress–strength model was shown to depend solely on the model parameters. The method of MLE was employed to estimate the parameters on simulated data, yielding satisfactory results. The applicability of the BDFWD was demonstrated through analyses of three real-life datasets: football score records, recurrence times to infection for kidney dialysis patients, and student marks from two internal examination statistical papers. Comparisons with five competitive bivariate distributions, including the bivariate discrete Rayleigh, bivariate discrete Weibull, bivariate discrete generalized Weibull, bivariate discrete inverse Weibull, and bivariate discrete generalized inverse Weibull, revealed the superiority of the BDFWD, as supported by goodness-of-fit tests. Ultimately, the study confirmed that the BDFWD outperformed the other five distributions in terms of efficiency across various discrete data applications.

Author Contributions

Conceptualization, D.D.; methodology, K.A.R. and M.S.E.; software, B.D., P.J.H. and M.S.E.; validation, T.S.A. and B.D.; formal analysis, P.J.H. and M.S.E.; investigation, D.D. and T.S.A.; resources, T.S.A., K.A.R. and P.J.H.; data curation, D.D., K.A.R. and B.D.; writing—original draft, D.D. and P.J.H.; writing—review & editing, T.S.A., B.D. and M.S.E.; visualization, D.D.; project administration, M.S.E. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Acknowledgments

The researchers would like to thank the Deanship of Graduate Studies and Scientific Research at Qassim University for financial support (QU-APC-2024-9/1).

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. The scatter plots of the joint pmf of BDFWD for different parameter values [Forumla omitted. See PDF.].

Figure 2. The scatter plots of the joint sf of BDFWD for different parameter values [Forumla omitted. See PDF.].

Figure 3. The scatter plots of the BHRF of BDFWD for different parameter values [Forumla omitted. See PDF.].

Figure 4. The scaled TTT plot of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for dataset I.

Figure 5. The estimated pmf of the marginals [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and min([Forumla omitted. See PDF.]) using dataset I.

Figure 6. The observed and estimated joint frequency distribution plots of dataset I.

Figure 7. The scaled TTT plot of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for dataset II.

Figure 8. The observed and expected frequency of the marginals [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and min([Forumla omitted. See PDF.]) for dataset II.

Figure 9. The observed and estimated joint frequency distribution plots of dataset II.

Figure 10. The scaled TTT plot of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for dataset III.

Figure 11. The observed and expected frequency of the marginals [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and min([Forumla omitted. See PDF.]) for dataset III.

Figure 12. The observed and estimated joint frequency distribution plot of dataset III.

Table 1

Numerical values for both $\mathrm{JSk}(Y_1,Y_2)$ and $\mathrm{JKu}(Y_1,Y_2)$.

Measure	$\mathit{\lambda }$	${\mathit{\theta }}_2=0.25,{\mathit{\theta }}_3=0.30$			${\mathit{\theta }}_1=0.35,{\mathit{\theta }}_3=0.75$
		${\mathit{\theta }}_{\mathbf{1}}$			${\mathit{\theta }}_{\mathbf{2}}$
		$\mathbf{0.15}$	$\mathbf{0.25}$	$\mathbf{0.55}$	$\mathbf{0.15}$	$\mathbf{0.25}$	$\mathbf{0.55}$
	$0.35$	$0.639$	$0.783$	$0.899$	$2.968$	$1.367$	$1.355$
$JSk(Y_1,Y_2)$	$0.55$	$0.784$	$0.772$	$0.793$	$3.057$	$1.591$	$1.429$
	$0.85$	$0.799$	$0.876$	$0.984$	$3.179$	$1.884$	$1.798$
	$0.35$	$1.221$	$1.287$	$1.367$	$2.971$	$2.998$	$3.568$
$JKu(Y_1,Y_2)$	$0.55$	$1.381$	$1.530$	$1.698$	$3.157$	$3.214$	$3.785$
	$0.85$	$1.967$	$2.061$	$2.102$	$3.264$	$3.288$	$3.859$

Table 2

Some descriptive statistics for BDFWD as ${\theta }_1,{\theta }_2,{\theta }_3$ and $\lambda$ increases.

${\mathit{\theta }}_1$	${\mathit{\theta }}_2$	${\mathit{\theta }}_3$	$\mathit{\lambda }$	Cov(${\mathit{Y}}_1,{\mathit{Y}}_2$)	Cor(${\mathit{Y}}_1,{\mathit{Y}}_2$)	$\mathit{\rho }$(${\mathit{Y}}_1,{\mathit{Y}}_2$)	$\mathit{\tau }$(${\mathit{Y}}_1,{\mathit{Y}}_2$)
0.3	0.2	0.6	5.0	0.214141	0.839211	0.833419	0.829402
0.5				0.139394	0.545766	0.531479	0.524273
0.7				0.123232	0.497665	0.475403	0.468685
0.9				0.123232	0.497665	0.475403	0.468685
0.4	0.2	0.5	4.0	0.211313	0.709154	0.675142	0.662239
	0.5			0.271717	0.892231	0.871426	0.871374
	0.7			0.271717	0.892231	0.871426	0.871374
	0.8			0.271717	0.892231	0.871426	0.871374
0.7	0.5	0.5	6.0	0.271111	1.00	1.00	1.00
		0.6		0.224465	0.831719	0.823383	0.823313
		0.8		0.126263	0.557209	0.553591	0.550924
		0.9		0.126263	0.557209	0.553591	0.550924
0.5	0.7	0.8	5.0	0.141010	0.584678	0.558551	0.558353
			7.0	0.118687	0.543689	0.543689	0.543689
			8.0	0.118687	0.543689	0.543689	0.543689
			10.0	0.118687	0.543689	0.543689	0.543689
0.2	0.5	0.7	3.0	0.410101	0.779594	0.670892	0.652277
			5.0	0.154950	0.622983	0.597717	0.597535
			7.0	0.127878	0.579494	0.582258	0.579497
			8.0	0.123131	0.580374	0.580374	0.580374
			9.0	0.123131	0.580374	0.580374	0.580374

Table 3

The average estimates, biases and MSEs for case I and case II.

		Case I:				Case II:
		(${\mathit{\theta }}_{\mathbf{1}},{\mathit{\theta }}_{\mathbf{2}},{\mathit{\theta }}_{\mathbf{3}},\mathit{\lambda }$) =				(${\mathit{\theta }}_{\mathbf{1}},{\mathit{\theta }}_{\mathbf{2}},{\mathit{\theta }}_{\mathbf{3}},\mathit{\lambda }$) =
		(0.12 , 0.74 , 0.09 , 2.21)				(0.2 , 0.5 , 0.6 , 3.0)
	n	$\widehat{{\mathit{\theta }}_{\mathbf{1}}}$	$\widehat{{\mathit{\theta }}_{\mathbf{2}}}$	$\widehat{{\mathit{\theta }}_{\mathbf{3}}}$	$\widehat{\mathit{\lambda }}$	$\widehat{{\mathit{\theta }}_{\mathbf{1}}}$	$\widehat{{\mathit{\theta }}_{\mathbf{2}}}$	$\widehat{{\mathit{\theta }}_{\mathbf{3}}}$	$\widehat{\mathit{\lambda }}$
Avg Est	20	0.0412	0.7787	0.2701	1.9745	0.3500	0.6221	0.4013	4.2101
	40	0.728	0.7101	0.1090	1.9905	0.2810	0.6009	0.4522	3.7722
	60	0.0851	0.7221	0.0551	2.044	0.2311	0.5742	0.5071	3.4432
	80	0.0905	0.7521	0.0808	2.1981	0.2051	0.5524	0.5743	3.1014
	100	0.1082	0.7413	0.0821	2.2212	0.2020	0.5005	0.5909	3.0015
Bias	20	−0.0788	0.0387	0.1801	−0.2355	0.1500	0.1221	−0.1987	1.2101
	40	−0.0472	−0.0299	0.01907	−0.2195	0.0810	0.1009	−0.1478	0.7722
	60	−0.0349	−0.0179	−0.0349	−0.1660	0.0311	0.0742	−0.0929	0.4432
	80	−0.0295	0.0121	−0.0092	−0.0119	0.0051	0.0524	−0.0257	0.1014
	100	−0.0098	0.0013	0.0029	0.0012	0.0020	0.0005	−0.0091	0.0015
MSE	20	0.5214	0.2887	0.6122	0.2041	0.4311	0.5951	0.2940	0.5424
	40	0.3551	0.1637	0.5162	0.1041	0.3520	0.4044	0.1952	0.3217
	60	0.2601	0.1083	0.3417	0.0247	0.2914	0.2315	0.1547	0.1128
	80	0.1434	0.1005	0.2107	0.0105	0.1059	0.1009	0.1220	0.0719
	100	0.1008	0.0097	0.1011	0.0097	0.0981	0.0077	0.1021	0.0052

Table 4

The football score data between ACF Fiorentina and Juventus.

Obs.	Match Date	${\mathit{Y}}_1$	${\mathit{Y}}_2$	Obs.	Match Date	${\mathit{Y}}_1$	${\mathit{Y}}_2$	Obs.	Match Date	${\mathit{Y}}_1$	${\mathit{Y}}_2$
1	25 October 2011	1	2	10	9 April 2006	1	1	19	18 December 1999	1	1
2	17 April 2011	0	0	11	4 December 2005	1	2	20	24 April 1999	1	2
3	27 November 2010	1	1	12	9 April 2005	3	3	21	12 December 1998	1	0
4	6 March 2010	1	2	13	10 November 2004	0	1	22	21 February 1998	3	0
5	17 October 2009	1	1	14	16 February 2002	1	2	23	4 October 1997	1	2
6	24 January 2009	0	1	15	19 December 2001	1	1	24	22 February 1997	1	1
7	31 August 2008	1	1	16	12 May 2001	1	3	25	28 September 1996	0	1
8	2 March 2008	3	2	17	6 January 2001	3	3	26	23 March 1996	0	1
9	7 October 2007	1	1	18	21 April 2000	0	1

Table 5

MLEs, AIC, BIC, CAIC, HQIC, ${\chi }^2$, and p-value of the marginals using DFWD for dataset I.

Marginals	Parameter Estimates		−log l	AIC	BIC	CAIC	HQIC	${\mathit{\chi }}^2$	p-Value
$Y_1$	$\widehat{{\theta }_1}=$ 0.238	$\widehat{\lambda }=$ 2.799	30.87	65.73	68.25	66.99	66.45	3.62	0.057
$Y_2$	$\widehat{{\theta }_2}=$ 0.096	$\widehat{\lambda }=$ 2.601	33.73	71.46	73.97	72.72	72.18	1.49	0.6854
min($Y_1,Y_2$)	$\widehat{{\theta }_3}=$ 0.234	$\widehat{\lambda }=$ 2.941	29.59	63.17	65.69	64.43	63.89	3.09	0.2138

Table 6

The MLEs, AIC, BIC, CAIC, and HQIC for dataset I.

Model	Parameter Estimates				−log l	AIC	BIC	CAIC	HQIC
BDR	$\widehat{p_1}$ = 0.623	$\widehat{p_2}$ = 0.817	$\widehat{p_3}$ = 0.926		74.05	154.1	157.87	155.99	155.18
BDW	$\widehat{p_1}$ = 0.759	$\widehat{p_2}$ = 0.819	$\widehat{p_3}$ = 0.943	$\widehat{\alpha }$ = 2.129	73.35	154.7	159.74	157.22	156.15
BDGW	$\widehat{q_1}$ = 0.752	$\widehat{q_2}$ = 0.788	$\widehat{q_3}$ = 0.992	$\widehat{\alpha }$ = 1.97	66.17	124.34	119.31	121.83	122.89
BDIW	$\widehat{{\alpha }_1}$ = 0.422	$\widehat{{\alpha }_2}$ = 0.047	$\widehat{{\alpha }_3}$ = 0.599	$\widehat{\beta }$ = 2.701	60.19	128.38	133.42	130.9	129.83
BDGIW	$\widehat{q_1}$ = 0.822	$\widehat{q_2}$ = 0.337	$\widehat{q_3}$ = 0.267	$\widehat{\alpha }$ = 2.608	58.53	125.06	130.1	127.58	126.51
BDFWD	$\widehat{{\theta }_1}$ = 0.822	$\widehat{{\theta }_2}$ = 0.327	$\widehat{{\theta }_3}$ = 0.257	$\widehat{\lambda }$ = 2.608	56.07	120.15	125.18	122.67	121.6

Table 7

The recurrence time to infection of 38 kidney dialysis patients.

Patient No.	${\mathit{Y}}_1$	${\mathit{Y}}_2$	Patient No.	${\mathit{Y}}_1$	${\mathit{Y}}_2$	Patient No.	${\mathit{Y}}_1$	${\mathit{Y}}_2$	Patient No.	${\mathit{Y}}_1$	${\mathit{Y}}_2$
1	1	1	11	1	1	21	1	1	31	1	1
2	1	0	12	1	0	22	1	0	32	0	1
3	1	1	13	1	1	23	1	1	33	1	1
4	1	1	14	0	0	24	1	0	34	1	0
5	1	1	15	1	0	25	1	1	35	1	1
6	1	1	16	1	0	26	0	1	36	0	0
7	1	1	17	1	1	27	1	1	37	0	1
8	1	1	18	1	1	28	1	1	38	1	0
9	1	1	19	0	0	29	1	1
10	1	1	20	1	0	30	1	1

Table 8

MLEs, AIC, BIC, CAIC, HQIC, ${\chi }^2$, and p-value of the marginals using DFWD for dataset II.

Marginals	Parameter Estimates		−log l	AIC	BIC	CAIC	HQIC	${\mathit{\chi }}^2$	p-Value
$Y_1$	$\widehat{{\theta }_1}$ = 0.158	$\widehat{\lambda }$ = 22.86	16.57	37.15	40.42	38.79	38.31	3.2 × ${10}^{-6}$	0.9986
$Y_2$	$\widehat{{\theta }_2}$ = 0.316	$\widehat{\lambda }$ = 24.934	23.70	51.40	54.67	53.04	52.56	1.4 × ${10}^{-6}$	0.9991
min($Y_1,Y_2$)	$\widehat{{\theta }_3}$ = 0.395	$\widehat{\lambda }$ = 23.656	25.49	54.98	58.26	56.62	56.15	6 × ${10}^{-7}$	0.9994

Table 9

The MLEs, AIC, BIC, CAIC and HQIC for dataset II.

Model	Parameter Estimates				−log l	AIC	BIC	CAIC	HQIC
BDR	$\widehat{p_1}$ = 0.585	$\widehat{p_2}$ = 0.524	$\widehat{p_3}$ = 0.996		59.46	124.92	129.83	127.38	126.67
BDW	$\widehat{p_1}$ = 0.642	$\widehat{p_2}$ = 0.814	$\widehat{p_3}$ = 0.951	$\widehat{\alpha }$ = 4.308	49.49	106.98	113.53	110.26	109.31
BDGW	$\widehat{q_1}$ = 0.654	$\widehat{q_2}$ = 0.694	$\widehat{q_3}$ = 0.992	$\widehat{\alpha }$ = 5.721	43.97	95.94	102.49	99.22	98.27
BDIW	$\widehat{{\alpha }_1}$ = 0.601	$\widehat{{\alpha }_2}$ = 0.424	$\widehat{{\alpha }_3}$ = 0.579	$\widehat{\beta }$ = 8.305	41.02	90.04	96.59	93.32	92.37
BDGIW	$\widehat{q_1}$ = 0.154	$\widehat{q_2}$ = 0.502	$\widehat{q_3}$ = 0.632	$\widehat{\alpha }$ = 16.559	37.84	83.68	90.23	86.96	86.01
BDFWD	$\widehat{{\theta }_1}$ = 0.35	$\widehat{{\theta }_2}$ = 0.69	$\widehat{{\theta }_3}$ = 0.532	$\widehat{\lambda }$ = 23.655	35.29	78.58	85.13	81.86	80.91

Table 10

The marks scored by 24-semester students in two papers.

Student No.	1	2	3	4	5	6	7	8	9	10	11	12
$Y_1$	3	2	1	0	3	4	1	3	3	0	2	3
$Y_2$	3	1	2	3	3	3	2	2	2	0	1	2
Student No.	13	14	15	16	17	18	19	20	21	22	23	24
$Y_1$	1	3	1	2	2	4	2	2	2	1	2	0
$Y_2$	2	2	2	1	1	4	2	3	2	0	1	1

Table 11

MLEs, AIC, BIC, CAIC, HQIC, ${\chi }^2$, and p-value of the marginals using DFWD for dataset III.

Marginals	Parameter Estimates		−log l	AIC	BIC	CAIC	HQIC	${\mathit{\chi }}^2$	p-Value
$Y_1$	$\widehat{{\theta }_1}=$ 0.071	$\widehat{\lambda }=$ 1.876	43.27	90.53	92.89	91.71	91.16	1.89	0.3874
$Y_2$	$\widehat{{\theta }_2}=$ 0.037	$\widehat{\lambda }=$ 2.183	38.90	81.81	84.16	82.98	82.43	0.14	0.712
min($Y_1,Y_2$)	$\widehat{{\theta }_3}=$ 0.137	$\widehat{\lambda }=$ 2.132	42.24	88.49	90.84	89.66	89.11	1.43	0.2325

Table 12

The MLEs, AIC, BIC, CAIC, and HQIC for dataset III.

Model	Parameter Estimates				−log l	AIC	BIC	CAIC	HQIC
BDR	$\widehat{p_1}$ = 0.842	$\widehat{p_2}$ = 0.754	$\widehat{p_3}$ = 0.998		78.57	163.14	166.67	164.91	164.08
BDW	$\widehat{p_1}$ = 0.791	$\widehat{p_2}$ = 0.84	$\widehat{p_3}$ = 0.995	$\widehat{\alpha }$ = 1.845	76.33	160.66	165.37	163.02	161.91
BDGW	$\widehat{q_1}$ = 0.852	$\widehat{q_2}$ = 0.835	$\widehat{q_3}$ = 0.999	$\widehat{\alpha }$ = 1.99	73.24	154.48	159.19	156.84	155.73
BDIW	$\widehat{{\alpha }_1}$ = 0.473	$\widehat{{\alpha }_2}$ = 0.501	$\widehat{{\alpha }_3}$ = 0.168	$\widehat{\beta }$ = 1.247	64.84	137.68	142.39	140.04	138.93
BDGIW	$\widehat{q_1}$ = 0.314	$\widehat{q_2}$ = 0.404	$\widehat{q_3}$ = 0.314	$\widehat{\alpha }$ = 1.876	61.37	130.74	135.45	133.1	131.99
BDFWD	$\widehat{{\theta }_1}$ = 0.585	$\widehat{{\theta }_2}$ = 0.786	$\widehat{{\theta }_3}$ = 0.074	$\widehat{\lambda }$ = 1.962	58.48	124.96	129.67	127.32	126.21

Appendix A

Proof to Proposition 2.

Given, $(X_1,X_2)\sim BFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$. Then, the marginals are $$X_1\sim FWD({\theta }_1{\theta }_3,\lambda )\mathrm{and}{X}_2\sim FWD({\theta }_2{\theta }_3,\lambda ).$$

Also, given that $Y_1=\lfloor X_1\rfloor$ and $Y_2=\lfloor X_2\rfloor$, where $\lfloor .\rfloor$ is the floor function.

Now by the definition of $DFWD(\theta ,\lambda )$, we have $$Y_1\sim DFWD({\theta }_1{\theta }_3,\lambda )\mathrm{and}{Y}_2\sim DFWD({\theta }_2{\theta }_3,\lambda ).$$

Thus, the pair $(Y_1,Y_2)\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$. Hence the result. □

Proof to Proposition 3.

If $(Y_1,Y_2)\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$, then by definition, we must have $$Y_1=max\{U_1,U_3\}\mathrm{and}{Y}_2=max\{U_2,U_3\},$$ where $U_1\sim DFWD({\theta }_1,\lambda )$, $U_2\sim DFWD({\theta }_2,\lambda )$ and $U_3\sim DFWD({\theta }_3,\lambda )$ and they are independently distributed. $$\begin{array}{ccc}\hfill P[Z\le z]& =& P[max\{Y_1,Y_2\}\le z]\hfill \\ & =& P[max\{U_1,U_2,U_3\}\le z]\hfill \\ & =& P[U_1\le z,U_2\le z,U_3\le z,]\hfill \\ & =& P[U_1\le z].P[U_2\le z].P[U_3\le z]\hfill \\ & =& {\theta }_1^{{(y+1)}^{-\lambda }}.{\theta }_2^{{(y+1)}^{-\lambda }}.{\theta }_3^{{(y+1)}^{-\lambda }}\hfill \\ & =& {\theta }^{{(y+1)}^{-\lambda }},\hfill \end{array}$$ which is the form of cdf of $DFWD(\theta ,\lambda )$, where $\theta ={\theta }_1{\theta }_2{\theta }_3$. Hence the result. □

Proof to Proposition 4.

Given, $(Y_{i1},Y_{i2})\sim BDFWD({\theta }_{i1},{\theta }_{i2},{\theta }_{i3},\lambda )$ for $i=1,2,\dots ,n$, and they are independently distributed. Then by definition, for each of $Y_{ij}(i=1,2,\dots ,n;j=1,2)$ we must have $$U_{i1}\sim DFWD({\theta }_{i1},\lambda ),U_{i2}\sim DFWD({\theta }_{i2},\lambda )\mathrm{and}{U}_{i3}\sim DFWD({\theta }_{i3},\lambda )$$ such that $Y_{i1}=max\{U_{i1},U_{i3}\}$ and $Y_{i2}=max\{U_{i2},U_{i3}\}$; $\forall i=1,2,\dots ,n$.

Now, let us define, $U_1=max\left\{U_{i1}\right\}=max\{U_{11},U_{21},U_{31},\dots ,U_{n1}\}$, then $$\begin{array}{ccc}\hfill P[U_1\le u]& =& P[max\{U_{11},U_{21},U_{31},\dots ,U_{n1}\}\le u]\hfill \\ & =& \prod _{i=1}^n{\theta }_{i1}^{{(u+1)}^{-\lambda }}={\left[\prod _{i=1}^n{\theta }_{i1}\right]}^{{(u+1)}^{-\lambda }},\hfill \end{array}$$ which is the form of cdf of DFWD. Thus, $U_1\sim DFWD\left(\left[{\prod }_{i=1}^n{\theta }_{i1}\right],\lambda \right)$. Similarly, $U_2=max\left\{U_{i2}\right\}\sim DFWD\left(\left[{\prod }_{i=1}^n{\theta }_{i2}\right],\lambda \right)$ and $U_3=max\left\{U_{i3}\right\}\sim DFWD\left(\left[{\prod }_{i=1}^n{\theta }_{i3}\right],\lambda \right);\forall i=1,2,\dots ,n.$

Now, we are given that $$\begin{array}{ccc}\hfill Z_1& =& max\{Y_{11},Y_{21},Y_{31},\dots ,Y_{n1}\}\hfill \\ & =& max\{(U_{11},U_{13}),(U_{21},U_{23}),(U_{31},U_{33}),\dots ,(U_{n1},U_{n3})\}\hfill \\ & =& max\{U_1,U_3\}.\hfill \end{array}$$

Similarly, $Z_2=max\{Y_{12},Y_{22},Y_{32},\dots ,Y_{n2}\}=max\{U_2,U_3\}.$

Thus, by the definition of BDFWD, we can write $$(Z_1,Z_2)\sim BDFWD\left(\prod _{i=1}^n{\theta }_{i1},\prod _{i=1}^n{\theta }_{i2},\prod _{i=1}^n{\theta }_{i3},\lambda \right).$$

Hence the result. □

Proof to Proposition 5.

If $(Y_1,Y_2)\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$, then by the Proposition 1, the cdf of the marginals $Y_1$ and $Y_2$ are given by $$F_{Y_i}\left(y\right)=F_{DFWD}(y;{\theta }_i{\theta }_3,\lambda )={\left({\theta }_i{\theta }_3\right)}^{{(y+1)}^{-\lambda }};i=1,2.$$

Thus the stress-strength reliability is given by $$\begin{array}{ccc}\hfill P[Y_1<Y_2]& =& \sum _{y=0}^{\infty }P[Y_2=y+1].P[Y_1\le y]\hfill \\ & =& \sum _{y=0}^{\infty }{\left({\theta }_1{\theta }_3\right)}^{{(y+1)}^{-\lambda }}\left[{\left({\theta }_2{\theta }_3\right)}^{{(y+2)}^{-\lambda }}-{\left({\theta }_2{\theta }_3\right)}^{{(y+1)}^{-\lambda }}\right].\hfill \end{array}$$

Hence the result. □

Proof to Proposition 6.

Let us recall that the function $\mathsf{{\rm Y}}(p,q):\mathcal{R}\times \mathcal{R}\to \mathcal{R}$ is said to have the total positivity of order two $\left(T{P}_2\right)$ property if $\mathsf{{\rm Y}}(p,q)$ satisfies the condition $$\mathsf{{\rm Y}}(p_1,q_1)\times \mathsf{{\rm Y}}(p_2,q_2)\ge \mathsf{{\rm Y}}(p_2,q_1)\times \mathsf{{\rm Y}}(p_1,q_2),$$ $\forall p_1,q_1,p_2,q_2\in \mathcal{R}$ (see [18]). Let us assume that $y_{11},y_{21},y_{12},y_{22}\in {\mathbf{Z}}_{+}$ are from $BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$, such that $y_{11}<y_{21}<y_{12}<y_{22}$. Also, let us denote the cdf of $(Y_1,Y_2)\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$ for $(y_{11},y_{21})$ as $$F_{BDFWD}(y_{11},y_{21};\Psi )=F_{Y_1,Y_2}(y_{11},y_{21};{\theta }_1,{\theta }_2,{\theta }_3,\lambda ).$$

Then, it can be shown that the joint cdf of $(Y_1,Y_2)\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$ satisfies the $\left(T{P}_2\right)$ property. $$\begin{array}{ccc}\hfill {\displaystyle \frac{F_{BDFWD}(y_{11},y_{21};\Psi )F_{BDFWD}(y_{12},y_{22};\Psi )}{F_{BDFWD}(y_{12},y_{21};\Psi )F_{BDFWD}(y_{11},y_{22};\Psi )}}& =& {\displaystyle \frac{F_{DFWD}(y_{21};{\theta }_2,\lambda )F_{DFWD}(y_{12};{\theta }_1{\theta }_3,\lambda )}{F_{DFWD}(y_{12};{\theta }_1,\lambda )F_{DFWD}(y_{21};{\theta }_2{\theta }_3,\lambda )}}\hfill \\ & =& {\theta }_3^{{(y_{12}+1)}^{-\lambda }-{(y_{21}+1)}^{-\lambda }}\ge 1.\hfill \end{array}$$

Thus, for $y_{11}<y_{21}<y_{12}<y_{22}$, we have $$F_{BDFWD}(y_{11},y_{21};\Psi )\times F_{BDFWD}(y_{12},y_{22};\Psi )\ge F_{BDFWD}(y_{12},y_{21};\Psi )\times F_{BDFWD}(y_{11},y_{22};\Psi ).$$

Similarly, considering all other cases such as $y_{11}=y_{21}<y_{12}<y_{22},y_{21}<y_{11}<y_{12}<y_{22}$, etc., it can be shown that the joint cdf of BDFWD satisfies the $\left(T{P}_2\right)$ property.

Hence the result. □

Proof to Proposition 7.

To prove that $Y_1$ and $Y_2$ are positive quadrant dependent $\left(PQD\right)$, we need to show that for $y_1$ and $y_2$ $$P(Y_1\le y_1,Y_2\le y_2)\ge P(Y_1\le y_1).P(Y_2\le y_2).$$

Case 1: Let us consider $y_1\ge y_2$, then $$\begin{array}{ccc}\hfill P(Y_1\le y_1,Y_2\le y_2)& =& F_{DFWD}(y_1;{\theta }_1,\lambda ).F_{DFWD}(y_2;{\theta }_2,\lambda ).F_{DFWD}(y_3;{\theta }_3,\lambda )\hfill \\ & =& {\theta }_1^{{(y_1+1)}^{-\lambda }}.{\theta }_2^{{(y_2+1)}^{-\lambda }}.{\theta }_3^{{(z+1)}^{-\lambda }};z=min(y_1,y_2),\hfill \end{array}$$ and $$\begin{array}{ccc}\hfill P(Y_1\le y_1).P(Y_2\le y_2)& =& F_{DFWD}(y_1;{\theta }_1{\theta }_3,\lambda ).F_{DFWD}(y_2;{\theta }_2{\theta }_3,\lambda )\hfill \\ & =& {\left({\theta }_1{\theta }_3\right)}^{{(y_1+1)}^{-\lambda }}.{\left({\theta }_2{\theta }_3\right)}^{{(y_2+1)}^{-\lambda }}.\hfill \end{array}$$

Clearly, it can be seen that $$P(Y_1\le y_1,Y_2\le y_2)\ge P(Y_1\le y_1).P(Y_2\le y_2);\forall y_1\ge y_2.$$

Case 2: Let us consider $y_1<y_2$, then $$\begin{array}{ccc}\hfill P(Y_1\le y_1,Y_2\le y_2)& =& F_{DFWD}(y_1;{\theta }_1{\theta }_3,\lambda ).F_{DFWD}(y_2;{\theta }_2,\lambda )\hfill \\ & =& F_{DFWD}(y_1;{\theta }_1,\lambda ).F_{DFWD}(y_1;{\theta }_3,\lambda ).F_{DFWD}(y_2;{\theta }_2,\lambda )\hfill \\ & =& {\theta }_1^{{(y_1+1)}^{-\lambda }}.{\theta }_3^{{(y_1+1)}^{-\lambda }}.{\theta }_2^{{(y_2+1)}^{-\lambda }}.\hfill \end{array}$$

For this case also it can be clearly seen that $$P(Y_1\le y_1,Y_2\le y_2)\ge P(Y_1\le y_1).P(Y_2\le y_2);\forall y_1<y_2.$$

Hence the result. □

Proof to Proposition 8.

Given, $(Y_1,Y_2)\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$, then the marginals $Y_i\sim DFWD({\theta }_i{\theta }_3,\lambda );(i=1,2)$. Now, using the inverse transformation method of cdf of DFWD, the quantile function for $Y_i$’s is given by $$q_u=-{\displaystyle \frac{1}{\lambda }}\left[{\displaystyle \frac{log\left(u\right)}{log\left({\theta }_i{\theta }_3\right)}}\right]-1;i=1,2,$$ where u has $Uniform(0,1)$ distribution. Thus, the median of $Y_1$ and $Y_2$ can be obtained by putting $u=0.5$ in the quantile function of DFWD, we get $$M_{Y_i}={\displaystyle \frac{1}{\lambda }}\left[{\displaystyle \frac{log\left(2\right)}{log\left({\theta }_i{\theta }_3\right)}}\right]-1;i=1,2.$$

Hence the result. □

Proof to Proposition 9.

Given $U_i={\lfloor X_i\rfloor }^{-{\displaystyle \frac{1}{\lambda }}}-1$, where $X_i\sim Exponential\left(p_i\right)(i=1,2,3)$ and the parameter $p_i=log\left({\displaystyle \frac{1}{{\theta }_i}}\right)$ such that $0<p_i<1;\forall {\theta }_i>0$. Let us denote the cdf of $U_i$ by $F_{U_i}\left(u\right)$, then for $u\in {\mathbf{Z}}_{+}$ we have $$\begin{array}{ccc}\hfill F_{U_i}\left(u\right)& =& P[U_i\le u]\hfill \\ & =& P[{\lfloor X_i\rfloor }^{-{\displaystyle \frac{1}{\lambda }}}-1\le u]\hfill \\ & =& P[\lfloor X_i\rfloor \ge {(u+1)}^{-\lambda }]\hfill \\ & =& exp\{-{p_i}^{{(\lfloor u\rfloor +1)}^{-\lambda }}\}\hfill \\ & =& {{\theta }_i}^{{(u+1)}^{-\lambda }};u\in {\mathbf{Z}}_{+}.\hfill \end{array}$$

Since $X_i$’s are independent, then $U_i$’s are also independently distributed. Let us consider $Y_1=max\{U_1,U_3\}$ and $Y_2=max\{U_2,U_3\}$, then the joint cdf of $(Y_1,Y_2)$ for $(y_1,y_2)\in {\mathbf{Z}}_{+}$, is given by $$\begin{array}{ccc}\hfill F_{Y_1,Y_2}(y_1,y_2)& =& P[Y_1\le y_1,Y_2\le y_2]\hfill \\ & =& P[max\{U_1,U_3\}\le y_1,max\{U_2,U_3\}\le y_2]\hfill \\ & =& P[U_1\le y_1,U_2\le y_2,U_3\le z]\hfill \\ & =& {\theta }_1^{{(y_1+1)}^{-\lambda }}.{\theta }_2^{{(y_2+1)}^{-\lambda }}.{\theta }_3^{{(z+1)}^{-\lambda }},z=min(y_1,y_2),\hfill \end{array}$$ which is the joint cdf of $BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$. Thus, $(Y_1,Y_2)=\left(max\{U_1,U_3\},max\{U_2,U_3\}\right)\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$.

Hence the result. □

Proof to Proposition 10.

Given $U_i=\left(X_i^{-{\displaystyle \frac{1}{\lambda }}}\right)-1$ where $X_i\sim Geometric(1-{\theta }_i);0<{\theta }_i<1,\lambda >0;i=1,2,3$.

Let us denote the cdf of $U_i$ by $F_{U_i}\left(u\right)$, then for $u\in {\mathbf{Z}}_{+}$ we have $$\begin{array}{ccc}\hfill F_{U_i}\left(u\right)& =& P[U_i\le u]\hfill \\ & =& P[X_i^{-{\displaystyle \frac{1}{\lambda }}}\le u]\hfill \\ & =& P[X_i\ge {(u+1)}^{-\lambda }]\hfill \\ & =& {{\theta }_i}^{{(u+1)}^{-\lambda }};(i=1,2,3).\hfill \end{array}$$

Since $X_i$’s are independent, then $U_i$’s are also independently distributed. Let us consider $Y_1=max\{U_1,U_3\}$ and $Y_2=max\{U_2,U_3\}$, then the joint cdf of $(Y_1,Y_2)$ for $(y_1,y_2)\in {\mathbf{Z}}_{+}$, is given by $$\begin{array}{ccc}\hfill F_{Y_1,Y_2}(y_1,y_2)& =& P[Y_1\le y_1,Y_2\le y_2]\hfill \\ & =& P[max\{U_1,U_3\}\le y_1,max\{U_2,U_3\}\le y_2]\hfill \\ & =& P[U_1\le y_1,U_2\le y_2,U_3\le z]\hfill \\ & =& {\theta }_1^{{(y_1+1)}^{-\lambda }}.{\theta }_2^{{(y_2+1)}^{-\lambda }}.{\theta }_3^{{(z+1)}^{-\lambda }},z=min(y_1,y_2),\hfill \end{array}$$ which is the joint cdf of $BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$. Thus, $(Y_1,Y_2)=\left(max\{U_1,U_3\},max\{U_2,U_3\}\right)$$\sim BDFWD({\theta }_1,{\theta }_2,{\theta }_3,\lambda )$. Hence the result. □

References

1. Jose, K.K.; Ristić, M.M.; Joseph, A. Marshall–Olkin bivariate Weibull distributions and processes. Stat. Pap.; 2011; 52, pp. 789-798. [DOI: https://dx.doi.org/10.1007/s00362-009-0287-8]

2. Kundu, D.; Gupta, R.D. Bivariate generalized exponential distribution. J. Multivar. Anal.; 2009; 100, pp. 581-593. [DOI: https://dx.doi.org/10.1016/j.jmva.2008.06.012]

3. Sarhan, A.; Hamilton, D.C.; Smith, B.; Kundu, D. The bivariate generalized linear failure rate distribution and its multivariate extension. Comput. Stat. Data Anal.; 2011; 55, pp. 644-654. [DOI: https://dx.doi.org/10.1016/j.csda.2010.06.006]

4. El-Sherpieny, E.A.; Ibrahim, S.A.; Bedar, R.E. A New Bivariate Distribution with Generalized Gompertz Marginals. Asian J. Appl. Sci.; 2013; 1, pp. 141-150.

5. Wagner, B.S.; Artur, J.L. Bivariate Kumaraswamy distribution: Properties and a new method to generate bivariate classes. Statistics; 2013; 47, pp. 1321-1342.

6. Balakrishnan, N.; Shiji, K. On a class of bivariate exponential distributions. Stat. Probab. Lett.; 2014; 85, pp. 153-160. [DOI: https://dx.doi.org/10.1016/j.spl.2013.11.009]

7. El-Bassiouny, A.H.; EL-Damcese, M.; Abdelfattah, M.; Eliwa, M.S. Bivariate exponentaited generalized Weibull–Gompertz distribution. J. Appl. Probab. Stat.; 2016; 11, pp. 25-46.

8. El-Gohary, A.; El-Bassiouny, A.H.; El-Morshedy, M. Bivariate exponentiated modified Weibull extension distribution. J. Appl. Probab. Stat.; 2016; 5, pp. 67-78. [DOI: https://dx.doi.org/10.18576/jsap/050106]

9. Hiba, Z.M. Bivariate inverse Weibull distribution. J. Stat. Comput. Simul.; 2016; 86, pp. 2335-2345.

10. Rasool, R.; Akbar, A.J. On bivariate exponentiated extended Weibull family of distributions. Ciência Nat. Santa Maria; 2016; 38, pp. 564-576.

11. Eliwa, M.S.; El-Morshedy, M. Bivariate Gumbel-G family of distributions: Statistical properties, Bayesian and non-Bayesian estimation with application. Ann. Data Sci.; 2019; 6, pp. 39-60. [DOI: https://dx.doi.org/10.1007/s40745-018-00190-4]

12. Eliwa, M.S.; El-Morshedy, M. Bivariate odd Weibull-G family of distributions: Properties, Bayesian and non-Bayesian estimation with bootstrap confidence intervals and application. J. Taibah Univ. Sci.; 2020; 14, pp. 331-345. [DOI: https://dx.doi.org/10.1080/16583655.2020.1741919]

13. Kocherlakota, S.; Kocherlakota, K. Bivariate Discrete Distributions; Marcel Dekker: New York, NY, USA, 1992.

14. Basu, A.P.; Dhar, S.K. Bivariate geometric distribution. J. Appl. Stat. Sci.; 1995; 2, pp. 33-44.

15. Kumar, C.S. A unified approach to bivariate discrete distributions. Metrika; 2008; 67, pp. 113-123. [DOI: https://dx.doi.org/10.1007/s00184-007-0125-8]

16. Kemp, A.W. New discrete Appell and Humbert distributions with relevance to bivariate accident data. J. Multivar. Anal.; 2013; 113, pp. 2-6. [DOI: https://dx.doi.org/10.1016/j.jmva.2011.08.011]

17. Lee, H.; Cha, J.H. On Two General Classes of Discrete Bivariate Distributions. Am. Stat.; 2015; 69, pp. 221-230. [DOI: https://dx.doi.org/10.1080/00031305.2015.1044564]

18. Nekoukhou, V.; Kundu, D. Bivariate discrete generalized exponential distribution. Statistics; 2017; 51, pp. 1143-1158. [DOI: https://dx.doi.org/10.1080/02331888.2017.1289534]

19. Eliwa, M.S.; El-Morshedy, M. Bayesian and non-Bayesian estimation of four-parameter of bivariate discrete inverse Weibull distribution with applications to model failure times, football and biological data. Filomat; 2020; 34, pp. 2511-2531. [DOI: https://dx.doi.org/10.2298/FIL2008511E]

20. El-Morshedy, M.; Eliwa, M.S.; El-Gohary, A.; Khalil, A.A. Bivariate exponentiated discrete Weibull distribution: Statistical properties, estimation, simulation and applications. Math. Sci.; 2020; 14, pp. 29-42. [DOI: https://dx.doi.org/10.1007/s40096-019-00313-9]

21. Ali, S.; Shafqat, M.; Shah, I.; Dey, S. Bivariate discrete Nadarajah and Haghighi distribution: Properties and different methods of estimation. Filomat; 2019; 33, pp. 5589-5610. [DOI: https://dx.doi.org/10.2298/FIL1917589A]

22. Kundu, D.; Nekoukhou, V. On bivariate discrete Weibull distribution. Commun. Stat.-Methods; 2019; 48, pp. 3464-3481. [DOI: https://dx.doi.org/10.1080/03610926.2018.1476712]

23. Shibu, D.S.; Beegum, N. Bivariate Discrete Modified Weibull Distribution. Statistica; 2021; 81, pp. 231-264.

24. Johnson, N.L.; Kotz, S.; Balakrishnan, N. Discrete Multivariate Distributions; John Wiley & Sons: New York, NY, USA, 1997.

25. Balakrishnan, N.; Lai, C. Continuous Bivariate Distributions; 2nd ed. Springer: New York, NY, USA, 2009.

26. Chakraborty, S. Generating Discrete Analogues of Continuous Probability Distributions—A Survey of Methods and Constructions. J. Stat. Distrib. Appl.; 2015; 2, pp. 1-35. [DOI: https://dx.doi.org/10.1186/s40488-015-0028-6]

27. Das, D.; Das, B. Discretized Fréchet–Weibull Distribution: Properties and Application. J. Indian Soc. Probab. Stat.; 2023; 24, pp. 1-40. [DOI: https://dx.doi.org/10.1007/s41096-023-00153-w]

28. Domma, F. Some properties of the bivariate Burr type III distribution. Statistics; 2010; 44, pp. 203-215. [DOI: https://dx.doi.org/10.1080/02331880902986547]

29. Basu, A.P. Bivariate failure rate. J. Am. Stat. Assoc.; 1971; 66, pp. 103-104. [DOI: https://dx.doi.org/10.1080/01621459.1971.10482228]

30. Cox, D.R. Regression Models and Life-Tables. J. R. Stat. Soc. Ser. B (Methodol.); 1972; 34, pp. 187-220. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1972.tb00899.x]

31. Glaser, R.E. Bathtub and related failure rate characterizations. J. Am. Stat. Assoc.; 1980; 75, pp. 667-672. [DOI: https://dx.doi.org/10.1080/01621459.1980.10477530]

32. Vollan, A.M. How to Identify a Bathtub Hazard Rate. IEEE Trans. Reliab.; 1987; 36, pp. 106-108.

33. Alessandro, B. Discrete analogues of continuous bivariate probability distributions. Ann. Oper. Res.; 2022; 312, pp. 23-43.

34. McGilchrist, C.A.; Aisbett, C.W. Regression with Frailty in Survival Analysis. Biometrics; 1991; 47, pp. 461-466. [DOI: https://dx.doi.org/10.2307/2532138] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/1912255]