A new flexible bivariate discrete Rayleigh distribution generated by the Marshall-Olkin family

Abstract

In the statistical analysis of bivariate data, it is possible to have discrete observations instead of continuous data, as observed in many studies with survival data. In this study it is introduced a new bivariate discrete distribution derived from two Rayleigh distributions using a method proposed by Marshall and Olkin (1997) where an additional parameter is is introduced to a family of distributions related to the dependence structure of two discrete random variables $X_{1}$ and $X_{2}$ . The study results show that this new bivariate distribution has good statistical properties and a simple mathematical expression for its correlation coefficient. The usual classical and Bayesian estimators for the parameters of the new distribution are also presented. A simulation study is carried out in order to evaluate some frequentist properties of the proposed model. The usefulness of the proposed model is illustrated using a real medical dataset introduced in the literature in presence of censoring and covariates.

Keywords

Bayesian inference correlation coefficient Marshall and Olkin approach maximum likelihood medical data Rayleigh distribution

1. Introduction

The introduction of new survival distributions has been the objective of many studies introduced in the literature since this class of models is widely used in many areas of application, especially in medical or engineering studies where the lifetimes usually are measured in a continuous scale. However, in many applications, it is possible to have the observed responses measured in a discrete scale where it is not appropriate to assume the data as continuous data. Nevertheless, in many applications despite the discrete data, it is common to use continuous models due to model analytical tractability or facilities in obtaining parameter estimators using existing statistical software where the most popular continuous lifetime distributions usually are implemented. However, this procedure may not be appropriate in many applications as it may be impossible to measure the lifetime of a component on a continuous scale, as in the on/off switching devices, number of cycles until failure, or the number of accidents in a road (see for example, Kundu, 2014, 2017; Kundu et al., 2010; Kundu & Gupta, 2009; Kundu & Nekoukhou, 2018).

As more specific examples of lifetime applications, we could consider the times of deterioration levels, the times of infection or the times to reaction for a treatment in pairs of lungs, kidneys, eyes or ears of humans that could be given as count of the number of days, weeks or months. In this case, the use of a continuous distribution could lead to inaccurate inference results which is a great motivation to introduce new and more flexible discrete bivariate lifetime distributions together with their mathematical properties and inference results. In this direction special attention has been given on bivariate geometric distributions and bivariate Poisson distributions (see for example Kocherlakota S & Kocherlakota K, 1992; Kocherlakota, 1995; Arnold, 1975; Basu & Dhar, 1995; Kumar, 2008; Kemp, 2013; Lee & Cha, 2014; Nekoukhou & Kundu, 2017; Kundu & Nekoukhou, 2018) as alternatives to many bivariate continuous models introduced in the literature (see for example, Block & Basu, 1974; Marshall & Olkin, 1967a,b; Downton, 1970; Freund, 1961; Sarkar, 1987; Arnold & Strauss, 1988; Gumbel, 1960; Hanagal, 2006; Hanagal & Ahmadi, 2008; Hawkes, 1972; Hougaard, 1986; Balakrishnan & Lai, 2009; Popović & Genç, 2018).

In some cases, the bivariate distributions of survival data introduced in the literature may present inconvenient forms to be used in applications, as the two classes of discrete bivariate distributions introduced by Lee and Cha (2014), although the motivation for this model to be quite simple based on the minimum and maximum of two different independent random variables. In this model there are great difficulties to calculate the estimators of the unknown parameters and to determine their properties (see also, Nekoukhou & Kundu, 2017).

From these considerations, there is a great motivation to introduce new bivariate lifetime distributions with simple mathematical properties and simplifications to get the inferences of interest, the main goal of this study. In this direction, this paper introduces a new discrete bivariate Rayleigh distribution obtained by using a procedure introduced by Marshall and Olkin (1997) where the marginal distributions are discrete generalized Rayleigh distributions implying in great model flexibility for the hazard function, which could be given by bathtub, increasing and increasing-decreasing-increasing shapes depending on the parameter values. This paper also introduces a multivariate extension of the proposed bivariate model. An additional justification for the use of bivariate survival distributions is the possible simplification of the mathematical expression for the likelihood function in the presence of censored data, a common situation in lifetime data applications as in medical or engineering studies, when compared to some existing continuous bivariate lifetime distributions where the likelihood function in presence of censored data usually depends on joint survival functions given in non-analytical mathematical forms.

In the construction of bivariate probability distributions, especially for the continuous case, the literature presents many different techniques such as: the use of copula functions, mixing and compounding; the use of trivariate reduction; the specification of a conditional and a marginal distribution; the use of a conditionally method; the construction of discrete bivariate distributions with specified marginals and correlation; the use of sums and limits of Bernoulli random variables; the use of clusters; the construction of finite bivariate distributions via extreme points using convex sets; the use of generalized distributions methods; the use of canonical correlation coefficients and semi-groups; the use of bivariate distributions generated from weight functions; the use of marginal transformation methods; the use of truncation or the use of two-stage failure risks of a two-component parallel system method (see Kocherlakota S & Kocherlakota K, 1992; Lai, 2006; Balakrishnan and Lai, 2009). However, unlike their continuous analogues, discrete bivariate distributions usually are harder to be constructed. Kemp and Papageorgiou (1982) pointed out that the main problem in the construction of bivariate distributions is the impossibility to have a standard set of criteria that can always be applied to produce a unique bivariate distribution obtained from an univariate distribution which could unequivocally be called the bivariate version.

Meanwhile, Marshall and Olkin (1997) introduced a method to obtain an extended family of distributions including one additional parameter denoted as univariate Marshall-Olkin family having cumulative distribution function $G(x)$ and survival function $\bar{G}(x)$ given, respectively, by,

$\displaystyle G(x)=\frac{F(x)}{\alpha+\bar{\alpha}F(x)}\ \text{and}\ \bar{G}(x% )=\frac{\alpha S(x)}{1-\bar{\alpha}S(x)}$ (1)

where $S(x)=\Pr(X>x)$ , $F(x)=1-S(x)$ , $\alpha>0,\bar{\alpha}=1-\alpha$ and $-\infty<x<\infty$ . This new family of distributions has an additional parameter $\alpha$ which generalizes the baseline distribution and is related to the dependence structure of two random variables. This family could also be extended to the multivariate case. In this way, let $\mathbf{X}=(X_{1},\ldots,X_{n})$ be a random vector with multivariate cumulative and survival functions given, respectively, by $F(x_{1},\ldots,x_{n})$ and $S(x_{1},\ldots,x_{n})$ ; the multivariate Marshall-Olkin family has cumulative distribution function $G(x_{1},\ldots,x_{n})$ and survival function $\bar{G}(x_{1},\ldots,x_{n})$ given, respectively, by,

$\displaystyle G(x_{1},\ldots,x_{n})=\frac{F(x_{1},\ldots,x_{n})}{\alpha+\bar{% \alpha}F(x_{1},\ldots,x_{n})}\ \text{and}\ \bar{G}(x_{1},\ldots,x_{n})=\frac{% \alpha S(x_{1},\ldots,x_{n})}{1-\bar{\alpha}S(x_{1},\ldots,x_{n})}$ (2)

where $\alpha>0,\bar{\alpha}=1-\alpha$ and $-\infty<x_{i}<\infty$ . Notice that the functions obtained in Eq. (2) are more flexible than the functions obtained from the product of independent distributions providing a great flexibility in the modeling of the dependence structure.

The main goal of this paper is to introduce a new bivariate discrete distribution with great flexibility for the correlation coefficient range and to avoid not very accurate inferential results usually obtained using other existing bivariate lifetime models. In this way, the Marshall and Olkin (1997) method is considered in this study in the construction of a bivariate discrete distribution due to its flexibility to introduce a new parameter related to the dependence of the two random variables $X_{1}$ and $X_{2}$ . The Rayleigh distribution is also used to compose the bivariate model due to its great applicability to reliability and survival medical data. In this way, the resulting distribution could be a good alternative to other existing bivariate discrete distributions introduced in the literature.

This paper is organized as follows: Section 2 presents the bivariate discrete generalized Rayleigh distribution and its basic properties. The structural properties including the mean, variance, the dependence structure, a formula to compute bivariate probabilities and its implementation are also presented in Section 2. Maximum likelihood and Bayesian estimators for the parameters of the proposed model are presented in Section 3. Section 4 presents the results of a simulation study; Section 5 presents two applications considering bivariate medical survival data. Finally, Section 6 closes the paper with some concluding remarks.

2. Bivariate discrete generalized rayleigh distribution

2.1 Model description

Let $X_{i},i=1,2$ be two independent discrete random variables having Rayleigh distributions (see Roy, 2004) with parameters $0<\lambda_{i}<1,i=1,2$ . Since $X_{i},i=1,2$ are independent random variables, the joint survival function of the bivariate random variables $X_{1}$ and $X_{2}$ is given by,

$\displaystyle\text{Pr}_{(X_{1},X_{2})}(X_{1}>x_{1},X_{2}>x_{2})=\lambda_{1}^{x% _{1}^{2}}\lambda_{2}^{x_{2}^{2}}.$ (3)

Observe that the joint survival function in Eq. (3) is restricted to independent lifetimes and cannot be applied directly assuming dependence structures in bivariate data. In this way, using the Marshall-Olkin survival function given in Eq. (2), the new proposed joint survival function is given by,

$\displaystyle\textrm{Pr}(X_{1}>x_{1},X_{2}>x_{2})=\frac{\alpha\lambda_{1}^{x_{% 1}^{2}}\lambda_{2}^{x_{2}^{2}}}{1-\bar{\alpha}\lambda_{1}^{x_{1}^{2}}\lambda_{% 2}^{x_{2}^{2}}}$ (4)

where $\bar{\alpha}=1-\alpha$ and $\alpha>0$ . Observe that Eq. (4) is more flexible than Eq. (3) and can be applied directly to model the dependence structure for correlated bivariate lifetimes. The joint survival function defined by Eq. (4) is denoted as a bivariate discrete generalized Rayleigh (BDGR) distribution and special cases of its contour plot are illustrated in Fig. 1, assuming different values for the parameters $\lambda_{1},\lambda_{2}$ and $\alpha$ .

Figure 1.

Discrete contour plots of the joint survival function for the BDGR model assuming different parameter values (Upper-panels: fixed values given by $\lambda_{1}=\lambda_{2}=$ 0.95, and $\alpha=$ 0.50 $\rightarrow$ 1.00 $\rightarrow$ 1.50s. Lower-panels: fixed values given by $\lambda_{1}=\lambda_{2}=$ 0.99, and $\alpha=$ 0.50 $\rightarrow$ 1.00 $\rightarrow$ 1.50).

.

Let $\mathbf{X}=(X_{1},X_{2})$ be a discrete random vector following the joint survival function given by Eq. (4) with parameters $0<\lambda_{1},\lambda_{2}<1$ and $\alpha>0$ . Defining $h(z_{1},z_{2})=\lambda_{1}^{z_{1}^{2}}\lambda_{2}^{z_{2}^{2}}$ , $z_{1},z_{2}\in\mathbb{R}_{+}$ , the joint probability mass function (pmf) of $\mathbf{X}$ is given by,

$\displaystyle\textrm{Pr}(X_{1}=x_{1},X_{2}=x_{2})=\frac{\alpha[1-\bar{\alpha}^% {2}h(k_{1},k_{2})]\left[1-\bar{\alpha}h(x_{1},x_{2})\right]^{-1}\textit{DR}(x_% {1},\lambda_{1})\textit{DR}(x_{2},\lambda_{2})}{\left[1-\bar{\alpha}\,h(x_{1}+% 1,x_{2})\right]\left[1-\bar{\alpha}\,h(x_{1},x_{2}+1)\right]\left[1-\bar{% \alpha}\,h(x_{1}+1,x_{2}+1)\right]}$ (5)

where $k_{i}=\sqrt{x_{i}^{2}+(x_{i}+1)^{2}},i=1,2$ ; $0<\lambda_{i}<1,i=1,2$ ; $\alpha>0$ ; $\bar{\alpha}=1-\alpha$ and $\textit{DR}(x_{i},\lambda_{i})=\lambda_{i}^{x_{i}^{2}}-\lambda_{i}^{(x_{i}+1)^% {2}},\,i=1,2$ denotes a univariate discrete Rayleigh distribution. Observe that Eq. (5) is a proper joint pmf since the series expressed by ${\sum}_{x_{1}=0}^{\infty}{\sum}_{x_{2}=0}^{\infty}h(x_{1},x_{2})$ converges to $\frac{1}{4}[\vartheta_{3}(0,\lambda_{1})+1][\vartheta_{3}(0,\lambda_{2})+1]$ for $z_{1},z_{2}\in\mathbb{R}_{+}$ where $\vartheta_{a}(x,q)$ is the Jacobi theta function.

2.2 Marginal and conditional distributions

For the Bivariate Discrete Generalized Rayleigh Distribution, the marginal distributions of the random variables $X_{1}$ and $X_{2}$ are given by discrete generalized Rayleigh (DGR) distributions with corresponding parameters $(\lambda_{1},\alpha)$ and $(\lambda_{2},\alpha)$ , respectively. Many properties of the continuous generalized Rayleigh distribution were introduced by MirMostafaee et al. (2017). Since the DGR inherits most of the properties of the continuous model in the univariate case, these marginals have great flexibility for the hazard function, given by bathtub, increasing and increasing-decreasing-increasing shapes depending on the parameter values. In this way, this bivariate distribution could be very useful especially to model positive data. The marginal survival and marginal pmf functions can be expressed, respectively, by,

$\displaystyle\text{Pr}(X_{i}>x_{i})=\frac{\alpha\lambda_{i}^{x_{i}^{2}}}{1-% \bar{\alpha}\lambda_{i}^{x_{i}^{2}}}\ \text{and}\ \text{Pr}(X_{i}=x_{i})=\frac% {\alpha\left(\lambda_{i}^{x_{i}^{2}}-\lambda_{i}^{(x_{i}+1)^{2}}\right)}{\left% (1-\bar{\alpha}\lambda_{i}^{(x_{i}+1)^{2}}\right)\left(1-\bar{\alpha}\lambda_{% i}^{x_{i}^{2}}\right)},\,\,i=1,2.$

The mean expected value and variance for the marginal probability distributions do not have closed forms, however they could be computed using numerical methods from the mean and variance around the origin given by,

$\displaystyle\mathbb{E}[X_{i}]=\alpha\left(\sum_{x_{i}=1}^{\infty}\frac{% \lambda_{i}^{x_{i}^{2}}}{1-\bar{\alpha}\lambda_{i}^{x_{i}^{2}}}\right)\quad% \text{and}\quad\text{Var}(X_{i})=\sum_{x_{i}=1}^{\infty}\frac{(2x_{i}^{2}-1)% \alpha\lambda_{i}^{x_{i}^{2}}}{1-\bar{\alpha}\lambda_{i}^{x_{i}^{2}}}-\left(% \sum_{x_{i}=1}^{\infty}\frac{\alpha\lambda_{i}^{x_{i}^{2}}}{1-\bar{\alpha}% \lambda_{i}^{x_{i}^{2}}}\right)^{2},i=1,2$

which could be approximated using the finite series given by:

$\displaystyle\sum_{x=1}^{\infty}\frac{\lambda_{i}^{x^{2}}}{1-\bar{\alpha}% \lambda_{i}^{x^{2}}}=\sum_{x=1}^{\infty}\frac{d}{d\alpha}\log(1-\bar{\alpha}% \lambda_{i}^{x^{2}})\approx\sum_{j=1}^{M}(-1)^{j+1}(\alpha-1)^{(j-1)^{2}}\frac% {\lambda_{i}^{j^{2}}}{1-\lambda_{i}^{j^{2}}}.$ (6)

Now, let us assume the transformation of the random variables $X_{1}$ and $X_{2}$ given by $W=\min(X_{1},X_{2})$ . In this case, the cumulative probability function of the random variable W is given by,

$\displaystyle\text{Pr}(W<w)=1-\text{Pr}(W>w)=1-\frac{\alpha\lambda_{1}^{w^{2}}% \lambda_{2}^{w^{2}}}{1-\bar{\alpha}\lambda_{1}^{w^{2}}\lambda_{2}^{w^{2}}}=1-% \frac{\alpha(\lambda_{1}\lambda_{2})^{w^{2}}}{1-\bar{\alpha}(\lambda_{1}% \lambda_{2})^{w^{2}}}$

which implies that the distribution of $W$ is a discrete generalized Rayleigh (DGR) distribution with parameters $(\lambda_{1}\lambda_{2},\alpha)$ . The mean expected value and variance could also be approximated by the finite series given in Eq. (6).

For the BDGR model, the conditional distribution of $X_{j}$ given $X_{i},i,j=1,2$ and $i\neq j$ is given by,

$\displaystyle\text{P}(X_{j}\mid X_{i}=x_{i})=\frac{[1-\bar{\alpha}^{2}h(k_{i},% k_{j})]\left[1-\bar{\alpha}h(x_{i},x_{j})\right]^{-1}[1-\bar{\alpha}(\textit{% DR}(x_{i},\lambda_{i})+\bar{\alpha}\lambda_{i}^{(x_{i}+1)^{2}}\lambda_{i}^{x_{% i}^{2}})]}{\left[1-\bar{\alpha}\,h(x_{i}+1,x_{j})\right]\left[1-\bar{\alpha}\,% h(x_{i},x_{j}+1)\right]}\times\frac{[\textit{DR}(x_{j},\lambda_{j})]}{\left[1-% \bar{\alpha}\,h(x_{i}+1,x_{j}+1)\right]}$

2.3 Covariance and correlation coefficient

For the BDGR distribution given by Eq. (5), the cross factorial moment between $X_{1}$ and $X_{2}$ is given by,

$\displaystyle\mu_{X_{1},X_{2}}=\mathbb{E}[X_{1}X_{2}]=\sum_{x_{1}=1}^{\infty}% \sum_{x_{2}=1}^{\infty}\frac{\alpha\lambda_{1}^{x_{1}^{2}}\lambda_{2}^{x_{2}^{% 2}}}{1-\bar{\alpha}\lambda_{1}^{x_{1}^{2}}\lambda_{2}^{x_{2}^{2}}}.$ (7)

which is a monotonic increasing function of $\lambda_{1},\lambda_{2}$ and $\alpha$ since,

$\displaystyle\frac{\partial\mu_{X_{1},X_{2}}}{\partial\alpha}=\alpha\sum_{x_{1% }=1}^{\infty}\sum_{x_{2}=1}^{\infty}\frac{(1-\lambda_{1}^{x_{1}^{2}}\lambda_{2% }^{x_{2}^{2}})\lambda_{1}^{x_{1}^{2}}\lambda_{2}^{x_{2}^{2}}}{(1-\bar{\alpha}% \lambda_{1}^{x_{1}^{2}}\lambda_{2}^{x_{2}^{2}})^{2}}>0\ \text{and}\ \frac{% \partial\mu_{X_{1},X_{2}}}{\partial\lambda_{i}}=\alpha\sum_{x_{1}=1}^{\infty}% \sum_{x_{2}=1}^{\infty}\frac{x_{i}^{2}\lambda_{i}^{x_{i}^{2}-1}\lambda_{j}^{x_% {j}^{2}}}{(1-\bar{\alpha}\lambda_{i}^{x_{i}^{2}}\lambda_{j}^{x_{j}^{2}})^{2}}>0$

for $i,j=1,2;i\neq j$ . However, the cross factorial also could be approximated using series representation of the logarithmic function. That is, suppose that $M$ is an integer sufficiently large and $|(\alpha-1)\lambda_{i}x_{i}^{2}|<1,i=1,2$ ; thus, we have,

$\displaystyle\mu_{X_{1},X_{2}}=\sum_{x_{1}=1}^{\infty}\sum_{x_{2}=1}^{\infty}% \frac{d}{d\alpha}\log(1-\bar{\alpha}\lambda_{1}^{x_{1}^{2}}\lambda_{2}^{x_{2}^% {2}})\approx\sum_{j=1}^{M}(\alpha-1)^{(j-1)^{2}}\frac{\lambda_{1}^{j^{2}}% \lambda_{2}^{j^{2}}}{(1-\lambda_{1}^{j^{2}})(1-\lambda_{2}^{j^{2}})}$

which is a finite series and can be determined when the parameters of the distribution are estimated from a dataset.

.

(Covariance signal) Let the covariance signal be defined by the function,

$\displaystyle\Psi(\alpha)=\left(1-\bar{\alpha}\lambda_{1}^{x_{1}^{2}}\right)% \left(1-\bar{\alpha}\lambda_{2}^{x_{2}^{2}}\right)-\alpha\left(1-\bar{\alpha}% \lambda_{1}^{x_{1}^{2}}\lambda_{2}^{x_{2}^{2}}\right),$ (8)

thus $\Psi(\alpha)>0$ if $0<\alpha<1$ , $\Psi(\alpha)=1$ if $\alpha=0$ and $\Psi(\alpha)<0$ if $\alpha>1$ . That is, the covariance signal only depends on the parameter $\alpha$ .

Proof..

Note that the function $\Psi(\alpha)$ is a continuous function on $\alpha$ and its first derivative is given by,

$\displaystyle\Psi^{\prime}(\alpha)=-\left(1-\lambda_{1}^{x_{1}^{2}}\right)% \left(1-\lambda_{2}^{x_{2}^{2}}\right)<0,\ \text{for all}\,\,\,\alpha>0.$

Note that, since $\Psi(0)=\left(1-\lambda_{1}^{x_{1}^{2}}\right)\left(1-\lambda_{2}^{x_{2}^{2}}% \right)>0$ and $\Psi(1)=0$ , then $\Psi(\alpha)>0$ if $0<\alpha<1$ . On the other hand, if $\alpha=1$ then $\Psi(\alpha)=0$ ; if $\alpha>1$ then $\Psi(\alpha)<0$ and the proof is complete. ∎

From the function $\Psi(\alpha)$ given in Eq. (8), the proposed model admits a very flexible behavior for the correlation coefficient $\rho$ of any sign. In fact it is observed that if $\Psi(\alpha)>0$ , $\rho>0$ and if $\Psi(\alpha)<0,\rho<0$ . If $\Psi(\alpha)=0$ the correlation coefficient is equal to zero. Although there is no closed form for the expressions of the covariance and correlation coefficient these expressions could be computed by taking a large number of terms in the series defined in Eq. (6). Table 1 presents the approximated results for the covariance and correlation coefficients for some fixed values of $\alpha,\lambda_{1}$ and $\lambda_{2}$ from which can be seen the flexibility of the correlation coefficient $\rho$ .

Table 1

Theoretical results under a BDGR distribution

$(\alpha,\lambda_{1},\lambda_{2})$	$\mathbb{E}[X_{1}X_{2}]$	$\text{cov}(X_{1},X_{2})$	$\rho$
(0.2, 0.2, 0.2)	0.0084	0.0061	0.1281
(0.8, 0.2, 0.2)	0.0327	0.0046	0.0310
(0.5, 0.5, 0.5)	0.1777	0.0433	0.0980
(1.0, 0.5, 0.5)	0.3186	0.0000	0.0000
(1.2, 0.5, 0.5)	0.3676	$-$ 0.0191	$-$ 0.0267
(1.5, 0.5, 0.5)	0.4348	$-$ 0.0466	$-$ 0.0579
(1.8, 0.8, 0.8)	2.6964	$-$ 0.2225	$-$ 0.0298
(2.0, 0.8, 0.8)	2.8597	$-$ 0.2721	$-$ 0.0342

2.4 Special cases and multivariate extension

.

(Special cases) Some especial cases of the BDGR are given by,

i.)
If $\lambda_{1}=\lambda_{2}$ , the pmf of BDGR is symmetric in its arguments, that is, $\text{Pr}(X_{1}=x_{1},X_{2}=x_{2})=\text{Pr}(X_{2}=x_{2},X_{1}=x_{1})$ for all $x_{1},x_{2}\in\mathbb{N}$ .
ii.)
If $\lambda_{1}=\lambda_{2}=\alpha=\lambda$ , the pmf of BDGR is also symmetric in its arguments and its reduced to a one parameter bivariate discrete distribution with probability mass function given by,

$\displaystyle\frac{\lambda^{x_{1}+x_{2}+1}[(1-\lambda)^{2}-(1-\lambda)^{4}% \lambda^{x_{1}^{2}+x_{2}^{2}+(x_{1}+1)^{2}+(x_{2}+1)^{2}}]}{(1-\lambda^{x_{1}+% x_{2}+1})(1-\lambda^{x_{1}+x_{2}+2})(1-\lambda^{x_{1}+x_{2}+3})}$ (9)
iii.)
If $0<\alpha<1$ and $\lambda_{1}=\lambda_{2}=\lambda$ , the pmf of BDGR can be rewritten as an infinite mixture of the product of two Rayleigh distributions.

Proof..

(i) and (ii) are trivial. For (iii), the result is obtained by considering the series representation,

$\displaystyle(1-b)^{-k}=\sum_{j=0}^{\infty}\frac{\Gamma(k+j)}{\Gamma(j+1)% \Gamma(k)}b^{j},\,|b|<1,\,k>0$ (10)

and noticing that $\text{Pr}(X_{1}>x_{1},X_{2}>x_{2})=\frac{\alpha\lambda^{x_{1}^{2}+x_{2}^{2}}}{% 1-\bar{\alpha}\lambda^{x_{1}^{2}+x_{2}^{2}}}=\alpha\lambda^{x_{1}^{2}+x_{2}^{2% }}{\sum}_{j=0}^{\infty}(\bar{\alpha}\lambda^{x_{1}^{2}+x_{2}^{2}})^{j}$ . ∎

.

(Multivariate extension)Let $X_{i},i=1,\ldots,n$ be $n$ independent discrete random variables having the Rayleigh distribution with parameters $0<\lambda_{i}<1,i=1,\ldots,n$ . Since $X_{i},i=1,\ldots,n$ are independent, the multivariate survival and multivariate pmf functions are given, respectively by,

$\displaystyle P(X_{1}>x_{1},\ldots,X_{n}>x_{n})=\frac{\alpha\lambda_{1}^{x_{1}% ^{2}}\ldots\lambda_{n}^{x_{n}^{2}}}{1-\bar{\alpha}\lambda_{1}^{x_{1}^{2}}% \ldots\lambda_{n}^{x_{n}^{2}}}$ (11)

and,

$\displaystyle P(X_{1}=x_{1},\ldots,X_{n}=x_{n})=\frac{\alpha{\prod}_{i=1}^{n}% \lambda_{i}^{x_{i}^{2}}}{1-\bar{\alpha}{\prod}_{i=1}^{n}\lambda_{i}^{x_{i}^{2}% }}-\frac{\alpha\lambda_{1}^{x_{1}^{2}}{\prod}_{i=2}^{n}\lambda_{i}^{(x_{i}+1)^% {2}}}{1-\bar{\alpha}\lambda_{1}^{x_{1}^{2}}{\prod}_{i=2}^{n}\lambda_{i}^{(x_{i% }+1)^{2}}}-\ldots+(-1)^{n}\frac{\alpha{\prod}_{i=1}^{n}\lambda_{i}^{(x_{i}+1)^% {2}}}{1-\bar{\alpha}{\prod}_{i=1}^{n}\lambda_{i}^{(x_{i}+1)^{2}}}$ (12)

For the multivariate model, the dependence among the n lifetimes is specified by the parameter $\alpha$ , where if $\alpha=1$ there is independence among the n lifetimes. The model adequacy could be checked by comparisons of the fitted marginal survival functions with the empirical estimates of the marginal survival distributions since the marginal survival functions are discrete generalized Rayleigh distributions as well. Moreover, the correlation coefficient for the multivariate case has the same properties of the correlation coefficient of the bivariate model, that is, it could be negative, positive or equal to zero. Since extensive calculations are required for the multivariate case, their properties are not derived here.
3. Inference methods

3.1 Maximum likelihood approach

This section introduces the maximum likelihood estimation (MLE) method in two situations: the situation assuming censored data and the situation with complete data. In both cases, the maximum likelihood estimators do not have closed forms requiring the use of numerical methods like Newton-Raphson or Nelder-Mead to get the estimators for each parameter.

3.1.1 Complete data

Let $(X_{11},X_{21}),(X_{12},X_{22}),\ldots,(X_{1n},X_{2n})$ be a random sample of size $n$ from a BDGR distribution. The log-likelihood function $\ell(\lambda_{1},\lambda_{2},\alpha)$ is given by:

$\displaystyle\ell(\lambda_{1},\lambda_{2},\alpha)=(n_{1}+n_{2})\log(\alpha)+% \sum_{s=1}^{n_{1}}\sum_{t=1}^{n_{2}}\log\left(1-\bar{\alpha}^{2}\lambda_{1}^{x% _{1s}^{2}}\lambda_{2}^{x_{2t}^{2}}\lambda_{1}^{(x_{1s}+1)^{2}}\lambda_{2}^{(x_% {2t}+1)^{2}}\right)+\sum_{s=1}^{n_{1}}\log\left(\lambda_{1}^{x_{1s}^{2}}-% \lambda_{1}^{(x_{1s}+1)^{2}}\right)+\sum_{t=1}^{n_{2}}\log\left(\lambda_{2}^{x% _{2t}^{2}}-\lambda_{2}^{(x_{2t}+1)^{2}}\right)-\sum_{s=1}^{n_{1}}\sum_{t=1}^{n% _{2}}\log\left(1-\bar{\alpha}\lambda_{1}^{(x_{1s}+1)^{2}}\lambda_{2}^{x_{2t}^{% 2}}\right)-\sum_{s=1}^{n_{1}}\sum_{t=1}^{n_{2}}\log\left(1-\bar{\alpha}\lambda% _{1}^{x_{1s}^{2}}\lambda_{2}^{(x_{2t}+1)^{2}}\right)-\sum_{s=1}^{n_{1}}\sum_{t% =1}^{n_{2}}\log\left(1-\bar{\alpha}\lambda_{1}^{x_{1s}^{2}}\lambda_{2}^{x_{2t}% ^{2}}\right)-\sum_{s=1}^{n_{1}}\sum_{t=1}^{n_{2}}\log\left(1-\bar{\alpha}% \lambda_{1}^{(x_{1s}+1)^{2}}\lambda_{2}^{(x_{2t}+1)^{2}}\right).$ (13)

The normal equations (the partial derivatives of the log-likelihood function with respect to the parameters are equal to zero) obtained from Eq. (13), are not reproduced here as they cannot be solved explicitly. They must be solved either by numerical methods as, for example, the Newton-Raphson optimization method or by directly maximization of the log-likelihood function. Since the global maximum of the logâ€“likelihood surface is not guaranteed, different initial values in the parameter space sould be considered as seed points. From the log-likelihood function, the first derivatives of $\ell(\lambda_{1},\lambda_{2},\alpha)$ with respect to $\lambda_{1},\lambda_{2}$ and $\alpha$ are given, respectively by

i.)

$\displaystyle\frac{\partial\ell}{\partial\lambda_{1}}=\mathop{\sum}_{s=1}^{n_{% 1}}\mathop{\sum}_{t=1}^{n_{2}}\displaystyle\frac{\left(\bar{\alpha}(x_{1s}+1)^% {2}\lambda_{1}^{(x_{1s}+1)^{2}-1}\lambda_{2}^{(x_{2t}+1)^{2}}\right)}{\left(1-% \bar{\alpha}\lambda_{1}^{(x_{1s}+1)^{2}}\lambda_{2}^{(x_{2t}+1)^{2}}\right)}+% \mathop{\sum}_{s=1}^{n_{1}}\mathop{\sum}_{t=1}^{n_{2}}\frac{\left(\bar{\alpha}% x_{1s}^{2}\lambda_{1}^{x_{1s}^{2}-1}\lambda_{2}^{(x_{2t}+1)^{2}}\right)}{\left% (1-\bar{\alpha}\lambda_{1}^{x_{1s}^{2}}\lambda_{2}^{(x_{2t}+1)^{2}}\right)}\\ \hskip 34.143307pt{+}\mathop{\sum}_{s=1}^{n_{1}}\displaystyle\frac{\left(x_{1s% }^{2}\lambda_{1}^{x_{1s}^{2}-1}-(x_{1s}+1)^{2}\lambda_{1}^{(x_{1s}+1)^{2}-1}% \right)}{\left(\lambda_{1}^{x_{1s}^{2}}-\lambda_{1}^{(x_{1s}+1)^{2}}\right)}+% \mathop{\sum}_{s=1}^{n_{1}}\mathop{\sum}_{t=1}^{n_{2}}\frac{\left(\bar{\alpha}% (x_{1s}+1)^{2}\lambda_{1}^{(x_{1s}+1)^{2}-1}\lambda_{2}^{x_{2t}^{2}}\right)}{% \left(1-\bar{\alpha}\lambda_{1}^{(x_{1s}+1)^{2}}\lambda_{2}^{x_{2t}^{2}}\right% )}\\ \hskip 34.143307pt+\mathop{\sum}_{s=1}^{n_{1}}\mathop{\sum}_{t=1}^{n_{2}}% \displaystyle\frac{\left(\bar{\alpha}x_{1s}^{2}\lambda_{1}^{x_{1s}^{2}-1}% \lambda_{2}^{x_{2t}^{2}}\right)}{\left(1-\bar{\alpha}\lambda_{1}^{x_{1s}^{2}}% \lambda_{2}^{x_{2t}^{2}}\right)}-\mathop{\sum}_{s=1}^{n_{1}}\mathop{\sum}_{t=1% }^{n_{2}}\displaystyle\frac{\bar{\alpha}^{2}x_{1s}^{2}\lambda_{1}^{x_{1s}^{2}-% 1}\lambda_{2}^{x_{2t}^{2}}(x_{1s}+1)^{2}\lambda_{1}^{(x_{1s}+1)^{2}-1}\lambda_% {2}^{(x_{2t}+1)^{2}}}{\left(1-\bar{\alpha}^{2}\lambda_{1}^{x_{1s}^{2}}\lambda_% {2}^{x_{2t}^{2}}\lambda_{1}^{(x_{1s}+1)^{2}}\lambda_{2}^{(x_{2t}+1)^{2}}\right)}$

ii.)

$\displaystyle\frac{\partial\ell}{\partial\lambda_{2}}=\mathop{\sum}_{t=1}^{n_{% 2}}\frac{\left(x_{2t}^{2}\lambda_{2}^{x_{2t}^{2}-1}-(x_{2t}+1)^{2}\lambda_{2}^% {(x_{2t}+1)^{2}-1}\right)}{\left(\lambda_{2}^{x_{2t}^{2}}-\lambda_{2}^{(x_{2t}% +1)^{2}}\right)}+\mathop{\sum}_{s=1}^{n_{1}}\mathop{\sum}_{t=1}^{n_{2}}% \displaystyle\frac{\left(\bar{\alpha}x_{2t}^{2}\lambda_{1}^{(x_{1s}+1)^{2}}% \lambda_{2}^{x_{2t}^{2}-1}\right)}{\left(1-\bar{\alpha}\lambda_{1}^{(x_{1s}+1)% ^{2}}\lambda_{2}^{x_{2t}^{2}}\right)}\\ \hskip 34.143307pt+\mathop{\sum}_{s=1}^{n_{1}}\mathop{\sum}_{t=1}^{n_{2}}% \displaystyle\frac{\left(\bar{\alpha}(x_{2t}+1)^{2}\lambda_{1}^{(x_{1s}+1)^{2}% }\lambda_{2}^{(x_{2t}+1)^{2}-1}\right)}{\left(1-\bar{\alpha}\lambda_{1}^{(x_{1% s}+1)^{2}}\lambda_{2}^{(x_{2t}+1)^{2}}\right)}+\mathop{\sum}_{s=1}^{n_{1}}% \mathop{\sum}_{t=1}^{n_{2}}\displaystyle\frac{\left(\bar{\alpha}(x_{2t}+1)^{2}% \lambda_{1}^{x_{1s}^{2}}\lambda_{2}^{(x_{2t}+1)^{2}-1}\right)}{\left(1-\bar{% \alpha}\lambda_{1}^{x_{1s}^{2}}\lambda_{2}^{(x_{2t}+1)^{2}}\right)}\\ \hskip 34.143307pt+\mathop{\sum}_{s=1}^{n_{1}}\mathop{\sum}_{t=1}^{n_{2}}% \displaystyle\frac{\left(\bar{\alpha}x_{2t}^{2}\lambda_{1}^{x_{1s}^{2}}\lambda% _{2}^{x_{2t}^{2}-1}\right)}{\left(1-\bar{\alpha}\lambda_{1}^{x_{1s}^{2}}% \lambda_{2}^{x_{2t}^{2}}\right)}-\mathop{\sum}_{s=1}^{n_{1}}\mathop{\sum}_{t=1% }^{n_{2}}\displaystyle\frac{\bar{\alpha}^{2}x_{2t}^{2}\lambda_{1}^{x_{1s}^{2}}% \lambda_{2}^{x_{2t}^{2}-1}(x_{2t}+1)^{2}\lambda_{1}^{(x_{1s}+1)^{2}}\lambda_{2% }^{(x_{2t}+1)^{2}-1}}{\left(1-\bar{\alpha}^{2}\lambda_{1}^{x_{1s}^{2}}\lambda_% {2}^{x_{2t}^{2}}\lambda_{1}^{(x_{1s}+1)^{2}}\lambda_{2}^{(x_{2t}+1)^{2}}\right)}$

iii.)

$\displaystyle\frac{\partial\ell}{\partial\alpha}=\displaystyle\frac{n_{1}+n_{2% }}{\alpha}+\mathop{\sum}_{s=1}^{n_{1}}\mathop{\sum}_{t=1}^{n_{2}}\displaystyle% \frac{2\bar{\alpha}\lambda_{1}^{x_{1s}^{2}}\lambda_{2}^{x_{2t}^{2}}\lambda_{1}% ^{(x_{1s}+1)^{2}}\lambda_{2}^{(x_{2t}+1)^{2}}}{\left(1-\bar{\alpha}^{2}\lambda% _{1}^{x_{1s}^{2}}\lambda_{2}^{x_{2t}^{2}}\lambda_{1}^{(x_{1s}+1)^{2}}\lambda_{% 2}^{(x_{2t}+1)^{2}}\right)}\\ \hskip 34.143307pt+\mathop{\sum}_{s=1}^{n_{1}}\mathop{\sum}_{t=1}^{n_{2}}% \displaystyle\frac{\lambda_{1}^{(x_{1s}+1)^{2}}\lambda_{2}^{(x_{2t}+1)^{2}}}{% \left(1-\bar{\alpha}\lambda_{1}^{(x_{1s}+1)^{2}}\lambda_{2}^{(x_{2t}+1)^{2}}% \right)}+\mathop{\sum}_{s=1}^{n_{1}}\mathop{\sum}_{t=1}^{n_{2}}\frac{\lambda_{% 2}^{(x_{2t}+1)^{2}}\lambda_{1}^{x_{1s}^{2}}}{\left(1-\bar{\alpha}\lambda_{2}^{% (x_{2t}+1)^{2}}\lambda_{1}^{x_{1s}^{2}}\right)}\\ \hskip 34.143307pt\mathop{\sum}_{s=1}^{n_{1}}\mathop{\sum}_{t=1}^{n_{2}}\frac{% \lambda_{1}^{(x_{1s}+1)^{2}}\lambda_{2}^{x_{2t}^{2}}}{\left(1-\bar{\alpha}% \lambda_{1}^{(x_{1s}+1)^{2}}\lambda_{2}^{x_{2t}^{2}}\right)}+\mathop{\sum}_{s=% 1}^{n_{1}}\mathop{\sum}_{t=1}^{n_{2}}\displaystyle\frac{\lambda_{1}^{x_{1s}^{2% }}\lambda_{2}^{x_{2t}^{2}}}{\left(1-\bar{\alpha}\lambda_{1}^{x_{1s}^{2}}% \lambda_{2}^{x_{2t}^{2}}\right)}$

Under standard asymptotic maximum likelihood theory, a consistent estimator for the covariance matrix of $(\widehat{\lambda_{1}},\widehat{\lambda_{2}},\widehat{\alpha})$ is obtained by the inverse of the expected Fisher information matrix of $(\lambda_{1},\lambda_{2},\alpha)$ , evaluated at $(\lambda_{1},\lambda_{2},\alpha)=(\widehat{\lambda_{1}},\widehat{\lambda_{2}},% \widehat{\alpha})$ . In this case, the Fisher information could be approximated by the observed Fisher information matrix where its elements are given by the second derivatives of the log-likelihood function with respect to $\lambda_{1},\lambda_{2}$ and $\alpha$ locally at the obtained MLE’s. Hypothesis testing and confidence intervals for $\lambda_{1},\lambda_{2}$ and $\alpha$ could be obtained by using the asymptotical normality (see Lawless, 1982) of the MLEs $\widehat{\lambda_{1}},\widehat{\lambda_{2}}$ and $\widehat{\alpha}$ , that is,

$\displaystyle(\widehat{\lambda_{1}},\widehat{\lambda_{2}},\widehat{\alpha})% \sim N[(\widehat{\lambda_{1}},\widehat{\lambda_{2}},\widehat{\alpha}),I_{0}^{-% 1}]$ (14)

where $I_{0}$ is the observed Fisher information matrix.

3.1.2 Censored data

In many applications related to lifetime data there is the presence of censored data, that could be right, left or interval censoring. In this section, let us assume the presence of right censored data, that is, associated to each lifetime $X_{j},j=1,2$ , there is a fixed censoring time $C_{j}$ and the data are given by $T_{1}=\min(X_{1},C_{1})$ and $T_{2}=\min(X_{2},C_{2})$ . The likelihood function for the parameters of the BDGR distribution based on a sample of size n $(i=1,\ldots,n)$ has the dataset classified in four regions:

•
$R_{1}$ : Both, $X_{1i}$ and $X_{2i}$ , are complete observations;
•
$R_{2}$ : $X_{1i}$ are complete and $X_{2i}$ are censored;
•
$R_{3}$ : $X_{1i}$ are censored and $X_{2i}$ are complete;
•
$R_{4}$ : Both, $X_{1i}$ and $X_{2i}$ , are censored observations.

Thus, the likelihood function for $\lambda_{1},\lambda_{2}$ and $\alpha$ based on $n$ bivariate observations $\mathbf{t}_{i}=(t_{1i},t_{2i}),i=1,2,\ldots,n$ is given by,

$\displaystyle L(\lambda_{1},\lambda_{2},\alpha)=\prod_{i\in R_{1}}\text{Pr}(T_% {1i}=t_{1i},T_{2i}=t_{2i})\times\prod_{i\in R_{2}}\text{Pr}(T_{1i}=t_{1i},T_{2% i}>t_{2i})\times\prod_{i\in R_{3}}\text{Pr}(T_{1i}>t_{1i},T_{2i}=t_{2i})\times% \prod_{i\in R_{4}}\text{Pr}(T_{1i}>t_{1i},T_{2i}>t_{2i})$ (15)

where $\text{Pr}(T_{1i}=t_{1i},T_{2i}>t_{2i})=\text{Pr}(T_{1i}=t_{1i},T_{2i}>t_{2i}+1% )-\text{Pr}(T_{1i}=t_{1i},T_{2i}=t_{2i})$ and $\text{Pr}(T_{1i}>t_{1i},T_{2i}=t_{2i})=\text{Pr}(T_{1i}>t_{1i},T_{2i}=t_{2i})-% \text{Pr}(T_{1i}=t_{1i},T_{2i}=t_{2i})$ . Let us also define the following indicator variables of censoring,

$\displaystyle\begin{cases}\delta_{1i}=1\ \text{if}\ X_{1i}<C_{1i}\ \text{and 0% , for other part.}\\ \delta_{2i}=1\ \text{if}\ X_{2i}<C_{2i}\ \text{and 0, for other part.}\end{cases}$ (16)

where $i=1,2,\ldots,n;(C_{1i},C_{2i})$ are the right censoring times. From Eq. (16) are obtained the following results needed for the likelihood function:

(a)
For $\delta_{1i}=\delta_{2i}=1$ , we have,

$\displaystyle\text{Pr}(T_{1i}=t_{1i},T_{2i}=t_{2i})=\Bigg{\{}\frac{\left[1-% \bar{\alpha}^{2}h\left(\sqrt{t_{1i}^{2}+(t_{1i}+1)^{2}},\sqrt{t_{2i}^{2}+(t_{2% i}+1)^{2}}\right)\right]}{\left[1-\bar{\alpha}\,h(t_{1i}+1,t_{2i})\right]\left% [1-\bar{\alpha}\,h(t_{1i},t_{2i}+1)\right]}\times\frac{\alpha DR(t_{1i},% \lambda_{1})DR(t_{2i},\lambda_{2})}{\left[1-\bar{\alpha}\,h(t_{1i}+1,t_{2i}+1)% \right]\left[1-\bar{\alpha}h(t_{1i},t_{2i})\right]}\Bigg{\}}^{\delta_{1i}% \delta_{2i}}$
(b)
For $\delta_{1i}=1,\delta_{2i}=0$ , we have,

$\displaystyle\text{Pr}(T_{1i}=t_{1i},T_{2i}>t_{2i})=\Bigg{\{}\frac{\alpha h(t_% {1i},t_{2i}+1)}{1-\bar{\alpha}h(t_{1i},t_{2i}+1)}-\frac{\alpha h(t_{1i},t_{2i}% )}{1-\bar{\alpha}h(t_{1i},t_{2i})}\Bigg{\}}^{\delta_{1i}(1-\delta_{2i})}$
(c)
For $\delta_{1i}=0,\delta_{2i}=1$ , we have,

$\displaystyle\text{Pr}(T_{1i}=t_{1i},T_{2i}>t_{2i})=\Bigg{\{}\frac{\alpha h(t_% {1i}+1,t_{2i})}{1-\bar{\alpha}h(t_{1i}+1,t_{2i})}-\frac{\alpha h(t_{1i},t_{2i}% )}{1-\bar{\alpha}h(t_{1i},t_{2i})}\Bigg{\}}^{\delta_{2i}(1-\delta_{1i})}$
(d)
For $\delta_{1i}=0,\delta_{2i}=0$ , we have,

$\displaystyle\text{Pr}(T_{1i}=t_{1i},T_{2i}>t_{2i})=\Bigg{\{}\frac{\alpha h(t_% {1i},t_{2i})}{1-\bar{\alpha}h(t_{1i},t_{2i})}\Bigg{\}}^{(1-\delta_{2i})(1-% \delta_{1i})}$

Remark 1 The dependence between the two lifetimes is specified by the parameter $\alpha$ , where if $\alpha=1$ there is independence between the two lifetimes. The adequacy of the marginal distributions in presence of censored data could be verified by plots of the fitted survival functions and non-parametric Kaplan-Meier estimates for the survival function. In situations not considering censored data is used standard empirical estimates of the survival function.
3.2 A Bayesian approach

The Bayesian paradigm is based on specifying a probability distribution for the observed data $D$ given a vector of unknown parameters $\bm{\eta}$ (assuming $\bm{\eta}$ as a vector of random variables) providing a rational method for updating the new information using the Bayes’ rule given the prior distribution specifying the uncertainty about the parameter (see Ibrahim et al., 2005).

In the determination of the Bayes estimators for the unknown parameters of the BDGR model based on the squared error loss function, $K(\eta,a)=(\eta-a)^{2}$ , suppose that the parameters $\lambda_{i},i=1,2$ and $\alpha$ have independent Beta( $\alpha_{i},\beta_{i}$ ) and flat (improper) prior distributions given respectively, by,

$\displaystyle\pi_{i}(\lambda_{i})=\frac{\Gamma(\alpha_{i})\Gamma(\beta_{i})}{% \Gamma(\alpha_{i}+\beta_{i})}\lambda_{i}^{(\alpha_{i}-1)}(1-\lambda_{i})^{% \beta_{i}-1},\,\,0<\lambda_{i}<1,i=1,2$ $\displaystyle\pi_{3}(\alpha)\propto 1,\,\,\alpha>0$ (17)

The joint posterior density function for the parameters $\lambda_{1},\lambda_{2}$ and $\alpha$ of the proposed BDGR distribution is obtained directly from the Bayes formula, that is,

$\displaystyle\pi(\lambda_{1},\lambda_{2},\alpha\mid\text{data})=\frac{L(% \lambda_{1},\lambda_{2},\alpha)\prod_{i=1}^{2}\pi_{i}(\lambda_{i})\pi_{3}(% \alpha)}{\int_{0}^{1}\int_{0}^{1}\int_{0}^{\infty}L(\lambda_{1},\lambda_{2},% \alpha)\prod_{i=1}^{2}\pi_{i}(\lambda_{i})\pi_{3}(\alpha)\,d\lambda_{1}d% \lambda_{2}d\alpha}.$ (18)

Therefore, the Bayes estimator of any function of $\lambda_{1}$ , $\lambda_{2}$ and $\alpha$ , say $\mu(\lambda_{1},\lambda_{2},\alpha)$ , assuming the squared error loss function is given by,

$\displaystyle\widehat{\mu}_{B}=\frac{\int_{0}^{1}\int_{0}^{1}\int_{0}^{\infty}% \mu(\lambda_{1},\lambda_{2},\alpha)L(\lambda_{1},\lambda_{2},\alpha)\prod_{i=1% }^{2}\pi_{i}(\lambda_{i})\pi_{3}(\alpha)\,d\lambda_{1}d\lambda_{2}d\alpha}{% \int_{0}^{1}\int_{0}^{1}\int_{0}^{\infty}L(\lambda_{1},\lambda_{2},\alpha)% \prod_{i=1}^{2}\pi_{i}(\lambda_{i})\pi_{3}(\alpha)\,d\lambda_{1}d\lambda_{2}d% \alpha}.$ (19)

Given the difficulties in analytically finding the Bayesian estimators of interest given by expression Eq. (19), MCMC (Markov Chain Monte Carlo) methods are used to get the posterior summaries of interest (see, for example, Gelfand & Smith, 1990; Chib & Greenberg, 1995; Achcar & Leandro, 1998). In this way, assuming the BDGR distribution for the bivariate lifetime data, Monte Carlo estimators for the parameters $\lambda_{1}$ , $\lambda_{2}$ and $\alpha$ under the squared error loss function are obtained from the simulated samples for the joint posterior distribution $\pi(\lambda_{1},\lambda_{2},\alpha\mid\text{data})$ using MCMC methods, as for example, the Gibbs sampling or the Metropolis-Hastings algorithms. The steps for the Gibbs sampling algorithm are given by,

•

Step 1: Choose initial values $\lambda_{1}^{(0)},\lambda_{2}^{(0)}$ and $\alpha^{(0)}$ for $\lambda_{1},\lambda_{2}$ and $\alpha$ , respectively. Denote the values of $\lambda_{1},\lambda_{2}$ and $\alpha$ at the ith step by $\lambda_{1}^{(i)},\lambda_{2}^{(i)}$ and $\alpha^{(i)}$ ;

•

Step 2: Generate $\lambda_{1}^{(i+i)}$ , $\lambda_{2}^{(i+1)}$ and $\alpha^{(i+1)}$ from the conditional posterior distributions $\pi(\lambda_{1}\mid\text{data},\lambda_{2}^{(i)},\alpha^{(i)})$ , $\pi(\lambda_{2}\mid\text{data},\lambda_{1}^{(i+1)},\alpha^{(i)})$ and $\pi(\alpha\mid\text{data},\lambda_{1}^{(i+1)},\lambda_{2}^{(i+1)})$ respectively; repeat this procedure N times;

•

Step 3: Calculate the Monte Carlo Bayesian estimate of $\mu(\lambda_{1},\lambda_{2},\alpha)$ by ${\sum}_{i=B+1}^{N}\frac{\mu\left(\lambda_{1}^{(i)},\lambda_{2}^{(i)},\alpha^{(% i)}\right)}{(N-B)}$ where $B$ is the burn-in period.

The posterior summaries of interest are computed using the package R2jags (Su & Yajima, 2012) from the R software (R Core Team, 2016) considering a “burn-in sample” of size 10,000 to eliminate the effect of the initial values and a final Gibbs sample of size 2,000 taking every 100th sample from 200,000 simulated Gibbs samples. In addition, the convergence of the Gibbs Sampling algorithm was monitored using standard graphical methods, as the trace plots of the simulated samples.

4. A simulation study

This section reports the results of a simulation study carried out to assess the performance of the MLEs of the BDGR model assuming complete data. The simulation study was performed using the library maxLik of the R software and considering the BFGS optimization method. To simulate observations from the BDGR model, the marginal distribution of $X_{1}$ and the conditional distribution of $X_{2}$ given $X_{1}$ were used following the steps:

•
Step 1: Generate $U_{1}\sim\textit{Uniform(0,1)}$ and $U_{2}\sim\textit{Uniform(0,1)}$ ;
•
Step 2: Generate a value $x_{1}$ of $X_{1}$ from the marginal distribution of $X_{1}$ using the inverse transformation method;
•
Step 3: Generate a value $x_{2}$ of $X_{2}$ using the inverse transformation method again based on the conditional distribution of $X_{2}$ given $X_{1}=x_{1}$ ;

The conditional distribution of $X_{j}$ given $X_{i},i,j=1,2$ and $i\neq j$ is given by,

$\displaystyle\text{Pr}(X_{j}\mid X_{i}=x_{i})=\frac{[1-\bar{\alpha}^{2}h(k_{i}% ,k_{j})]\left[1-\bar{\alpha}h(x_{i},x_{j})\right]^{-1}}{\left[1-\bar{\alpha}\,% h(x_{i}+1,x_{j})\right][DR(x_{j},\lambda_{j})]^{-1}}\times\frac{[1-\bar{\alpha% }(DR(x_{i},\lambda_{i})+\bar{\alpha}\lambda_{i}^{(x_{i}+1)^{2}}\lambda_{i}^{x_% {i}^{2}})]}{\left[1-\bar{\alpha}\,h(x_{i},x_{j}+1)\right]\left[1-\bar{\alpha}% \,h(x_{i}+1,x_{j}+1)\right]}$

where $k_{a}=\sqrt{x_{a}^{2}+(x_{a}+1)^{2}}$ .
•
Step 4: Return $X=(X_{1},X_{2})$ .

We performed the simulation study under three scenarios considering the following parameter values assumed for better computational stability: ( $\lambda_{1},\lambda_{2},\alpha$ ) $=$ ( $\lambda_{1}$ , 0.90, 2.00) where $\lambda_{1}=$ 0.40, 0.50, 0.70, 0.80 for the first scenario; ( $\lambda_{1},\lambda_{2},\alpha$ ) $=$ (0.95, $\lambda_{2}$ , 2.50) where $\lambda_{2}=$ 0.40, 0.50, 0.70, 0.80 for the second scenario; and ( $\lambda_{1},\lambda_{2},\alpha$ ) $=$ (0.95, 0.97, $\alpha$ ) where $\alpha=$ 0.50, 1.00, 1.50, 2.00 for the third scenario. We also considered the sample sizes $n=$ 10, $\ldots$ ,100, each one with 10,000 Monte Carlo replications.

For each scenario, based on the average of the 10,000 simulated parameter components of the vector of parameters ( $\lambda_{1},\lambda_{2},\alpha$ ), the biases and the RMSE were computed using the expressions:

$\displaystyle\textit{BIAS}(\widehat{\bm{\Psi}})=\frac{1}{N}\sum_{i=1}^{N}(% \widehat{\bm{\Psi}}_{i}-\bm{\Psi}),\,\textit{RMSE}(\widehat{\bm{\Psi}})=\sqrt{% \frac{1}{N}\sum_{i=1}^{N}(\widehat{\bm{\Psi}}_{i}-\bm{\Psi})^{2}}$

where $N=$ 10,000 is the number of simulations and $\bm{\Psi}$ denotes each parameter $\lambda_{1},\lambda_{2}$ or $\alpha$ . The obtained simulation results for each scenario are illustrated, respectively, in Figs 2–4.

Figure 2.
The biases (upper panels) and RMSEs (lower panels) for the BDGR distribution assuming ( $\lambda_{1},\lambda_{2},\alpha$ ) $=$ ( $\lambda_{1}$ , 0.90, 2.00) where $\circ:\,\lambda_{1}=$ 0.40; $\triangle:\,\lambda_{1}=$ 0.50; $+$ : $\lambda_{1}=$ 0.60; $\times:\,\lambda_{1}=$ 0.70.

Figure 3.
The biases (upper panels) and RMSEs (lower panels) for the BDGR distribution assuming ( $\lambda_{1},\lambda_{2},\alpha$ ) $=$ (0.95, $\lambda_{2}$ , 2.50) where $\circ:\,\lambda_{2}=$ 0.35; $\triangle:\,\lambda_{2}=$ 0.40; $+$ : $\lambda_{2}=$ 0.45; $\times:\,\lambda_{2}=$ 0.50.

From the simulation results illustrated in Figs 2–4, it is possible to conclude that,

i.)
For all considered scenarios, the optimization method converged successful in the determination of the MLEs estimates with no instability for the determination of the biases and RMSEs for all parameters;
ii.)
For all considered scenarios, the biases and RMSEs tends to zero when the sample size increases. The convergence to zero is much faster in the third scenario ( $\lambda_{1},\lambda_{2}$ fixed);
iii.)
The parameters $\lambda_{1}$ and $\lambda_{2}$ have negative biases for all considered scenarios. The parameter $\alpha$ has a positive bias for all scenarios;
iv.)
The parameters $\lambda_{1}$ and $\lambda_{2}$ have small values for the biases and RMSEs; however the parameter $\alpha$ has high values for the biases and RMSEs;
v.)
The smallest values for the biases were obtained for the third scenario; the smallest values for the RMSE were also obtained for the third scenario;
vi.)
From the simulation results, we concluded that the BDGR distribution has better asymptotically non-biased estimation in the third scenario since $E(\lambda_{i})\approx\lambda_{i},i=$ 1, 2 and $E(\alpha)\approx\alpha$ in this scenario; for the other scenarios, we need a sample size n larger than 100 to get the results $E(\lambda_{i})\approx\lambda_{i}$ and $E(\alpha)\approx\alpha$ ;
vii.)
It is important to point out that the simulation study also could be made using a Bayesian approach assuming different prior distributions for the parameters of the BDGR distribution.
viii.)
Based on the results of these simulation studies we conclude that the BDGR distribution could be used as an alternative to other existing discrete bivariate distributions (such as the Basu-Dhar bivariate geometric distribution introduced by Basu & Dhar, 1995) to describe bivariate lifetimes with good accuracy in applications.

5. A medical example with censored observations and covariates

To illustrate the usefulness of the proposed methodology, we present in this section, the analysis of a real medical dataset related to lung cancer assuming the bivariate discrete generalized Rayleigh distribution. The dataset is introduced by Ding et al. (2017) and corresponds to the lifetimes of Chinese patients with pathologically confirmed lung cancer who received EGFR, KRAS, and BARF mutation tests at the Thoracic Cancer Institute, Tongji University from January 2012 to April 2016.

For the statistical analysis, we assumed as lifetimes the overall survival times (calculated from the date of lung cancer diagnosis to death from any reason or censored at the last follow-up date), and the progression-free survival times (the times from the treatment start time until the date of systemic progression or death). The data set consists of $n=$ 28 patients with not censored observations for the overall survival times and 4 censored observations for the progression-free survival times. The non-parametric estimators for the means obtained from the kaplan and Meier (1958) non-parametric estimators for the survival functions are given, respectively, by 16.57 months for the overall survival time ( $T_{1}$ ) and 5.71 months for the progression-free survival time ( $T_{2}$ ). We also considered in the statistical analysis of the data, a prognostic factor (smoking status) and an oncological factor (cancer stage) as two covariates possibly affecting the survival times of the patients.

5.1 A Bayesian analysis not considering the presence covariates

As a first analysis, the presence of the prognostic and oncology factors (smoking status and cancer stage) were not considered. The fit of BDGR distribution was compared to the fit of a Basu-Dhar bivariate geometric distribution (BDBG) (Basu & Dhar, 1995); to the fit of a Arnold bivariate geometric distribution (ABG) (Arnold, 1975); and to the fit of a generalized bivariate geometric distribution (GBG) (Gómez-Déniz et al., 2017).

The parameters of the models were estimated under a Bayesian approach assuming independent $\text{Beta}(1,1)$ prior distributions for the parameters $\lambda_{i},i=$ 1, 2, 12 and a gamma prior distribution with hyperparameter values $\left(\overline{T_{1}}+\overline{T_{2}}\right)/2$ and $\left[\text{Var}(T_{1})+\text{Var}(T_{2})\right]/2$ where $\overline{T_{1}}=$ 5.0714, $\overline{T_{2}}=$ 15.3214, $\text{Var}{(T_{1})}=$ 10.2910 and $\text{Var}{(T_{2})}=$ 36.5966 for the parameter $\alpha$ , that is, we used empirical Bayesian methods in the elicitation of the prior distribution for the dependence parameter $\alpha$ (see Carlin & Louis, 2000). The DIC criteria was considered to discriminate the best model among all proposed models (smaller values indicate better models, see Spiegelhalter et al., 2014). From the obtained results presented in Table 2, we conclude that the BDGR distribution is the best model fitted by the dataset.

Table 2
Bayesian estimates assuming bivariate discrete models for the lung cancer dataset

Model	Param.	Mean	S.D.	95% Conf. Int.	Model	Param.	Mean	S.D.	95% Cred. Int.
BDGR	$\alpha$	0.4309	0.1139	(0.2386, 0.6909)	GBG	$\alpha$	0.8805	0.1040	(0.8255, 0.9973)
	$\lambda_{1}$	0.9848	0.0041	(0.9756, 0.9918)		$\lambda_{1}$	0.8606	0.0289	(0.8417, 0.9134)
	$\lambda_{2}$	0.9978	0.0006	(0.9963, 0.9988)		$\lambda_{2}$	0.9413	0.0115	(0.9338, 0.9626)
BDBG	$\lambda_{1}$	0.7982	0.0335	(0.7768, 0.8570)	ABG	$\lambda_{1}$	0.1897	0.0311	(0.1675, 0.2533)
	$\lambda_{2}$	0.9697	0.0188	(0.9553, 0.9988)		$\lambda_{2}$	0.0583	0.0114	(0.0506, 0.0817)
	$\lambda_{12}$	0.9692	0.0187	(0.9551, 0.9984)

DIC ${}_{\textit{BDGR}}=$ 313.75; DIC ${}_{\textit{ABG}}=$ 329.70; DIC ${}_{\textit{BDBG}}=$ 330.20; DIC ${}_{\textit{GBG}}=$ 354.60.

Figure 4.

The biases (upper panels) and RMSEs (lower panels) for the BDGR distribution assuming ( $\lambda_{1},\lambda_{2},\alpha$ ) $=$ (0.95, 0.97, $\alpha$ ) where $\circ:\,\alpha=$ 0.50; $\triangle:\,\alpha=$ 0.75; $+$ : $\alpha=$ 1.00; $\times:\,\alpha=$ 1.25.

Figure 5 presents the Kaplan-Meier marginal survival plots for both lifetimes and also the estimated survival BDGR functions since this model was verified to be the best fitted model for the dataset. From the graphs of Fig. 5, we conclude that the proposed model has a good accuracy to predict both survival times.

The suitability of the BDGR distribution can also be verified from the marginal distributions which must also be generalized Rayleigh distributions (see Balakrishnan & Ristić, 2016; Ristić et al., 2018). In our application, this is observed from graphs of the Kaplan and Meier non-parametric estimators of the survival functions and the Bayesian estimators of the marginal generalized Rayleigh survival functions for both lifetime times.

In addition, the estimator of the correlation coefficient for the proposed model was compared to the empirical correlation coefficient obtained using the Iterative Multiple Imputation (IMI) method (for more details, see Schemper et al., 2013) available in the package SurvCorr from R software. The empirical correlation coefficient was computed using the function survcorr setting the argument MCMCSteps $=$ 210,000. Assuming the BDGR under a Bayesian approach, the Monte Carlo estimator for the correlation coefficient based on the simulated Gibbs samples is given by $\hat{\rho}_{\textit{BDGR}}=$ 0.8726 which is very close to the IMI correlation coefficient given by $\hat{\rho}_{\textit{IMI}}=$ 0.9328, that is, the proposed model have a good accurate inference result for the correlation coefficient even it is based on a infinite series.

5.2 A Bayesian analysis in the presence of the covariates

In the statistical analysis in presence of the prognostic and oncology factors assuming the BDGR distribution, we assumed logistic regression models for the parameters $\lambda_{1}$ and $\lambda_{2}$ given by,

$\displaystyle\text{logit}(\lambda_{1i})=\beta_{10}+\beta_{11}\,\text{smoker}_{% i}+\beta_{12}\,\text{stage}_{i}$ $\displaystyle\text{logit}(\lambda_{2i})=\beta_{20}+\beta_{21}\,\text{smoker}_{% i}+\beta_{22}\,\text{stage}_{i}$

where $i=$ 1, $\ldots$ , 28; $\text{logit}(p)=\log\left(p/(1-p)\right)$ ; $\text{smoker}_{i}=$ 0 (never-smoker) and $\text{smoker}_{i}=$ 1 (former/current smoker); $\text{stage}_{i}=$ 0 (stage at diagnosis IIIB) and $\text{stage}_{i}=$ 1 (stage at diagnosis IV). For a Bayesian analysis of this regression model we assumed normal prior distributions for the regression parameters with zero mean and large variance (approximately non-informative priors) and the same gamma prior distribution considered in Section 5.1 for the parameter $\alpha$ . The inference results are reported in Table 3.

Table 3
Maximum likelihood and Bayesian estimates in presence of prognostic and oncology factors for the lung cancer data set

Param.	Classical approach			Bayesian approach
	MLE	S.E.	95% Conf. Int.	Mean	S.D.	95% Cred. Int.
$\alpha$	0.6758	0.3063	(0.0755, 1.2761)	0.4342	0.1195	(0.2285, 0.6971)
$\beta_{10}$	2.0580	0.7146	(0.6575, 3.4585)	$-$ 0.1184	1.7907	( $-$ 3.7383, 3.4465)
$\beta_{11}$	1.6732	0.5933	(0.5104, 2.8360)	1.6135	0.6642	(0.3050, 2.9042)
$\beta_{12}$	0.8725	0.5111	( $-$ 0.1291, 1.8742)	0.9004	0.5579	( $-$ 0.1979, 2.0299)
$\beta_{20}$	4.8502	0.6927	(3.4924, 6.2080)	3.8752	1.6987	(0.4833, 7.1144)
$\beta_{21}$	0.9195	0.5638	( $-$ 0.1856, 2.0247)	0.8446	0.6434	( $-$ 0.4747, 2.1223)
$\beta_{22}$	0.4400	0.4753	( $-$ 0.4917, 1.3716)	0.4860	0.5083	( $-$ 0.4628, 1.4760)
Criteria	AIC $=$ 314.20			DIC $=$ 313.18

Figure 5.

Kaplan-Meier estimators versus Bayesian estimated survival functions (including the 95% credible interval) for the marginal survival functions $T_{1}$ and $T_{2}$ (blue line: fitted model; black line: Kaplan-Meier).

From the results of Table 3, we observed that the Bayesian estimates are very close to the MLEs. Furthermore, we concluded that only the covariate smoker has significant effect on the overall survival time, since the 95% confidence and credible intervals for the regression parameter $\beta_{11}$ does not contain the zero value. From the results of Table 3 it is also observed that the parameter $\beta_{10}$ has very different classical and Bayesian inference results possible due to the effect of a very non-informative prior for this regression parameter (large 95% credible interval for the parameter). Under the classical approach the obtained inference is obtained under asymptotical results, although there is a small sample size (inference possible not very accurate).

6. Concluding remarks

The use of existing bivariate discrete distributions could be a good alternative to analyze bivariate lifetime data in presence of censored data and covariates. However, there are few discrete bivariate lifetime distributions introduced in the literature. In general, bivariate lifetime datasets usually are analyzed using standard continuous bivariate lifetime distributions as the popular Block and Basu bivariate exponential distribution (see Block & Basu, 1974) or the Marshall-Olkin bivariate exponential distribution (see Marshall & Olkin, 1967a, 1967b).

In this study, we introduced a new bivariate distribution, denoted as bivariate discrete generalized Rayleigh (BDGR) distribution, obtained using the Marshall and Olkin (1997) method to add a new parameter to the survival function of the discrete Rayleigh distribution proposed by Roy (2004) in order to propose a more flexible joint survival function as an alternative to existing discrete models as the popular Arnold (see Arnold, 1975) and Basu-Dhar (see Basu & Dhar, 1995) bivariate geometric distributions to analyze bivariate discrete lifetime data in presence of censored data and covariates. Some properties of this new distribution were also discussed in this study and an extension to multivariate case was provided.

An extensive simulation study was performed to verify the effectiveness of the maximum likelihood method assuming different fixed values for the parameters of the model and different sample sizes. The results obtained from Monte Carlo studies showed that the biases and RMSEs of the BDGR distribution are asymptotically non-biased and tends to zero when the sample size increases even assuming negative values for $\lambda_{i},i=$ 1, 2 in some scenarios.

In the application with real lifetime data presented in this study, we observed that, with the use of the BDGR distribution, it is possible to obtain in a simple way the inferences of interest for bivariate lifetime datasets in presence or not of covariates and censored data with small computational costs as compared to many existing bivariate parametric lifetime distributions introduced in the literature or bivariate models derived from copula functions for continuous bivariate lifetime data (see, for example, Achcar et al., 2016). The identification of important covariates was also easily obtained assuming the BDGR distribution even using non-informative priors for the parameters of the model, under a Bayesian approach. These results could be of great interest for the search of appropriate bivariate lifetime distributions especially in engineering, medical studies or other areas of interest.

Footnotes

Acknowledgments

The authors are grateful to the referees for their valuable suggestions that significantly improved the manuscript.

References

Achcar

J. A.

, & Leandro

R. A

. (1998). Use of Markov Chain Monte Carlo methods in a bayesian analysis of the block and basu bivariate exponential distribution. Annals of the Institute of Statistical Mathematics, 50(3), 403-416.

Achcar

J. A.

Martinez

E. Z.

, & Cuevas

J. R. T

. (2016). Bivariate lifetime modeling using copula functions in presence of mixture and non-mixture cure fraction models, censored data and covariates. Model Assisted Statistics and Applications, 11(4), 261-276.

Arnold

B. C

. (1975). A characterization of the exponential distribution by multivariate geometric compounding. Sankhyā: The Indian Journal of Statistics, Series A, pp. 164-173.

Arnold

B. C.

, & Strauss

. (1988). Bivariate distributions with exponential conditionals. Journal of the American Statistical Association, 83(402), 522-527.

Balakrishnan

, & Lai

. (2009). Continuous Bivariate Distributions. Springer, New York.

Balakrishnan

, & Ristić

M. M

. (2016). Multivariate families of gamma-generated distributions with finite or infinite support above or below the diagonal. Journal of Multivariate Analysis, 143, 194-207.

Basu

A. P.

, & Dhar

. (1995). Bivariate geometric distribution. Journal of Applied Statistical Science, 2(1), 33-44.

Block

H. W.

, & Basu

. (1974). A continuous, bivariate exponential extension. Journal of the American Statistical Association, 69(348), 1031-1037.

Carlin

B. P.

, & Louis

T. A

. (2000). Empirical Bayes: Past, present and future. Journal of the American Statistical Association, 95(452), 1286-1289.

10.

Chib

, & Greenberg

. (1995). Understanding the Metropolis-Hastings algorithm. The american statistician, 49(4), 327-335.

11.

Ding

Zhang

Jiang

Zhao

, & Zhou

. (2017). Clinicopathologic characteristics and outcomes of Chinese patients with non-small-cell lung cancer and BRAF mutation. Cancer Medicine, 6(3), 555-562.

12.

Downton

. (1970). Bivariate exponential distributions in reliability theory. Journal of the Royal Statistical Society. Series B (Methodological), pp. 408-417.

13.

Freund

J. E

. (1961). A bivariate extension of the exponential distribution. Journal of the American Statistical Association, 56(296), 971-977.

14.

Gelfand

A. E.

, & Smith

A. F

. (1990). Sampling-based approaches to calculating marginal densities. Journal of the American statistical association, 85(410), 398-409.

15.

Gómez-Déniz

Ghitany

, & Gupta

R. C

. (2017). A bivariate generalized geometric distribution with applications. Communications in Statistics-Theory and Methods, 46(11), 5453-5465.

16.

Gumbel, E. J . (1960). Bivariate exponential distributions. Journal of the American Statistical Association, 55(292), 698-707.

17.

Hanagal

, & Ahmadi

. (2008). Estimation of parameters by EM algorithm in bivariate exponential distribution based on censored samples. Econ. Quality Control, 23(2), 257-66.

18.

Hanagal

D. D

. (2006). Bivariate Weibull regression model based on censored samples. Statistical Papers, 47(1), 137-147.

19.

Hawkes

A. G

. (1972). A bivariate exponential distribution with applications to reliability. Journal of the Royal Statistical Society. Series B (Methodological), pp. 129-131.

20.

Hougaard

. (1986). A class of multivariate failure time distributions. Biometrika, pp. 671-678.

21.

Ibrahim

J. G.

Chen

M.-H.

, & Sinha

. (2005). Bayesian survival analysis. Wiley Online Library.

22.

Kaplan

E. L.

, & Meier

. (1958). Nonparametric estimation from incomplete observations. Journal of the American statistical association, 53(282), 457-481.

23.

Kemp

A. W

. (2013). New discrete appell and humbert distributions with relevance to bivariate accident data. Journal of Multivariate Analysis, 113, 2-6. Kemp and Papageorgiou, 1982 Kemp, C., & Papageorgiou, H. (1982). Bivariate Hermite distributions. Sankhyā: The Indian Journal of Statistics, Series A, pp. 269-280. Kocherlakota, 1995 Kocherlakota, S. (1995). Discrete bivariate weighted distributions under multiplicative weight function. Communications in statistics-theory and methods, 24(2), 533-551. Kocherlakota and Kocherlakota, 1992 Kocherlakota, S., & Kocherlakota, K. (1992). Bivariate discrete distributions. Wiley Online Library.

24.

Kumar

C. S

. (2008). A unified approach to bivariate discrete distributions. Metrika, 67(1), 113-123.

25.

Kundu,

. (2014). Geometric skew normal distribution. Sankhya B, 76(2), 167-189.

26.

Kundu

. (2017). Multivariate geometric skew-normal distribution. Statistics, 51(6), 1377-1397.

27.

Kundu

Balakrishnan

, & Jamalizadeh

. (2010). Bivariate birnbaum–saunders distribution and associated inference. Journal of Multivariate Analysis, 101(1), 113-125.

28.

Kundu

, & Gupta

R. D

. (2009). Bivariate generalized exponential distribution. Journal of Multivariate Analysis, 100(4), 581-593.

29.

Kundu

, & Nekoukhou

. (2018). Univariate and bivariate geometric discrete generalized exponential distributions. Journal of Statistical Theory and Practice, 12(3), 595-614.

30.

Lai

C.-D

. (2006). Constructions of discrete bivariate distributions. In Advances in Distribution Theory, Order Statistics, and Inference, pp. 29-58. Springer.

31.

Lawless

J. F

. (1982). Statistical models and methods for lifetime data. John Wiley & Sons.

32.

Lee

, & Cha

J. H

. (2014). On construction of general classes of bivariate distributions. Journal of Multivariate Analysis, 127, 151-159.

33.

Marshall

A. W.

, & Olkin

. (1967a). A generalized bivariate exponential distribution. Journal of Applied Probability, 4(2), 291-302.

34.

Marshall

A. W.

, & Olkin

. (1967b). A multivariate exponential distribution. Journal of the American Statistical Association, 62(317), 30-44.

35.

Marshall

A. W.

, & Olkin

. (1997). A new method for adding a parameter to a family of distributions with application to the exponential and Weibull families. Biometrika, 84(3), 641-652.

36.

MirMostafaee

Mahdizadeh

, & Lemonte

A. J

. (2017). The Marshall – Olkin extended generalized Rayleigh distribution: Properties and applications. Communications in Statistics-Theory and Methods, 46(2), 653-671.

37.

Nekoukhou

, & Kundu

. (2017). Bivariate discrete generalized exponential distribution. Statistics, 51(5), 1143-1158.

38.

Popović

B. V.

, & Genç

A. İ

. (2018). On extremes of two-dimensional Student-t distribution of the Marshall–Olkin type. Mediterranean Journal of Mathematics, 15(4), 153.

39.

R Core Team (2016). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.

40.

Ristić

M. M.

Popović

B. V.

Zografos

, & Balakrishnan

(2018). Discrimination among bivariate beta-generated distributions. Statistics, 52(2), 303-320.

41.

Roy

. (2004). Discrete Rayleigh distribution. IEEE Transactions on Reliability, 53(2), 255-260.

42.

Sarkar

S. K

. (1987). A continuous bivariate exponential distribution. Journal of the American Statistical Association, 82(398), 667-675.

43.

Schemper

Kaider

Wakounig

, & Heinze

. (2013). Estimating the correlation of bivariate failure times under censoring. Statistics in medicine, 32(27), 4781-4790.

44.

Spiegelhalter

D. J.

Best

N. G.

Carlin

B. P.

, & Van Der Linde

. (2014). The deviance information criterion: 12 years on. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(3), 485-493.

45.

Y.-S.

, & Yajima

. (2012). R2jags: A package for running jags from r. R package version 0.03-08, URL http://CRAN. R-project. org/package= R2jags.

A new flexible bivariate discrete Rayleigh distribution generated by the Marshall-Olkin family

Abstract

Keywords

1. Introduction

2.1 Model description

.

.

Proof..

.

Proof..

.

3.1 Maximum likelihood approach

3.1.1 Complete data

5.1 A Bayesian analysis not considering the presence covariates

Table 2 Bayesian estimates assuming bivariate discrete models for the lung cancer dataset

Table 3 Maximum likelihood and Bayesian estimates in presence of prognostic and oncology factors for the lung cancer data set

Footnotes

Acknowledgments

References

Table 2
Bayesian estimates assuming bivariate discrete models for the lung cancer dataset

Table 3
Maximum likelihood and Bayesian estimates in presence of prognostic and oncology factors for the lung cancer data set