Gumbel's Type I Bivariate Distribution Family Function

Estimate the association parameter of Gumbel's Type I bivariate distribution by maximum likelihood estimation.

bigumbelIexp(lapar = "identitylink", iapar = NULL, imethod = 1)

Arguments

lapar

Link function applied to the association parameter $\alpha$. See Links for more choices.

iapar

Numeric. Optional initial value for $\alpha$. By default, an initial value is chosen internally. If a convergence failure occurs try assigning a different value. Assigning a value will override the argument imethod.

imethod

An integer with value 1 or 2 which specifies the initialization method. If failure to converge occurs try the other value, or else specify a value for ia.

Details

The cumulative distribution function is $$P(Y_1 \leq y_1, Y_2 \leq y_2) = e^{-y_1-y_2+\alpha y_1 y_2} + 1 - e^{-y_1} - e^{-y_2} $$ for real $\alpha$. The support of the function is for $y_1>0$ and $y_2>0$. The marginal distributions are an exponential distribution with unit mean.

A variant of Newton-Raphson is used, which only seems to work for an intercept model. It is a very good idea to set trace=TRUE.

Value

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

References

Gumbel, E. J. (1960). Bivariate Exponential Distributions. Journal of the American Statistical Association, 55, 698–707.

Author

T. W. Yee

Note

The response must be a two-column matrix. Currently, the fitted value is a matrix with two columns and values equal to 1. This is because each marginal distribution corresponds to a exponential distribution with unit mean.

This VGAM family function should be used with caution.

Examples

nn <- 1000
gdata <- data.frame(y1 = rexp(nn), y2 = rexp(nn))
if (FALSE)  with(gdata, plot(cbind(y1, y2)))  # \dontrun{}
fit <- vglm(cbind(y1, y2) ~ 1, bigumbelIexp, gdata, trace = TRUE)
#> Iteration 1: loglikelihood = -2043.3139
#> Iteration 2: loglikelihood = -2042.999
#> Iteration 3: loglikelihood = -2042.9984
#> Iteration 4: loglikelihood = -2042.9984
coef(fit, matrix = TRUE)
#>                    apar
#> (Intercept) -0.04082615
Coef(fit)
#>        apar 
#> -0.04082615 
head(fitted(fit))
#>   y1 y2
#> 1  1  1
#> 2  1  1
#> 3  1  1
#> 4  1  1
#> 5  1  1
#> 6  1  1