EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 17, No. 2, 2024, 790-809 ISSN 1307-5543 – ejpam.com Published by New York Business Global On the Bivariate Extension of the extended Standard U-quadratic Distribution Idzhar A. Lakibul1,2,∗, Daisy Lou L. Polestico1,2, Arnulfo P. Supe1,2 1 Department of Mathematics and Statistics, Mindanao State University - Iligan Institute of Technology, Iligan City, Philippines 2 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 bivariate version of the extended standard U-quadratic (eSU) distribution using the compounding method or pseudo-family of distributions. The joint probabil- ity and cumulative distribution functions of the derived distribution are obtained and it is observed that the said distribution can generate bivariate shape distributions with the following properties: (i) X and Y have bathtub shapes; (ii) X has a constant distribution and Y has a bathtub shape; and (iii) X has an inverted bathtub and Y has a bathtub shape. Moreover, some properties of this proposed distribution are derived such as the marginal distribution, conditional distribution, conditional moments, product and ratio moments, Pearson correlation coefficient, joint moment generating function, and the stress - strength parameter. Further, the maximum likelihood esti- mation is performed to estimate the parameters of the derived distribution. A simulation study is carried out to evaluate the behavior of the estimates of the parameters. A new bivariate Ku- maraswamy distribution is derived and used to simulate bivariate data with X and Y having bathtub shapes. A new bivariate version of the Cubic Transmuted Uniform (CTU) distribution is also derived. Finally, the proposed Bivariate eSU distribution is applied to simulated data and compared with the Bivarite Cubic Transmuted Uniform distribution. The results show that the proposed Bivariate eSU distribution provides a better fit on the simulated dataset compared with the Bivariate CTU distribution. 2020 Mathematics Subject Classifications: 60E05, 62E10, 65C10 Key Words and Phrases: Standard U-quadratic distribution, Kumaraswamy distribution, Bi- variate distribution, Bivariate Pseudo Family, bathtub shape distribution ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v17i2.5136 Email addresses: idzhar.lakibul@g.msuiit.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 790 © 2024 EJPAM All rights reserved. I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 791 1. Introduction Nowadays, one of the developing research areas in the field of distribution theory is the generalization of existing univariate distribution into a bivariate case. A bivariate distribution is useful for modeling two related random variables. Filus and Filus [1] proposed a method of generating bivariate or multivariate distri- butions as the linear combinations of two or more random variables. In particular, they derived the pseudo-Weibull and pseudo-gamma distributions as the linear combinations of the Weibull and gamma random variables, respectively. Shahbaz et al. [10] formed a bivariate exponential distribution as the compound distribution of two exponential ran- dom variables. Aside from the bivariate exponential distribution, many derived bivariate distributions were formulated using this idea such as the bivariate pseudo-Weibull distri- bution [9], bivariate pseudo-Rayleigh distribution [7], bivariate pseudo-inverse Rayleigh distribution [5], bivariate pseudo-Gumbel distribution [8], bivariate inverse exponential distribution [11], among others. Lakibul and Tubo [4] developed a probability distribution called the extended Standard U-quadratic distribution with support on [0, 1] 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 observed that this distribution can generate three different types of shape for its pdf, namely, inverted bathtub for λ ∈ [−0.5, 0), constant for λ = 0 and bathtub for λ ∈ (0, 1]. In addition, some properties of this distribution can be found in the paper of Lakibul and Tubo [3]. This distribution can be used as an alternative to the Beta distribution and Kumaraswamy [2] distribution for modeling data with support on [0, 1], particularly, those data follow bathtub, inverted bathtub, and constant behavior. However, if we want to model simultaneously two related random variables where each of the variables has a bathtub shape then the eSU distribution cannot be used. Thus, there is a need to extend this distribution into the bivariate case. In this paper, we will follow the idea of Shahbaz et al. [10] to expand the extended Stan- dard U-quadratic distribution into a bivariate case. We will also derive some properties of the proposed distribution such as the marginal distribution, conditional distribution, conditional moments, product and ratio moments, Pearson correlation coefficient, joint moment generating function, and the stress - strength parameter. Investigation on the performance of the proposed distribution is done by applying it on a simulated dataset. The rest of the paper is structured as follows: Section 2 presents the construction of the Bivariate extended standard U-quadratic (BeSU) distribution . Section 3 provides derivations of some properties of the proposed BeSU distribution. Section 4 discusses the maximum likelihood estimation for estimating the proposed bivariate distribution’s parameter. 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 I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 792 7 presents the application of the proposed bivariate distribution on a simulated dataset. Section 8 gives some concluding remarks about the paper and recommendations for future studies. 2. The Bivariate extended Standard U-quadratic distribution This section presents the derivation of the new bivariate distribution with support on [0, 1]× [0, 1]. 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], (2) where v(x) ∈ [−0.5, 1]. Following the idea of Shahbaz [10] and by the definition of the conditional probability, the joint probability distribution function of X and Y is given by f(x, y) = [ 1− v(x) + 3v(x)(2y − 1)2 ] [ 1− λ+ 3λ(2x− 1)2 ] . Note that, v(x) ∈ [−0.5, 1], then v(x) can be defined in many ways. For example, for x ∈ [0, 1], v(x) can be defined as v(x) = 1.5xτ−0.5, τ > 0 or v(x) = 1.5(1−(1−x)a)b−0.5, a > 0, b > 0. In this paper, we focus on the special case where we have no additional parameter in the model to have a simple bivariate model with a single parameter. A simple bivariate model with a single parameter is usually preferred for modeling some bivariate data than a bivariate model with more parameters because sometimes a model with more parameters will over parameterize the data. Hence, we define v(x) as v(x) = 1.5x − 0.5, where x ∈ [0, 1]. Thus, 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 ] . 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 ] , (3) where 0 ≤ (x, y) ≤ 1 and λ ∈ [−0.5, 1]. Theorem 1. Let (X,Y ) be the bivariate random vector with joint pdf given in Equation (3), then the joint cdf of (X,Y ) is given by F (x, y) = 3y ( 2y2 − 3y + 1 ) M(x)− (2y − 3)y2F (x), (4) where M(x) = (1 + 2λ) 2 x2 − 4λx3 + 3λx4, and F (x) = (1 + 2λ)x− 6λx2 + 4λx3. I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 793 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 = ( 6y3 − 9y2 + 3y ) u− (2y3 − 3y2). Thus, F (x, y) = ∫ x 0 ( 1− λ+ 3λ(2u− 1)2 ) [( 6y3 − 9y2 + 3y ) u− (2y3 − 3y2) ] du =3y ( 2y2 − 3y + 1 ) M(x)− (2y − 3)y2F (x), where M(x) = (1 + 2λ)x2 2 − 4λx3 + 3λx4, and F (x) =(1 + 2λ)x− 6λx2 + 4λx3. This BeSU distribution can be used to model the lifetimes of two related electronic devices where the distributions of the lifetimes of two related devices follow the bathtub shape. An example of the electronic dataset that follows a bathtub shape is the dataset used by Rahman et al. [6] in their study and that data is about the lifetimes of 30 electronic devices. 3. Some properties of the Bivariate extended Standard U-quadratic distribution Theorem 2. Let (X,Y ) be a bivariate random vector with joint probability density func- tion given in Equation (3), then the marginal density function of Y follows an extended Standard U-quadratic (eSU) distribution with parameter λ∗ = 0.25. 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 I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 794 f(y) =1.5 [ 1− (2y − 1)2 ] ∫ 1 0 f(x)dx− 1.5 [ 1− 3(2y − 1)2 ] ∫ 1 0 xf(x)dx. (5) Observe that the first part of Equation (5) 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 xf(x)dx = ∫ 1 0 x [ (1− λ) + 3λ(2x− 1)2 ] dx = 1 2 . It follows that Equation (5) simplifies to f(y) =1− λ∗ + 3λ∗(2y − 1)2, where y ∈ [0, 1] and λ∗ = 0.25. Thus, the marginal distribution of Y is eSU distributed with parameter λ∗ = 0.25. Figure 1: joint pdf plot of the Bivariate extended Standard U-quadratic distribution for λ = −0.5. Figure 1 presents the joint pdf of the Bivariate extended Standard U-quadratic dis- tribution and it shows that this distribution can generate a combination of the inverted bathtub shape for X and bathtub shape for Y . Figure 2 presents the joint pdf of the Bivariate extended Standard U-quadratic distribution and it shows that this distribution can generate a combination of the constant shape for X and bathtub shape for Y . Figure 3 presents the joint pdf of the Bivariate extended Standard U-quadratic distribution and it shows that this distribution can generate a combination of a bathtub shape for X and a bathtub shape for Y . I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 795 Figure 2: joint pdf plot of the Bivariate extended Standard U-quadratic distribution for λ = 0. Theorem 3. Let X and Y be any two random variables with joint pdf given in Equa- tion (3). 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 ] 0.75 + 0.75(2y − 1)2 , (6) where λ ∈ [−0.5, 1]. The proof follows directly from the definition of the conditional distribution of X given Y = y. Theorem 4. Let X and Y be any two random variables with joint pdf given in Equation (3). 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 {[ 1− (2y − 1)2 ] E [Xr]− [ 1− 3(2y − 1)2 ] E [ Xr+1 ]} 0.75 + 0.75(2y − 1)2 , (7) where E [Xr] = (1 + 2λ)r2 + (5− 2λ)r + 6 (r + 1)(r + 2)(r + 3) , and E [ Xr+1 ] = (1 + 2λ)(r + 1)2 + (5− 2λ)(r + 1) + 6 (r + 2)(r + 3)(r + 4) . Proof. The rth conditional moment of X given Y = y is defined as E [Xr|y] = ∫ 1 0 xrf(X|Y = y)dx I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 796 Figure 3: joint pdf plot of the Bivariate extended Standard U-quadratic distribution for λ = 0.5. 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− (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 + 1)2 + (5− 2λ)(r + 1) + 6 (r + 2)(r + 3)(r + 4) ] = 1.5 f(y) {[ 1− (2y − 1)2 ] E [Xr]− [ 1− 3(2y − 1)2 ] E [ Xr+1 ]} , where E [Xr] = (1 + 2λ)r2 + (5− 2λ)r + 6 (r + 1)(r + 2)(r + 3) , and E [ Xr+1 ] = (1 + 2λ)(r + 1)2 + (5− 2λ)(r + 1) + 6 (r + 2)(r + 3)(r + 4) . Remark 1. The conditional mean of X given Y is given by E[X|y] = 1.5{0.5[1− (2y − 1)2]− ( 5+λ 15 ) [1− 3(2y − 1)2]} 0.75 + 0.75(2y − 1)2 . (8) I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 797 Remark 2. The conditional variance of X given Y is given by V ar(X|y) = E[X2|y]− (E[X|y])2 , (9) where E[X|y] is the conditional mean of X given Y , and E[X2|y] = 1.5{ ( 5+λ 15 ) [1− (2y − 1)2]− ( 5+2λ 20 ) [1− 3(2y − 1)2]} 0.75 + 0.75(2y − 1)2 . (10) Theorem 5. Let X and Y be any two random variables with joint pdf given in Equa- tion (3). If the conditional distribution of Y given X = x has an eSU distribution with parameter v(x) = 1.5x− 0.5, then the sth conditional moment of Y given X = x is E[Y s|X = x] = 6(s+ 1) + 3s(s− 1)x (s+ 1)(s+ 2)(s+ 3) . (11) Proof. The sth conditional moment of Y given X = x is defined as E[Y s|X = x] = ∫ 1 0 ysf(y|x)dy. Recall from [3] that the rth moment of an eSU distribution is given by E[Xr] = (1 + 2λ)r2 + (5− 2λ)r + 6 (r + 1)(r + 2)(r + 3) , (12) where λ ∈ [−0.5, 1]. Using Equation (12) for the sth moment of Y |X, we have E[Y s|X = x] = (1 + 2v(x))s2 + (5− 2v(x))s+ 6 (s+ 1)(s+ 2)(s+ 3) = 6(s+ 1) + 3s(s− 1)x (s+ 1)(s+ 2)(s+ 3) . Theorem 6. Let X and Y be any two random variables with joint pdf given in Equation (3), 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+1] ] , (13) and E [ XrY −s ] = 1 (1− s)(2− s)(3− s) [ 6(1− s)E[Xr] + 3s(s+ 1)E[Xr+1] ] , (14) where E[Xr] = (1 + 2λ)r2 + (5− 2λ)r + 6 (r + 1)(r + 2)(r + 3) , and E[Xr+1] = (1 + 2λ)(r + 1)2 + (5− 2λ)(r + 1) + 6 (r + 2)(r + 2)(r + 4) . I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 798 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+1f(x)dx ] = 1 (s+ 1)(s+ 2)(s+ 3) [ 6(s+ 1)E[Xr] + 3s(s− 1)E[Xr+1] ] , where E[Xr] = (1 + 2λ)r2 + (5− 2λ)r + 6 (r + 1)(r + 2)(r + 3) , and E[Xr+1] = (1 + 2λ)(r + 1)2 + (5− 2λ)(r + 1) + 6 (r + 2)(r + 3)(r + 4) . 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+1] ] . Remark 3. 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 4. 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 7. 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+1] (s+ 1) ] , (15) I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 799 where E[Xr] = (1 + 2λ)r2 + (5− 2λ)r + 6 (r + 1)(r + 2)(r + 3) , and E[Xr+1] = (1 + 2λ)(r + 1)2 + (5− 2λ)(r + 1) + 6 (r + 2)(r + 3)(r + 4) . 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 et1xet2y [ 1.5− 1.5x+ 3(1.5x− 0.5)(2y − 1)2 ] 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! ys [ 1.5− 1.5x+ 3(1.5x− 0.5)(2y − 1)2 ] f(x)dydx = ∫ 1 0 ∞∑ r1=0 tr1 r! xr {∫ 1 0 ∞∑ s=0 ts2 s! ys [ 1.5− 1.5x+ 3(1.5x− 0.5)(2y − 1)2 ] dy } f(x)dx = 3 (s+ 2)(s+ 3) ∞∑ r=0 tr1 r! ∞∑ s=0 ts2 s! [ 2E[Xr] + s(s− 1)E[Xr+1] (s+ 1) ] , where E[Xr] = (1 + 2λ)r2 + (5− 2λ)r + 6 (r + 1)(r + 2)(r + 3) , and E[Xr+1] = (1 + 2λ)(r + 1)2 + (5− 2λ)(r + 1) + 6 (r + 2)(r + 3)(r + 4) . Theorem 8. 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) = 1 70 ( 63 2 − λ ) , (16) where λ ∈ [−0.5, 1]. I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 800 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 =6E[X2]− 11E[X3] + 6E[X4] = 1 70 ( 63 2 − λ ) . Table 1: Some values of the stress - strength parameter for different values of λ. λ P(Y < X) λ P(Y < X) −0.5 0.4571429 0.3 0.4457143 −0.4 0.4557143 0.4 0.4442857 −0.3 0.4542857 0.5 0.4428571 −0.2 0.4528571 0.6 0.4414286 −0.1 0.4514286 0.7 0.4400000 0.0 0.4500000 0.8 0.4385714 0.1 0.4485714 0.9 0.4371429 0.2 0.4471429 1.0 0.4357143 From Table 1 we can see that as λ increases, the stress - strength seems to decrease. 4. Maximum Likelihood Estimation Let (X1, Y1),(X2, Y2),...,(Xn, Yn) be a random sample of size n from a Bivariate ex- tended Standard U-quadratic (BeSU) Distribution. Then the likelihood function is defined by L = n∏ i [ 1.5− 1.5xi + 3(1.5xi − 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.5xi + 3(1.5xi − 0.5)(2yi − 1)2 ] + n∑ i log [ 1− λ+ 3λ(2xi − 1)2 ] . The partial derivative of logL with respect to the parameter λ is given by ∂ logL ∂λ = n∑ i 3(2xi − 1)2 − 1 1− λ+ 3λ(2xi − 1)2 . I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 801 The maximum likelihood estimate of the parameter λ of the BeSU distribution is computed by solving n∑ i 3(2xi − 1)2 − 1 1− λ+ 3λ(2xi − 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. Consider the algorithm of Lakibul and Tubo [4] for generating random numbers from a T- extended Standard U-quadratic (TeSU) - G family of distributions. 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, (17) where λ ∈ [−0.5, 1] and G(x) is any baseline cumulative distribution function. The algo- rithm 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 ) . 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 I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 802 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 ) . 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: Step 1.** Draw a random sample X of size n from an extended Standard U-quadratic distribution with parameter λ. Step 2.** For each observation of X, draw a sample of size 1 from an extended Standard U-quadratic distribution with parameter 1.5x − 0.5. Repeat this process for all observations of X. Denote this sample as Y . Step 3.** Finally, the desired random sample is (x, y). I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 803 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: Step 1.*** Draw sample of size n, n = 50, 100, 200, 500, 1000, from a bivariate extended Standard U-quadratic distribution with parameter λ using the algorithm given in Section 5. Step 2.*** Using the bivariate sample (x, y) obtained in Step 1 above, compute the maximum likelihood estimate of λ. Step 3.*** Repeat the preceding Steps 1-2 N = 1000 times to get 1000 estimates of λ. Step 4.*** Compute the mean, bias, and mean squared error (MSE) of the 1000 esti- mates 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 . Table 2: Results of the Simulation Study for λ = −0.3, 0 and 0.5. λ = −0.3 λ = 0 λ = 0.5 n AE Bias MSE AE Bias MSE AE Bias MSE 50 −0.302 −0.002 0.016 −0.004 −0.004 0.025 0.494 −0.006 0.025 100 −0.297 0.003 0.008 0.002 0.002 0.011 0.498 −0.002 0.011 200 −0.299 0.001 0.004 0.001 0.001 0.006 0.498 −0.002 0.006 500 −0.300 0.000 0.002 −0.001 −0.001 0.002 0.498 −0.002 0.002 1000 −0.300 0.000 0.001 −0.001 −0.001 0.001 0.499 −0.001 0.001 Table 3: Results of the Simulation Study for λ = −0.5, 0.2 and 0.8. λ = −0.5 λ = 0.2 λ = 0.8 n AE Bias MSE AE Bias MSE AE Bias MSE 50 −0.471 0.029 0.004 0.195 −0.005 0.027 0.794 −0.006 0.015 100 −0.481 0.019 0.001 0.200 0.000 0.012 0.797 −0.003 0.007 200 −0.485 0.015 0.001 0.200 0.000 0.007 0.797 −0.003 0.004 500 −0.492 0.008 0.000 0.198 −0.002 0.003 0.797 −0.003 0.001 1000 −0.495 0.005 0.000 0.199 −0.001 0.001 0.799 −0.001 0.001 Considering nine different values of λ, Tables 2 to 4 show that as n becomes large, the average estimate (AE) of λ goes closer to the true value while the bias and the MSE I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 804 Table 4: Results of the Simulation Study for λ = −0.1, 0.3 and 1. λ = −0.1 λ = 0.3 λ = 1 n AE Bias MSE AE Bias MSE AE Bias MSE 50 −0.104 −0.004 0.023 0.294 −0.006 0.027 0.979 −0.021 0.002 100 −0.097 0.003 0.010 0.299 −0.001 0.012 0.987 −0.013 0.001 200 −0.099 0.001 0.006 0.299 −0.001 0.007 0.992 −0.008 0.000 500 −0.101 −0.001 0.002 0.298 −0.002 0.003 0.996 −0.004 0.000 1000 −0.100 0.000 0.001 0.299 −0.001 0.001 0.997 −0.003 0.000 diminish to zero. Thus, the maximum likelihood estimate of the proposed distribution parameter is consistent. 7. Application In this section, we derive a bivariate version of the Kumaraswamy distribution and use it to simulate a bivariate dataset. Since there is a limited to none existing bivariate distributions following a bathtub shape with support particularly on (0,1), we derive a bivariate extension of the Cubic Transmuted Uniform (CTU) distribution to compare it with the proposed Bivariate eSU distribution using the simulated data. Consider the Kumaraswamy (Km) distribution of Kumaraswamy [2], that is, for a given random variable K, the probability density function (pdf) of a Km distribution is given by f(k) = abka−1(1− ka)b−1, (18) with corresponding cumulative distribution function (cdf) F (k) = 1− (1− ka)b, (19) where k ∈ (0, 1), a > 0 and b > 0. It was observed that the Km distribution can generate shapes for a bathtub, constant, inverted bathtub, increasing and decreasing distribution. To generate a random sample from a Km distribution, we consider K = [ 1− (1− u) 1 b ] 1 a , (20) where u follows a uniform distribution defined on (0, 1). Let K1 be a random variable that follows a Km distribution with parameters a and b. Let K2 be another random variable such that the conditional distribution of K2 given K1 follows a Km distribution with parameters h(k1) > 0 and b > 0. From the work of Shahbaz et al. [10], the joint pdf of K1 and K2 is given by f(k1, k2) = ab2h(k1)k a−1 1 k h(k1)−1 2 [ (1− ka1)(1− k h(k1) 2 ) ]b−1 . I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 805 Further, we can set h(k1) = − log(1−k1) since 0 < − log(1−k1) < ∞. Thus, the joint pdf of K1 and K2 reduces to f(k1, k2) = ab2ka−1 1 k − log(1−k1)−1 2 [ (1− ka1)(1− k − log(1−k1) 2 ) ]b−1 log ( 1 1− k1 ) . (21) Equation (21) gives the joint pdf of the Bivariate Kumaraswamy (BKm) distribution. The following algorithm presents how to generate a random sample from a BKm dis- tribution: Step 1.∗∗∗∗ Draw sample of size n from the Kumaraswamy distribution with parameters a and b using Equation (20). Denote this sample as X∗. Step 2.∗∗∗∗ For each observation of X∗, draw a sample of size 1 from the Kumaraswamy distribution with parameters − log(1 − x) and b. Repeat this process for all obser- vations of X∗. Denote this sample as Y ∗. Step 3.∗∗∗∗ Finally, the desired random sample from the BKm distribution is (x∗, y∗). Here, we generate a bivariate set of data such that X and Y have bathtub shapes from the derived Bivariate Kumaraswamy distribution with a = 0.5 and b = 0.5. The following figures show the joint and marginal distributions of X and Y based on the simulated dataset. Figure 4: histogram plot of the simulated data for variable X. Next, consider the Cubic Transmuted Uniform (CTU) distribution of Rahman et al. [6], that is, for a given random variable C, the probability density function (pdf) of a CTU distribution is given by f(c) = 1− δ + 6δc− 6δc2, (22) with corresponding cumulative distribution function (cdf) F (x) = (1− δ)c+ 3δc2 − 2δc3, (23) I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 806 Figure 5: histogram plot of the simulated data for variable Y . where c ∈ [0, 1] and δ ∈ [−1, 1]. The CTU distribution can generate the following shapes: (i) bathtub for δ ∈ [−1, 0); (ii) constant for δ = 0; and inverted bathtub for δ ∈ (0, 1]. In the work of Rahman et al. [6], the CTU distribution was applied to model the lifetimes of 30 electronic devices where the lifetimes of the devices follow a bathtub shape. They compared the CTU distribution with the Beta, Kumaraswamy and skew uniform distri- butions using the said data. They found out that the CTU distribution provided better fit for the said dataset as compared with the said competing distributions. Let C1 be a random variable that follows a CTU distribution with parameter δ. Let C2 be another random variable such that the conditional distribution of C2 given C1 follows a CTU distribution with parameters Q(c1). Following the results of Shahbaz et al. [10], the joint pdf of C1 and C2 is given by f(c1, c2) = [ 1− δ + 6δc1 − 6δc21 ] [ 1−Q(c1) + 6Q(c1)c2 − 6Q(c1)c 2 2 ] . Further, we can set Q(c1) = 2c1 − 1 since −1 ≤ 2c1 − 1 ≤ 1. Thus, the joint pdf of C1 and C2 reduces to f(c1, c2) = [ 1− δ + 6δc1 − 6δc21 ] [ 1− (2c1 − 1) + 6(2c1 − 1)c2 − 6(2c1 − 1)c22 ] . (24) We name the joint pdf in (24) as the joint pdf of the Bivariate Cubic Transmuted Uniform (BCTU) distribution. In the analysis, we use the R-package ”bbmle” to compute the maximum likelihood estimates of the parameters of the proposed BeSU and BCTU distributions. In addition, the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) are used to assess and compare the performance of the proposed bivariate distributions. Table 5: Estimates and some diagnostic values of the fitted models for the simulated dataset. Distribution Estimate Std.Error −2logLik AIC BIC BeSU λ̂ = 0.489109 0.014249 2243.733 2245.733 2252.25 BCTU δ̂ = −0.978219 0.028498 3127.089 3129.089 3135.606 Table 5 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 I. A. Lakibul, D. L. Polestico, A. P. Supe / Eur. J. Pure Appl. Math, 17 (2) (2024), 790-809 807 Figure 6: Bivariate Histogram of the simulated data. values of the AIC and BIC than the BCTU distribution. Thus, the BeSU distribution provides a better fit for this simulated data than the BCTU distribution. 8. Conclusions and Recommendations In this communication, the bivariate extension of the extended Standard U-quadratic distribution called the bivariate extended Standard U-quadratic (BeSU) distribution has been derived. Some properties of the proposed BeSU distribution such as the marginal dis- tribution, conditional distribution, conditional moments, the product and ratio moments, Pearson correlation coefficient, joint moment generating function and stress - strength parameter were computed. Maximum likelihood estimation was implemented to estimate the parameters of the BeSU distribution and a simulation study was carried out to assess the behavior of the estimates of the parameters of this derived distribution. It was ob- served that the maximum likelihood estimate of the parameter of the BeSU distribution is consistent. Moreover, bivariate version of the Kumaraswamy (BKm) distribution has been derived. The new bivariate version of the Cubic Transmuted Uniform distribution is also obtained and compared with the proposed BeSU distribution using simulated data. It reveals that the proposed BeSU distribution provides a better fit for the dataset generated from the BKm distribution as compared with the derived bivariate Cubic Transmuted Uniform distribution. For future studies in this field, it is recommended to use another bi- variate family of distributions like the bivariate Marshall-Olkin family and Farlie - Gumbel - Morgenstern (FGM) family for generalizing the extended Standard U-quadratic distribu- REFERENCES 808 tion into another form of the bivariate extended Standard U-quadratic distributions and to compare it with the proposed BeSU distribution presented in this paper. In addition, it is also suggested to use the BeSU distribution to model the failure rates of two related components in a system or to model the lifetimes of two related electronic devices where the failure rates or lifetimes follow the bathtub shapes on the interval [0, 1]. 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 and the MSU-Iligan Institute of Technology. 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] P. Kumaraswamy. A generalized probability density function for double-bounded random processes. Journal of Hydrology, 46:79–88, 1980. [3] 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. [4] 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. [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. M. Rahman, B. Al-Zahrani, S. H. Shahbaz, and M. Q. Shahbaz. Cubic Transmuted Uniform Distribution: An Alternative to Beta and Kumaraswamy Distributions. Eu- ropean Journal of Pure and Applied Mathematics, 12:1106–1121, 2019. [7] M. Q. Shahbaz and S. Shahbaz. Order Statistics and Concomitants of Bivariate Pseudo- Rayleigh Distribution. World Applied Sciences Journal, 7:826–828, 2009. [8] M. Q. Shahbaz, S. Shahbaz, and A. Rafiq. A New Bivariate Gumbel distribution. Nonlinear Analysis Forum, 16:133–136, 2011. [9] S. Shahbaz and M. Ahmad. Concomitants of Order Statistics for Bivariate Pseudo- Weibull Distribution. World Applied Sciences Journal, 6:1409–1412, 2009. REFERENCES 809 [10] 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. [11] 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.