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 = "loge", lshape1 = "loge", lshape2 = "loge",
              iscale = NULL, ishape1 = NULL, ishape2 = NULL,
              imethod = 1, zero = 2:3)

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) / [Γ (p) Γ (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.

Examples

Run this code

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)

Run the code above in your browser using DataLab

Get 50% off unlimited learning