VGAM (version 1.0-5)

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 = "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)/[\mathrm{\Gamma }(p)\mathrm{\Gamma }(q)]$$ for $a>0$, $p>0$, $q>0$ and $0<y_1<y_2$ (Mckay, 1934). Here, $\mathrm{\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)
# }