VGAM (version 0.8-1)

bivgamma.mckay: Bivariate Gamma: McKay's Distribution

Description

Estimate the three parameters of McKay's bivariate gamma distribution by maximum likelihood estimation.

Usage

bivgamma.mckay(lscale = "loge", lshape1 = "loge", lshape2 = "loge",
               iscale = NULL, ishape1 = NULL, ishape2 = NULL,
               method.init=1, zero = 1)

Arguments

lscale, lshape1, lshape2

Link functions applied to the (positive) parameters $a$, $p$ and $q$ respectively. See Links for more choices.

iscale, ishape1, ishape2

Optional initial values for $a$, $p$ and $q$ respectively. The default is to compute them internally.

method.init, zero

See CommonVGAMffArguments.

Value

• An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm and vgam.

Details

One of the earliest forms of the bivariate gamma distribution has a joint probability density function given by $$f(y_1,y_2;a,p,q) = (1/a)^{p+q} y_1^{p-1} (y_2-y_1)^{q-1} \exp(-y_2 / a) / [\Gamma(p) \Gamma(q)]$$ for $a > 0$, $p > 0$, $q > 0$ and $0 < y_1 < y_2$ (Mckay, 1934). Here, $\Gamma$ is the gamma function, as in gamma. By default, the linear/additive predictors are $\eta_1=\log(a)$, $\eta_2=\log(p)$, $\eta_3=\log(q)$.

The marginal distributions are gamma, with shape parameters $p$ and $p+q$ respectively, but they have a common scale parameter $a$. Pearson's product-moment correlation coefficient of $y_1$ and $y_2$ is $\sqrt{p/(p+q)}$. This distribution is also known as the bivariate Pearson type III distribution. Also, $Y_2 - y_1$, conditional on $Y_1=y_1$, has a gamma distribution with shape parameter $q$.

References

McKay, A. T. (1934) Sampling from batches. Journal of the Royal Statistical Society---Supplement, 1, 207--216.

Kotz, S. and Balakrishnan, N. and Johnson, N. L. (2000) Continuous Multivariate Distributions Volume 1: Models and Applications, 2nd edition, New York: Wiley.

Balakrishnan, N. and Lai, C.-D. (2009) Continuous Bivariate Distributions, 2nd edition. New York: Springer.

See Also

gamma2.

Examples

shape1 = exp(1); shape2 = exp(2); scalepar = exp(3)
mdata = data.frame(y1 = rgamma(nn <- 1000, shape=shape1, scale=scalepar))
mdata = transform(mdata, zedd = rgamma(nn, shape=shape2, scale=scalepar))
mdata = transform(mdata, y2 = y1 + zedd) # Z is defined as Y2-y1|Y1=y1
fit = vglm(cbind(y1, y2) ~ 1, bivgamma.mckay, mdata, trace = TRUE)
coef(fit, matrix = TRUE)
Coef(fit)
vcov(fit)

colMeans(fit@y)    # Check moments
head(fitted(fit), 1)