Inferences of a Mixture Bivariate Alpha Power Exponential Model with Engineering Application

: The univariate alpha power exponential (APE) distribution has several appealing charac-teristics. It behaves similarly to Weibull, Gamma, and generalized exponential distributions with two parameters. In this paper, we consider different bivariate mixture models starting with two independent univariate APE models, and, in the latter case, starting from two dependent APE models. Several useful structural properties of such a mixture model (under the assumption of two independent APE distribution) are discussed. Bivariate APE (BAPE), in short, modelled under the dependent set up are also discussed in the context of a copula-based construction. Inferential aspects under the classical and under the Bayesian paradigm are considered to estimate the model parameters, and a simulation study is conducted for this purpose. For illustrative purposes, a well-known motor data is re-analyzed to exhibit the ﬂexibility of the proposed bivariate mixture model.


Introduction
Univariate continuous alpha power exponential distribution with two parameters has received a considerable amount of attention among researchers working in the area of distribution theory. Since the introduction of the model, an appreciable amount of work has taken place on the development of bivariate APE distribution and its different generalizations; for more information, see Ref.
[1] and the references cited therein. In this paper, we discuss several different bivariate APE models with various mixing mechanisms. Noticeably, mixture distribution plays an important role in various data-analysis purposes; see [2]. Precisely, if the data are coming from different subpopulations, and their identifications are unknown, then the mixture distribution can be used quite effectively to analyze the data set. Unlike any other symmetric distribution (for example, normal) APE distribution is a skewed distribution. Consequently, if the subpopulations do not have symmetric distribution, the mixture of APE resulting in a bivariate APE distribution may be used to analyze these data. Furthermore, a multi-model distribution can be well approximated by a mixture distribution. Not much work has been carried out on a bivariate mixture of APE distribution. This is an effort towards that direction. We consider two different mixing strategies to construct the bivariate mixture APE (BMAPE type I). In both cases, we start with the following: let X ∼ APE(α 1 , θ 1 ) and Y ∼ APE(α 2 , θ 2 ), and let them be independent. In continuous probability model arising from two independent univariate APE distributions, which we call BMAPE, in short. In this paper, we consider various mechanisms of mixing univariate APE distributions to construct a bivariate mixture of APE models, which are listed below: (i) In the first case, we start with two independent and univariate APE with parameters (α 1 , α 2 , θ) and (µ 1 , µ 2 , λ) with α 1 = α 2 , µ 1 = µ 2 . Further, we assume that α 1 , α 2 , µ 1 , µ 2 are random variables. Then, using the methodology of compounding, we develop the bivariate APE distribution, and, for the sake of notational simplicity, call it the BMAPE (type-I) distribution. (ii) In the second scenario, we consider a copula-based construction approach. We consider a bivariate Gaussian copula in which the margins are two univariate and independent APE distributions with the coupling parameter ρ. For notational simplicity, we call this the BMAPE (type-II) distribution, for more details of copula based on construction and associated properties, see [25].
The remainder of the paper is laid out as follows. In Section 2, we describe the BMAPE (type-I) distribution and discuss its properties. In Section 3, we provide some useful properties of the BMAPE (type-I) distribution, including the survival function, hazard rate function and conditional moments. In Section 4, we look at how the expectation-maximization (EM) technique can be used to estimate model parameters for BMAPE (type-I) distribution. The estimation approach for the BMAPE (type-II) distribution formed using the bivariate Gaussian copula is discussed in Section 5. In Section 6, we investigate how to estimate model parameters under the Bayesian paradigm. In Section 7, the simulation results are reported. A well-known Motor data set was re-analyzed in Section 8 to demonstrate the efficacy of the suggested BMAPE type models. Finally, in Section 9, we wrap up the paper with some concluding remarks.

The APE Model with a Bivariate Mixture
We will start with two independent univariate APE distributions with parameters of (α, θ) and (µ, λ), respectively. We employ here the strategy of compounding by assuming that α and µ are random variables; then, the marginal joint distribution of X in and Y in may be derived from the joint distribution of α and µ respectively, which is given as follows h1 x,y (x, y) = h1(x, y, α, µ)dα dµ, where X in and Y in are independent APE distributions. Next, we consider two different scenarios in the context of the nature of the distributions of X in and Y in . From (1), one may obtain the desired BMAPE (type-I) distribution. We generate a bivariate mixture of the APE distribution in two scenarios. In the first scenario, X in and Y in are independent APE distributions with a generalized bivariate Bernoulli distribution for the shape parameters as µ 1 µ 2 P = α 1 α 2 P α 1 µ 1 P α 1 µ 2 P α 2 µ 1 P α 2 µ 2 .
A random variable with an APE distribution has a CDF and a PDF for X in > 0, given by, respectively, f (x in , α, θ) = θlog(α)e −θx in α 1−e −θx in α − 1 , x in > 0, where α and θ are the shape and scale parameters.
Next, in order to construct a finite mixture of bivariate alpha power exponential model, we proceed as follows: (i) We consider X i ∼ APE (α i , θ) ∀i = 1, 2, where α 1 and α 2 are random and θ is fixed, and each are independent for fixed choices of α 1 and α 2 . Likewise, Y i ∼ APE (µ i , λ) ∀i = 1, 2, where µ 1 and µ 2 are random and λ is fixed. (ii) Next, we further assume that (X 1 , X 2 ) and (Y 1 , Y 2 ) are independent but (α i , µ i ), ∀i = 1, 2, are correlated through their bivariate generalized Bernoulli distribution with the shape parameter (α i , µ i ), assuming the following probability matrix.
where P is the mixing component with the condition ∑ 2 i=1 P α i µ i = 1. Let h1(x, y) be the joint PDF created from the marginal densities of (x i , y i ), i = 1, 2, given earlier.
To make things easier in terms of writing and visualization, let a = P α 1 µ 1 , b = P α 1 µ 2 , c = P α 2 µ 1 and d = P α 2 µ 2 . Depending on a variety of choices for the mixing parameters a, b, c, d, θ, λ, α i and µ i for i = 1, 2, several shapes of the joint PDF in (5) are provided in Figure 1. Moreover, several contour plots for various values of a, b, c and d when (α 1 , α 2 , µ 1 , µ 2 , θ, λ) = (9, 7, 6, 5, 1.2, 0.6) are displayed in Figure 2. As one might expect, not all four components may be required in applications to real-world data sets. Consequently, some constraints can be applied, such as b = c = 0, a = d = 0 or a = b = c = 0, as appropriate. These constraints result in correlation values of +1, −1, and 0 among the shape parameter, respectively.     The marginal densities of X and Y, respectively, from the joint PDF in (5) are not available in closed form; one may consider some numerical methods to obtain an approximate expression. The CDF for the joint model will be For (6), to be a valid CDF, we must have Next, one may obtain the marginal survival functions for the smallest order statistics for a random sample of size two. For example, let w 1 = min{X, Y}, then, From (7), one can obtain the associated PDF for W 2 . We discuss some useful structural properties in the next section for the BMAPE (type-I) distribution.
The BMAPE (type-I) distribution's survival function will be The marginal distribution functions of X and Y are F X (x) and F Y (y), respectively. The following is the hazard rate function (hrf): For each fixed Y = y, the conditional PDF of X given Y as well as Y given X are given, respectively, by Then, (X,Y) has a positive association with a total positivity of order 2 (TP 2 ), provided (α 1 , α 2 , µ 1 , µ 2 ) > 1.
Observe that an absolute continuous bivariate random vector, such as (U 1 , U 2 ), has the TP 2 property if and only if for each u 11 , u 12 , u 21 , and u 22 , whenever u 11 < u 12 and u 21 > u 22 , is the joint PDF. Take distinct ordered u 11 , u 12 , u 21 and u 22 such that u 11 < u 12 and u 21 > u 22 , and our result follows instantly. Consequently, the positive quadrant dependence property (also known as the TP 2 property) would reveal a number of non-increasing features relating to conditional survival, the conditional CDF of X given Y, and Y given X, as well as the following finding: g 1 (.) and g 2 (.), Cov(g 1 (X), g 2 (Y) ) ≥ 0.

Proposition 2.
Shape of the distribution A critical point in a two-variable function is where the first-order partial derivatives are equal to zero. The following are the two most compelling reasons to investigate the critical points of a bivariate distribution: (1) A real-life data set (or sets) can assume a variety of shapes. A study such as this can help determine the adaptability of this bivariate distribution. (2) When working with bivariate distributions, it is essential to examine the joint PDF's tails as well as the point of inflection. This study is useful in the context of data fitting, as data sets having heavier tails and/or which are small can be adequately modeled by a distribution exhibiting one of these features, as appropriate. Critical-point(s) knowledge will aid in a better understanding of these features, as will be illustrated in Appendix A.

Inference
In this section, we consider the estimation of the model parameters of the BMAPE (type-I) distribution using the EM algorithm. The EM algorithm, originally proposed by [2], is an iterative approach for determining the maximum likelihood estimator of a parameter in a parametric probability distribution. To use the EM technique, we combine the data (x k , y k ), k = 1, . . . , n (to be obtained from a random sample of size n obtained from the associated joint density) with the group membership variables Φ k = (a k , b k , c k ), k = 1, . . . , n, where a k is an indicator variable that is one if the k th observation is in f (x, α 1 , µ 1 ) and zero otherwise. Similarly, we have four groups f ij , i, j = 1, 2, for b k , c k , with the following densities of where P( f 11 ) = a, P( f 12 ) = b, P( f 21 ) = c, and P( f 22 ) = 1 − a − b − c are the equivalent mixing proportions.
The associated log likelihood function for the complete sample ij (x, y) = log f ij (x, y), is symbolized by: The EM method has two steps: Step-E (expectation) and Step-M (maximization) in each iteration. As the group membership variables represented as Φ k are linear, we replace the related expected values in (11) with the current estimates (α 1 ,μ 1 ,α 2 ,μ 2 ,λ,θ,â,b,ĉ) of the parameters determined and, subsequently, an update occurs. For more information on this topic, see [26] as: Similarly, for b k and c k , apply the same technique, as follows: Following that, in the M-step, we maximize (11) across (λ, θ, α 1 , α 2 , µ 1 , µ 2 ) for fixed values of Φ k . The conditional independence of X and Y given the group membership achieves maximization by differentiating (11), which are: ∂ι and ∂ι ∂λ The M-step is then completed by settinĝ The model parameters must have a starting value to initiate the process, which is denoted by the symbol Φ (0) . To ensure that the proposed EM algorithm's rate of convergence does not become significantly sluggish, careful selection of these initial values is required. Another factor to consider is that the maximum likelihood equation may have several solutions corresponding to local maxima, necessitating careful selection of starting values. In the literature, some discussion on the identification of starting value(s) have been discussed and the references cited therein. The "coda" and "maxLik" of R package were used to quantitatively solve these equations. After obtaining the maximum likelihood estimators (MLEs) for α 1 , µ 1 , α 2 , µ 2 , θ, λ, δ, ρ 11 , ρ 12 , ρ 21 and ρ 22 , we replace these estimates in (a k , b k , c k ). Then, the M-step is completed by settingâ = n −1 ∑ n k=1 a k , and the process continues. To produce initial values for the mixing proportions individually, the matching moment's method for the marginal univariate APE parameter is employed. The BMAPE (type-I) parameter estimates are used as the starting points for the EM technique. Then, assuming that the two variables X and Y are independent, we combine the marginal mixing parameter moment estimators to obtain initial values for the bivariate mixing parameters. Specific details are provided in Section 7, and, more specifically, in Tables 2-4. By assuming that the two variables X and Y are independent, as in (3) and (4), we combine the marginal mixing parameter of moment estimators to obtain initial values for the bivariate mixing parameters. We equate (12)- (17) to zero to obtain the estimates of the BMAPE (type-I) distribution's parameter. Approximate confidence intervals (ACIs) of the unknown parameters are also obtained. This part is not investigated here but its results are presented numerically in Sections 7 and 8.

BMAPE (Type-II) Model Gaussian Copula Distribution
In this scenario, the proposed BMAPE model uses APE marginals and provides a reasonable flexibility in terms of marginals and the associated correlation structure. A copula is a useful tool for constructing bivariate and multivariate distributions, see [24,25] and the references cited therein. Copulas are extensively employed and play a vital part in the study of dependency. Several bivariate and multivariate lifetime distributions have been proposed by Refs. [27,28] on using a wide variety of couples. In this paper, we discuss the construction of a BMAPE (type-II) distribution via bivariate Gaussian copula with two APE marginals.

BMAPE (Type-II) Distribution Based on a Gaussian Copula
Any multivariate distribution may be decomposed into a copula and its continuous marginal, according to Ref. [29]. Copulas are useful tools in a bivariate scenario to combine two marginal distributions such that, for every bivariate distribution function F(x,y), with continuous marginal F(x) and F(y), one can write The density function associated with the bivariate distribution will be where is the density function of a copula; see [30]. As explained, there are a variety of options for building BMAPE (type-II) distributions using copulas and APE marginals (1). A bivariate Gaussian copula is defined as where Φ 2 is the joint CDF of a bivariate normal vector with zero means, unit variates, and correlation coefficient ρ.
Here, we denote the joint distribution of the two variables as X d and Y d where d symbolizes the phrase dependence. The joint PDF of X d and Y d based on Gaussian copula becomes where Z 1 and Z 2 are latent random variable as well as f (x d ) and f (y d ) being the density functions of APE distributions provided in (3) and (4), and ρ ∈ [−1, 1] the dependence parameter. The BAPE (type-II) distribution PDF is generated if the marginals are from an APE distribution.
and ( α, µ, θ, λ) > 0. Assuming that X d and Y d are dependent APE distributions with a generalized bivariate Bernoulli distribution for the shape parameters α i , µ i and, hence, the PDF of the BMAPE (type-II) distribution is defined as where the mixing proportions are m i , and it must satisfy ∑ 2 i=1 m i = 1, and m i ≥ 0, all of which are unknown. The first component of APE's density is provided by (3). The PDF of the first component of APE is given by where θ is a fixed-scale parameter, θ > 0, and α is a random-shape parameter, α > 0, which takes two distinct values of α 1 and α 2 . Let Y d have an APE mixing density and the PDF of the second component APE be determined by Y d for a fixed-scale parameter λ with the values µ 1 and µ 2 being a random-shape parameter. We suppose that X d and Y d are dependent for given values (α, µ) and that α and µ are associated by their generalized bivariate distribution with the following probability matrix provided by The joint PDF in (25) can take on a variety of shapes, similar to the joint PDF in (4). For simplicity, let x d = x, y d = y, a = m α 1 µ 1 , b = m α 1 µ 2 , c = m α 2 µ 1 and m α 2 µ 2 ; as a result, the BMAPE (type-II) distribution PDF will be

Gaussian Copula and the EM Algorithm
Here again, the EM algorithm is introduced. To use the EM technique, we supplement the data (x k , y k ); k = 1, . . . , n, with the group membership variables (a dk , b dk , c dk ), k = 1, . . . , n, where a dk is one if the k th observation is in f ij (x, y, α 1 , µ 1 , θ, λ, ρ 11 ), and zero otherwise. Similarly, we have four groups f ij , i, j = 1, 2 for b dk , c dk , with densities as where and then the EM algorithm as a technique of estimate by finding the complete log likelihood, say d , is given as follows: As the group membership variables (a dk , b dk , c dk ) are linear, we enter their predicted values into (28) in the E-step, given the current estimates (α 1 ,μ 1 ,α 2 ,μ 2 ,θ,θ,λ,λ,ρ 11 ,ρ 12 , ρ 21 ,ρ 22 ,â d ,b d ,ĉ d ) of the parameter, computed as: Similarly, we use the same approach for b dk and c dk . It is worth noting that the above algebraic simplification may be required to obtain numerical solution. For fixed values of (a dk , b dk , c dk ), we maximize (28) over (α 1 , µ 1 , α 2 , µ 2 , θ, λ, ρ 11 , ρ 12 , ρ 21 , ρ 22 ). The univariates and the Gaussian copula parameter can be dealt with separately.
The method of matching moments, which is generated from the marginal univariate APE and the Gaussian copula parameter separately, are used to obtain the initial values of the parameters for the mixing proportions. The BAPE parameter estimators obtained are used as starting values for the EM method. Then, assuming that the dependency between two variables X and Y is zero, we combine the moment estimators of the marginal mixing parameters to generate initial values for the bivariate mixing parameters. Ref. [31] for more information. Following that, we present the process for estimating the unknown parameters for the density in (28). We use two approaches in copula-based estimation: parametric and semiparametric.

Maximum Likelihood Estimation
We look at how to estimate the unknown parameters of BMAPE (type-II) distributions using the maximum likelihood technique and two-step estimation. This entails a two-step process in which the marginal and copula functions are individually estimated; see [32]. The log-likelihood function is defined as follows: The log-likelihood function in (40) may be rewritten as The MLEs are separately used to estimate the parameters of the marginal distribution F 1 and F 2 , as follows: Then, by maximizing the copula density, the copula parameters may be estimated.
The maximum likelihood estimation will estimate the parameters of each marginal distribution by evaluating the first step with APE distributions. If ( x 1 , . . . , x n ) is a random sample from APE(α, θ) and ( y 1 , . . . , y n ) is a random sample from APE(µ, λ), the loglikelihood functions are provided by and As a result, the maximum likelihood equations are as follows: ∂l ∂α = n αlog(α) and ∂l ∂µ = n µlog(µ) The MLEs of α, µ, θ and λ are obtained by solving the nonlinear (46)-(49) system. The density of copulas is then calculated as follows: whereF 1 (x) and ,F 2 (y) are the maximum likelihood estimates of the first-step parameters. Therefore, the maximum likelihood estimate of γ can be obtained by solving the nonlinear (50).

Estimation via Semiparametric Methods
The inversion Kendall's and inversion of Spearman's semiparametric methods for estimating the copula parameter in copula models are compared to the two methods of moments approaches, inversion Kendall's τ and inversion of Spearman's ρ, respectively.

Moments Method
We present a brief description of moment's method of inversion of Kendall's τ and the inversion of Spearman's ρ as indicated in Ref. [33]. Let c be a bivariate random sample from a CDF, C γ [F 1 (x ), F 2 (y)], where F 1 and F 2 are continuous CDF and C γ is an absolutely continuous copula such that γ ∈ O. Moreover, unless otherwise specified, let R 1 , . . . , R n be the vectors of ranks associated with (x 1 , . . . , x n ). All vectors in the following are row vectors. The inversion of a consistent estimator of a moment of the copula C γ is used in moments methods. The two most well-known moments, Spearman's ρ and Kendall's τ, are given by Spearman and Kendall, respectively, as and These two moments' consistent estimators can be stated as and Consistent estimators of ρ and τ will be if and are one-to-one γ n,ρ = ρ −1 (ρ n ), γ n,τ = τ −1 (τ n ), respectively.
Inversion of Kendall's τ and inversion of Spearman's ρ are two terms for the same thing. Ref. [33] and the references referenced therein for more information. As previously stated, the moment's method of τ and ρ estimation for copula can be classified as semiparametric approach estimation.

Copula Fit Tests
We want to compare the empirical copula with the parametric estimator generated under the null hypothesis; for more information, see [34,35]. According to the theory, a test to see if C is well-represented by a given copula C γ can be represented as There are several well-known ways in the literature, such as Ref. [36] or the rapid multiplier approach, as described by Refs. [37,38]. The empirically based goodness of fit test is written as where C n (u, v) is the empirical copula of the X and Y data.
, the pseudo-observations U i,n and V i,n from C, are calculated as: U i,n = R 1i n+1 , V i,n = R 2i n+1 ; R 1i and R 2i are, respectively, the ranks of X i , Y i . C n (u, v) is an estimator obtained using the pseudo observations and C n (u, v) is a consistent estimator. The relevant test statistic, is the Cramer-von Miss, which is defined as [39] as well as the references within.

Bayesian Estimation
The Bayes estimates of the joint posterior density (54) model parameters are computed in this part under the assumption that the random variable Φ = (α 1 , α 2 , µ 1 , µ 2 , θ, λ) have independent gamma prior distributions with hyperparameter w k and m k and k = 1, 2, 3, 4, 5, 6, provided by and ρ 11 , ρ 12 , ρ 21 , ρ 22 theses have a non-informative prior. The joint posterior density for the vector Φ given the data is obtained by multiplying (9) or (28) as Integrating out the (nuisance) hyperparameters yields marginal distributions of Φ. As a result, under the square error loss function, the Bayesian estimators of the parameters can be derived as follows:

Bayes MCMC Estimates
As the integrals in (56) are not in a closed form, the Markov-chain Monte-Carlo (MCMC) method is used. The posterior distribution and intractable integrals are derived using simulated samples from the posterior distribution in MCMC algorithms. Gibbs sampling and the Metropolis-Hastings (M-H) algorithm are also employed as MCMC techniques, see [40,41].
The M-H algorithm assumes that a candidate value can be generated from a proposed distribution for each iteration of the process. As a result, the applicant value is allowed if the approval chance is high enough. This method ensures that the Markov chain for the target density will converge. Finally, we can see that the advantage of the MCMC technique over the MLE method is that by creating probability intervals based on empirical posterior distribution, we can always obtain an acceptable interval estimate of the parameters. In MLE, this is frequently unavailable.

Bayes Credible Intervals
The following expression can be used to obtain a symmetric 100(1 − ε)% percent two-sided Bayes probability interval estimate of Φ, denoted by [L Φ , We need appropriate numerical approaches to solve this non-linear equation because finding the interval L Φ and U Φ analytically is difficult.

Monte-Carlo Simulation
To assess the adaptability and flexibility of the proposed MBAPE model for different combinations of the true parameter value and sample size, an extensive Monte-Carlo simulation is conducted based on both independent and dependent cases. The point (MLE and MCMC) estimates and interval (ACI and BCI) estimates are evaluated based on 2,000 random samples of sizes n(= 50, 100, 150, 200). Extensive computations were performed using R statistical software via two useful statistical packages, namely, (i) the "coda" package proposed by Ref. [42], and (ii) the "maxLik" package, which uses the Newton-Raphson method of maximization in the likelihood computations, proposed by Ref. [43].
Comparison between point estimators is made in terms of their mean absolute bias (MAB) and root mean squared-error (RMSE) values. In addition, the behavior of 95% asymptotic/credible intervals are compared by their average confidence lengths (ACLs). In Table 1, different actual values of the model parameters and the associated hyperparameters are reported. Clearly, the values of hyperparameters are specified in such a way that the prior mean becomes the expected value of the target parameter. In this study, for both sets III and IV, we take (ρ 11 , ρ 12 , ρ 21 , ρ 22 ) as (0.2,0.15,0.05,0.1) and (0.5,0.4,0.2,0.3) as well as also assuming they have improper gamma priors. The MCMC estimates under squarederror loss are computed by running the M-H algorithm 12,000 times and ignoring the first 2,000 MCMC iterations. The calculated MLE of each unknown parameter is taken as an initial guess to obtain the corresponding Bayes estimate; for more information in this topic see [44][45][46]. The simulation outcomes are provided in Tables 2-6. In these tables, the Bayes MCMC results under Prior 1 (as an example) are denoted as "MCMC.1"; similarly, the BCI estimates under Prior 1 are also denoted as "BCI.1", for brevity. From Tables 2-6, it can be seen that the proposed estimates of the unknown MBAPE parameters are sufficiently efficient with respect to their minimum simulated values of RMSEs, MABs and ACLs. For all true parametric combinations, when n increases under independent and dependent assumptions, the RMSEs, MABs and ACLs of the classical and Bayesian estimates decrease, as expected. Therefore, one can maximize the efficiency associated with all estimates by increasing the overall sample size. As α, θ, µ, λ and ρ increase, it is observed that their RMSEs, MABs and ACLs increase. A similar pattern is observed in the case of α 1 , α 2 , µ 1 , µ 2 and ρ ij for i, j = 1, 2 increases, a dependent case. Since the Bayesian estimates include the flexible gamma information, the Bayesian estimates of all unknown parameters have overall lower RMSE and MAB values, so they performed better compared to the classical estimates in both independent and dependent cases. Consequently, the credible intervals estimates are also better than the asymptotic confidence intervals in terms of shortest ACLs. In particular, as the calculated variance in Prior 2 is less than Prior 1, it is also seen that the Bayes (point and interval) estimates have performed better based on Prior 2 than those obtained based on Prior 1. All accuracy criteria developed in this work support this statement.

Motor Data Analysis
To demonstrate the practical utility of the proposed methodologies in an engineering phenomenon, one real data set reported by Ref. [47] is analyzed. This data represents the failure times (in days) of a parallel system constituted by two identical motors, namely, (X) and (Y), see Table 7. We first check the suitability of the APE model to the complete motor data sets. For this purpose, for each X and Y data point, the Kolmogorov-Smirnov (KS) statistics and its p-value, as well as the MLEs with their standard errors (SEs), are obtained and also provided in Table 7. It shows that the APE distribution fits the motor data sets quite well. Moreover, some graphical methods for goodness of fit for the motor data sets are developed. Plots of fitted density, empirical cumulative distribution, probability-probability (PP) and quantilequantile (QQ) for X and Y data sets are displayed in Figures 3 and 4, respectively. It is clear that the two variables X and Y are fitted for the marginal APE distribution.  Using the complete X and Y data sets, the MLEs and 95% ACIs of α, θ, µ, λ and ρ are calculated. Furthermore, since no prior information is available about the MBAPE parameters, the Bayesian (point/credible) estimates of α, θ, µ, λ and ρ are developed using the non-informative priors. Taking the maximum likelihood estimates of the same unknown quantities as initial values, the MCMC sampler is replicated 30,000 times and then the first 5,000 iterations removed as a burn-in. The point (with their SEs) estimates and interval (with their lengths) estimates of the unknown parameters based on X and Y data sets are computed and listed in Tables 8 and 9, respectively. This indicates that the estimates of the unknown parameters obtained by maximum likelihood and Bayesian methods are very similar. Moreover, trace plots (in Figure 5) for MCMC draws of the posterior generated samples of α, θ, µ, λ and ρ showed that the MCMC simulated variates converge very well. It represents the sample mean and two bounds of 95% BCIs with solid (-) and dashed (---) lines, respectively. Moreover, the MCMC frequency estimates of α, θ, µ, λ and ρ using the Gaussian kernel are represented in Figure 6. It represents the sample mean with a vertical dash-dotted (:) line. It is evident that the generated posteriors of all unknown parameters are quite symmetrical. To sum up, the numerical findings of the suggested estimation methods based on motor data provide a good demonstration of the proposed model.

Conclusions
In this paper, we propose and study a new class of bivariate alpha power distributions via finite mixtures involving APE distributions as marginals, and also consider the copula-based construction of bivariate alpha power distributions. The developed class of distribution is made up of two different types of mixture, namely, type I, which starts with two independent APE distributions, and type II, which starts with a bivariate Gaussian copula. The model parameters are estimated using both frequentist (for example, method of moments and the method of maximum likelihood) and under the Bayesian (using

Conclusions
In this paper, we propose and study a new class of bivariate alpha power dis tions via finite mixtures involving APE distributions as marginals, and also consid copula-based construction of bivariate alpha power distributions. The developed cl distribution is made up of two different types of mixture, namely, type I, which start two independent APE distributions, and type II, which starts with a bivariate Gau copula. The model parameters are estimated using both frequentist (for example, m of moments and the method of maximum likelihood) and under the Bayesian (

Conclusions
In this paper, we propose and study a new class of bivariate alpha power distributions via finite mixtures involving APE distributions as marginals, and also consider the copula-based construction of bivariate alpha power distributions. The developed class of distribution is made up of two different types of mixture, namely, type I, which starts with two independent APE distributions, and type II, which starts with a bivariate Gaussian copula. The model parameters are estimated using both frequentist (for example, method of moments and the method of maximum likelihood) and under the Bayesian (using independent gamma priors) paradigm. This distribution can be employed in practice for non-negative and positively correlated random variables because the joint distribution function and joint density function are in closed forms; in addition, it can be applied to model censored data as well. Under the Bayesian paradigm, we adopted the EM approach, which works quite well and may be efficiently employed to compute the MLEs, because the MLEs of the unknown parameters cannot be computed in closed form. To compare the behavior/efficiency of the proposed methods and to see the applicability of the different estimators, a small simulation study is considered, and one real data set is investigated to examine the applicability of the proposed model. From the simulation study, it appears that the point and interval estimates of all unknown parameters using the MCMC algorithm are performing better. Therefore, the Bayes M-H method used to estimate the unknown MBAPE parameters is recommended. However, whether this is a general trend or not is subject to a separate study with a wide range of possible values of the model parameters.

Data Availability Statement:
The data used to support the findings of this study are included within the article.