bigamma.mckay: Bivariate Gamma: McKay's Distribution

View source: R/family.bivariate.R

bigamma.mckayR Documentation

Bivariate Gamma: McKay's Distribution


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


bigamma.mckay(lscale = "loglink", lshape1 = "loglink",
              lshape2 = "loglink", iscale = NULL, ishape1 = NULL,
              ishape2 = NULL, imethod = 1, zero = "shape")


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.

imethod, zero

See CommonVGAMffArguments.


One of the earliest forms of the bivariate gamma distribution has a joint probability density function given by

f(y1,y2;a,p,q) = (1/a)^(p+q) y1^(p-1) (y2-y1)^(q-1) exp(-y2/a) / [gamma(p) gamma(q)]

for a > 0, p > 0, q > 0 and 0<y1<y2 (Mckay, 1934). Here, gamma is the gamma function, as in gamma. By default, the linear/additive predictors are eta1=log(a), eta2=log(p), eta3=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 y1 and y2 is sqrt(p/(p+q)). This distribution is also known as the bivariate Pearson type III distribution. Also, Y2 - y1, conditional on Y1=y1, has a gamma distribution with shape parameter q.


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


The response must be a two column matrix where the first column is y1 and the second y2. It is necessary that each element of the vectors y1 and y2-y1 be positive. Currently, the fitted value is a matrix with two columns; the first column has values ap for the marginal mean of y1, while the second column has values a(p+q) for the marginal mean of y2 (all evaluated at the final iteration).


T. W. Yee


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



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

colMeans(depvar(fit))  # Check moments
head(fitted(fit), 1)

VGAM documentation built on July 6, 2022, 5:05 p.m.