A bivariate geometric distribution is an extension to a univariate geo­metric distribution where the occurrence of three different types of events is considered. Many statisticians have studied and given different forms of a bivariate geometric distribution. In this paper, we considered the form given by Phatak & Sreehari.¹ We estimated the parameters of this distribution under three different models using maximum likelihood estimation (mle) and derived deviances as the goodness of fit statistics for testing the parameters and deviance difference for comparing two models. Using simulated data we found that the deviance measure works well to test a reduced model against a full model.

Keywords: bivariate geometric distribution, deviance, deviance difference

Many situations in real world cannot be described by a single variable. Simul­taneous occurrence of multiple events warrants multivariate distributions. For instance, univariate geometric distribution can represent occurrence of failure of one component of a system. However, to study systems with several com­ponents that may have different types of failures, such as twin engines of an airplane or the paired organ in a human body, bivariate geometric distributions are suitable. Bivariate geometric distribution has increasingly important roles in various fields, including reliability and survival analysis. There are different forms of a bivariate geometric distribution. Phatak & Sreehari¹ pro­vided a form of the bivariate geometric distribution which is considered here. They introduced a form of probability mass function which take into considera­tion of three different types of events. There are other forms which can be seen in Nair & Nair,² Hawkes,³ Arnold et al.⁴ and Sreehari & Vasudeva.⁵ Basu & Dhar⁶ proposed a bivariate geometric model which is analogous to bivariate exponential model developed by Marshal & Olkin.⁷ Characterization results are developed by Sun & Basu,⁸ Sreehari,⁹ and Sreehari & Vasudeva.⁵

Omey & Minkova¹⁰ considered the bivariate geometric distribution with negative correlation coefficient and analyzed some properties, probability generating function, probability mass function, moments and tail probabilities. Krishna & Pundir,¹¹ studied the plausibility of a bivariate geometric distribution as a reliability model. They derived the maximum likelihood esti­mators and Bayes estimators of the parameters and various reliability charac­teristics. They also compared these estimators using Monte-Carlo simulation.

In this paper, the parameters of a saturated model, reduced model and gen­eralized linear model (glm) for a bivariate geometric distribution are estimated using the maximum likelihood method. We also derived deviances as the good­ness of fit statistics for testing parameters corresponding to these models and deviance difference to compare two related models in order to determine which model fits the data well. Rest of the paper is organized as follows: section 2 describes the univariate geometric distribution, section 3 presents the bivariate geometric distribution, section 4 presents hypothesis testing, section 5 discusses a numerical example with simulated data and section 6 has the conclusion.

The probability mass function (pmf) of a random variable Y which follows a geometric distribution with probability of success p can be written as,

$P\left(Y=y\right)=p{\left(1-p\right)}^y,y = 0,1,2,\mathrm{...}; 0<p<1, 0<q=1-p<1$ .

The moment generating function can be given by,

$M_Y\left(t\right)=\frac{p}{1-q{e}^t}$

the mean and the variance of this distribution are

$E\left(Y\right)={\mu }_Y=\frac{1-p}{p}=\frac{q}{p} and V ar\left(Y\right)=\frac{1-p}{p^2}=\frac{q}{p^2}$

An extension to the univariate geometric distribution is the bivariate geo­metric distribution which is discussed in the next section.

The joint probability mass function of a bivariate geometric distribution can be obtained by the product of a marginal and a conditional distribution, introduced by Phatak & Sreehari.¹ They considered a process from which the units could be classified as good, marginal and bad with probabilities $q_1$ , $q_2$  and $q_3=\left(1-q_1-q_2\right)$  respectively. They proposed that the probability mass function of observing the first bad unit after several good and marginal units are passed as follows:

$P\left(Y_1=y_1,Y_2=y_2\right)=\left(\begin{array}{l}{y}_1+y_2\\ y_1\end{array}\right)q_1{}^{y_1}{q}_2{}^{y_2}\left(1-q_1-q_2\right),y_1,y_2=0,1,2,\mathrm{...};$  (1)

$0 < q_1+q_2< 1$ .

Here $Y_1$  and $Y_2$  denote the number of good and marginal units respectively before the first bad unit is observed.

The marginal distribution of $Y_1$ is a geometric distributions with probability of success $\left(\frac{1-q_1-q_2}{1-q_2}\right)$ , and can be written as follows,

$P\left(Y_1=y_1\right)=\left(\frac{1-q_1-q_2}{1-q_2}\right){\left(\frac{q_1}{1-q_2}\right)}^{y_1}, y_1=0,1,2,\mathrm{...}$

The conditional distribution of $Y_2$ given $Y_1$  is

$P\left(Y_2=y_2/Y_1=y_2\right)= \left(\begin{array}{l}{y}_1+y_2\\ y_2\end{array}\right)q_{{}_2}{}^{y_2}{\left(1-q_2\right)}^{y_1+1}, y_1,y_2= 0,1,2,\mathrm{...}$  (3)

The product of the marginal distribution of $Y_1$  in equation (2) and the con­ditional distribution of $Y_2$  given $Y_1$  in equation (3) gives the mass function of bivariate geometric distribution in equation (1).

Maximum likelihood estimation

Estimation of parameters in the absence of regressors

In order to find the maximum likelihood estimators (mle)s from a saturated model (parameters are different for each pair of observations), it suffices to con­sider the likelihood functions based on the marginal and conditional mass func­tions. Let $Y_1,\mathrm{....},Y_n$  be independent random vectors each having bivariate geo­metric distribution with different pairs of parameters $\left(q_{1i},q_{2i}\right)$ for $i=1,2,\mathrm{....},n.$ .

The log likelihood function based on the conditional distribution of $Y_2$  given $Y_1$  can be written as follows using (3):

$l= {\displaystyle \sum _{i=1}^n\left[\ln \left(1-q_{1i}-q_{2i}\right)-\ln \left(1-q_{2i}\right)+y_{1i}\ln q_{1i}-y_{1i}\ln \left(1-q_{2i}\right)\right]}$  (4)

Differentiating (4) with respect to $q_{2i}$ and setting it equal to zero, we get the mle of $q_{2i}$  as,

${\widehat{q}}_{2i}=\frac{y_{2i}}{y_{1i}+y_{2i}+1}$  (5)

The log likelihood function based on the marginal distribution of $Y_1$  from (2) is,

$l= {\displaystyle \sum _{i=1}^n\left[\ln \left(1-q_{1i}-q_{2i}\right)-\ln \left(1-q_{2i}\right)+y_{1i}\ln q_{1i}-y_{1i}\ln \left(1-q_{2i}\right)\right]}$

Differentiating (6) with respect to $q_{1i}$ and setting it equal to zero, the mle of $q_{1i}$ can be derived as,

${\widehat{q}}_{1i}=\frac{y_{1i}}{y_{1i}+y_{2i}+1}$  (7)

Here, ${\widehat{q}}_{1i}$ and $q_{1i}$ are the maximum likelihood estimators of $q_{1i}$ and $q_{2i}$ , $i=1,\mathrm{...},n$ respectively under the saturated model.

Similarly the maximum likelihood estimators (mle)s from a reduced model (parameters are the same for each pair of observations) can be obtained as:

${\widehat{q}}_2=\frac{{\stackrel{\leftharpoonup }{y}}_2}{{\stackrel{\rightharpoonup }{y}}_1+{\stackrel{\leftharpoonup }{y}}_2+1}$  (8)

${\widehat{q}}_1=\frac{{\stackrel{\leftharpoonup }{y}}_1}{{\stackrel{\rightharpoonup }{y}}_1+{\stackrel{\leftharpoonup }{y}}_2+1}$  (9)

Where ${\widehat{q}}_1$ and ${\widehat{q}}_2$  are the maximum likelihood estimators of $q_1$  and $q_2$ respectively under the reduced model.

Estimation of parameters in the presence of regressors: In the presence of regressors, one can employ a generalized linear model and hence estimate the parameters in terms of the estimated model parameters. The conditional distribution of $Y_2$  given $Y_1$  in (3) can be set as exponential family representation as follows,

$P\left(Y_2=y_2/Y_1=y_1\right) = \exp \left[y_2\ln q_2-\left\{-\left(y_1+1\right)\ln \left(1-q_2\right)\right\}+\ln \frac{\left(y_1+y_2\right)!}{y_1!y_2!}\right]$

Here the natural parameter and the function of the natural parameter respec­tively are,

$\theta = \ln q_2$

$b\left(\theta \right) = - \left(y_1+1\right)\ln \left(1-q_2\right)$

Thus the mean of the conditional distribution of $Y_2$  given $Y_1$ is

${\mu }_i= E\left[Y_2/Y_1=y_1\right]/= b\text{'}\left(\theta \right)=\frac{y_1+1}{1-q_2}$

A generalized linear model based on the conditional distribution of $Y_2$ given $Y_1$ can be written as,

$g\left({\mu }_i\right) = ln{\mu }_i=/ n ={\displaystyle \sum _{j=1}^p{x}_{2ij}{\beta }_{2i}; i = } 1,2,\mathrm{....},n; n>p$

Since, $Y_2$ represents the number of trials before a certain event can occur it is considered as count response, the linear predictor can be written as the log­arithm of the mean ${\mu }_i$ . Thus the conditional link function can be expressed as,

$g\left({\mu }_i\right)=\ln \frac{y_{1i}+1}{1-q_{2i}}={\displaystyle \sum _{j=1}^p{x}_{2ij}{\beta }_{2i}}$

$\Rightarrow$ ${\widehat{q}}_{2i}= 1- \frac{y_{1i}+1}{\exp \left\{{\displaystyle {\sum }_{j=1}^p{x}_{2ij}{\beta }_{2j}}\right\}}$  (10)

Here, $\beta$  is an element of the matrix $\beta$ corresponding to the covariate $x_{2ij}$ which represents the effect of covariate to the mean responses through the link function $g\left({\mu }_i\right)=n_i$ .

Differentiating (6) again with respect to $q_{1i}$ , setting it to zero and using (10) we get,

${\widehat{q}}_{1i}=\frac{y_{1i}}{\exp \left\{{\displaystyle \sum _{j=1}^p{x}_{2ij}{\beta }_{2j}}\right\}}$ (11)

In order to test the identical parameter assumption across each pair of observed data, we derived deviance as a goodness of fit statistics. Additional deviance statistics are derived for generalized linear model (glm) to compare two nested glms.

Deviance for reduced model with identical parameter assumption

The log likelihood function for the saturated model can be written using (1) and the maximum likelihood estimates of the parameters $q_{1i}$ and $q_{2i}$ from equations (5) and (7) respectively as follows,

$l\left(b_{\max };y\right)={\displaystyle \sum _{i=1}^n\left[y_{1i}\ln y_{1i}-y_{1i}\ln \left(y_{1i}+y_{2i}+1\right)+y_{2i}\ln y_{2i}-y_{2i}\ln \left(y_{1i}+y_{2i}+1\right)-\ln \left(y_{1i}+y_{2i}+1\right)+\ln \left(y_{1i}+y_{2i}\right)!-\ln y_{1i}!-\ln y_{2i}!\right]}$ (12)

Similarly, the log likelihood function of the reduced model can be written using (1) and the maximum likelihood estimates of $q_1$ and $q_2$ from equations (8) and (9) respectively as follows,

$l\left(b;y\right)={\displaystyle \sum _{i=1}^{i=n}\left[y_{1i}\ln \widehat{y}-y_{1i}\ln \left({\widehat{y}}_1+{\widehat{y}}_2+1\right)+y_{2i}\ln {\widehat{y}}_2-y_{2i}\ln \left({\widehat{y}}_{{}_1}+{\widehat{y}}_2+1\right)-\ln \left({\widehat{y}}_{{}_1}+{\widehat{y}}_2+1\right)+\ln \left(y_{1i}+y_{2i}\right)!-\ln y_{1i}!-\ln y_{2i}!\right]}$  (13)

Thus the deviance statistic for testing the identical parameter for each observed pair of data can be expressed as follows,

$D_I=2\left[l\left(b_{\max };y\right)-l\left(b;y\right)\right]=2{\displaystyle \sum _{i=1}^{i=n}\left[y_{1i}\ln \frac{y_{1i}}{{\widehat{y}}_{{}_1}}-y_{1i}\ln \frac{y_{1i}+y_{2i}+1}{{\widehat{y}}_1+{\widehat{y}}_2+1}+y_{2i}\ln \frac{y_{2i}}{{\widehat{y}}_2}-y_{2i}\ln \frac{y_{1i}+y_{2i}+1}{{\widehat{y}}_1+{\widehat{y}}_2+1}-\ln \frac{y_{1i}+y_{2i}+1}{{\widehat{y}}_1+{\widehat{y}}_2+1}\right]}$  (14)

According to Dobson (2001), $D_I$ follows a $X^2$ distribution with $\left(2n-2\right)$ degrees of freedom.

Deviance for a GLM

The deviance statistic for the glm of interest can be written using (1) and the maximum likelihood estimates of $q_{1i}$ and $q_{2i}$ based on the glm from equations (10) and (11) respectively as follows,

$l\left(b;y\right)={\displaystyle \sum _{i=1}^{i=n}\left[y_{1i}\ln \frac{y_{1i}}{\exp \left\{{\displaystyle {\sum }_{j=1}^{j=p}{x}_{2ij}}{\beta }_{2j}\right\}}+y_{2i}\ln \left(1-\frac{y_{{}_{1i}+1}}{exp\left\{{\displaystyle {\sum }_{j=1}^{j=p}{x}_{2ij}}{\beta }_{2j}\right\}}\right)+\ln \frac{1}{exp\left\{{\displaystyle {\sum }_{j=1}^{j=p}{x}_{2ij}}{\beta }_{2j}\right\}}+\ln \left(y_{1i}+y_{2i}\right)!-\ln y_{1i}!-\ln y_{2i}!\right]}$ (15)

Thus the deviance can be expressed as follows,

$D_{II}=2\left[l\left(b_{\max };y\right)-l\left(b;y\right)\right]=2{\displaystyle \sum _{i=1}^{i=n}\left[y_{2i}\ln y_{2i}-\left(y_{1i}+y_{2i}+1\right)\ln \left(y_{1i}+y_{2i}+1\right)+\left\{{\displaystyle \sum _{j=1}^{j=p}{x}_{2ij}{\beta }_{2j}}\right\}\left(y_{1i}+y_{2i}+1\right)-y_{2i}\ln \left(\exp \left\{{\displaystyle \sum _{j=1}^{j=p}{x}_{2ij}{\beta }_{2j}}\right\}-y_{1i}-1\right)\right]}$ (16)

According to Dobson (2001), $D_{II}$ follows ${\chi }^2$ distribution with $\left(2n-p\right)$ degrees of freedom.

Comparison between two GLMs

In order to compare two nested generalized linear models, we consider the fol­lowing hypotheses. The null hypothesis corresponding to a smaller model $\left(M_0\right)$ in terms of number of regression parameters is

$H_0:\beta ={\beta }_0=\left[\begin{array}{l}{\beta }_1\\ {\beta }_2\\ .\\ .\\ .\\ {\beta }_q\end{array}\right]$

The alternative hypothesis corresponding to a bigger model (M1 with q < p < n) within which the smaller model is nested can be written as,

$H_1:\beta ={\beta }_1=\left[\begin{array}{l}{\beta }_1\\ {\beta }_2\\ .\\ .\\ .\\ {\beta }_p\end{array}\right]$

We can test $H_0$ against $H_1$ using the difference of the deviance statistics. Here, $l\left(b_0;y\right)$ is used to denote the likelihood function corresponding to the model $M_0$ and $l\left(b_1;y\right)$  to denote the likelihood function corresponding to the model $M_1$ . Hence the deviance difference can be written as,

$\begin{array}{l}△D=D_0-D_1=2\left[l\left(b_{\max };y\right)-l\left(b_0;y\right)\right]-2\left[l\left(b_{\max };y\right)-l\left(b_1;y\right)\right]=2\left[l\left(b_1;y\right)-l\left(b_0;y\right)\right]\\ =2{\displaystyle \sum _{i=1}^{i=n}\left[\left\{{\displaystyle \sum _{j=1}^{j=p}{x}_{2ij}{\beta }_{2j}}\right\}\left(y_{1i}+y_{2i}+1\right)-y_{2i}\ln \left(\exp \left\{{\displaystyle \sum _{j=1}^{j=p}{x}_{2ij}{\beta }_{2j}}\right\}-y_{1i}-1\right)-\left\{{\displaystyle \sum _{j=1}^{j=q}{x}_{2ij}{\beta }_{2j}}\right\}\left(y_{1i}+y_{2i}+1\right)+y_{2i}\ln \left(\exp \left\{{\displaystyle \sum _{j=1}^{j=q}{x}_{2ij}{\beta }_{2j}}\right\}-y_{1i}-1\right)-\left\{{\displaystyle \sum _{j=1}^{j=q}{x}_{2ij}{\beta }_{2j}}\right\}\left(y_{1i}+y_{2i}+1\right)+y_{2i}\ln \left(\exp \left\{{\displaystyle \sum _{j=1}^{j=q}{x}_{2ij}{\beta }_{2j}}\right\}-y_{1i}-1\right)\right]}\end{array}$

According to Dobson¹² this $\Delta D$ follows ${\chi }^2$ distribution with $p-q$ degrees of freedom.

If the value of $\Delta D$ is consistent with the ${\chi }_{\left(p-q\right)}^2$ distribution we would generally choose the $M_0$ corresponding to $H_0$ because it is simpler. On the other hand, if the value of $\Delta D$ is in the critical region i.e., greater than the upper tail $100\times \alpha \%$ point of the ${\chi }_{\left(p-q\right)}^2$  distribution then would reject $H_0$ in favor of $H_1$ on the grounds that model $M_1$ provides a significantly better description of the data.

To determine the efficiency of our derived deviances we need to have data with known parameters. However, we cannot generate data directly from bivariate geometric distribution using the available computer software packages. Kr­ishna and Pundir suggested an algorithm based on a theorem given by Hogg et al.¹³ to generate random numbers from bivariate geometric distribution. According to this, paired values can be generated from a bivariate geometric distribution using the following steps,

1. Step 1: Generate krandom numbers from univariate geometric distribu­tion with probability of success $\left(\frac{1-q_1-q_2}{1-q_2}\right)$ .
2. Step 2: Suppose that our generated random numbers from the geometric distribution are $x_1,x_2,\mathrm{...},x_k$ .
3. Step 3: Generate k random numbers $y_{ij}$ , k times each from a negative binomial distribution with parameters $x_i+1$ and $\left(1-q_2\right)$ .
4. Step 4: These generated pairs are from the bivariate geometric distribu­tion with parameters $q_1$ and $q_2$ .

In this subsection, we use the following steps to check our derived deviance for the reduced model with identical values of parameters $\left(q_1,q_2\right)$  for each observed pair of data.

1. Step 1: Assume some fixed values of $q_1$ and $q_2$ .
2. Step 2: Generate k random numbers from univariate geometric distribu­tion with probability of success $\left(\frac{1-q_1-q_2}{1-q_2}\right)$ using the assumed values of $q_1$ and $q_2$ from Step 1.
3. Step 3: Suppose that our generated random numbers from the geometric distribution are $x_1,x_2,\mathrm{...},x_k$ .
4. Step 4: Generate k random numbers $y_{ij}$ , k times each from the negative binomial distribution with parameters $x_i+1$ and $\left(1-q_2\right)$ .
5. Step 5: The generated pairs are from the bivariate geometric distribution with parameters $q_1$ and $q_2$ .
6. Step 6: Estimate deviance which is derived in (14).

We take the values of $q_1$ and $q_2$ ranging from 0.10 to 0.90 and satisfying the constraint $q_1+q_2<1$ . We considered several values for the pair $\left(q_1,q_2\right)$ and generate random pairs to observe the efficiency of our derived deviance under different parametric values. For each specified pairs of parameters $\left(q_1,q_2\right)$ , we ran this experiment twice to see whether there is a change in our decision due to randomness. The values of the pair of parameters and the corresponding deviance values are tabulated as follows.

The deviance we derived to test the parameters of the reduced model works well as we see that all, but four of the values of the deviances are smaller than ­ $X_{198}^2\left(0.95\right)$ . However, among these four values of the deviances three are greater than ­ $X_{198}^2\left(0.95\right)$ , but less than $X_{198}^2\left(0.99\right)$ . So, it can be concluded that our derived deviance works well. On the other hand, if most of the values of the deviances had a larger value than our desired ­ ${\chi }^2$ value, then we had to conclude that our derived deviance does not work in testing hypothesis regarding the parameters of the reduced model.

In this paper, we addressed an important problem of inference regarding bivariate geometric distribution and developed testing procedure for the parameters of this distribution with and without covariate information. Our method depends on deriving the deviance statistics using maximum likelihood estimators (mle) of parameters. Our mles of the parameters of the bivariate geometric distribution are obtained using the conditional and the marginal distributions.

We conducted a numerical analysis based on simulated data for the testing the 11 identical parameter assumption for each pair of observed data. Our numerical example did not consider any covariate information. We found that without covariate information our derived deviance worked well in most cases.

None.

Author declares that there are no conflicts of interest.

1. Phatak AG, Sreehari M. Some characterizations of a bivariate geometric distributions. Journal of Indian Statistical Association. 1981;19:141‒146.
2. Nair KRM, Nair NU. On characterizing the bivariate exponential and geometric distributions. Ann Inst Statist Math. 1988;40(2):267‒271.
3. Hawkes AG. On characterizing the bivariate exponential and geometric distributions. Journal of the Royal Statistical Society, Series. 1972;B34:1293.
4. Arnold BC, Castillo E, Erbia J. Conditionally specified distributions. Springer, New York, USA. 1992.
5. Sreehar M, Vasudev R. Characterizations of multivariate geometric distributions in terms of conditional distributions. Metrica. 2012;75(2):271‒286.
6. Basu AP, Dhar SK. Bivariate geometric distribution. Journal of Applied Statistical Sciences. 1995;2(1):12.
7. Marshal AW, Olkin IA. A multivariate exponential distribution. Journal of the American Statistical Association. 1967;62(317):30‒44.
8. Sun K, Basu AP. A characterization of a bivariate geometric distribution. Statistics and Probability Letters. 1995;23(4):307‒311.
9. Sreehari G. Characterization via conditional distributions. Journal of Indian Statistical Association. 2005;43:77‒93.
10. Omey Edward, Minkova DL. Bivariate geometric distributions. Hub Research Paper 2013=02 Economics & Science. 2013.
11. Krishna Hare, Pundir PS. A bivariate geometric distribution with applications to reliability. Communications in Statistics‒ Theory and Methods. 2009;38(7):1079‒1093.
12. Dobson JA. An introduction to generalized linear models. (2nd edn), Chapman & Hall, CRC press, UK. 2001.
13. Hogg RV, McKean JW, Craig AT. Introduction to mathematical statistics. (6th edn), New Delhi: Pearson Education, India. 2005.