VGAM (version 1.1-2)

bigamma.mckay: Bivariate Gamma: McKay's Distribution

Description

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

Usage

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

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.

imethod, 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

# NOT RUN {
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, bigamma.mckay, data = mdata, trace = TRUE)
coef(fit, matrix = TRUE)
Coef(fit)
vcov(fit)

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