A Class of Bivariate Semiparametric Families of Distributions

The study of semiparametric families is useful because it provides methods of extending families for adding flexibility in fitting data. The main aim of this paper is to introduce a class of bivariate semiparametric families of distributions. One especial bivariate family of the introduced semiparametric families is discussed in details with its sub-models and different properties. In most of the cases the joint probability distribution, joint distribution and joint hazard functions can be expressed in compact forms. The maximum likelihood and Bayesian estimation are considered for the vector of the unknown parameters. For illustrative purposes a data set has been re-analyzed and the performances are quite satisfactory. A simulation study is performed to see the performances of the estimators.


Introduction
To mathematically describe any family of distributions, various alternative functions are in common use. These functions include distribution functions, survival functions, densities; hazard functions, reversed hazard functions, cumulative hazard and cumulative reversed hazard functions. When they exist, any of these functions can be obtained from any other.
The distribution function F (·) and the survival function S(·) defined on (−∞, ∞) as where f (x) is the probability density function for a continuous random variable. The hazard function h(·) and cumulative hazard function H(·) are defined on (−∞, ∞) respectively, as The reversed hazard function r(·) and cumulative reversed hazard R(·) function are defined on (−∞, ∞) respectively, as R (x) = log F (x) .
it is possible to find a three parameter family that includes both the gamma and Weibull families as special cases . So, it provides methods of extending families for added flexibility in fitting data. In this section some semiparametric families that introduced by Olkin (2007) will be discussed in one dimension and its bivariate extension will be obtained in the next section.
It is important to see that the main criteria for the semiparametric families is that the underlying distribution is a member of the parametric family. And the second criteria, is that once the semiparametric family is used to add a parameter, its reuse may reparameterize the family, but it should fail to again add a new parameter. this is a kind of stability property (Olkin, 2007, p.609).

Univariate Power Parameter Family (UPP)
Let F B be a baseline cdf. Suppose that that F UPP (·; α) is defined in terms of F B by the formula Then α is called a power parameter and {F UPP (·; α) α > 0} is an univariate power parameter family with underling distribution F B . The corresponding probability density function (pdf) and hazard function are, respectively, where f B and h B are a baseline pdf and hazard functions respectively. The family of Weibull distributions is the prime example of a Power parameter family. For this example the underlying distribution is an exponential distribution. Also Ghitany et al. (2013) introduced a power Lindley distribution, the underlying distribution is a Lindley distribution.

Univariate Frailty Parameter Family (UFP)
Let S B be a baseline survival function with cumulative hazard function H B = − log S B . Suppose that S UFP (·; α) is defined in terms of S B by In this case α is called a frailty parameter and S UFP (·; α), α > 0, is a frailty parameter family, or alternatively, a proportional hazard family with underlying distribution S B . The corresponding pdf and hazard function is given respectively, as where f B and h B are a baseline pdf and hazard functions respectively. For exponential and Weibull distributions, introducing powers of the survival function does not introduce a new parameter because these families are already proportional hazards families. For a number of other families, however, a new parameter is introduced such as Pareto type I (Olkin, 2007), exponentiated Fréchet (Nadarajah and Kotz, 2003), exponentiated Gumble (Nadarajah, 2006), and extended Lindley Distributions (Bakouch et al., 2012).

Univariate Resilience Parameter Family (URP)
Let F B be a baseline distribution function with cumulative reversed hazard function In this case α is called a resilience parameter and F URP (·; α) , α > 0 is a resilience parameter family, or alternatively, a proportional reversed hazard family with underlying distribution F B .
The corresponding pdf and hazard function is given respectively, as where f B and r B are a baseline pdf and reversed hazard functions respectively. There exist number of distributions that produced by adding a resilience power parameter to some failure functions such as generalized exponential distribution (Gupta and Kundu, 1999), Generalized linear failure rate distribution (Sarhan and Kundu, 2009), exponetiated gamma (Nadarajah and Kotz, 2006), exponetiated Kumaraswamy (Lemonte et al., 2013), exponetiated Weibull (Mudholkar and Srivastava, 1993), exponetiated inverted Weibull (Flaih et al., 2012), and so on.
Resilience and frailty parameter families have the stability property, i.e., once a resilience or frailty parameter has been introduced; the reintroducing of the same kind of parameter does not extend the family.

Univariate Hazard Power Parameter Family (UHPP)
According to Equation (7) the survival function and its corresponding cumulative hazard function are related via the formula S (x) = e −H(x) , ∀x.
It follows that if H is a cumulative hazard function, then H α is accumulative hazard function. Thus, Defines a survival function for all α > 0, and S UHPPF (x; α), α > 0 is a semiparametric family. The parameter α is called hazard powerparameter family.The corresponding pdf is given by differentiating (13) as where h B (·) and H B (·) are the baseline hazard and cumulative hazard functions respectively. Accordingly the hazard function for UHPP family is given as It is follows if h B increasing and α ≥ 1, then h UHPPF is increasing; if h B decreasing and 0 < α < 1, then h UHPPF is decreasing.

Examples: 1. Univariate Exponential Distribution
Suppose that the underlying distribution is univariate exponential (UE) with scale parameter λ with survival function S B (x) = e −λx and H B (x) = λx. Then according to (13)- (14) and after adding a hazard power parameter α the survival, density and hazard functions respectively, are For x > 0 and α > 0. (15) is the survival function of an univariate Weibull (UW) distribution with hazard power parameter α and scale parameter λ. it is noted that α can be regarded not only as a hazard power parameter but also as a power parameter as mention above in Section 2.1.

Univariate Gompertz Distribution
Suppose that the underlying distribution is univariate Gompertz (UG) distribution with survival function S B (x) = e −ξ(e λx −1) , ξ, λ > 0, x > 0, according to (13) the survival function of the UG distribution with hazard power parameter α is given as and the corresponding pdf and hazard function are given respectively, as Its observed by (Olkin, 2007) that the hazard function (19) is convex. It is increasing when α ≥ 1 and when α < 1 the hazard function has minimum at x = [− log α]/λ.

Univariate Pareto Type I Distribution
Suppose the underlying distribution is univariate Pareto (UP) distribution with survival function S B (x) = [1 + λx] −1 , λ > 0, x > 0. After adding the hazard power parameter α, the survival function for the new distribution is given as Accordingly, It is noted from (22) that hazard function is decreasing for α ≤ 1. For α > 1 the hazard rate is unimodal with mode at e α−1 /λ.

Univariate Uniform Distribution
Suppose the underlying distribution is univariate uniform (UU) distribution with S B (x) = 1 − x and H B (x) = − log(1 − x). After adding the hazard power parameter α then, the survival function for the new distribution is given as Consequently,

Univariate Reversed Hazard Power Parameter Family (URPP)
It's possible to define a reversed hazard power parameter models via the formula F (x) = e R(x) , as following where α > 0 and R (x) = log F is a cumulative reversed hazard function.
where r B (·) and R B (·) are the baseline reversed hazard and cumulative reversed hazard functions respectively. Accordingly the reversed hazard function for URPP family is given as

Bivariate Hazard Power Parameter (BHPP) Family
Assume the univariate hazard power parameter model is denoted by UHPP(α, Θ) where α is the hazard power parameter and Θ may be a vector of parameters for an underlying distribution. Now suppose that U i ∼ UHPP (α i , Θ), i = 1, 2, 3 such that U i s are mutually independent random variables and define X j = min (U j , U 3 ), j = 1, 2. Such that X j s are dependent random variables. Hence the joint survival function of the vector (X 1 , X 2 ) denoted by S BHPP (x 1 , x 2 ) is given as where The joint survival function of BHPP model can be stretching in the following form where Accordingly, the joint pdf of BHPP model can be obtained as where The joint distribution function of (X 1 , X 2 ) is given by The joint hazard function of the dependent variables (X 1 , X 2 ) is obtained as follows where The marginal survival and densities of X 1 and X 2 is given respectively, as follows Further, for the BHPP family the conditional density of X 1i given X 2j = x 2j is given by
The log-likelihood function of the sample of size n from BHPP (α 1 , α 2 , α 3 ) is given by l (α) = n 1 log α 1 + n 2 log α 2 + n 3 log α 3 Accordingly, the likelihood equations can be written as The second derivatives are given as follows The asymptotic variance-covariance matrix ofα 1 ,α 2 andα 3 is obtained by inverting the Fisher information matrix with elements that are negatives of expected values of the second order derivatives of logarithms of the likelihood function. In the present situation, it seems appropriate to approximate the expected values by their maximum likelihood estimates. Accordingly; the asymptotic variance-covariance matrix can be written as follows Now, the asymptotic normality results will be stated to obtain the asymptotic confidence intervals ofα 1 , α 2 and α 3 . Under particular regularity conditions it can be stated as follow where F −1 is the variance-covariance matrix,α = (α 1 ,α 2 ,α 3 ) and α = (α 1 , α 2 , α 3 ). Since α is unknown, then F −1 (α) is estimated by F −1 (α); the asymptotic variance-covariance matrix that defined above and this can be used to obtain the asymptotic confidence intervals of α 1 , α 2 and α 3 .

Bayesian Estimation
In this section, the Bayesian estimation for the BHPP models parameters is considered under the assumption that the random variables α 1 , α 2 and α 3 are independently distributed with gamma prior distributions as follows The joint prior density of α 1 , α 2 and α 3 can be written as

Hyper-Parameters Determination
The hyper-parameters involved in priors can be easily evaluated, if the prior mean and variance are considered to be known. The prior mean and prior variance will be obtained from the maximum likelihood estimates of (α 1 , α 2 , α 3 ) by equating the mean and variance of (α j 1 ,α j 2 ,α j 3 ) with the mean and variance of the considered priors (gamma prior), where j = 1, 2, . . . , k and k is the number of random samples generated from the model. Thus, on equating the mean and variance of (α j 1 ,α j 2 ,α j 3 ) with the mean and variance of gamma priors, gets Now on solving the above two equations, the estimated hyper-parameters can be written as Similar procedure for determining the hyper parameters (a 2 , b 2 , a 3 , b 3 ) can be used for α 2 and α 3 . Since hence, the joint posterior density function of θ = (α 1 , α 2 , α 3 ) given the data D, denoted by π 1 (θ | D) can be written as where D = {(x 11 , x 21 ) , (x 12 , x 22 ) , . . . , (x 1n , x 2n ) } and L(D | θ) is the likelihood function. Therefore, the Bayes estimates of the unknown parameters θ = (α 1 , α 2 , α 3 ) under square error loss function (SEL) can be calculated through the following equations as follows Obviously, the 3 integrals given by (29) cannot be obtained in a closed form. In this case, the MCMC technique to generate samples from the posterior distributions is used and then compute the Bayes estimators for the individual parameters. A wide variety of MCMC schemes are available, and it can be difficult to choose among them. An important sub-class of MCMC methods is Gibbs sampling and more general Metropolis within Gibbs samplers. The advantage of using the MCMC method over the MLE method is that we can always obtain a reasonable interval estimate of the parameters by constructing the probability intervals based on empirical posterior distribution. To generate samples from the proposed family, the Metropolis-Hastings (M-H) method (Metropolis et al. (1953) with normal proposal distribution) is used. Thus, the following steps of M-H algorithm are performed to draw samples from the posterior density and in turn compute the Bayes estimates (BEs) of θ = (α 1 , α 2 , α 3 ).

The approximate BEs of θ
, i = 1, . . . , N with respect to squared error loss (SEL) function is given byθ whereθ is BEs under SEL and M is the burn-in-period (that is, a number of iterations before the stationary distribution is achieved).

Bivariate Weibull Distribution
Using Equations (15)- (17) in Equations (26)-(28), a new bivariate Weibull distribution denoted by BW(α 1 , α 2 , α 3 , λ) can be defined by the joint survival function The joint survival function of BW model can be stretching in the following form Accordingly, the joint pdf of BW model can be obtained as where

Application of BHPP Models to Real Data Set
To see how the BHPP models work in practices, one data set will be reanalyzed in this section.
The data set has been obtained from Meintanis (2007). The data represent the football (soccer) data where at least one goal scored by the home team and at least one goal scored directly from a penalty kick, foul kick or any other direct kick (all of them together will be called as kick goal) by any team have been considered. Here X 1 represents the time in minutes of the first kick goal scored by any team and X 2 represents the first goal of any type scored by the home team. In this case all possibilities are open, for example X 1 < X 2 or X 1 > X 2 or X 1 = X 2 = X. These data were analyzed by Meintanis (2007), who considered the Marshall-Olkin bivariate exponential distribution, and by many authors such as Kundu and Dey (2009), Gupta (2009), Muhammed (2016), Muhammed (2017), Muhammed (2019), and Muhammed (2020). Here, these data will be fitted to three BHPP models namely: (i) bivariate Weibull (BW), (ii) bivariate hazard power parameter Gompertz (BHPG) and (iii) bivariate hazard power parameter Pareto (BHPP) distributions. Note that both BW and BHPP are four parameters models but BHPG is a five parameter model.The main aim is to see, how the different BHPP models and the MLE works in practice.
Before trying to analyze the data using the BHPP models, first fit the marginals (UHPP) models to X 1 , X 2 separately. The UHPP models are (i) univariate Weibull(UW), (ii) univariate hazard power parameter Gompertz (UHPG) and (iii) univariate hazard power parameter Pareto (UHPP). Table 1 shows the MLEs, the Kolmogorov-Smirnov ((KS) distances between the fitted distribution and the empirical distribution function for X 1 and X 2 with correspondence p-value, the Akaike information criterion (AIC), Bayesian information criterion (BIC), the consistent Akaike information criterion (CAIC) and Hannan-Quinn information criterion (HQIC). That gives an indication that the BHPP models may be used to analyze this data set. Now, the data will fitted to the three BHPP models defined above, the MLEs and the standard error (SE) will be calculated for each bivariate model. To compare these model with each other or with any other bivariate models that represent this data the AIC, BIC, HQIC and CAIC are calculated for the three PHPP models as shown in Table 2. The BW model provides a better fit than the other tested models, because it has the smallest value among AIC, BIC, HQIC and CAIC and standard error.

Simulation Study
In this section, the results of a Monte Carlo simulation study testing the performance of MLE and Bayesian estimation for the BHPP model Parameters will be introduced in general and   (30) and denoted by BW(α 1 , α 2 , α 3 , λ) and belongs to the BHPP models The evaluation of the MLE and the Bayes estimation was performed based on the following quantities for each sample size: the Average Estimates (AE), the Mean Squared Error (MSE), Bias and confidence interval length (CL) are estimated from R = 10000 replications forα 1 ,α 2 ,α 3 andλ the sample size has been considered at n = 20, 30, 40, 50, 70 and 100, and some values for the parameters α 1 , α 2 , α 3 and λ have been considered.
The algorithm to generate from BHPP Models goes as follows.
Step 4. Define the indicator functions as Step 5. The corresponding sample size n must satisfy n = n 1 + n 2 + n 3 such that For different choices of sample sizes, 10000 data sets were generated. The results are summarized in Table 3 and 4. The estimates are work well and MSE and RAB decreases as the sample Assuming that U 1 , U 2 and U 3 are mutually independent random variables such that U 1 ∼UPP(α 1 ), U 2 ∼UPP(α 2 ) and U 3 ∼UPP(α 3 ). Define X 1 = max(U 1 , U 3 ) and X 2 = max(U 2 , U 3 ) then by using Equations (8) and (9), the bivariate Power parameter family of distributions denoted by BPP(α 1 , α 2 , α 3 ) is defined by the joint cdf as follows where x 3 = min(x 1 , x 2 ).
Then, (X 1 , X 2 ) constitute a BRPP class of distribution with the following cdf and pdf where x 3 = min(x 1 , x 2 ).
The joint cdf of BHPP model can be stretching in the following form And it will be considered in details in a future work.

Conclusion
In this paper, a review of some univariate semi parametric families of distributions such as power parameter, frailty parameter, resilience parameter, hazard power parameter and reversed hazard power parameter is discussed in details. Moreover, proposed bivariate extensions for these families are introduced based on Marshal-Olkin idea. The bivariate hazard power parameter family is discussed with its main properties and the MLE and Bayesian estimation are also considered for the shape parameters. A simulation study and a real data set are considered. As a future work, other bivariate extensions for these families will be introduced based on other copulas soon as possible.