Freund's (1961) Bivariate Extension of the Exponential Distribution
Estimate the four parameters of the Freund (1961) bivariate extension of the exponential distribution by maximum likelihood estimation.
freund61(la = "loglink", lap = "loglink", lb = "loglink", lbp = "loglink",
ia = NULL, iap = NULL, ib = NULL, ibp = NULL,
independent = FALSE, zero = NULL)la, lap, lb, lbp |
Link functions applied to the (positive)
parameters alpha, alpha',
beta and beta', respectively
(the “ |
ia, iap, ib, ibp |
Initial value for the four parameters respectively. The default is to estimate them all internally. |
independent |
Logical. If |
zero |
A vector specifying which
linear/additive predictors are modelled as intercepts only.
The values can be from the set {1,2,3,4}.
The default is none of them.
See |
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.
An object of class "vglmff" (see vglmff-class).
The object is used by modelling functions such as vglm
and vgam.
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].
T. W. Yee
Freund, J. E. (1961). A bivariate extension of the exponential distribution. Journal of the American Statistical Association, 56, 971–977.
fdata <- data.frame(y1 = rexp(nn <- 1000, rate = exp(1)))
fdata <- transform(fdata, y2 = rexp(nn, rate = exp(2)))
fit1 <- vglm(cbind(y1, y2) ~ 1, fam = freund61, data = fdata, trace = TRUE)
coef(fit1, matrix = TRUE)
Coef(fit1)
vcov(fit1)
head(fitted(fit1))
summary(fit1)
# y1 and y2 are independent, so fit an independence model
fit2 <- vglm(cbind(y1, y2) ~ 1, freund61(indep = TRUE),
data = fdata, trace = TRUE)
coef(fit2, matrix = TRUE)
constraints(fit2)
pchisq(2 * (logLik(fit1) - logLik(fit2)), # p-value
df = df.residual(fit2) - df.residual(fit1), lower.tail = FALSE)
lrtest(fit1, fit2) # Better alternativePlease choose more modern alternatives, such as Google Chrome or Mozilla Firefox.