EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 3, Article Number 5921 ISSN 1307-5543 – ejpam.com Published by New York Business Global On the Generalization of the Bivariate extended Standard U-quadratic Distribution Idzhar A. Lakibul1,∗, Daisy Lou L. Polestico2,3, Arnulfo P. Supe2,3 1 Department of Mathematics, College of Science and Mathematics, Mindanao State University - Sulu, Jolo, Sulu, Philippines 2 Department of Mathematics and Statistics, Mindanao State University - Iligan Institute of Technology, Iligan City, Philippines 3 Premier Research Institute of Science and Mathematics - Center for Computational Analytics and Modeling, Mindanao State University - Iligan Institute of Technology, Iligan City, Philippines Abstract. This paper derives the generalized version of the Bivariate extended standard U- quadratic distribution using the compounding method. The joint probability and cumulative distribution functions of the derived distribution are obtained and it is observed that the said dis- tribution can generate bivariate shape distributions with the following properties: (i)X and Y have bathtub shapes; (ii) X and Y have inverted bathtub shapes; (iii) X and Y have constant shapes; (iv) X has a constant distribution and Y has a bathtub shape; (v) X has a constant distribution and Y has inverted bathtub shape; and (vi) X has an inverted bathtub distribution and Y has a bathtub shape. Moreover, two special cases of the generalized BeSU distribution are constructed, and these are called the Special Bivariate extended Standard U-Quadratic Type I (SBeSU-Type I) and Special Bivariate extended Standard U-Quadratic Type II (SBeSU-Type II) distributions. Fur- ther, some properties of this proposed distribution are derived such as the marginal distribution, conditional distribution, conditional moments, conditional mean, conditional variance, product and ratio moments, Pearson correlation coefficient, joint moment generating function, Kendall’s tau coefficient, Spearman’s rho, and the stress - strength parameter. In addition, maximum likeli- hood estimation is performed to estimate the parameters of the derived distribution. A simulation study is carried out to evaluate the behavior of the parameter estimates. Finally, the proposed generalization of the Bivariate eSU distribution is applied to simulated data and compared with the Bivariate extended Standard U-quadratic distribution. The results show that the proposed generalization of the bivariate eSU distribution provides a better fit on the simulated data set than the bivariate eSU distribution. 2020 Mathematics Subject Classifications: 60E05, 62E10, 65C10 Key Words and Phrases: Standard U-quadratic distribution, Kumaraswamy distribution, Bi- variate distribution, Bivariate extended Standard U-quadratic Distribution, bathtub shape distri- bution ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.5921 Email addresses: idzhar.lakibul@msusulu.edu.ph (I. A. Lakibul), daisylou.polestico@g.msuiit.edu.ph (D. L. Polestico), arnulfo.supe@g.msuiit.edu.ph (A. P. Supe) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 2 of 24 1. Introduction A prominent area of research in distribution theory is the generalization of univariate distributions to their corresponding bivariate cases. A bivariate distribution is useful for modeling the relationship between two random variables. Filus and Filus [1] introduced a method for generating bivariate or multivariate dis- tributions using linear combinations of random variables. Specifically, they derived the pseudo-Weibull and pseudo-gamma distributions as linear combinations of Weibull and gamma random variables, respectively. Shahbaz et al. [2] developed a bivariate expo- nential distribution as the compound distribution of two exponential random variables. In addition to the bivariate exponential distribution, many known bivariate distributions were derived using this, such as the bivariate pseudo-Weibull distribution [3], bivariate pseudo-Rayleigh distribution [4], bivariate pseudo-inverse Rayleigh distribution [5], bivari- ate pseudo-Gumbel distribution [6], bivariate inverse exponential distribution [7], among others. Lakibul and Tubo [8] formed a probability distribution in the interval [0, 1] called the extended standard U-quadratic distribution, and it is given in the following definition. A random variable X is said to have an extended Standard U-quadratic distribution denoted by ”eSU” if the probability density function (pdf) of X is given by f(x) = 1− λ+ 3λ(2x− 1)2, x ∈ [0, 1], (1) where λ ∈ [−0.5, 1]. It was investigated that this distribution can form three different types of shapes namely, the inverted bathtub for λ ∈ [−0.5, 0), constant for λ = 0 and bathtub for λ ∈ (0, 1]. Moreover, several properties of this distribution are detailed in the paper of Lakibul and Tubo [9]. Furthermore, the generalized version of this distribution is given in the paper of Lakibul, Polestico and Supe [10]. This distribution can be used as an option to the Beta distribution and the Kumaraswamy [11] distribution to model data with support on [0, 1], particularly those data that follow the bathtub, the inverted bathtub, and constant behavior. Lakibul, Polestico and Supe [12] expanded the eSU distribution into the Bivariate extended Standard U-quadratic Distribution (BeSU), and it is define in the following statement. A bivariate random vector (X,Y ) is said to have a Bivariate extended Standard U- quadratic (BeSU) distribution if the joint pdf of X and Y is given by f(x, y) = [ 1.5− 1.5x+ 3(1.5x− 0.5)(2y − 1)2 ] [ 1− λ+ 3λ(2x− 1)2 ] , (2) where 0 ≤ (x, y) ≤ 1 and λ ∈ [−0.5, 1]. It was observed that this BeSU distribution can describe three different types of bivariate shapes, namely, X and Y have bathtub shapes, X has a constant distribution and Y has a bathtub shape, and X has an inverted bathtub and Y has a bathtub shape. However, this distribution cannot model the bivariate shape distribution with the following properties: (i) X and Y have inverted bathtub shapes; (ii) X and Y have constant shapes; and (iii) X has a constant distribution and Y has inverted I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 3 of 24 shape. In this paper, we will use the idea of Shahbaz et al. [2] to generalize the extended standard bivariate U-quadratic distribution by adding an additional parameter to it to accommodate other combinations of the bivariate shape distribution. We will also derive some properties of the proposed generalized distribution such as the marginal distribu- tion, conditional distribution, conditional moments, conditional mean, conditional vari- ance, product and ratio moments, Pearson correlation coefficient, joint moment generating function, Kendall’s tau coefficient, Spearman’s rho, and the stress - strength parameter. Observation on the performance of the proposed generalized distribution is done by ap- plying it on a simulated dataset. The rest of the paper is structured as follows: Section 2 presents the construction of the proposed generalized Bivariate extended standard U-quadratic distribution. Section 3 provides derivations of some properties of the proposed generalized BeSU distribution. Section 4 discusses the maximum likelihood estimation for estimating the parameter of the proposed bivariate distribution. Section 5 deals with the random number generation of the proposed bivariate distribution. Section 6 presents the simulation results for assessing the behavior of the maximum likelihood estimate of the proposed bivariate distribution’s parameter. Section 7 presents the application of the proposed bivariate distribution on a simulated dataset. Section 8 gives some concluding remarks about the paper and recom- mendations for future studies. 2. The ρ - Bivariate extended Standard U-quadratic distribution This section presents the derivation of the generalized bivariate eSU distribution called as the ρ - Bivariate extended Standard U-quadratic (ρ-BeSU) Distribution. LetX be a random variable that follows an extended Standard U-quadratic distribution with pdf given in Equation (1). Let Y be another random variable such that the conditional probability density function of Y given X = x follows an extended Standard U-quadratic distribution, that is, f(y|x) = 1− v(x) + 3v(x)(2y − 1)2, y ∈ [0, 1], (3) where v(x) ∈ [−0.5, 1]. Following the idea of Shahbaz [2] and by the definition of the conditional probability, the joint probability distribution function of X and Y is given by f(x, y) =f(y|x)f(x) = [ 1− v(x) + 3v(x)(2y − 1)2 ] [ 1− λ+ 3λ(2x− 1)2 ] . (4) Note that v(x) ∈ [−0.5, 1] can be defined in many ways. In this paper, we define v(x) = 1.5xρ − 0.5, ρ ≥ 0 for x ∈ [0, 1], which reduces to v(x) defined in the construction of the simple BeSU distribution, when ρ = 1. Thus, the joint pdf of X and Y is given by f(x, y) = [ 1− (1.5xρ − 0.5) + 3(1.5xρ − 0.5)(2y − 1)2 ] [ 1− λ+ 3λ(2x− 1)2 ] = [ 1.5− 1.5xρ + 3(1.5xρ − 0.5)(2y − 1)2 ] [ 1− λ+ 3λ(2x− 1)2 ] . (5) I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 4 of 24 Definition 1. A bivariate random vector (X,Y ) is said to have a ρ - Bivariate extended Standard U-quadratic (ρ - BeSU) distribution if the joint pdf of X and Y is given by f(x, y) = [ 1.5− 1.5xρ + 3(1.5xρ − 0.5)(2y − 1)2 ] [ 1− λ+ 3λ(2x− 1)2 ] , (6) where 0 ≤ (x, y) ≤ 1, ρ ≥ 0 and λ ∈ [−0.5, 1]. Remark 1. If ρ = 1, then ρ - BeSU distribution reduces to the BeSU distribution. Theorem 1. Let (X,Y ) be the bivariate random vector with joint pdf given in Equation (6), then the joint cdf of (X,Y ) is given by F (x, y) = 3y ( 2y2 − 3y + 1 ) Mρ(x)− (2y − 3)y2F (x), (7) where Mρ(x) = (1 + 2λ)xρ+1 ρ+ 1 − 12λxρ+2 ρ+ 2 + 12λxρ+3 ρ+ 3 , and F (x) = (1 + 2λ)x− 6λx2 + 4λx3. Proof. The joint cumulative distribution function of X and Y is defined as F (x, y) = ∫ x 0 ∫ y 0 f(u, v)dvdu = ∫ x 0 f(u) [∫ y 0 f(v|u)dv ] du. Observe that,∫ y 0 f(v|u)dv = ∫ y 0 [ 1.5− 1.5uρ + 3(1.5uρ − 0.5)(2v − 1)2 ] dv =(1.5− 1.5uρ) ∫ y 0 dv + 3 (1.5uρ − 0.5) ∫ y 0 ( 4v2 − 4v + 1 ) dv =(1.5− 1.5uρ) ( v ∣∣∣∣y 0 ) + 3 (1.5uρ − 0.5) ( 4 v3 3 ∣∣∣∣y 0 − 4 v2 2 ∣∣∣∣y 0 + v ∣∣∣∣y 0 ) =(1.5− 1.5uρ) y + (1.5uρ − 0.5) ( 4y3 − 6y2 + 3y ) = [1.5− 1.5uρ + 3 (1.5uρ − 0.5)] y + (1.5uρ − 0.5) ( 4y3 − 6y2 ) =3uρy + (1.5uρ − 0.5) ( 4y3 − 6y2 ) =3uρy + 1.5uρ ( 4y3 − 6y2 ) − 0.5 ( 4y3 − 6y2 ) = ( 6y3 − 9y2 + 3y ) uρ − (2y3 − 3y2). Thus, F (x, y) = ∫ x 0 ( 1− λ+ 3λ(2u− 1)2 ) [( 6y3 − 9y2 + 3y ) uρ − (2y3 − 3y2) ] du I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 5 of 24 F (x, y) = ( 6y3 − 9y2 + 3y ) ∫ x 0 uρf(u)du− (2y3 − 3y2) ∫ x 0 f(u)du =3y ( 2y2 − 3y + 1 ) Mρ(x)− (2y − 3)y2F (x), where Mρ(x) = (1 + 2λ)xρ+1 ρ+ 1 − 12λxρ+2 ρ+ 2 + 12λxρ+3 ρ+ 3 , and F (x) =(1 + 2λ)x− 6λx2 + 4λx3. Corollary 1. Let (X,Y ) be a bivariate random variable with ρ−BeSU distribution joint CDF given in Equation (7). If λ = 0, then the joint CDF of (X,Y ) simplifies to F (x, y) = 3y ( 2y2 − 3y + 1 ) xρ+1 ρ+ 1 − (2y − 3)y2x, (8) where 0 ≤ (x, y) ≤ 1 and ρ ≥ 0. The proof follows easily from Theorem 1 by setting λ = 0 in Equation (7). Corollary 2. Let (X,Y ) be a bivariate random variable with ρ−BeSU distribution joint CDF given in Equation (7). If ρ = 0, then the joint CDF of (X,Y ) is given by F (x, y) = ( 4y2 − 6y + 3 ) ( 1 + 2λ− 6λx+ 4λx2 ) xy, (9) where 0 ≤ (x, y) ≤ 1 and λ ∈ [−0.5, 1]. The proof is straightforward by inserting ρ = 0 into Equation (7) of Theorem 1. I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 6 of 24 (a) (b) (c) (d) (e) (f) Figure 1: PDF plots of ρ-BeSU distribution for ρ = 2 and different values of λ: (a) λ = −0.5; (b) λ = −0.25; (c) λ = 0; (d) λ = 0.25; (e) λ = 0.5; and (f) λ = 1. (a) (b) (c) (d) (e) (f) Figure 2: PDF plots of ρ - BeSU distribution for ρ = 0.5 and different values of λ: (a) λ = −0.5; (b) λ = −0.25; (c) λ = 0; (d) λ = 0.25; (e) λ = 0.5; and (f) λ = 1. I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 7 of 24 Figures 1 - 2 present the bivariate plots of the joint PDF of the ρ-BeSU distribution. It is observed from the said figures that the ρ-BeSU distribution can generate different bivariate behaviors such as combinations of bathtub and inverted bathtub shapes, bathtub and bathtub shapes, inverted bathtub and inverted bathtub shapes, among others. 2.1. Special Cases of the ρ - BeSU distribution This section presents two new special cases of the proposed generalized BeSU distri- bution. 1. If λ = 0, then the ρ−BeSU distribution in Equation (6) reduces to f(x, y) = [ 1.5− 1.5xρ + 3(1.5xρ − 0.5)(2y − 1)2 ] . (10) We refer to the PDF in Equation (10) the special bivariate extended standard U-quadratic - Type I (SBeSU-Type I) distribution. (a) (b) (c) (d) (e) (f) Figure 3: PDF plots of SBeSU-Type I distribution for different values of ρ: (a) ρ = 0; (b) ρ = 0.5; (c) ρ = 1; (d) ρ = 1.5; (e) ρ = 2; and (f) ρ = 2.5. Figure 3 presents the plots of the joint PDF of the SBeSU-Type I distribution for vary- ing values of ρ. It is observed that the SBeSU-Type I distribution can represent bivariate behaviors with the following combination: (i) X has constant and Y has bathtub shapes; and (ii) X has constant and Y has inverted bathtub behaviors. I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 8 of 24 2. If ρ = 0, then the ρ−BeSU distribution reduces to f(x, y) = 3(2y − 1)2 [ 1− λ+ 3λ(2x− 1)2 ] . (11) Equation (11) is the PDF of the Special Bivariate extended Standard U-quadratic Type - II (SBeSU - Type II) distribution. (a) (b) (c) (d) (e) (f) Figure 4: PDF plots of SBeSU-Type II distribution for different values of λ: (a) λ = −0.5; (b) λ = −0.25; (c) λ = 0; (d) λ = 0.25; (e) λ = 0.5; and (f) λ = 1. Figure 4 presents the plots of the joint PDF of the SBeSU-Type II distribution for varying values of λ. It is observed that the SBeSU-Type II distribution can generate bivariate behaviors with the following combination: (i) X and Y have bathtub shapes; (ii) X has bathtub and Y has inverted bathtub behaviors; (iii) X has bathtub and Y has constant shapes. 3. Some properties of the ρ - Bivariate extended Standard U-quadratic distribution This section presents some properties of the proposed generalized BeSU distribution such as the marginal distributions, conditional distribution, conditional moment, condi- tional mean, conditional variance, product and ratio moments, covariance, Pearson corre- lation, joint moment generating function, Kendall’s tau coefficient, Spearman’s rho coef- ficient, and the Stress - Strength parameter . I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 9 of 24 Theorem 2. Let (X,Y ) be a bivariate random vector with joint probability density func- tion given in Equation (6). Then the marginal density function of Y follows an ex- tended Standard U-quadratic (eSU) distribution with parameter λ∗ = 1.5δ − 0.5, where δ = (1+2λ)ρ2+(5−2λ)ρ+6 (ρ+1)(ρ+2)(ρ+3) , λ ∈ [−0.5, 1] and ρ ≥ 0. Proof. The marginal distribution of Y is defined as f(y) = ∫ 1 0 f(x, y)dx = ∫ 1 0 [ 1.5− 1.5xρ + 3(1.5xρ − 0.5)(2y − 1)2 ] f(x)dx =1.5 [ 1− (2y − 1)2 ] ∫ 1 0 f(x)dx− 1.5 [ 1− 3(2y − 1)2 ] ∫ 1 0 xρf(x)dx. (12) Observe that the first part of Equation (12) simplifies to 1.5 [ 1− (2y − 1)2 ] since f(x) is a pdf of an eSU distribution for a random variable X. Furthermore,∫ 1 0 xρf(x)dx = ∫ 1 0 xρ [ (1− λ) + 3λ(2x− 1)2 ] dx = (1 + 2λ)ρ2 + (5− 2λ)ρ+ 6 (ρ+ 1)(ρ+ 2)(ρ+ 3) = δ. It follows that Equation (12) simplifies to f(y) =1− λ∗ + 3λ∗(2y − 1)2, where y ∈ [0, 1] and λ∗ = 1.5δ − 0.5. Thus, the marginal distribution of Y is eSU with parameter λ∗. Corollary 3. Let Y be a random variable that follows an eSU distribution with parameter λ∗, where λ∗ and δ are given in Theorem 2. If λ = 0, then the PDF of Y follows an eSU distribution with parameter λ∗ = 1−0.5ρ ρ+1 , where ρ ≥ 0. Corollary 4. Let Y be a random variable that follows an eSU distribution with parameter λ∗, where λ∗ and δ are given in Theorem 2. If ρ = 0, then the PDF of Y follows an eSU distribution with parameter λ∗ = 1. Theorem 3. Let X and Y be any two random variables with joint pdf given in Equa- tion (6). If the marginal distribution of Y is given in Theorem 2, then the conditional distribution of X given Y = y is f(X|Y = y) = [ 1.5− 1.5xρ + 3(1.5xρ − 0.5)(2y − 1)2 ] [ 1− λ+ 3λ(2x− 1)2 ] 1− (1.5δ − 0.5) + 3(1.5δ − 0.5)(2y − 1)2 , (13) where δ = (1+2λ)ρ2+(5−2λ)ρ+6 (ρ+1)(ρ+2)(ρ+3) , λ ∈ [−0.5, 1] and ρ ≥ 0. The proof follows directly from the definition of the conditional distribution of X given Y = y. I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 10 of 24 Theorem 4. Let X and Y be any two random variables with joint pdf given in Equation (6). If the marginal distribution of Y is given in Theorem 2 then the rth conditional moment of X given Y = y is E [Xr|y] = 1.5 f(y) {[ 1− (2y − 1)2 ] E [Xr]− [ 1− 3(2y − 1)2 ] E [ Xr+ρ ]} , (14) where f(y) is the marginal distribution of Y of the ρ−BeSU distribution, E [Xr] = (1 + 2λ)r2 + (5− 2λ)r + 6 (r + 1)(r + 2)(r + 3) , and E [ Xr+ρ ] = (1 + 2λ)(r + ρ)2 + (5− 2λ)(r + ρ) + 6 (r + ρ+ 1)(r + ρ+ 2)(r + ρ+ 3) . Proof. The rth conditional moment of X given Y = y is defined as E [Xr|y] = ∫ 1 0 xrf(X|Y = y)dx E [Xr|y] = ∫ 1 0 xr f(x, y) f(y) dx. It follows that the rth conditional moment of X given Y = y becomes E [Xr|y] = 1 f(y) ∫ 1 0 xr [ 1.5 ( 1− (2y − 1)2 ) − 1.5x ( 1− 3(2y − 1)2 )] f(x|λ)dx = 1.5 ( 1− 3(2y − 1)2 ) f(y) [ (1 + 2λ)r2 + (5− 2λ)r + 6 (r + 1)(r + 2)(r + 3) ] − 1.5 ( 1− 3(2y − 1)2 ) f(y) [ (1 + 2λ)(r + ρ)2 + (5− 2λ)(r + ρ) + 6 (r + ρ+ 1)(r + ρ+ 2)(r + ρ+ 3) ] = 1.5 f(y) {[ 1− (2y − 1)2 ] E [Xr]− [ 1− 3(2y − 1)2 ] E [ Xr+ρ ]} , where f(y) is the marginal distribution of Y of the ρ−BeSU distribution, E [Xr] = (1 + 2λ)r2 + (5− 2λ)r + 6 (r + 1)(r + 2)(r + 3) , and E [ Xr+ρ ] = (1 + 2λ)(r + ρ)2 + (5− 2λ)(r + ρ) + 6 (r + ρ+ 1)(r + ρ+ 2)(r + ρ+ 3) . I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 11 of 24 Remark 2. The conditional mean of X given Y is given by E [X|y] = 1.5 f(y) { 0.5 [ 1− (2y − 1)2 ] − [ 1− 3(2y − 1)2 ] E [ Xρ+1 ]} , (15) where ρ ≥ 0, f(y) is the marginal distribution of Y of the ρ - BeSU distribution, and E [ Xρ+1 ] = (1 + 2λ)(ρ+ 1)2 + (5− 2λ)(ρ+ 1) + 6 (ρ+ 2)(ρ+ 3)(ρ+ 4) . Remark 3. The conditional variance of X given Y is given by V ar(X|y) = E[X2|y]− (E[X|y])2 , (16) where E[X|y] is the conditional mean of X given Y , E [ X2|y ] = 1.5 f(y) {[ 1− (2y − 1)2 ](5 + λ 15 ) − [ 1− 3(2y − 1)2 ] E [ Xρ+2 ]} , (17) ρ ≥ 0, λ ∈ [−0.5, 1], f(y) is the marginal distribution of Y of the ρ - BeSU distribution, and E [ Xρ+2 ] = (1 + 2λ)(ρ+ 2)2 + (5− 2λ)(ρ+ 2) + 6 (ρ+ 3)(ρ+ 4)(ρ+ 5) . Theorem 5. Let X and Y be any two random variables with joint pdf given in Equation (6), then the product and ratio moments are given by E [XrY s] = 1 (s+ 1)(s+ 2)(s+ 3) [ 6(s+ 1)E[Xr] + 3s(s− 1)E[Xr+ρ] ] (18) and E [ XrY −s ] = 1 (1− s)(2− s)(3− s) [ 6(1− s)E[Xr] + 3s(s+ 1)E[Xr+ρ] ] , (19) where E[Xr] = (1+2λ)r2+(5−2λ)r+6 (r+1)(r+2)(r+3) and E[Xr+ρ] = (1+2λ)(r+ρ)2+(5−2λ)(r+ρ)+6 (r+ρ+1)(r+ρ+2)(r+ρ+3) . Proof. The product moment of X and Y is defined as E [XrY s] = ∫ 1 0 ∫ 1 0 xrysf(x, y)dxdy = 1 (s+ 1)(s+ 2)(s+ 3) [ 6(s+ 1) ∫ 1 0 xrf(x)dx + 3s(s− 1) ∫ 1 0 xr+ρf(x)dx ] = 1 (s+ 1)(s+ 2)(s+ 3) [ 6(s+ 1)E[Xr] + 3s(s− 1)E[Xr+ρ] ] , I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 12 of 24 where E[Xr] = (1 + 2λ)r2 + (5− 2λ)r + 6 (r + 1)(r + 2)(r + 3) , and E[Xr+ρ] = (1 + 2λ)(r + ρ)2 + (5− 2λ)(r + ρ) + 6 (r + ρ+ 1)(r + ρ+ 2)(r + ρ+ 3) . Next, the ratio moment of X and Y is defined as E [ Xr Y s ] =E [ XrY −s ] = 1 (1− s)(2− s)(3− s) [ 6(1− s)E[Xr] + 3s(s+ 1)E[Xr+ρ] ] . Corollary 5. Let X and Y be random variables with joint PDF of a ρ−BeSU distribution with product moment given in Equation (18). If r = s = 1, then E[XY ] = 1 4 . (20) Proof. Substituting r = s = 1 to Equation (18) gives us the E[XY ] = 1 4 since E[X] = 1 2 . Remark 4. The covariance of X and Y is given by Cov(X,Y ) = E[XY ]− E[X]E[Y ] = 1 4 − 1 2 ( 1 2 ) = 0. Remark 5. The Pearson correlation or correlation coefficient of X and Y is zero. Note that the correlation coefficient of X and Y is zero since the covariance of X and Y is zero. This implies that there is no linear relationship between X and Y . A nonlinear form may best describe the relationship between X and Y . Theorem 6. Let (X,Y ) be a bivariate random vector that follows a ρ - bivariate extended Standard U-quadratic distribution, then the joint moment generating function of X and Y is given by MX,Y (t1, t2) = 3 (s+ 2)(s+ 3) ∞∑ r=0 tr1 r! ∞∑ s=0 ts2 s! [ 2E[Xr] + s(s− 1)E[Xr+ρ] (s+ 1) ] , I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 13 of 24 where E[Xr] = (1 + 2λ)r2 + (5− 2λ)r + 6 (r + 1)(r + 2)(r + 3) , and E[Xr+ρ] = (1 + 2λ)(r + ρ)2 + (5− 2λ)(r + ρ) + 6 (r + ρ+ 1)(r + ρ+ 2)(r + ρ+ 3) . Proof. The joint moment generating function of X and Y is defined as MX,Y (t1, t2) = ∫ 1 0 ∫ 1 0 et1x+t2yf(x, y)dydx = ∫ 1 0 ∫ 1 0 et1xet2yf(y|x)f(x)dydx. Recall that etx = ∑∞ r=0 tr r!x r, then we have MX,Y (t1, t2) = ∫ 1 0 ∫ 1 0 ∞∑ r=0 tr1 r! xr ∞∑ s=0 ts2 s! ysf(y|x)f(x)dydx = ∞∑ r=0 tr1 r! ∞∑ s=0 ts2 s! ∫ 1 0 ∫ 1 0 xrysf(y|x)f(x)dydx = ∞∑ r=0 tr1 r! ∞∑ s=0 ts2 s! ∫ 1 0 xrf(x) (∫ 1 0 ysf(y|x)dy ) dx = ∞∑ r=0 tr1 r! ∞∑ s=0 ts2 s! ∫ 1 0 xrf(x) ( 6(s+ 1) + 3s(s− 1)xρ (s+ 1)(s+ 2)(s+ 3) ) dx = 1 (s+ 2)(s+ 3) ∞∑ r=0 tr1 r! ∞∑ s=0 ts2 s! [ 6 ∫ 1 0 xrf(x)dx+ 3s(s− 1) s+ 1 ∫ 1 0 xr+ρf(x)dx ] = 3 (s+ 2)(s+ 3) ∞∑ r=0 tr1 r! ∞∑ s=0 ts2 s! [ 2E[Xr] + s(s− 1) s+ 1 E[Xr+ρ] ] , where E[Xr] = (1 + 2λ)r2 + (5− 2λ)r + 6 (r + 1)(r + 2)(r + 3) , and E[Xr+ρ] = (1 + 2λ)(r + ρ)2 + (5− 2λ)(r + ρ) + 6 (r + ρ+ 1)(r + ρ+ 2)(r + ρ+ 3) . I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 14 of 24 Lemma 1. Let X be a random variable that follows an eSU distribution, then∫ 1 0 F (x)f(x)dx = 1 2 , where F (x) and f(x) are CDF and PDF of the eSU distribution. Proof. Let X be a random variable that follows an eSU distribution with CDF F (x) and PDF f(x), respectively. Then, using the rth moment of eSU distribution, we have∫ 1 0 F (x)f(x)dx = ∫ 1 0 [ (1 + 2λ)x− 6λx2 + 4λx3 ] f(x)dx =(1 + 2λ)E[X]− 6λE[X2] + 4λE[X3] = 1 2 . Theorem 7. Let (X,Y ) be the random vector that follows a ρ−BeSU distribution, then the Kendall’s tau coefficient is zero. Proof. Let (X,Y ) be the random vector that follows a ρ − BeSU distribution with joint PDF and CDF given in Equations (6) and (7), respectively. Then, the Kendall’s tau coefficient is defined as τ =4 ∫ 1 0 ∫ 1 0 FX,Y (x, y)fX,Y (x, y)dxdy − 1, (21) where FX,Y (x, y) and fX,Y (x, y) are joint CDF and PDF of the random variables X and Y , respectively. Let us first consider∫ 1 0 FX,Y (x, y)fX,Y (x, y)dv = ∫ 1 0 [ 3y ( 2y2 − 3y + 1 ) Mρ(x)− (2y − 3)y2F (x) ] {[ 1.5− 1.5xρ + 3(1.5xρ − 0.5)(2y − 1)2 ] f(x) } dy =Mρ(x)f(x) ∫ 1 0 3y ( 2y2 − 3y + 1 ) {[ 1.5− 1.5xρ + 3(1.5xρ − 0.5)(2y − 1)2 ]} dy − F (x)f(x) ∫ 1 0 (2y − 3)y2×[ 1.5− 1.5xρ + 3(1.5xρ − 0.5)(2y − 1)2 ] dy. Now set A ≡ 1.5− 1.5xρ+3(1.5xρ− 0.5)(2y− 1)2, which is the f(y|x), where Y |X follows an eSU distribution with parameter 1.5xρ − 0.5.∫ 1 0 FX,Y (x, y)fX,Y (x, y)dy =Mρ(x)f(x) ∫ 1 0 6y3Ady −Mρ(x)f(x) ∫ 1 0 9y2Ady I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 15 of 24 +Mρ(x)f(x) ∫ 1 0 3yAdy − F (x)f(x) ∫ 1 0 2y3Ady + F (x)f(x) ∫ 1 0 3y2Ady =Mρ(x)f(x) ( 6E[Y 3|X = x]− 9E[Y 2|X = x] + 3E[Y |X = x] ) − F (x)f(x) ( 2E[Y 3|X = x] + 3E[Y 2|X = x] ) . Observe from the construction of the ρ − BeSU distribution that the random variable Y |X follows an eSU distribution with parameter 1.5xρ − 0.5. Now, using the rth moment of eSU distribution for Y |X, we have∫ 1 0 FX,Y (x, y)fX,Y (x, y)dv =Mρ(x)f(x) [ 6 ( 4 + 3xρ 20 ) − 9 ( 3 + xρ 10 ) + 3 ( 1 2 )] − F (x)f(x) [ 2 ( 4 + 3xρ 20 ) + 3 ( 3 + xρ 10 )] = 1 2 F (x)f(x). Now, we have τ =4 ∫ 1 0 1 2 f(x)F (x)dx− 1 =2 ∫ 1 0 f(x)F (x)dx− 1. From Lemma (1), we have ∫ 1 0 F (x)f(x)du = 1 2 . Thus, the Kendall’s tau coefficient becomes τ = 0. Theorem 8. Let (X,Y ) be the random vector that follows a ρ−BeSU distribution, then the Spearman’s rho coefficient is zero. Proof. Let (X,Y ) be the random vector that follows a ρ − BeSU distribution with joint PDF and CDF given in Equations (6) and (7), respectively. Then, the Spearman’s rho coefficient is defined as ρ∗ =12 ∫ 1 0 ∫ 1 0 [F (x, y)− F (x)F (y)] f(x)f(y)dxdy I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 16 of 24 ρ∗ =12 ∫ 1 0 ∫ 1 0 F (x, y)f(x)f(y)dxdy − 12 ∫ 1 0 ∫ 1 0 F (x)F (y)f(x)f(y)dxdy. Set S = ∫ 1 0 F (x, y)f(x)f(y)dy. It follows that S =f(x) ∫ 1 0 F (x, y)f(y)dy =f(x) ∫ 1 0 [ 3y ( 2y2 − 3y + 1 ) Mρ(x)− (2y − 3)y2F (x) ] f(y)dy =f(x)Mρ(x) ∫ 1 0 ( 6y3 − 9y2 + 3y ) f(y)dy − f(x)F (x) ∫ 1 0 ( 2y3 − 3y2 ) f(y)dy = 1 2 f(x)F (x). Also, from Lemma (1), we have ∫ 1 0 F (x)f(x)dx = 1 2 . Hence, the Spearman’s rho is simplified to ρ∗ =12 ∫ 1 0 (∫ 1 0 F (x, y)f(x)f(y)dy ) dx− 12 ∫ 1 0 f(x)F (x) (∫ 1 0 F (y)f(y)dy ) dx = 0. In the following theorem we derive the Stress-Strength parameter for the ρ − BeSU distribution. The Stress-Strength parameter is a measure used in reliability engineering and Statistics to evaluate the performance and reliability systems under stress. The vari- able X represents the strength of a system or component, while variable X represents the applied stress. The parameter P (Y < X) measures the probability that the system’s strength exceeds the applied stress, which is critical in assessing the reliability of materials and components. Theorem 9. Let (X,Y ) be a bivariate random vector that follows a ρ - bivariate extended Standard U-quadratic distribution, then the Stress - Strength parameter of X and Y is given by P (Y < X) =3 ( (1 + 2λ)(ρ+ 1)2 + (5− 2λ)(ρ+ 1) + 6 (ρ+ 2)(ρ+ 3)(ρ+ 4) ) − 9 ( (1 + 2λ)(ρ+ 2)2 + (5− 2λ)(ρ+ 2) + 6 (ρ+ 3)(ρ+ 4)(ρ+ 5) ) + 6 ( (1 + 2λ)(ρ+ 3)2 + (5− 2λ)(ρ+ 3) + 6 (ρ+ 4)(ρ+ 5)(ρ+ 6) ) + 5 + 4λ 10 , (22) where ρ ≥ 0 and λ ∈ [−0.5, 1]. I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 17 of 24 Proof. The stress - strength parameter of random variables X and Y is defined by P(Y < X) = ∫ 1 0 ∫ x 0 f(x, y)dydx = ∫ 1 0 ∫ x 0 [ 1.5− 1.5xρ + 3(1.5xρ − 0.5)(2y − 1)2 ] [ 1− λ+ 3λ(2x− 1)2 ] dydx =3E[Xρ+1]− 9E[Xρ+2] + 3E[X2] + 6E[Xρ+3]− 2E[X3] =3 ( (1 + 2λ)(ρ+ 1)2 + (5− 2λ)(ρ+ 1) + 6 (ρ+ 2)(ρ+ 3)(ρ+ 4) ) − 9 ( (1 + 2λ)(ρ+ 2)2 + (5− 2λ)(ρ+ 2) + 6 (ρ+ 3)(ρ+ 4)(ρ+ 5) ) + 6 ( (1 + 2λ)(ρ+ 3)2 + (5− 2λ)(ρ+ 3) + 6 (ρ+ 4)(ρ+ 5)(ρ+ 6) ) + 5 + 4λ 10 , where ρ ≥ 0 and λ ∈ [−0.5, 1]. 4. Maximum Likelihood Estimation Let (X1, Y1),(X2, Y2),...,(Xn, Yn) be a random sample of size n from a ρ - Bivariate extended Standard U-quadratic (ρ - BeSU) Distribution. Then the likelihood function is defined by L = n∏ i [ 1.5− 1.5xρi + 3(1.5xρi − 0.5)(2yi − 1)2 ] [ 1− λ+ 3λ(2xi − 1)2 ] , with its log-likelihood function is given by logL = n∑ i log [ 1.5− 1.5xρi + 3(1.5xρi − 0.5)(2yi − 1)2 ] + n∑ i log [ 1− λ+ 3λ(2xi − 1)2 ] . The partial derivative of logL with respect to the parameters λ and ρ are respectively, given by ∂ logL ∂λ = n∑ i 3(2xi − 1)2 − 1 1− λ+ 3λ(2xi − 1)2 , and ∂ logL ∂ρ = n∑ i 4.5xρi log xi − 1.5xρi log xi(2yi − 1)2 1.5− 1.5xρi + 3(1.5xρi − 0.5)(2yi − 1)2 . I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 18 of 24 The maximum likelihood estimates of the parameters λ and ρ of the BeSU distribution are computed by solving the following system of non-linear equations: n∑ i 3(2xi − 1)2 − 1 1− λ+ 3λ(2xi − 1)2 = 0, and n∑ i 4.5xρi log xi − 1.5xρi log xi(2yi − 1)2 1.5− 1.5xρi + 3(1.5xρi − 0.5)(2yi − 1)2 = 0. 5. Random Number Generation This section presents the algorithm for the generation of bivariate random numbers from the ρ - Bivariate extended Standard U-quadratic (ρ - BeSU) distribution. Let us first consider the random number generation from the T-extended Standard U-quadratic (TeSU)-G family of distributions. To generate random samples from TeSU-G family of distributions, we follow the algorithm proposed by Lakibul and Tubo [8]. The cumulative distribution function of the TeSU-G family is given by F (x) = (1 + 2λ)G(x)− 6λ(G(x))2 + 4λ(G(x))3, (23) where λ ∈ [−0.5, 1] and G(x) is any baseline cumulative distribution function. Due to the nonlinear nature of the CDF, it is not possible to obtain a simple closed-form expression for its inverse. Therefore, a numerical algorithm is employed to approximate the inverse CDF, enabling efficient random number generation. The algorithm to generate random numbers from TeSU-G family is given as follows. Let v follow a uniform distribution (0, 1). Step 1. Compute Q = 1 2 ( 1− 1 λ ) , λ ̸= 0; R = 1− 2v 16λ . Step 2. If R2 > Q3, then compute A = −sign(R) ( |R|+ √ R2 −Q3 ) 1 3 ; B = A, ifA = 0 Q A , otherwise ; x = G−1 ( A+B + 1 2 ) . I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 19 of 24 Otherwise, θ = arccos ( R√ Q3 ) ; x = G−1 ( 1 2 − 2 √ Q cos ( θ − 2π 3 )) , where G−1(x) is the inverse function of any baseline distribution function G(x). If λ = 0, then x = G−1(v). Note that, the extended Standard U-quadratic distribution is derived from the T-extended Standard U-quadratic - G family of distributions by taking G(x) = x. Thus, we have the following modified algorithm to generate random numbers from an ex- tended Standard U-quadratic distribution. Let v follow a uniform distribution (0, 1). If λ = 0, then x = v. Otherwise, it is given as follows: Step 1.* Compute Q = 1 2 ( 1− 1 λ ) , λ ̸= 0; R = 1− 2v 16λ . Step 2.* If R2 > Q3, then A = −sign(R) ( |R|+ √ R2 −Q3 ) 1 3 ; B = A, ifA = 0 Q A , otherwise ; x = A+B + 1 2 . Otherwise, θ = arccos ( R√ Q3 ) ; x = 1 2 − 2 √ Q cos ( θ − 2π 3 ) . The inverse CDF of the TeSU-G family is not readily available in a simple closed form due to its nonlinear and complex structure. As a result, the numerical algorithm presented above is used to approximate the inverse CDF, enabling random number generation in a computationally feasible manner. This approach is a standard method for generating random variables from distributions where the inverse CDF is difficult to compute directly. I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 20 of 24 In addition, to generate a random sample from the ρ - bivariate extended Standard U-quadratic distribution, we use the conditional approach given in the following algorithm: Steps Description 1 Draw a random sample X of size n from an extended Standard U-quadratic distribution with parameter λ. 2 For each observation X, draw a sample of size 1 from an extended Standard Standard U-quadratic distribution with parameter 1.5xρ − 0.5. Repeat this process for all observations of X. Denote this sample as Y . 3 Finally, the desired random sample is (x, y). 6. Simulation Study This section presents the simulation results to assess the behavior of the maximum likelihood estimate of the parameter of the proposed bivariate distribution. The simulation algorithm is given as follows: Steps Description 1 Draw sample of size n, n = 50, 100, 200, 500, 1000 from a ρ - bivariate extended Standard U-quadratic distribution with parameter λ and ρ using the algorithm given in the previous section. 2 Using the bivariate sample (x, y) obtained in Step 1 above, compute the maximum likelihood estimate of λ and ρ. 3 Repeat the preceding Steps 1-2 N = 1000 times to get 1000 estimates of λ and ρ. 4 Compute the mean, bias, and mean squared error (MSE) of the 1000 estimates obtained in Step 3 to get the desired results The mean (AE), Bias, and MSE are, respectively, defined by AE = N∑ i=1 λi N , Bias = AE − λ and MSE = N∑ i=1 (λi − λ)2 N . Considering two different sets of values of the parameters λ and ρ, the following table shows that as n becomes large, the average estimate (AE) of the parameters λ and ρ are going closer to the true values, respectively, while the bias and the MSE diminish to zero. Thus, the maximum likelihood estimates of the proposed distribution parameters are consistent. I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 21 of 24 Table 1: Results of the simulation study for the following set of values of the parameters of ρ - BeSU distribution: (a.) λ = −0.3 and ρ = 0.1 (Left); and (b.) λ = 0.5 and ρ = 1.5 (Right) n MLE AE Bias MSE 50 λ̂ −0.304 −0.004 0.017 100 λ̂ −0.298 0.002 0.009 200 λ̂ −0.300 0.000 0.004 500 λ̂ −0.301 −0.001 0.002 1000 λ̂ −0.301 −0.001 0.001 50 ρ̂ 0.106 0.006 0.008 100 ρ̂ 0.102 0.002 0.004 200 ρ̂ 0.102 0.002 0.002 500 ρ̂ 0.102 0.002 0.001 1000 ρ̂ 0.101 0.001 0.000 n MLE AE Bias MSE 50 λ̂ 0.495 −0.005 0.024 100 λ̂ 0.504 0.004 0.011 200 λ̂ 0.499 −0.001 0.006 500 λ̂ 0.499 −0.001 0.002 1000 λ̂ 0.498 −0.002 0.001 50 ρ̂ 1.808 0.308 1.720 100 ρ̂ 1.598 0.098 0.271 200 ρ̂ 1.541 0.041 0.122 500 ρ̂ 1.525 0.025 0.045 1000 ρ̂ 1.515 0.015 0.023 7. Application In this section, we apply the proposed generalized bivariate distribution and compare with the BeSU distribution. Here, we use the simulated bivariate data from the study of Lakibul, Polestico and Supe [12]. This bivariate data has the following property: X and Y have bathtub shapes. It was generated from the Bivariate Kumaraswamy distribution of Lakibul, Polestico and Supe [12], with a = 0.5 and b = 0.5. The following figures show the joint and marginal distributions of X and Y of the said bivariate simulated data. (a) (b) (c) Figure 5: histogram plot of the simulated data : (a) for variable X; (b) for variables X and Y ; and (c) for variable Y . In the analysis, we use the R-package ”bbmle” to compute the maximum likelihood I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 22 of 24 estimates of the parameters of the proposed ρ - BeSU and BeSU distributions. In addition, the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) are used to assess and compare the performance of the proposed generalized bivariate distributions. Table 2: Estimates and some diagnostic values of the fitted models for the simulated dataset. Distribution Estimate Std.Error −2logLik AIC BIC ρ−BeSU λ̂ = 0.4891090 0.0142492 -2724.289 -2720.289 -2707.255 ρ̂ = 0.0953629 0.0055758 BeSU λ̂ = 0.489109 0.014249 2243.733 2245.733 2252.25 Table 2 shows some diagnostic statistics and maximum likelihood estimates of the fitted models for the simulated dataset. It is observed that the ρ - BeSU distribution has smaller values of the AIC and BIC than the BeSU distribution. Thus, the ρ - BeSU distribution provides a better fit for this simulated data than the BeSU distribution. 8. Conclusions and Recommendations In this paper, we have derived the generalized version of the Bivariate extended Standard U-quadratic distribution, referred to as the ρ-Bivariate extended Standard U- quadratic (ρ-BeSU) distribution. Several important properties of the proposed ρ-BeSU distribution were computed, including the marginal and conditional distributions, con- ditional moments, conditional mean, conditional variance, product and ratio moments, Pearson correlation coefficient, joint moment generating function, Kendall’s tau coeffi- cient, Spearman’s rho coefficient, and the stress-strength parameter. Maximum likelihood estimation was applied to estimate the parameters of the ρ-BeSU distribution, and a simulation study was conducted to assess the behavior of the parameter estimates. The results of the simulation study demonstrated that the maximum likelihood estimate of the parameter of the ρ-BeSU distribution is consistent. Furthermore, it was observed that the proposed ρ-BeSU distribution provides a better fit for the simulated bivariate dataset compared to the standard BeSU distribution, highlighting its potential as a more flexible model for bivariate data. While the ρ-BeSU distribution offers notable improvements, it is important to acknowledge some potential limitations. For instance, the current study focuses on a specific class of bivariate distributions, and there may be alternative distribu- tions or estimation techniques that could provide additional benefits, particularly in cases with more complex dependence structures. Future research could explore extensions of the ρ-BeSU distribution to higher dimensions, which would enhance its applicability to multivariate data. Furthermore, alternative estimation methods, such as Bayesian infer- ence or the use of copulas, could be considered to better account for varying dependence structures in different applications. In addition to these theoretical advancements, we recommend further investigations using the ρ-BeSU distribution to model the failure rates of two related components in a system or to model the lifetimes of two related electronic I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 23 of 24 devices, particularly when the failure rates or lifetimes exhibit bathtub, inverted, or con- stant shapes on the interval [0, 1]. Additionally, future studies could consider extending the ρ-BeSU distribution by comparing it to other bivariate families of distributions, such as the bivariate Marshall-Olkin family or copulas, to evaluate its performance in diverse practical contexts. Acknowledgements The authors are grateful to the anonymous referees for their valuable comments and suggestions. Moreover, Idzhar A. Lakibul is also grateful to the Department of Science and Technology - Accelerated Science and Technology Human Resource Development Program (DOST-ASTHRDP) for giving him financial support to study at Mindanao State Univer- sity - Iligan Institute of Technology (MSU-IIT). In addition, this paper is also supported by the DOST-ASTHRDP, MSU-Iligan Institute of Technology and the Mindanao State University - Sulu. References [1] J. K. Filus and L. Z. Filus. On some new classes of Multivariate Probability Distri- butions. Pakistan Journal of Statistics, 22:21–42, 2006. [2] S. Shahbaz, M. Q. Shahbaz, and M. Mohsin. On Concomitants of Order Statistics for Bivariate Pseudo Exponential distribution. World Applied Sciences Journal, 6:1151– 1156, 2009. [3] S. Shahbaz and M. Ahmad. Concomitants of Order Statistics for Bivariate Pseudo- Weibull Distribution. World Applied Sciences Journal, 6:1409–1412, 2009. [4] M. Q. Shahbaz and S. Shahbaz. Order Statistics and Concomitants of Bivariate Pseudo- Rayleigh Distribution. World Applied Sciences Journal, 7:826–828, 2009. [5] M. Mohsin, M. Ahmad, S. Shahbaz, and M. Q. Shahbaz. Concomitants of Lower Records for Bivariate Pseudo – Inverse Rayleigh Distribution. Sci. Int. (Lahore), 21:21–23, 2009. [6] M. Q. Shahbaz, S. Shahbaz, and A. Rafiq. A New Bivariate Gumbel distribution. Nonlinear Analysis Forum, 16:133–136, 2011. [7] S. H. Shahbaz and M. Q. Shahbaz. On Concomitants of Dual Generalized Order Statistics for a Bivariate Inverse Exponential Distribution. European Journal of Pure and Applied Mathematics, 11:929–936, 2018. [8] I. A. Lakibul and B. F. Tubo. On the TeSU-G family of distributions applied to life data analysis. Reliability:Theory and Applications, 18:24–38, 2023. [9] I. A. Lakibul and B. F. Tubo. On the Four-Parameter T-extended Standard U- quadratic Exponentiated Weibull distribution. The Mindanawan Journal of Mathe- matics, 5:17–33, 2023. [10] I. A. Lakibul, D. L. Polestico, and A. P. Supe. On the Generalized Version of the extended Standard U-quadratic Distribution. European Journal of Pure and Applied Mathematics, 18(2), 2025. I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 18 (3) (2025), 5921 24 of 24 [11] P. Kumaraswamy. A generalized probability density function for double-bounded random processes. Journal of Hydrology, 46:79–88, 1980. [12] I. A. Lakibul, D. L. Polestico, and A. P. Supe. On the Bivariate Extension of the extended Standard U-quadratic Distribution. European Journal of Pure and Applied Mathematics, 17:790–809, 2024.