Estimate the parameter of a gamma hyperbola bivariate distribution by maximum likelihood estimation.

gammahyperbola(ltheta = "loglink", itheta = NULL, expected = FALSE)

Arguments

ltheta

Link function applied to the (positive) parameter \(\theta\). See Links for more choices.

itheta

Initial value for the parameter. The default is to estimate it internally.

expected

Logical. FALSE means the Newton-Raphson (using the observed information matrix) algorithm, otherwise the expected information matrix is used (Fisher scoring algorithm).

Details

The joint probability density function is given by $$f(y_1,y_2) = \exp( -e^{-\theta} y_1 / \theta - \theta y_2 )$$ for \(\theta > 0\), \(y_1 > 0\), \(y_2 > 1\). The random variables \(Y_1\) and \(Y_2\) are independent. The marginal distribution of \(Y_1\) is an exponential distribution with rate parameter \(\exp(-\theta)/\theta\). The marginal distribution of \(Y_2\) is an exponential distribution that has been shifted to the right by 1 and with rate parameter \(\theta\). The fitted values are stored in a two-column matrix with the marginal means, which are \(\theta \exp(\theta)\) and \(1 + 1/\theta\).

The default algorithm is Newton-Raphson because Fisher scoring tends to be much slower for this distribution.

Value

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

References

Reid, N. (2003). Asymptotics and the theory of inference. Annals of Statistics, 31, 1695–1731.

Author

T. W. Yee

Note

The response must be a two-column matrix.

See also

Examples

gdata <- data.frame(x2 = runif(nn <- 1000))
gdata <- transform(gdata, theta = exp(-2 + x2))
gdata <- transform(gdata, y1 = rexp(nn, rate = exp(-theta)/theta),
                          y2 = rexp(nn, rate = theta) + 1)
fit <- vglm(cbind(y1, y2) ~ x2, gammahyperbola(expected = TRUE), data = gdata)
coef(fit, matrix = TRUE)
#>             loglink(theta)
#> (Intercept)     -1.9947417
#> x2               0.9763341
Coef(fit)
#> (Intercept)          x2 
#>  -1.9947417   0.9763341 
head(fitted(fit))
#>          y1       y2
#> 1 0.4051930 4.331844
#> 2 0.3051449 5.166170
#> 3 0.2954764 5.276041
#> 4 0.2081538 6.721628
#> 5 0.3512969 4.723570
#> 6 0.1791261 7.509636
summary(fit)
#> 
#> Call:
#> vglm(formula = cbind(y1, y2) ~ x2, family = gammahyperbola(expected = TRUE), 
#>     data = gdata)
#> 
#> Coefficients: 
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept) -1.99474    0.04302  -46.37   <2e-16 ***
#> x2           0.97633    0.07060   13.83   <2e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Name of linear predictor: loglink(theta) 
#> 
#> Log-likelihood: -2224.826 on 998 degrees of freedom
#> 
#> Number of Fisher scoring iterations: 17 
#> 
#> No Hauck-Donner effect found in any of the estimates
#>