freund61 {VGAM} | R Documentation

## Freund's (1961) Bivariate Extension of the Exponential Distribution

### Description

Estimate the four parameters of the Freund (1961) bivariate extension of the exponential distribution by maximum likelihood estimation.

### Usage

    freund61(la="loge", lap="loge", lb="loge", lbp="loge",
             ia=NULL, iap=NULL, ib=NULL, ibp=NULL,
             independent=FALSE, zero=NULL)

### Arguments

`la,lap,lb,lbp` | Link function applied to the (positive) parameters alpha, alpha', beta and beta', respectively (the ```p`'' stands for ``prime''). See `Links` for more choices.
`ia,iap,ib,ibp` | Initial value for the four parameters respectively. The default is to estimate them all internally.
`independent` | Logical. If `TRUE` then the parameters are constrained to satisfy alpha=alpha' and beta=beta', which implies that y1 and y2 are independent and each have an ordinary exponential distribution.
`zero` | An integer-valued vector specifying which linear/additive predictors are modelled as intercepts only. The values must be from the set {1,2,3,4}. The default is none of them.

### Details

This model represents one type of bivariate extension of the exponential distribution that is applicable to certain problems, in particular, to two-component systems which can function if one of the components has failed. For example, engine failures in two-engine planes, paired organs such as peoples' eyes, ears and kidneys. Suppose y1 and y2 are random variables representing the lifetimes of two components A and B in a two component system. The dependence between y1 and y2 is essentially such that the failure of the B component changes the parameter of the exponential life distribution of the A component from alpha to alpha', while the failure of the A component changes the parameter of the exponential life distribution of the B component from beta to beta'.

The joint probability density function is given by

f(y1,y2) = alpha * beta' * exp(-beta' * y2 - (alpha+beta-beta') * y1)

for 0 < y1 < y2, and

f(y1,y2) = beta * alpha' * exp(-alpha' * y1 - (alpha+beta-alpha') * y2)

for 0 < y2 < y1. Here, all four parameters are positive, as well as the responses y1 and y2. Under this model, the probability that component A is the first to fail is alpha/(alpha+beta). The time to the first failure is distributed as an exponential distribution with rate alpha+beta. Furthermore, the distribution of the time from first failure to failure of the other component is a mixture of Exponential(alpha') and Exponential(beta') with proportions beta/(alpha+beta) and alpha/(alpha+beta) respectively.

The marginal distributions are, in general, not exponential. By default, the linear/additive predictors are eta1=log(alpha), eta2=log(alpha'), eta3=log(beta), eta4=log(beta').

A special case is when alpha=alpha' and beta'=beta', which means that y1 and y2 are independent, and both have an ordinary exponential distribution with means 1/alpha and 1/beta respectively.

Fisher scoring is used, and the initial values correspond to the MLEs of an intercept model. Consequently, convergence may take only one iteration.

### Value

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

### Note

To estimate all four parameters, it is necessary to have some data where y1<y2 and y2<y1.

The response must be a two-column matrix, with columns y1 and y2. Currently, the fitted value is a matrix with two columns; the first column has values (alpha'+beta)/(alpha' * (alpha+beta)) for the mean of y1, while the second column has values (beta'+alpha)/(beta' * (alpha+beta)) for the mean of y2. The variance of y1 is

[(alpha')^2 + 2 * alpha * beta + beta^2]/ [(alpha')^2 * (alpha + beta)^2],

the variance of y2 is

[(beta')^2 + 2 * alpha * beta + alpha^2]/ [(beta')^2 * (alpha + beta)^2],

the covariance of y1 and y2 is

[alpha' * beta' - alpha * beta]/ [alpha' * beta' * (alpha + beta)^2].

### Author(s)

T. W. Yee

### References

Freund, J. E. (1961) A bivariate extension of the exponential distribution. Journal of the American Statistical Association, 56, 971–977.

### See Also

`exponential`.

### Examples

    y1 = rexp(n <- 200, rate=4)
    y2 = rexp(n, rate=8)
    ymat = cbind(y1,y2)
    fit =  vglm(ymat ~ 1, fam=freund61, trace=TRUE)
    coef(fit, matrix=TRUE)
    Coef(fit)
    vcov(fit)
    head(fitted(fit))
    summary(fit)

    # y1 and y2 are independent, so fit an independence model
    fit2 = vglm(ymat ~ 1, fam=freund61(indep=TRUE), trace=TRUE)
    coef(fit2, matrix=TRUE)
    constraints(fit2)
    1 - pchisq(2*(logLik(fit)-logLik(fit2)), df=2) # p-value

[Package VGAM version 0.7-9 Index]