- Email: [email protected]
Inference procedures for bivariate exponential model of Gumbel
Statistics & Probability North-Holland
Letters 12 (1991) 37-50
July 1991
Inference procedures for bivariate exponential model of Gumbel Jye-Chyi Department
Lu of Statistics,
North Carolina State University, Raleigh, NC 27695-8203,
USA
Gouri K. Bhattacharyya Department Received
of Statistics,
University of Wisconsin,
Madison,
WI 53706, USA
May 1990
Abstract: In this article, statistical properties of an absolutely continuous bivariate exponential model due to Gumbel (1960) are examined and the Fisher information matrix is derived. Simple estimators are motivated from a structural representation of the model and are compared both in terms of asymptotic efficiency and finite sample simulation. Besides their good performance, some of the simple estimators lead to exact or nearly exact confidence interval and hypotheses test procedures. The proposed methods are illustrated by reanalyzing a data set to which asymptotic methods were employed by Gross and Lam (1981) under a different bivariate exponential model. Keywords:
Dependent
lifetimes,
Fisher information,
simple estimators,
asymptotic
efficiency,
confidence
intervals,
simulation.
1. Introduction Many bivariate extensions of the exponential distribution have been proposed in the literature in order to model dependence of the lifetimes of paired organs or components in a system. The models derived from physical behavior of the failure process are generally more dependable in reliability studies than those that are just mathematical constructs, and are especially so when a large enough sample is not available for a test of fit. Tractability of statistical inferences is yet another criterion for a model to be useful in analyzing life test data. In both these regards, the bivariate exponential model of Marshall and Olkin (1967) has considerable advantages over the others. However, because of a singular part, its use is not appropriate in the situations where simultaneous failures of the components are quite unlikely to occur. Freund’s (1961) bivariate extension of the exponential is also well-grounded on a physical basis and is absolutely continuous but does not have exponential marginals. Two absolutely continuous modifications of the Marshall-Olkin model have been proposed: ACBVE by Block and Basu (1974) and ACBVE, by Sarkar (1987). The former is a sub-class of Freund’s model while the latter does not seem to have as strong a physical basis and is not convenient for statistical inferences. In this article, we investigate the properties and inference procedures for a rather neglected bivariate exponential model whose survival function and probability density function (p.d.f.) are respectively given by F(x,,
x2) = exp
I
-
{
(x,/0,)“’
+ (xa/Ws}
(l.la)
“1)
(x,/e,)'/"]"-2 f(x,, ~,)=(s,e,)~'~"(x,x~)~~~-'[(~~/e~)~~*+ + (x2/8Jis}'+
i/s-
i] exp[-
{ (x1/e,)1/6+(x2/e2)1/'}6]. (l.lb)
0167.7152/91/$03.50
0 1991 - Elsevier Science
Publishers
B.V. (North-Holland)
37
Volume
12, Number
1
STATISTICS
& PROBABILITY
LETTERS
July 1991
0 < x,, x2 < cc, 0 < 0,, 19~< cc, 0 i S < 1. It is a transformation of Gumbel’s type B bivariate extreme-value distribution (cf. Johnson and Kotz, 1972, p. 251) proposed by Gumbel (1960) and we will henceforth refer to it as GBVE(8,, 8,, S). Here, 8, and 8, are the scale parameters, S is the dependence parameter, and 6 = 1 corresponds to the case of independence. In spite of its long history, the GBVE has received little attention in modeling dependence of paired lifetimes seemingly for two reasons: (i) it was proposed solely as a mathematical construct and (ii) the complicated form of its p.d.f. poses considerable difficulties with statistical inference. That it is indeed a meaningful physical model surfaces from Hougaard (1986b) who provides a general construction of dependent life models from the consideration of random stress environment affecting both components (see Lu and Bhattacharyya, 1988, for other derivations and extensions of this model). Specifically, if conditionally given the stress S = s, the component lives X, and X, are independent Weibull with survival i = 1, 2, respectively, and if S has a positive stable functions exp[-H,(t)s] and H,(t) = (r/0,)“S, distribution (cf. Hougaard, 1986a), then the unconditional distribution of (Xi, X,) is GBVE(e,, 0,, 6). Thus, while the source of dependence in the Marshall-Olkin model is a common shock destroying both components, and that in the Freund model is a failure event internal to the system, it is the random mixing effect of an external stress that serves to explain the dependence in GBVE. Such stochastic failure rate models are long recognized in reliability modeling, see for instance, Mann, Schafer and Singpurwalla (1974, p. 145). With the recognition that GBVE has a sound physical basis in modeling bivariate life times, our major goal here is to develop procedures for parameter estimation, confidence statements and hypotheses tests. Some statistical properties of the model, including the derivation of the Fisher information matrix, are presented in Section 2. Despite the complicated form of the p.d.f. due to which the MLE’s could not be derived in closed forms, it has been possible to obtain closed form expressions for the elements of the Fisher information matrix. These are especially handy in evaluating the asymptotic efficiency (AE) of some simple estimators constructed in Section 3. A representation, due to Lee (1979) of (X,, X,) in terms of certain independence random variables proves instrumental for these calculations and also for the development of some simple inference procedures. Based on a random sample (Xi,, X,,), i = 1,. . . , n, from GBVE(B,, 8,, S), some simple estimators of the model parameters are constructed in Section 3. These are compared both in terms of AE and Monte Carlo simulation. For the scale parameters S,, the sample means x1 = n-‘C:=,X,;, j = 1, 2, are found to retain a remarkably high AE under the bivariate model and they also do very well in small samples. A simple estimator of the dependence parameter 6 is motivated from the exponential distribution of min( X,, X,), and its AE is evaluated. Other simple estimators are constructed from the order statistics of log( X,/X,,), i = 1,. _. , n, or by referring to their likelihood. These are found to be substantially better than the first and quite competitive with the MLE. This latter construction draws from the property that log( X,/X,) generates a location-scale logistic family, which itself is a consequence of the independence representation of GBVE. Problems of confidence intervals and hypotheses tests are taken up in Section 4. Pivotal quantities are identified for the basic model parameters as well as for (3,/e,, the ratio of the component mean lives. Exact, nearly exact and asymptotic procedures are discussed. Applications of the proposed methods are illustrated by means of the data set of a paired comparison experiment which was previously analyzed by by Gross and Lam (1981) under the ACBVE model.
2. Properties
of GBVE
and the Fisher
information
Several distributional properties of GBVE can be conveniently deduced from a structural representation due to Lee (1979) in terms of independent random variables. A slight variant of this result is given in the as”, exp(a) for an exponential following lemma. The notation - will be used to mean “is distributed 38
distribution with mean parameter h.
a, and logistic( a, b) for a logistic distribution
with location
parameter
a and scale
Lemma 2.1 (Lee, 1979). Consider independent random variables U, V,,, V,, and M6 such that U where uniform(O, I), V,, - exp(l), i = 1, 2, and M, = 0 or 1 with probability 1 - 6 and 6, respectively, 6 E (0, l] is a constant. Letting 8, and f& be positive constants, define V= V,, + MJl:2,
Then (X,,
x, = 6,U6V,
x, = f?,(l - U)%.
X,) - GBVE(8,,
6,, 6).
(2.1)
0
The moments of GBVE are readily obtained from (2.1) and the facts that b + 1) and E(V”“) = (1 + a6)F(a + l), a, b > 0. In particular, E( X;lX;‘)
= 0;l@F(
k,, k, > 0, and the correlation p(X,,
k,S + l)F( coefficient
k,S + l)F(
k, + k, + l),‘F(
E[ZY’(l - U)h] = B(a + 1,
k,S + k,6 + l),
(2.2)
of X, and X, is given by
x,)=2r2(6+1)/r(26+1)-1,
(2.3)
which ranges from 0 to 1 and is monotone decreasing in 6. Let T = min( X,, X,), W = log( X,/X,) and D = I( X, < X,) where Z( .) denotes the indicator function. Then, the following properties hold: (a) W - logistic(7, 6) with n = log(0,/8,), (b) Pr( D = 1) = and (d) T and D are independent. Property (a) e;‘/s/(e;l/s + ep”), (c) T - exp[(e, p’/s + e;‘/‘)-“1, readily follows from (2.1) and the fact that log[U/(l - U)] - logistic(0, 1). Also, since X, < X, is equivalent to log[U/(l - U)] <) - q/6, (b) follows from the logistic c.d.f. Verification of (c) and (d) involves some routine algebra which we omit for brevity. Of course, (c) is obvious from the survival function (l.la). Incidentally, the Marshall-Olkin model also has similar properties (c) and (d). It is commonly known that the logarithm of the ratio of two i.i.d. exp(1) random variables has the standard logistic distribution. As a dependent exponential model, GBVE has the interesting feature that the logistic family is preserved under the aforementioned transformation, and the dependence parameter 6 becomes a scale parameter. Consider the function
Q(X,,
s
X2>= [w4Y’~s + (x,/B:)“‘]
(2.4)
which appears in the exponent of the survival function (1.1). From the representation (2.1), we note that Q( X,, X,) = V which is a binary mixture of independent exp(1) and gamma(1, 2). In particular, for 6 = 1 (independence), it has the &i-square distribution with d.f. = 4. Expectations of some functions of Q will be needed in connection with the calculation of the Fisher information matrix for GBVE. Its identification with I’, which has a simple distribution, facilitates the work. A few key results are listed in the following lemma. Lemma
2.2. Let a(S) = l/6
- 1. For the random variable
(a>
E[{V+a(G)}-‘1
=E[V{V+a(G)}-‘]=a,
@I
E[{V+a(S)}-‘1
=aq0,,
cc>
E[V{v+a(G)}-*I
=6-
V defined in (2.1) we have
(1 -6)q0,,
where qol = exp[a(6)]E,[a(6)] and E,(t) = jp”x-’ exp( - x) dx is the classical exponential for which extensive tables are available (cf. Abramowitz and Stegun, 1970).
integral function
39
Proof. Let V= V’, + V’, where V” and V’, are independent exp(l), and for integers j > 0, k > 1, define the expectation of the random variable Rlk = V’[ V + q(j, k) = E[V/,{ V’, + ~(6>}~‘1. C onsidering a(6)lek where v= V’, + M,V’,, and first taking expectation over M,, we obtain
E(R,,)=(l-S)q(j,
k)+SE[V1{V+u(8)jPk]
= (1 - S)q(j,
k) + Sq( j + 1, k) = 6q( j, k - 1).
(2.5)
The second equality follows from the fact that V.- gamma(1, 2) so E[g( V)] = E[V”g( V,‘)] holds for any integrable function g. The last equality obtains from the identity (1 - S) + SV,, = S[T/,’ + a(6)]. Taking (j, k) = (1, 1) and (0, 2), the first two of the stated results follow from (2.5) once we note that q(1, 0) = 1 and qol = q(0, 1) = exp[4U4[4W For the last result we take (j, k) = (1, 2) and use the relation q(1, 1) = 1 - a(6)q(O, 1). q 2.1. Fisher information
for GBVE(B,,
0,, 6)
Because 8, and 8, are scale parameters in the model, it will be convenient to work with scaled variables Z, = X,/0,, j = 1, 2, whose joint p.d.f. will be denoted by g( z,, z2, 8). Define X = log g( z,, z2, 6) and Q = (Z:‘fi + z:‘“)s and referring to (l.l), we have X=a(G)(log The log-likelihood I=
Z, +log
Z,) + (l-2/6)
log Q+log[Q+a@)]
-Q.
(2.6)
I of ( X,, X,) is then given by
-loge,-log8,+A.
(2.7)
For the convenience of notation, we identify 6 as the argument 19’ in I and as Z, in X, denote 8 = (8’, 19,, I!?,)‘, Z= (Z’, Z,, Z,)’ and view X as a function of Z. Also, we denote the partial derivatives i,= algae,, I,k = a21/ae,ae,, x,= ax/az,, x,,= aQ/az,az,, j, k = 1, 2, 3. Theorem 2.1. The Fisher information 1” = E( -I’,)
matrix
= (38,?-‘[26
i,, = (3t3’0z)P’[6
Z = (i,,)
= E( - I,,)
of GBVE( 0,, f?,, S) is given by
+ K2 + (1 - S)‘S-‘qO,],
- K2 + (26))‘(1
i,, = (0’/82)2i’,,
- S)‘qo’], (2.X)
i ,3 = (20,)P’[SP’ ‘33
-S
- (26))‘(1
- 6){2-‘(1
- 6) + K’}qo’],
i,, = (t3’/B2)i13,
=6-‘(1/S-1)+2q,(S-6-~)+qo,[q’S-‘(l-6)2-S~~(1+S+S~‘)],
where qo, is as defined in Lemma 2.2 and q, = $ + f[{ q(2) - S(4)}2 - ‘k,(4)] = 0.2835533, where ‘P(t) = d log r(t)/dt and q’(t) = d \k(t)/dt, are the digamma and trigumma function, respectively. Proof.
Recalling
that aZ’/M’
I,, = e;2[1+ Differentiating
= -0;‘Z’
and aZ,/M,
= 0, j = 2, 3, relation
(2.7) yields
2Z,X’ + Z’Jh”].
(2.9)
(2.6) we obtain
h,=a@)Z;‘+(Z,/Q)“‘“‘[(1-2/S)QP’+{Q+a(8)}-l-l],
C-?/Q) 20(B’[(1-2/S)QP2+
A,’ = -@)Z,‘+ a(6){
40
Q’-‘/9;/S-2
_
Q’-2/9;/8-
{Q+a(s)}-‘1 2}[(1-2/6)Q-‘+
{Q+a(8)}-l-11.
(2.10)
Substituting (2.10) in (2.9) and referring and Z, = U&V we obtain @I,, = 1-t 2a(6) -a(6)
+ 24(1-
Use of Lemma 2.2 then leads provided in Appendix A. 0
3. Estimation
2/6)
cP[(l -
+a(s)u(l-
to the representation
2/S)
+ v{ V+ a(S)}
+ vZ{ v+
o’)[(l-2/6)+ to the stated
(2.1) and its consequences
a(S)}
-’ - v] -‘I
Y{V+a(G)}_‘expression
such as Q = V
v].
for i,,. The derivations
of i,,,
ii,
and
i,,
are
of parameters
In this section, we consider estimation of the parameters of GBVE( 8,, t?,, 6) from a random sample (X,,, X2,), i = 1,. . . , n. The density model (1.1) is not an exponential family, and thus no drastic sufficient reduction is possible. The maximum likelihood estimators have to be calculated through some numerical procedures. Our goal here is to construct simple estimators which are intuitively appealing, perform well in small or moderate samples and do not incur much loss of asymptotic efficiency compared to the MLE’s. Since the marginal distributions are exponential with means 6, and &, it is natural to estimate 0, by $ = 3 = nP’xy=, X,,, j = 1, 2. These are simple estimators with known exact properties such as: unbiased, Var(0,) = 8:,/n and 2n8;/8, - ~2,. The asymptotic efficiency, AE( 8,) = i”/Bf which depends only on the parameter S, is computed by inverting the Fisher information matrix given in Section 2. The numerical results show that AE(&,) has a U-shaped curve and attains its minimum 98.215% at S = 0.725. Besides their uniformly high AE and good performance in small samples (see Table l), the fact that these estimators are robust with respect to the dependence structure of the bivariate model makes them very useful for practical application. As for the dependence parameter (S), a moment-type estimator based on the sample correlation Y is a natural choice. In view of (2.3) this estimator s”,, is the solution of the equation l+r=2P(8+1)/r(26”+1),
(3.1)
if r > 0 and is defined to be 1 if r < 0. To obtain the asymptotic distribution of this moment recall that the asymptotic distribution of nil2 (p - p), where b is the sample product-moment coefficient is normal with mean 0 and variance crP2as P2(~22~~I’+4~‘[~40~2C,1+~22(~~2~2~)~’+~ir4~011]
estimator, we correlation
-[O~~~((Y:U(YI,)-‘+~~~((YO~(YII)-’]),
where (Y,, = E{[ X- E( X)]‘[Y - E(Y)]‘}, and (X, Y) has a general bivariate distribution. In the GBVE distribution, the moments aA, can be obtained using the expression (2.2). By an application of delta method it follows that ~“~(8, - S) is asymptotically N(0, cr$) where CT; = u;166’r2(26
+ 1)r-4(6
+ l)[ \k(S + 1) - 9(26
+ l)] -2.
The lack of a closed form expression s”, and the lowness (cf. Table 2) of its asymptotic efficiency are obvious disadvantages. Two approaches for the construction of simple estimators of S turn out to be fruitful. The first draws from the fact that the scaled pair (X,/8,, X,/e,) - GBVE(1, 1, 6) so 7;(o) = min( X,/8,, X2,/e,), i=l ,...> n, are i.i.d. exp(2-S). If 0 were known, 7;,(o) = n-‘Z3f=,7;(8) could be used to estimate 2-“. Substitution of 2 = (x1, x2)’ for 8 motivates the estimator 8,= -(log2)P’T,(x)=
-(log2))W
,g, min { X1,/X, ? x2,/%
}.
(3.2)
41
Table
1
Results of a Simulation n
6
10
1 .oo
from GBVE(l, 61
1, 6) e,
82
B,
0.9904
0.9870
1.0148
1.0174
0.8427
0.8434
0.8850
0.1077
0.1038
0.1106
0.0313
0.0277
0.0225
0.0216
Mse
0.1032
0.1079
0.1040
0.1109
0.0561
0.0522
0.0357
0.0357
39.4
1 .oo
1.0038
1.0038
1.0026
1 .oooo
0.8360
0.8322
0.8737
0.8730
0.1096
0.0951
0.1023
0.0327
0.0261
0.0240
0.0212
Mse
0.1071
0.1096
0.0951
0.1023
0.0298 41 .o
0.0272 37.2
Mean
0.9986
0.9909
0.9968
0.9887
0.6969
0.6762
0.7331
0.7141
0.0962
0.0949
0.0909
0.0916
0.;0390
0.0268
0.0282
0.0262
Mse
0.0962
0.0950
0.0909
0.0918
0.0418
0.0322
0.0285
0.0275
13.2
6.6
12.4
7.2
Mean
1.0002
0.9871
1.0050
1.0037
0.4662
0.4546
0.4976
0.4820
Var
0.1092
0.1058
0.1099
0.1065
0.0301
0.0166
0.0199
0.0174
Mse
0.1092
0.1059
0.1099
0.1065
0.0313
0.0187
0.0199
0.0177
1.2
0.0
0.2
0.0
Mean
1.0044
1.0049
1.0030
1.0025
0.2271
0.2229
0.2459
0.2393
Var
0.0935
0.0929
0.1026
0.0997
0.0073
0.0043
0.0054
0.0049
Mse
0.0935
0.0929
0.1026
0.0997
0.0079
0.0050
0.0054
0.0051
0.0
0.0
0.0
0.0
0.9238
0.9357
Mean
1.0063
1.0026
1.0027
0.9975
0.9015
0.9067
Var
0.0496
0.0515
0.0467
0.0471
0.0161
0.0127
0.0109
0.0086
Mse
0.0497
0.0516
0.0467
0.0471
0.0258
0.0214
0.0167
0.0127
47.6
41.0
50.4
48.6 0.9045
Mean
0.9864
0.9828
1.0122
1.0126
0.8760
0.8811
0.9011
Var
0.0447
0.0455
0.0484
0.0493
0.0182
0.0142
0.0125
0.0105
Mse
0.0449
0.0458
0.0485
0.0495
0.0236
0.0189
0.0149
0.0126
34.0
29.6
37.6
33.0
Mean
1.0047
1.0002
1.0137
1.0073
0.7506
0.7253
0.7507
0.7462
Var
0.0462
0.0458
0.0499
0.0512
0.0255
0.0172
0.0176
0.0176
Mse
0.0463
0.0458
0.0500
0.0513
0.0255
0.0178
0.0176
0.0176
4.8
7.4
5.8
12.0
Mean
0.9921
0.9859
1.0015
1.0029
0.4937
0.4856
0.5034
0.4975
Var
0.0552
0.0523
0.0550
0.0535
0.0173
0.0091
0.0098
0.0092
Mse
0.0552
0.0525
0.0550
0.0535
0.0174
0.0093
0.0098
0.0092
0.2
0.0
0.0
0.0
%Tr 0.25
0.0400 27.6
Var
%Tr 0.50
0.0457 41 .o
STr 0.75
40.2
0.1071
%Tr 0.95
45.8
Mean
%Tr 20
35.0
0.8812
VW
%Tr 0.25
s^
0.1031
‘%Tr 0.50
$x
Mean
%Tr 0.75
Lx_
Var %Tr 0.95
8,
Mean
0.9986
0.9980
1.0054
1.0065
0.2411
0.2416
0.2505
0.2507
Var
0.0466
0.0458
0.0491
0.0476
0.0045
0.0022
0.0024
0.0024
Mse
0.0466
0.0458
0.0492
0.0477
0.0046
0.0023
0.0024
0.0024
0.0
0.0
0.0
0.0
%Tr
By the strong law of large numbers, T,(e) + 2-’ w.p. 1, and it is not difficult to verify that the same property holds for T,(x) so & is strongly consistent for 6. Although the central limit theorem readily entails that as n + cc n’/‘[T (a) - 2-‘] is asymptotically N(0, 2P2s), the replacement of 8 by x brings some intricacies in the derivat;on of the limit distribution and also changes the end result, as the following theorem shows. 42
Table 2 Asymptotic
efficiencies
of the simple estimators
of the dependence
parameter
S
6
AU&)
AWk,)
A&&)
1 .oo 0.99 0.95 0.90 0.80 0.70 0.60 0.50 0.40 0.30 0.20 0.10 0.05
0.01 0.11 0.30 0.39 0.48 0.52 0.56 0.58 0.60 0.61 0.62 0.61 0.61
0.01 0.11 0.2x 0.32 0.26 0.17 0.10 0.05 0.02 0.00 0.00 0.00 0.00
0.02 0.17 0.46 0.61 0.71 0.87 0.96 0.99 0.99 0.99 0.99 0.99 0.99
Theorem
3.1. The asymptotic distribution o~Pz”~[T,(X) qf = 3(2P2s)
- 2-j’+’
+ (2-*‘-
- 2-‘]
2-‘+‘)T2(8
is N(0, u,‘) where
+ l)/I’(26
+ 1).
(3.3)
Proof. Since r,(x) is invariant with respect to scale changes in the marginals, we can w.1.o.g. assume that expansion of $ =2, = 1 i.e. (X,,, X2,), i = 1,. . , n, are i.i.d. GBVE(1, 1, 8). To arrive at an appropriate r,(X) we first consider a single pair (X,, X2) and define s= T(Y)
[X,
X/y,,
where Y, and y2 are arbitrary T=
X,Z(S)
+ X21(F) T(Y)
G X,/YZl>
s, = [X,/Y,
T=
X2/~2},
positive
T(1)
Y = (Y,? Y2L = min{ X,, X2},
We then have T(y) = (X,/y,)Z(
constants.
(3.4) S,) + ( X2/y2)I(S;),
and
so
- T= X,(1/y,
- 1)1(S)
+ X2(1/Y, = (l/u,
+ (X,/Y,)[
- l)I(S”)
- l)TI(S)
I&)
--I(S)]
+ (XIL/Y*)[Q)
+ (l/v,
-
- l)TI(S’)
W)]
+ W(Y),
(3.5)
where W(Y) = (X,/YJ[qs,q
-@;s)]
+ (x:/Y,)[qs:‘s)
-QS’)]
=(~,/Y,-x2/Y2)~(~,~“)+(x2/Y,-x,/Y,)~(~,’~).
Now, we refer to the random T(Y)
=min{X,,/~,,
S,=[Xlr~X2,19
sample X2,/~21,
x,
(X,,,
(3.6)
X2,), i = 1,.
, n, and denote
7; = 7; (1) = min{
= t JfI,/Y, G X2,/YZl>
i=l,
From (3.5) and (3.6), we then have T(X)
- T=
(l/X,
- l)TZ(S,)
+ (l/X2-
1)qqs:’
+ u’;(Q
(3.7) 43
where q(x) is given by (3.6). with X,, X,, S,, S replaced We define T, = n-‘C:=,q, w,(y) = nP’C:=,u/;(y>, Zln = K’ From
f: T,I(S,), I=1
z,,
= n-l
i:
by X,,, X,,, S,Y, S,, respectively,
and y by x.
T,z(S,c).
(3.8)
r=l
(3.7) it follows that 7;,(X)
- T, = (l/X,
- l)Zln
+ (l/X,
- l)Z,,
+ W,(X),
(3.9)
so Tn(b)
-2-s
= v, + W,(X),
where v,=
(T,-2-“)
+ (Z,,/X,)(l
-X,)
+ (Z2JX2)(’
-2,).
The stated result follows from the facts that n ‘j27/ n. is asymptotically - _ and n’/‘E,( x) = o,(l) which are proved in Appendix B. q Theorem Corollary
3.1 along with the relation 3.1. The asymptotic uf = (log 2)-2{
8, = -(log
distribution
3 - 2-‘+I
2)-r
log T,(x)
of n”*( 6, - 8) is N(0,
(3.10) normal
with mean 0 and variance
u(:
yields: CJ:), where q
+(1-2fi+‘)rZ(6+1)/r(26+1)}.
Remark. The representation (3.10) and Appendix Lemma B.l also yield the asymptotic joint normality of g,, & and e;. We simply record that the asymptotic covariance of RI/*( 6, - 6,) and n’j2(8, - 8) is B,(log 2))‘{(1 - 2-s) + (1 - 2”)P(26 + l)/r(2s + 1)). Our second approach to construct a simple estimator of 8 draws from the structural representation of GBVE. It was noted in Section 2 that log( X,/X,) - logistic(q, S) where n = log(0,/0,). Linear functions of order statistics or L-estimators are known to be effective for location-scale families especially in regard to their robustness properties. Although 6 is not a location or scale parameter in GBVE, because of the aforementioned property, an L-estimator of 6 can be taken to be of the form 8r = C:=,arnWCrj, where WC,, denotes the ith order statistic of v = log( X1,/X2,), i = 1,. . _, n, and the constants a!,, are the weights appropriate for an L-estimator of the scale for the logistic family. For finite sample optimality, one can use the weights that provide the best linear unbiased estimator (O-BLUE). These involve numerical integration and are tabulated by Gupta, Qureishi and Shah (1967) for n = 2, 5, 10, 15, 20, 25. We denote the corresponding 8, by gooL. Alternatively, one can use the asymptotically optimal weights a,,J[i/(n + l)] where J(u)=9(a2+3)-1[2u-1+2u(l-u)log{u/(l-u)}],
O
The corresponding 8, is denoted by gAoL. With either choice, n’12(8, - 8) is asymptotically N(0, 0:) where 0: = 9(7r* + 3) -I&2. Also, these are asymptotically equivalent to the estimator 8* which maximizes the (logistic) likelihood of W,, i = 1,. . . , n. The asymptotic efficiency AE( 8,) = i’3/aE is computed by inverting the Fisher information matrix given in Section 2, and as Table 2 shows, 8, has a considerably in some matrices inverse, the AE’s are higher AE than both 8, and 8,. Due to the numerical inaccuracy significant up to the second decimal. A Monte Carlo simulation was carried out to examine the behavior of the simple estimators and to compare with the MLE’s for sample sizes n = 10 and 20. Since the estimators are equivalent with respect 44
to scale changes of the marginals, there is no loss of generality in taking 8, = & = 1 in the entire simulation study. For each S and sample size n, four independent random samples were generated: two of these samples were taken from exp(l), one from uniform(O, 1) and one from binomial(1, S). By means of the representation (2.1), the simulated values were transformed to yield a sample of size n from GBVE(l, 1, 6). Repeating this sampling 500 times, we computed the mean, variance and mean squared error of the of 6 were estimates 4, j = 1, 2, &, s”AOL, goi_ as well as the MLE’s (8,, &, 8). The simple estimators truncated at 1 (ie. an estimated value was taken to be 1 whenever the estimator yielded a numerical value greater than 1). The percentage of truncated (%Tr) was recorded in each case. Because of the complexity of the likelihood equations and the restriction 0 < 6 < 1, the MLE’s were computed by three-dimensional grid search method using the error bounds of 0.005, 0.005 and 0.001 for 6,, 19, and 6, respectively. A different seed was used for each case of S, and each n. The results are summarized in Table 1. To get some idea of the accuracy of the simulated Mse’s of the sample means, we repeated the simulation 25 times (for several cells) and found that the results agreed up to three decimals. For estimation of the scale parameters, the simple estimators 8; = x,, j = 1, 2, perform very well. In fact, as S gets close to 1 (independence), the Mse’s of the sample means tend to be smaller than those of the MLE’s in some cases. Also, the simple estimators gAo,_ and gooL based on the order statistics of log( X,/X,) are competitive with the MLE, especially for smaller values of 6. The estimator go,, has a larger proportion of truncation. The other simple estimator 6, has generally a larger Mse than the others. These findings are in general agreement with the asymptotic efficiency values shown in Table 2. All simulations were performed on the University of Wisconsin-Madison VAX 750 computer under the Berkeley Unix operating system. IMSL subroutines GGAMR, GGUBS and GGBN were used to generate exponential, uniform and binomial random variables.
4. Confidence
intervals
and hypotheses
tests
Large-sample confidence interval and hypotheses test procedures based on the asymptotic theory of MLE’s and likelihood ratio tests can be formulated in the usual manner and we do not address them here. They are difficult to implement because of the computational complexity associated with the MLE’s. Also, the restricted range, 0 < 6 < 1, requires some additional care with the determination of the correct limit distribution (see Self and Liang 1987). Here again our focus is on the development of simple methods that have good properties, are intuitively appealing, and are applicable to small as well as large samples. With regard to the scale parameters, we have already observed the remarkably good performance of the simple estimators 8; = x, both from small sample simulation results as well as asymptotic efficiency. Since 2&‘0, - Xz,, the usual X2-based confidence intervals that are appropriate for a univariate exponential can be used to set exact confidence intervals for the parameters 6, and 0, of the GBVE model. As for the dependence parameters 6, nearly exact confidence interval procedure emerges from the fact that go,,/6 is a pivotal quantity because 60,, is the O-BLUE of the scale parameter of a location-scale family. The exact distribution of this pivotal is not tractable but, as highlighted by Bain (1978) one can simulate the desired percentage points as precisely as one wants by sampling from the standard logistic. One can also take advantage of the existing tables (see Bain 1978, Chapter 7) by working with the alternative pivotal 8*/S where s’* 1s the MLE based on the (logistic) likelihood of v = log( X,/X,,), i=l,...,n. Computationof s”* mvolves numerical iteration but is considerably easier than the computation of the MLE s^ from the original bivariate sample. Besides, no pivotal in terms of s^ or even in terms of the simple estimator &t is known to exist. By virtue of the duality between testing hypotheses and confidence statements, exact or nearly exact hypotheses test procedures for the parameters of GBVE are apparent from the above discussion. Besides testing the null hypothesis of independence (6 = l), the method allows testing an arbitrary dependence hypothesis 6 = 6,. Another problem of general interest is testing equality of the marginals, H: 0, = 8,, for 45
which two simple procedures are suggested: (a) Since P( X, < X,) = a,-‘/“[B,-‘I* + 8y1’6]-1 = (>)i if and only if 8, = ( > )0,, we have an exact test based on the binomial count D, = # { X,, < X2,, i = I,. . . , n }, (b) The null hypothesis is also equivalent to n = 0 where n = log(8,/8,) is the location parameter of the logistic distribution of log( X,/X,). Consequently, we have a nearly exact test procedure based on the statistic +*/g* where (?I*, S*) are the MLE’s computed from the logistic likelihood of y, i = 1,. . . , n. The first test is of course simpler but is not expected to be as effective as the second. In fact, the Pitman efficiency of D relative to the second test turns out to be i irrespective of the true 6.
5. An example
To illustrate the proposed methods, we use the data set analyzed by Gross and Lam (1981). Each of 10 patients were given a standard treatment (Tr 1) and a new treatment (Tr 2) for headache on separate occasions, and the lengths of time required to have relief from headache were recorded. The paired data of relief times (in minutes) are: (8.4, 6.9), (7.7, 6.8) (10.1, 10.3) (9.6, 9.4) (9.3, 8.0), (9.1, 8.8) (9.0, 6.1) (7.7, 7.4) (8.1, 8.0) and (5.3, 5.1). For the plausibility of exponential marginals, Gross and Lam transformed the data by subtracting a threshold value of 5.0 from each data point. Concerned with the problem of testing equality of the two treatments, they considered two models: (i) independent exponentials and (ii) ACBVE of Block and Basu (1974). Under each model, they derived the p.d.f. of the ratio R = X,/X,, calculated the likelihood ratio statistic X based on the n = 10 observations of R, and employed the asymptotic X2 test based on -2 log h. We analyze the same transformed data under the GBVE(8,, e,, f) model. The MLE’s computed by three-dimensional grid search method, are found to be 8, = 3.451, e2 = 2.659, 6 = 0.252, and the simple estimates are 8, = 3.43, & = 2.68, 8, = 0.207, Joe, = 0.263. g*,,,_ = 0.239 and (+*, s”*) = (0.354, 0.251) where n = log(B,/e,). Using the pivotal quantity 6*/6 and the table of its (simulated) percentage points provided by Bain (1978, p. 369) a 95% confidence interval for 6 is found to be (0.166, 0.505). Therefore, the null hypothesis of independence (6 = 1) is rejected. Use the pivotal quantity (?* - n)/6”* and Table 1 of Bain, page 367, yields the 95% confidence interval (0.032, 0.676) for n or equivalently, (1.033, 1.967) for 8,/e,; the null hypothesis H,: 8i = 8, is therefore rejected. Concerning the significance tests, the same conclusions were obtained by Gross and Lam. However, their use of the asymptotic critical values is questionable. Also, the present analysis has the added strength that it provides exact confidence intervals for the parameters of interest. Based on a preliminary simulation study, we conclude that the confidence intervals of 6 and n can also be obtained from the simulated percentage points of the pivotal quantities $r/6 and (?joL -
Remark.
n)/gooL, or &oL/a and (qAOL - n)/%oL.
6. Conclusions Many bivariate exponential models have been proposed in the literature but few are known to combine the advantages of (i) a physical basis, (ii) simplicity of the distribution, (iii) tractability of statistical properties, and (iv) availability of exact methods of hypotheses tests and confidence intervals. Ordinarily, (ii) goes hand-in-hand with (iii) and (iv) but, as the present article shows, the GBVE is an exception. Despite its complicated p.d.f., it has a clear edge over all other existing BVE models in regard to (iii) and (iv) with no shortcomings in (i). 46
Appendix A. Derivation of i,2, ii3 and i,, The following three lemmas list the moment results needed for the derivation of the elements and i,,. The proof of Lemma A.1 is quite similar to that of Lemma 2.2 and is omitted. Lemma
A.1. Let a(6) = l/S
- 1. For the random variable
E(V)=l-tS,
E[V{V+a(G)}-‘1
of it,,
i,,
V defined in (2.1) we have
=6
and E[Vz{V+a(6)}-*] Lemma A.2. E{ V(log V)[V+
=26-1+6P1(1-6)2exp[a(S)]E,[a(6)]. a(6)]-‘}
0
= 6[1 + ‘P(l)].
Proof. We denote E{ V(log V)[(V+ a(6)]-‘} = k, and we define r(j) = E[V,/,(log V){ V,, + a(S)}-‘], j = 0, 1. Since r(2) = E(V,, log V,,) - a(S)r(l), we have k = (1 - S)r(l) + 6r(2) = 6E(VI:, log VI,), and the stated result then follows from
1% v,,)
E(Kl
/,
Lemma A.3. For the random variable
V defined in (2.1) we have
(4
E(log
V) = 6 + P(l),
(b)
E[(log
(c)
E( v log V) = (1 + 26) + 9(1)(1
(4
E[ v(log
v)?]
Proof. Differentiating
I0
= (1 - S)[ \I;(l)
v)*] = (1 - S)[ q,(2)
+ \k2(1)]
+ 4 %(2)
+ \k2(2)],
+ S), + q*(2)]
both sides of the identity
r(t
+ 26[ \k,(3) + q2(3)]. + 1) = JgPx’ exp( -x)
dx with respect
to t, we obtain
mxr(logx)exp(-x)dx=!P(t+l)T(r+l),
~~x’(logx)2exp(-x)dx=[!Pl(t+l)+!P2(t+l)]T(t+l). 0
The stated results follow by substituting t = 0, 1, 2 in the above relations \k(t + 1) = ‘P(t) + l/t and ‘r;(t + 1) = q’(t) - l/t2. 0 TO evaluate 42
Taking
i,, we recall that aZ,/ae, =
and aZ,/a0,
= 0, j = 1, 3, so relation
M~2)-14z2~12.
the derivative 4,
= -tI;‘Z,
relations
(2.7) yields (A.11
of h, given in (2.10) with respect
= -(aQ/aZ,)(aQ/az,)[(L
and using the recursive
-
‘WQ-~
+(a2p/azlaz~)[(l-2/6)~-*+{~+~(6)}-’-i].
to Z,, we obtain
+ {Q + QW’] (A.2)
Using aQ/aZ, = Q1/‘/sZ:/s-l, j = 1, 2, and a2 Q/aZ,aZ, = (1 - 1/6)Q’-2’“(Z,Z2)““~‘, we substitute (A.2) in (A.l) and then referring to the representation (2.1) and its consequences Q = I/, Z, = U6V and Z, = (1 - U)“V, we obtain I9182112= - U(1 - o’)[(l +(1-l/S)U(l-
- 2/6) + P{
v+ o(s)}P2]
U)[(l-2/S)
-I- V{v+a(6)}P’-
v].
Use of Lemma A.1 then leads to the expression of i,, given in (2.8). To evaluate i,,, we note that I,, = -Z,X,,/t$ and h,, = -S-‘Z;’
+ S-‘(aQ/aZ,)[2Q-r
-(aQ/aZ,)(aQ/W[(L
-
+ {Q + a(S)}
ww*
-‘I
+ {Q + 4w2]
+(a2Q/az,as){(i-2p3)~-1+[~+a(s)]m’-i}, where aQ/M = a-‘{
Q log Q - Q’-‘/sB(S)}
a2Q/az,as
and {log Q - I!-’
= K’(Q/Z,)‘-“”
log z, - (1 - 1,,‘6)Q_“%(6)},
and B(6) = Z,“s log z, + z, r/6 log Z,. Rep resentation of Q and its derivatives as functions of U, V and use of Lemmas 2.2, A.1 and A.2 yield the stated expression of i,,. Finally, from (2.6) and (2.7) we obtain I,, = h,, = 2K3(log
Z, + log Z,) - 4K3
-(l-2,‘6){Q-*Q2-Q-IQ} +[Q+a(6)]P’[Q+26-3] where Q = aQ/M
log Q + 46-*Q-‘Q - [Q+a(6)]-2{Q--6-2}2
-Q,
and
Q = i3*Q/&S2 = -S-*Q
log Q + S-‘Q
- S-‘Q’-““B(6)[6P2
log Q + S-IQ + S-2Q’P”6B(S) log Q - { +)/Q}Q]
+ K3Q’-“‘C(6),
and C(6) = Zr’/“(log Z,)’ + Z:/‘(log Z2)2. Representing Z,, j = 1, 2, Q and its derivatives as functions of U and V, and referring to the expectations of functions of U and V given in Lemmas 2.2, A.l-A.3, we obtain i,, as stated in Theorem 2.1.
Appendix
B
n, be i.i.d. GBVE(l, 1, a), and V, as defined in (3.10). Then nl”Vn Lemma B.l. Let (Xl,, X2,), i= l,,.., is asymptotically normal with mean 0 and variance qf, where CT: is given in (3.3). Proof. If (X,, X2) - GBVE(l, 1, S) and T= min( XI, X,) we have Cov( X,, X2) = p and Cov( XI, T) = ulr 2T26 - 2-’ + T2(6 + l)/r(26 + l), so according to the multivariate CLT, the asymptotic distribution of J/*,(X, - l), (X2 - l), (T, - 2-6)] IS normal with mean 0 and covariance matrix Z = (a,,) and a,, = a,, for i, j = 1, 2, 3, u,, = u2* = 1, u33 = 2-*‘, ur2 = p and u,~ = u23 = ulr_ Using this, the stated result follows from the fact that 2,/x, + E(Z,) w.p. 1, and E(Z,)
=E[TZ(X,
=E(T)P(X,
< X2) =2-(‘+l).
Here, we have used the fact that T and Z( XI < X,) are independent. 48
0
Lemma B.2. Let W,(y) be as defined in (3.8). Then, n”2W,,(X) Proof. For arbitrary
0 < 7 < co, consider
N,(T) =
a sequence
3 0.
of neighborhoods
of 1 as
(y: Il(Y,‘>VT’) - (1~1)l(--“*J~
where I/. 11denotes the Euclidean norm (the choice of the reciprocal notational simplicity). For an arbitrary fixed E’ > 0, we define P,,(T)
Considering
= Pr
[
sup ln1’2W,(y)( VEUT)
the event [x E N,(T)]
(l/y,,
l/y,)
is only
for
1, Pin=
Pr[XeN,(r)].
and its complement,
> E’]
Pr[l n”%QX))
> E’
scale
we then have the inequality
+ P*,(T).
(B.1)
We are to show that for any E > 0, the right-hand side of (B.l) is smaller than E for all sufficiently large n. Since n’/‘(X - 1) is bounded in probability, we can choose 7 such that P2,,( r) < SE for all sufficiently large n. Henceforth 7 is fixed at this chosen value, and we proceed to show that pi,,(~) -+ 0 as n + cc. In order to bound n1/2W,(y) in the set y E N,(r) we refer to the definition of W(y) in (3.6). In the set S,,S’ we have (n/n)x,
GX*
T= X2,
G (Yl/Y2)X2?
which implies T(l/y, - l/y*) < X,/y1 - X,/y, < 0. Similarly, in the set S;S, we have T = X, and T(l/y, - l/y,) < X2/y2 - X,/y, d 0. We note that both these sets are empty if yi = y2. If y, # y2 one of the sets is empty. Since 1 - ~n-‘/~ < l/y, < 1 + rn-‘I*, i = 1, 2, we have 11/y, - l/y, 1 < 2rn-1/2, so IW(Y)I~2~n Hence,
-1/2T[Z(S,‘)
+Z(S,‘S)]
for YEN,(T)
for each i = 1,. . . , n, we have
(B.2) We now note that S,,S:: = 1x2,
< Xl, * X1,/y,
< X2,/y,]
c
[I
< &/X2,
< 1 + T’npl”],
for sufficiently large n and T’ = 37, and likewise, S&?, c [l - T’n-“‘* -c X1,/X2, < 11. Denoting A, = [l -r’np’/2, 1 + *‘n-1’2], and referring to (B.2) we then have,_for y E N,(T), 1w(y) I < T’n-“*T,Z( X1,/X2, E A,). This implies sup, E N?ol r?*W,, ( y) 1Q T ‘v,, where V, is the average of i.i.d. y = q.Z( Xu/X2, E n. Now, considering V = TZ( X,/X, E A,) and denoting p,, = Pr( X,/X, E A,), we have by A,), i= l,..., the Cauchy_Schwarz inequality, E(F,) = E(V) < { E(~)P,}“~. Since E(T*) < cc and p, -+ 0 as n + co, we have E( V,) ---f0 which implies that supyE Nn(7jln”‘W, ( y) 1% 0, and the proof is concluded. 0
Acknowledgement The research for this article AFOSR-87-0256.
was sponsored
by the Air Force
Office
of Scientific
Research
under
Grant
49
References Abramowitz, M. and I.E. Stegun, Handbook of Mathematical Funcrions (National Bureau of Standards, Washington, DC). Bain, L.J. (1978) Statistical Analysis of Reliability and LijeTesting Models (Dekker, New York). Block, H.W. and A.P. Basu (1974), A continuous bivariate exponential extension, J. Amer. Statist. Assoc. 69, 10311037. Freund, J.E. (1961) A bivariate extension of exponential distribution, J. Amer. Statist. Assoc. 56, 971-977. Gross, A.J. and C.F. Lam (1981) Paired observations from a survival distribution, Btometrics 37, 5044511. Gumbel, E.J. (1960) Bivariate exponential distribution, J. Amer. Statist. Assoc. 55, 698-707. Gupta, S’S,, A.S. Qureishi and B.K. Shah (1967) Best linear unbiased estimators of the parameters of the logistic distribution using order statistics, Technometrics 9, 43-56. Hougaard, P. (1986a), Survival models for heterogeneous populations derived from stable distributions, Biometrika 73, 387-396. Hougaard, P. (1986b), A class of multivariate failure time distributions, Biometrika 73, 671-678.
50
Johnson, N.L. and S. Kotz (1972) Distributions in Statrstics: Continuous Multivariate Distributions (Wiley, New York). Lee, L. (1979), Multivariate distributions having Weibull properties, J. Multivariate Anal. 9, 267-277. Lu, J.C. and G.K. Bhattacharyya (1988), Some bivariate extensions of the Weibull distribution, Tech. Rept. No. 821, Univ. of Wisconsin (Madison, WI). Mann, N.R., R.E. Schafer and N.D. Singpurwalla (1974) Methods for Statistical Analysis of Reliability and Lifetime Data (Wiley, New York). Marshall, A.W. and I. Olkin (1967) A multivariate exponential distribution, J. Amer. Statist. Assoc. 62, 30-44. Sarkar, S.K. (1987) A continuous bivariate exponential distribution, J. Amer. Statist. Assoc. 82, 667-675. Self, S.G. and K.Y. Liang (1987) Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions, J. Amer. Statist. Assoc. 82, 605-610.
Letters 12 (1991) 37-50
July 1991
Inference procedures for bivariate exponential model of Gumbel Jye-Chyi Department
Lu of Statistics,
North Carolina State University, Raleigh, NC 27695-8203,
USA
Gouri K. Bhattacharyya Department Received
of Statistics,
University of Wisconsin,
Madison,
WI 53706, USA
May 1990
Abstract: In this article, statistical properties of an absolutely continuous bivariate exponential model due to Gumbel (1960) are examined and the Fisher information matrix is derived. Simple estimators are motivated from a structural representation of the model and are compared both in terms of asymptotic efficiency and finite sample simulation. Besides their good performance, some of the simple estimators lead to exact or nearly exact confidence interval and hypotheses test procedures. The proposed methods are illustrated by reanalyzing a data set to which asymptotic methods were employed by Gross and Lam (1981) under a different bivariate exponential model. Keywords:
Dependent
lifetimes,
Fisher information,
simple estimators,
asymptotic
efficiency,
confidence
intervals,
simulation.
1. Introduction Many bivariate extensions of the exponential distribution have been proposed in the literature in order to model dependence of the lifetimes of paired organs or components in a system. The models derived from physical behavior of the failure process are generally more dependable in reliability studies than those that are just mathematical constructs, and are especially so when a large enough sample is not available for a test of fit. Tractability of statistical inferences is yet another criterion for a model to be useful in analyzing life test data. In both these regards, the bivariate exponential model of Marshall and Olkin (1967) has considerable advantages over the others. However, because of a singular part, its use is not appropriate in the situations where simultaneous failures of the components are quite unlikely to occur. Freund’s (1961) bivariate extension of the exponential is also well-grounded on a physical basis and is absolutely continuous but does not have exponential marginals. Two absolutely continuous modifications of the Marshall-Olkin model have been proposed: ACBVE by Block and Basu (1974) and ACBVE, by Sarkar (1987). The former is a sub-class of Freund’s model while the latter does not seem to have as strong a physical basis and is not convenient for statistical inferences. In this article, we investigate the properties and inference procedures for a rather neglected bivariate exponential model whose survival function and probability density function (p.d.f.) are respectively given by F(x,,
x2) = exp
I
-
{
(x,/0,)“’
+ (xa/Ws}
(l.la)
“1)
(x,/e,)'/"]"-2 f(x,, ~,)=(s,e,)~'~"(x,x~)~~~-'[(~~/e~)~~*+ + (x2/8Jis}'+
i/s-
i] exp[-
{ (x1/e,)1/6+(x2/e2)1/'}6]. (l.lb)
0167.7152/91/$03.50
0 1991 - Elsevier Science
Publishers
B.V. (North-Holland)
37
Volume
12, Number
1
STATISTICS
& PROBABILITY
LETTERS
July 1991
0 < x,, x2 < cc, 0 < 0,, 19~< cc, 0 i S < 1. It is a transformation of Gumbel’s type B bivariate extreme-value distribution (cf. Johnson and Kotz, 1972, p. 251) proposed by Gumbel (1960) and we will henceforth refer to it as GBVE(8,, 8,, S). Here, 8, and 8, are the scale parameters, S is the dependence parameter, and 6 = 1 corresponds to the case of independence. In spite of its long history, the GBVE has received little attention in modeling dependence of paired lifetimes seemingly for two reasons: (i) it was proposed solely as a mathematical construct and (ii) the complicated form of its p.d.f. poses considerable difficulties with statistical inference. That it is indeed a meaningful physical model surfaces from Hougaard (1986b) who provides a general construction of dependent life models from the consideration of random stress environment affecting both components (see Lu and Bhattacharyya, 1988, for other derivations and extensions of this model). Specifically, if conditionally given the stress S = s, the component lives X, and X, are independent Weibull with survival i = 1, 2, respectively, and if S has a positive stable functions exp[-H,(t)s] and H,(t) = (r/0,)“S, distribution (cf. Hougaard, 1986a), then the unconditional distribution of (Xi, X,) is GBVE(e,, 0,, 6). Thus, while the source of dependence in the Marshall-Olkin model is a common shock destroying both components, and that in the Freund model is a failure event internal to the system, it is the random mixing effect of an external stress that serves to explain the dependence in GBVE. Such stochastic failure rate models are long recognized in reliability modeling, see for instance, Mann, Schafer and Singpurwalla (1974, p. 145). With the recognition that GBVE has a sound physical basis in modeling bivariate life times, our major goal here is to develop procedures for parameter estimation, confidence statements and hypotheses tests. Some statistical properties of the model, including the derivation of the Fisher information matrix, are presented in Section 2. Despite the complicated form of the p.d.f. due to which the MLE’s could not be derived in closed forms, it has been possible to obtain closed form expressions for the elements of the Fisher information matrix. These are especially handy in evaluating the asymptotic efficiency (AE) of some simple estimators constructed in Section 3. A representation, due to Lee (1979) of (X,, X,) in terms of certain independence random variables proves instrumental for these calculations and also for the development of some simple inference procedures. Based on a random sample (Xi,, X,,), i = 1,. . . , n, from GBVE(B,, 8,, S), some simple estimators of the model parameters are constructed in Section 3. These are compared both in terms of AE and Monte Carlo simulation. For the scale parameters S,, the sample means x1 = n-‘C:=,X,;, j = 1, 2, are found to retain a remarkably high AE under the bivariate model and they also do very well in small samples. A simple estimator of the dependence parameter 6 is motivated from the exponential distribution of min( X,, X,), and its AE is evaluated. Other simple estimators are constructed from the order statistics of log( X,/X,,), i = 1,. _. , n, or by referring to their likelihood. These are found to be substantially better than the first and quite competitive with the MLE. This latter construction draws from the property that log( X,/X,) generates a location-scale logistic family, which itself is a consequence of the independence representation of GBVE. Problems of confidence intervals and hypotheses tests are taken up in Section 4. Pivotal quantities are identified for the basic model parameters as well as for (3,/e,, the ratio of the component mean lives. Exact, nearly exact and asymptotic procedures are discussed. Applications of the proposed methods are illustrated by means of the data set of a paired comparison experiment which was previously analyzed by by Gross and Lam (1981) under the ACBVE model.
2. Properties
of GBVE
and the Fisher
information
Several distributional properties of GBVE can be conveniently deduced from a structural representation due to Lee (1979) in terms of independent random variables. A slight variant of this result is given in the as”, exp(a) for an exponential following lemma. The notation - will be used to mean “is distributed 38
distribution with mean parameter h.
a, and logistic( a, b) for a logistic distribution
with location
parameter
a and scale
Lemma 2.1 (Lee, 1979). Consider independent random variables U, V,,, V,, and M6 such that U where uniform(O, I), V,, - exp(l), i = 1, 2, and M, = 0 or 1 with probability 1 - 6 and 6, respectively, 6 E (0, l] is a constant. Letting 8, and f& be positive constants, define V= V,, + MJl:2,
Then (X,,
x, = 6,U6V,
x, = f?,(l - U)%.
X,) - GBVE(8,,
6,, 6).
(2.1)
0
The moments of GBVE are readily obtained from (2.1) and the facts that b + 1) and E(V”“) = (1 + a6)F(a + l), a, b > 0. In particular, E( X;lX;‘)
= 0;l@F(
k,, k, > 0, and the correlation p(X,,
k,S + l)F( coefficient
k,S + l)F(
k, + k, + l),‘F(
E[ZY’(l - U)h] = B(a + 1,
k,S + k,6 + l),
(2.2)
of X, and X, is given by
x,)=2r2(6+1)/r(26+1)-1,
(2.3)
which ranges from 0 to 1 and is monotone decreasing in 6. Let T = min( X,, X,), W = log( X,/X,) and D = I( X, < X,) where Z( .) denotes the indicator function. Then, the following properties hold: (a) W - logistic(7, 6) with n = log(0,/8,), (b) Pr( D = 1) = and (d) T and D are independent. Property (a) e;‘/s/(e;l/s + ep”), (c) T - exp[(e, p’/s + e;‘/‘)-“1, readily follows from (2.1) and the fact that log[U/(l - U)] - logistic(0, 1). Also, since X, < X, is equivalent to log[U/(l - U)] <) - q/6, (b) follows from the logistic c.d.f. Verification of (c) and (d) involves some routine algebra which we omit for brevity. Of course, (c) is obvious from the survival function (l.la). Incidentally, the Marshall-Olkin model also has similar properties (c) and (d). It is commonly known that the logarithm of the ratio of two i.i.d. exp(1) random variables has the standard logistic distribution. As a dependent exponential model, GBVE has the interesting feature that the logistic family is preserved under the aforementioned transformation, and the dependence parameter 6 becomes a scale parameter. Consider the function
Q(X,,
s
X2>= [w4Y’~s + (x,/B:)“‘]
(2.4)
which appears in the exponent of the survival function (1.1). From the representation (2.1), we note that Q( X,, X,) = V which is a binary mixture of independent exp(1) and gamma(1, 2). In particular, for 6 = 1 (independence), it has the &i-square distribution with d.f. = 4. Expectations of some functions of Q will be needed in connection with the calculation of the Fisher information matrix for GBVE. Its identification with I’, which has a simple distribution, facilitates the work. A few key results are listed in the following lemma. Lemma
2.2. Let a(S) = l/6
- 1. For the random variable
(a>
E[{V+a(G)}-‘1
=E[V{V+a(G)}-‘]=a,
@I
E[{V+a(S)}-‘1
=aq0,,
cc>
E[V{v+a(G)}-*I
=6-
V defined in (2.1) we have
(1 -6)q0,,
where qol = exp[a(6)]E,[a(6)] and E,(t) = jp”x-’ exp( - x) dx is the classical exponential for which extensive tables are available (cf. Abramowitz and Stegun, 1970).
integral function
39
Proof. Let V= V’, + V’, where V” and V’, are independent exp(l), and for integers j > 0, k > 1, define the expectation of the random variable Rlk = V’[ V + q(j, k) = E[V/,{ V’, + ~(6>}~‘1. C onsidering a(6)lek where v= V’, + M,V’,, and first taking expectation over M,, we obtain
E(R,,)=(l-S)q(j,
k)+SE[V1{V+u(8)jPk]
= (1 - S)q(j,
k) + Sq( j + 1, k) = 6q( j, k - 1).
(2.5)
The second equality follows from the fact that V.- gamma(1, 2) so E[g( V)] = E[V”g( V,‘)] holds for any integrable function g. The last equality obtains from the identity (1 - S) + SV,, = S[T/,’ + a(6)]. Taking (j, k) = (1, 1) and (0, 2), the first two of the stated results follow from (2.5) once we note that q(1, 0) = 1 and qol = q(0, 1) = exp[4U4[4W For the last result we take (j, k) = (1, 2) and use the relation q(1, 1) = 1 - a(6)q(O, 1). q 2.1. Fisher information
for GBVE(B,,
0,, 6)
Because 8, and 8, are scale parameters in the model, it will be convenient to work with scaled variables Z, = X,/0,, j = 1, 2, whose joint p.d.f. will be denoted by g( z,, z2, 8). Define X = log g( z,, z2, 6) and Q = (Z:‘fi + z:‘“)s and referring to (l.l), we have X=a(G)(log The log-likelihood I=
Z, +log
Z,) + (l-2/6)
log Q+log[Q+a@)]
-Q.
(2.6)
I of ( X,, X,) is then given by
-loge,-log8,+A.
(2.7)
For the convenience of notation, we identify 6 as the argument 19’ in I and as Z, in X, denote 8 = (8’, 19,, I!?,)‘, Z= (Z’, Z,, Z,)’ and view X as a function of Z. Also, we denote the partial derivatives i,= algae,, I,k = a21/ae,ae,, x,= ax/az,, x,,= aQ/az,az,, j, k = 1, 2, 3. Theorem 2.1. The Fisher information 1” = E( -I’,)
matrix
= (38,?-‘[26
i,, = (3t3’0z)P’[6
Z = (i,,)
= E( - I,,)
of GBVE( 0,, f?,, S) is given by
+ K2 + (1 - S)‘S-‘qO,],
- K2 + (26))‘(1
i,, = (0’/82)2i’,,
- S)‘qo’], (2.X)
i ,3 = (20,)P’[SP’ ‘33
-S
- (26))‘(1
- 6){2-‘(1
- 6) + K’}qo’],
i,, = (t3’/B2)i13,
=6-‘(1/S-1)+2q,(S-6-~)+qo,[q’S-‘(l-6)2-S~~(1+S+S~‘)],
where qo, is as defined in Lemma 2.2 and q, = $ + f[{ q(2) - S(4)}2 - ‘k,(4)] = 0.2835533, where ‘P(t) = d log r(t)/dt and q’(t) = d \k(t)/dt, are the digamma and trigumma function, respectively. Proof.
Recalling
that aZ’/M’
I,, = e;2[1+ Differentiating
= -0;‘Z’
and aZ,/M,
= 0, j = 2, 3, relation
(2.7) yields
2Z,X’ + Z’Jh”].
(2.9)
(2.6) we obtain
h,=a@)Z;‘+(Z,/Q)“‘“‘[(1-2/S)QP’+{Q+a(8)}-l-l],
C-?/Q) 20(B’[(1-2/S)QP2+
A,’ = -@)Z,‘+ a(6){
40
Q’-‘/9;/S-2
_
Q’-2/9;/8-
{Q+a(s)}-‘1 2}[(1-2/6)Q-‘+
{Q+a(8)}-l-11.
(2.10)
Substituting (2.10) in (2.9) and referring and Z, = U&V we obtain @I,, = 1-t 2a(6) -a(6)
+ 24(1-
Use of Lemma 2.2 then leads provided in Appendix A. 0
3. Estimation
2/6)
cP[(l -
+a(s)u(l-
to the representation
2/S)
+ v{ V+ a(S)}
+ vZ{ v+
o’)[(l-2/6)+ to the stated
(2.1) and its consequences
a(S)}
-’ - v] -‘I
Y{V+a(G)}_‘expression
such as Q = V
v].
for i,,. The derivations
of i,,,
ii,
and
i,,
are
of parameters
In this section, we consider estimation of the parameters of GBVE( 8,, t?,, 6) from a random sample (X,,, X2,), i = 1,. . . , n. The density model (1.1) is not an exponential family, and thus no drastic sufficient reduction is possible. The maximum likelihood estimators have to be calculated through some numerical procedures. Our goal here is to construct simple estimators which are intuitively appealing, perform well in small or moderate samples and do not incur much loss of asymptotic efficiency compared to the MLE’s. Since the marginal distributions are exponential with means 6, and &, it is natural to estimate 0, by $ = 3 = nP’xy=, X,,, j = 1, 2. These are simple estimators with known exact properties such as: unbiased, Var(0,) = 8:,/n and 2n8;/8, - ~2,. The asymptotic efficiency, AE( 8,) = i”/Bf which depends only on the parameter S, is computed by inverting the Fisher information matrix given in Section 2. The numerical results show that AE(&,) has a U-shaped curve and attains its minimum 98.215% at S = 0.725. Besides their uniformly high AE and good performance in small samples (see Table l), the fact that these estimators are robust with respect to the dependence structure of the bivariate model makes them very useful for practical application. As for the dependence parameter (S), a moment-type estimator based on the sample correlation Y is a natural choice. In view of (2.3) this estimator s”,, is the solution of the equation l+r=2P(8+1)/r(26”+1),
(3.1)
if r > 0 and is defined to be 1 if r < 0. To obtain the asymptotic distribution of this moment recall that the asymptotic distribution of nil2 (p - p), where b is the sample product-moment coefficient is normal with mean 0 and variance crP2as P2(~22~~I’+4~‘[~40~2C,1+~22(~~2~2~)~’+~ir4~011]
estimator, we correlation
-[O~~~((Y:U(YI,)-‘+~~~((YO~(YII)-’]),
where (Y,, = E{[ X- E( X)]‘[Y - E(Y)]‘}, and (X, Y) has a general bivariate distribution. In the GBVE distribution, the moments aA, can be obtained using the expression (2.2). By an application of delta method it follows that ~“~(8, - S) is asymptotically N(0, cr$) where CT; = u;166’r2(26
+ 1)r-4(6
+ l)[ \k(S + 1) - 9(26
+ l)] -2.
The lack of a closed form expression s”, and the lowness (cf. Table 2) of its asymptotic efficiency are obvious disadvantages. Two approaches for the construction of simple estimators of S turn out to be fruitful. The first draws from the fact that the scaled pair (X,/8,, X,/e,) - GBVE(1, 1, 6) so 7;(o) = min( X,/8,, X2,/e,), i=l ,...> n, are i.i.d. exp(2-S). If 0 were known, 7;,(o) = n-‘Z3f=,7;(8) could be used to estimate 2-“. Substitution of 2 = (x1, x2)’ for 8 motivates the estimator 8,= -(log2)P’T,(x)=
-(log2))W
,g, min { X1,/X, ? x2,/%
}.
(3.2)
41
Table
1
Results of a Simulation n
6
10
1 .oo
from GBVE(l, 61
1, 6) e,
82
B,
0.9904
0.9870
1.0148
1.0174
0.8427
0.8434
0.8850
0.1077
0.1038
0.1106
0.0313
0.0277
0.0225
0.0216
Mse
0.1032
0.1079
0.1040
0.1109
0.0561
0.0522
0.0357
0.0357
39.4
1 .oo
1.0038
1.0038
1.0026
1 .oooo
0.8360
0.8322
0.8737
0.8730
0.1096
0.0951
0.1023
0.0327
0.0261
0.0240
0.0212
Mse
0.1071
0.1096
0.0951
0.1023
0.0298 41 .o
0.0272 37.2
Mean
0.9986
0.9909
0.9968
0.9887
0.6969
0.6762
0.7331
0.7141
0.0962
0.0949
0.0909
0.0916
0.;0390
0.0268
0.0282
0.0262
Mse
0.0962
0.0950
0.0909
0.0918
0.0418
0.0322
0.0285
0.0275
13.2
6.6
12.4
7.2
Mean
1.0002
0.9871
1.0050
1.0037
0.4662
0.4546
0.4976
0.4820
Var
0.1092
0.1058
0.1099
0.1065
0.0301
0.0166
0.0199
0.0174
Mse
0.1092
0.1059
0.1099
0.1065
0.0313
0.0187
0.0199
0.0177
1.2
0.0
0.2
0.0
Mean
1.0044
1.0049
1.0030
1.0025
0.2271
0.2229
0.2459
0.2393
Var
0.0935
0.0929
0.1026
0.0997
0.0073
0.0043
0.0054
0.0049
Mse
0.0935
0.0929
0.1026
0.0997
0.0079
0.0050
0.0054
0.0051
0.0
0.0
0.0
0.0
0.9238
0.9357
Mean
1.0063
1.0026
1.0027
0.9975
0.9015
0.9067
Var
0.0496
0.0515
0.0467
0.0471
0.0161
0.0127
0.0109
0.0086
Mse
0.0497
0.0516
0.0467
0.0471
0.0258
0.0214
0.0167
0.0127
47.6
41.0
50.4
48.6 0.9045
Mean
0.9864
0.9828
1.0122
1.0126
0.8760
0.8811
0.9011
Var
0.0447
0.0455
0.0484
0.0493
0.0182
0.0142
0.0125
0.0105
Mse
0.0449
0.0458
0.0485
0.0495
0.0236
0.0189
0.0149
0.0126
34.0
29.6
37.6
33.0
Mean
1.0047
1.0002
1.0137
1.0073
0.7506
0.7253
0.7507
0.7462
Var
0.0462
0.0458
0.0499
0.0512
0.0255
0.0172
0.0176
0.0176
Mse
0.0463
0.0458
0.0500
0.0513
0.0255
0.0178
0.0176
0.0176
4.8
7.4
5.8
12.0
Mean
0.9921
0.9859
1.0015
1.0029
0.4937
0.4856
0.5034
0.4975
Var
0.0552
0.0523
0.0550
0.0535
0.0173
0.0091
0.0098
0.0092
Mse
0.0552
0.0525
0.0550
0.0535
0.0174
0.0093
0.0098
0.0092
0.2
0.0
0.0
0.0
%Tr 0.25
0.0400 27.6
Var
%Tr 0.50
0.0457 41 .o
STr 0.75
40.2
0.1071
%Tr 0.95
45.8
Mean
%Tr 20
35.0
0.8812
VW
%Tr 0.25
s^
0.1031
‘%Tr 0.50
$x
Mean
%Tr 0.75
Lx_
Var %Tr 0.95
8,
Mean
0.9986
0.9980
1.0054
1.0065
0.2411
0.2416
0.2505
0.2507
Var
0.0466
0.0458
0.0491
0.0476
0.0045
0.0022
0.0024
0.0024
Mse
0.0466
0.0458
0.0492
0.0477
0.0046
0.0023
0.0024
0.0024
0.0
0.0
0.0
0.0
%Tr
By the strong law of large numbers, T,(e) + 2-’ w.p. 1, and it is not difficult to verify that the same property holds for T,(x) so & is strongly consistent for 6. Although the central limit theorem readily entails that as n + cc n’/‘[T (a) - 2-‘] is asymptotically N(0, 2P2s), the replacement of 8 by x brings some intricacies in the derivat;on of the limit distribution and also changes the end result, as the following theorem shows. 42
Table 2 Asymptotic
efficiencies
of the simple estimators
of the dependence
parameter
S
6
AU&)
AWk,)
A&&)
1 .oo 0.99 0.95 0.90 0.80 0.70 0.60 0.50 0.40 0.30 0.20 0.10 0.05
0.01 0.11 0.30 0.39 0.48 0.52 0.56 0.58 0.60 0.61 0.62 0.61 0.61
0.01 0.11 0.2x 0.32 0.26 0.17 0.10 0.05 0.02 0.00 0.00 0.00 0.00
0.02 0.17 0.46 0.61 0.71 0.87 0.96 0.99 0.99 0.99 0.99 0.99 0.99
Theorem
3.1. The asymptotic distribution o~Pz”~[T,(X) qf = 3(2P2s)
- 2-j’+’
+ (2-*‘-
- 2-‘]
2-‘+‘)T2(8
is N(0, u,‘) where
+ l)/I’(26
+ 1).
(3.3)
Proof. Since r,(x) is invariant with respect to scale changes in the marginals, we can w.1.o.g. assume that expansion of $ =2, = 1 i.e. (X,,, X2,), i = 1,. . , n, are i.i.d. GBVE(1, 1, 8). To arrive at an appropriate r,(X) we first consider a single pair (X,, X2) and define s= T(Y)
[X,
X/y,,
where Y, and y2 are arbitrary T=
X,Z(S)
+ X21(F) T(Y)
G X,/YZl>
s, = [X,/Y,
T=
X2/~2},
positive
T(1)
Y = (Y,? Y2L = min{ X,, X2},
We then have T(y) = (X,/y,)Z(
constants.
(3.4) S,) + ( X2/y2)I(S;),
and
so
- T= X,(1/y,
- 1)1(S)
+ X2(1/Y, = (l/u,
+ (X,/Y,)[
- l)I(S”)
- l)TI(S)
I&)
--I(S)]
+ (XIL/Y*)[Q)
+ (l/v,
-
- l)TI(S’)
W)]
+ W(Y),
(3.5)
where W(Y) = (X,/YJ[qs,q
-@;s)]
+ (x:/Y,)[qs:‘s)
-QS’)]
=(~,/Y,-x2/Y2)~(~,~“)+(x2/Y,-x,/Y,)~(~,’~).
Now, we refer to the random T(Y)
=min{X,,/~,,
S,=[Xlr~X2,19
sample X2,/~21,
x,
(X,,,
(3.6)
X2,), i = 1,.
, n, and denote
7; = 7; (1) = min{
= t JfI,/Y, G X2,/YZl>
i=l,
From (3.5) and (3.6), we then have T(X)
- T=
(l/X,
- l)TZ(S,)
+ (l/X2-
1)qqs:’
+ u’;(Q
(3.7) 43
where q(x) is given by (3.6). with X,, X,, S,, S replaced We define T, = n-‘C:=,q, w,(y) = nP’C:=,u/;(y>, Zln = K’ From
f: T,I(S,), I=1
z,,
= n-l
i:
by X,,, X,,, S,Y, S,, respectively,
and y by x.
T,z(S,c).
(3.8)
r=l
(3.7) it follows that 7;,(X)
- T, = (l/X,
- l)Zln
+ (l/X,
- l)Z,,
+ W,(X),
(3.9)
so Tn(b)
-2-s
= v, + W,(X),
where v,=
(T,-2-“)
+ (Z,,/X,)(l
-X,)
+ (Z2JX2)(’
-2,).
The stated result follows from the facts that n ‘j27/ n. is asymptotically - _ and n’/‘E,( x) = o,(l) which are proved in Appendix B. q Theorem Corollary
3.1 along with the relation 3.1. The asymptotic uf = (log 2)-2{
8, = -(log
distribution
3 - 2-‘+I
2)-r
log T,(x)
of n”*( 6, - 8) is N(0,
(3.10) normal
with mean 0 and variance
u(:
yields: CJ:), where q
+(1-2fi+‘)rZ(6+1)/r(26+1)}.
Remark. The representation (3.10) and Appendix Lemma B.l also yield the asymptotic joint normality of g,, & and e;. We simply record that the asymptotic covariance of RI/*( 6, - 6,) and n’j2(8, - 8) is B,(log 2))‘{(1 - 2-s) + (1 - 2”)P(26 + l)/r(2s + 1)). Our second approach to construct a simple estimator of 8 draws from the structural representation of GBVE. It was noted in Section 2 that log( X,/X,) - logistic(q, S) where n = log(0,/0,). Linear functions of order statistics or L-estimators are known to be effective for location-scale families especially in regard to their robustness properties. Although 6 is not a location or scale parameter in GBVE, because of the aforementioned property, an L-estimator of 6 can be taken to be of the form 8r = C:=,arnWCrj, where WC,, denotes the ith order statistic of v = log( X1,/X2,), i = 1,. . _, n, and the constants a!,, are the weights appropriate for an L-estimator of the scale for the logistic family. For finite sample optimality, one can use the weights that provide the best linear unbiased estimator (O-BLUE). These involve numerical integration and are tabulated by Gupta, Qureishi and Shah (1967) for n = 2, 5, 10, 15, 20, 25. We denote the corresponding 8, by gooL. Alternatively, one can use the asymptotically optimal weights a,,J[i/(n + l)] where J(u)=9(a2+3)-1[2u-1+2u(l-u)log{u/(l-u)}],
O
The corresponding 8, is denoted by gAoL. With either choice, n’12(8, - 8) is asymptotically N(0, 0:) where 0: = 9(7r* + 3) -I&2. Also, these are asymptotically equivalent to the estimator 8* which maximizes the (logistic) likelihood of W,, i = 1,. . . , n. The asymptotic efficiency AE( 8,) = i’3/aE is computed by inverting the Fisher information matrix given in Section 2, and as Table 2 shows, 8, has a considerably in some matrices inverse, the AE’s are higher AE than both 8, and 8,. Due to the numerical inaccuracy significant up to the second decimal. A Monte Carlo simulation was carried out to examine the behavior of the simple estimators and to compare with the MLE’s for sample sizes n = 10 and 20. Since the estimators are equivalent with respect 44
to scale changes of the marginals, there is no loss of generality in taking 8, = & = 1 in the entire simulation study. For each S and sample size n, four independent random samples were generated: two of these samples were taken from exp(l), one from uniform(O, 1) and one from binomial(1, S). By means of the representation (2.1), the simulated values were transformed to yield a sample of size n from GBVE(l, 1, 6). Repeating this sampling 500 times, we computed the mean, variance and mean squared error of the of 6 were estimates 4, j = 1, 2, &, s”AOL, goi_ as well as the MLE’s (8,, &, 8). The simple estimators truncated at 1 (ie. an estimated value was taken to be 1 whenever the estimator yielded a numerical value greater than 1). The percentage of truncated (%Tr) was recorded in each case. Because of the complexity of the likelihood equations and the restriction 0 < 6 < 1, the MLE’s were computed by three-dimensional grid search method using the error bounds of 0.005, 0.005 and 0.001 for 6,, 19, and 6, respectively. A different seed was used for each case of S, and each n. The results are summarized in Table 1. To get some idea of the accuracy of the simulated Mse’s of the sample means, we repeated the simulation 25 times (for several cells) and found that the results agreed up to three decimals. For estimation of the scale parameters, the simple estimators 8; = x,, j = 1, 2, perform very well. In fact, as S gets close to 1 (independence), the Mse’s of the sample means tend to be smaller than those of the MLE’s in some cases. Also, the simple estimators gAo,_ and gooL based on the order statistics of log( X,/X,) are competitive with the MLE, especially for smaller values of 6. The estimator go,, has a larger proportion of truncation. The other simple estimator 6, has generally a larger Mse than the others. These findings are in general agreement with the asymptotic efficiency values shown in Table 2. All simulations were performed on the University of Wisconsin-Madison VAX 750 computer under the Berkeley Unix operating system. IMSL subroutines GGAMR, GGUBS and GGBN were used to generate exponential, uniform and binomial random variables.
4. Confidence
intervals
and hypotheses
tests
Large-sample confidence interval and hypotheses test procedures based on the asymptotic theory of MLE’s and likelihood ratio tests can be formulated in the usual manner and we do not address them here. They are difficult to implement because of the computational complexity associated with the MLE’s. Also, the restricted range, 0 < 6 < 1, requires some additional care with the determination of the correct limit distribution (see Self and Liang 1987). Here again our focus is on the development of simple methods that have good properties, are intuitively appealing, and are applicable to small as well as large samples. With regard to the scale parameters, we have already observed the remarkably good performance of the simple estimators 8; = x, both from small sample simulation results as well as asymptotic efficiency. Since 2&‘0, - Xz,, the usual X2-based confidence intervals that are appropriate for a univariate exponential can be used to set exact confidence intervals for the parameters 6, and 0, of the GBVE model. As for the dependence parameters 6, nearly exact confidence interval procedure emerges from the fact that go,,/6 is a pivotal quantity because 60,, is the O-BLUE of the scale parameter of a location-scale family. The exact distribution of this pivotal is not tractable but, as highlighted by Bain (1978) one can simulate the desired percentage points as precisely as one wants by sampling from the standard logistic. One can also take advantage of the existing tables (see Bain 1978, Chapter 7) by working with the alternative pivotal 8*/S where s’* 1s the MLE based on the (logistic) likelihood of v = log( X,/X,,), i=l,...,n. Computationof s”* mvolves numerical iteration but is considerably easier than the computation of the MLE s^ from the original bivariate sample. Besides, no pivotal in terms of s^ or even in terms of the simple estimator &t is known to exist. By virtue of the duality between testing hypotheses and confidence statements, exact or nearly exact hypotheses test procedures for the parameters of GBVE are apparent from the above discussion. Besides testing the null hypothesis of independence (6 = l), the method allows testing an arbitrary dependence hypothesis 6 = 6,. Another problem of general interest is testing equality of the marginals, H: 0, = 8,, for 45
which two simple procedures are suggested: (a) Since P( X, < X,) = a,-‘/“[B,-‘I* + 8y1’6]-1 = (>)i if and only if 8, = ( > )0,, we have an exact test based on the binomial count D, = # { X,, < X2,, i = I,. . . , n }, (b) The null hypothesis is also equivalent to n = 0 where n = log(8,/8,) is the location parameter of the logistic distribution of log( X,/X,). Consequently, we have a nearly exact test procedure based on the statistic +*/g* where (?I*, S*) are the MLE’s computed from the logistic likelihood of y, i = 1,. . . , n. The first test is of course simpler but is not expected to be as effective as the second. In fact, the Pitman efficiency of D relative to the second test turns out to be i irrespective of the true 6.
5. An example
To illustrate the proposed methods, we use the data set analyzed by Gross and Lam (1981). Each of 10 patients were given a standard treatment (Tr 1) and a new treatment (Tr 2) for headache on separate occasions, and the lengths of time required to have relief from headache were recorded. The paired data of relief times (in minutes) are: (8.4, 6.9), (7.7, 6.8) (10.1, 10.3) (9.6, 9.4) (9.3, 8.0), (9.1, 8.8) (9.0, 6.1) (7.7, 7.4) (8.1, 8.0) and (5.3, 5.1). For the plausibility of exponential marginals, Gross and Lam transformed the data by subtracting a threshold value of 5.0 from each data point. Concerned with the problem of testing equality of the two treatments, they considered two models: (i) independent exponentials and (ii) ACBVE of Block and Basu (1974). Under each model, they derived the p.d.f. of the ratio R = X,/X,, calculated the likelihood ratio statistic X based on the n = 10 observations of R, and employed the asymptotic X2 test based on -2 log h. We analyze the same transformed data under the GBVE(8,, e,, f) model. The MLE’s computed by three-dimensional grid search method, are found to be 8, = 3.451, e2 = 2.659, 6 = 0.252, and the simple estimates are 8, = 3.43, & = 2.68, 8, = 0.207, Joe, = 0.263. g*,,,_ = 0.239 and (+*, s”*) = (0.354, 0.251) where n = log(B,/e,). Using the pivotal quantity 6*/6 and the table of its (simulated) percentage points provided by Bain (1978, p. 369) a 95% confidence interval for 6 is found to be (0.166, 0.505). Therefore, the null hypothesis of independence (6 = 1) is rejected. Use the pivotal quantity (?* - n)/6”* and Table 1 of Bain, page 367, yields the 95% confidence interval (0.032, 0.676) for n or equivalently, (1.033, 1.967) for 8,/e,; the null hypothesis H,: 8i = 8, is therefore rejected. Concerning the significance tests, the same conclusions were obtained by Gross and Lam. However, their use of the asymptotic critical values is questionable. Also, the present analysis has the added strength that it provides exact confidence intervals for the parameters of interest. Based on a preliminary simulation study, we conclude that the confidence intervals of 6 and n can also be obtained from the simulated percentage points of the pivotal quantities $r/6 and (?joL -
Remark.
n)/gooL, or &oL/a and (qAOL - n)/%oL.
6. Conclusions Many bivariate exponential models have been proposed in the literature but few are known to combine the advantages of (i) a physical basis, (ii) simplicity of the distribution, (iii) tractability of statistical properties, and (iv) availability of exact methods of hypotheses tests and confidence intervals. Ordinarily, (ii) goes hand-in-hand with (iii) and (iv) but, as the present article shows, the GBVE is an exception. Despite its complicated p.d.f., it has a clear edge over all other existing BVE models in regard to (iii) and (iv) with no shortcomings in (i). 46
Appendix A. Derivation of i,2, ii3 and i,, The following three lemmas list the moment results needed for the derivation of the elements and i,,. The proof of Lemma A.1 is quite similar to that of Lemma 2.2 and is omitted. Lemma
A.1. Let a(6) = l/S
- 1. For the random variable
E(V)=l-tS,
E[V{V+a(G)}-‘1
of it,,
i,,
V defined in (2.1) we have
=6
and E[Vz{V+a(6)}-*] Lemma A.2. E{ V(log V)[V+
=26-1+6P1(1-6)2exp[a(S)]E,[a(6)]. a(6)]-‘}
0
= 6[1 + ‘P(l)].
Proof. We denote E{ V(log V)[(V+ a(6)]-‘} = k, and we define r(j) = E[V,/,(log V){ V,, + a(S)}-‘], j = 0, 1. Since r(2) = E(V,, log V,,) - a(S)r(l), we have k = (1 - S)r(l) + 6r(2) = 6E(VI:, log VI,), and the stated result then follows from
1% v,,)
E(Kl
/,
Lemma A.3. For the random variable
V defined in (2.1) we have
(4
E(log
V) = 6 + P(l),
(b)
E[(log
(c)
E( v log V) = (1 + 26) + 9(1)(1
(4
E[ v(log
v)?]
Proof. Differentiating
I0
= (1 - S)[ \I;(l)
v)*] = (1 - S)[ q,(2)
+ \k2(1)]
+ 4 %(2)
+ \k2(2)],
+ S), + q*(2)]
both sides of the identity
r(t
+ 26[ \k,(3) + q2(3)]. + 1) = JgPx’ exp( -x)
dx with respect
to t, we obtain
mxr(logx)exp(-x)dx=!P(t+l)T(r+l),
~~x’(logx)2exp(-x)dx=[!Pl(t+l)+!P2(t+l)]T(t+l). 0
The stated results follow by substituting t = 0, 1, 2 in the above relations \k(t + 1) = ‘P(t) + l/t and ‘r;(t + 1) = q’(t) - l/t2. 0 TO evaluate 42
Taking
i,, we recall that aZ,/ae, =
and aZ,/a0,
= 0, j = 1, 3, so relation
M~2)-14z2~12.
the derivative 4,
= -tI;‘Z,
relations
(2.7) yields (A.11
of h, given in (2.10) with respect
= -(aQ/aZ,)(aQ/az,)[(L
and using the recursive
-
‘WQ-~
+(a2p/azlaz~)[(l-2/6)~-*+{~+~(6)}-’-i].
to Z,, we obtain
+ {Q + QW’] (A.2)
Using aQ/aZ, = Q1/‘/sZ:/s-l, j = 1, 2, and a2 Q/aZ,aZ, = (1 - 1/6)Q’-2’“(Z,Z2)““~‘, we substitute (A.2) in (A.l) and then referring to the representation (2.1) and its consequences Q = I/, Z, = U6V and Z, = (1 - U)“V, we obtain I9182112= - U(1 - o’)[(l +(1-l/S)U(l-
- 2/6) + P{
v+ o(s)}P2]
U)[(l-2/S)
-I- V{v+a(6)}P’-
v].
Use of Lemma A.1 then leads to the expression of i,, given in (2.8). To evaluate i,,, we note that I,, = -Z,X,,/t$ and h,, = -S-‘Z;’
+ S-‘(aQ/aZ,)[2Q-r
-(aQ/aZ,)(aQ/W[(L
-
+ {Q + a(S)}
ww*
-‘I
+ {Q + 4w2]
+(a2Q/az,as){(i-2p3)~-1+[~+a(s)]m’-i}, where aQ/M = a-‘{
Q log Q - Q’-‘/sB(S)}
a2Q/az,as
and {log Q - I!-’
= K’(Q/Z,)‘-“”
log z, - (1 - 1,,‘6)Q_“%(6)},
and B(6) = Z,“s log z, + z, r/6 log Z,. Rep resentation of Q and its derivatives as functions of U, V and use of Lemmas 2.2, A.1 and A.2 yield the stated expression of i,,. Finally, from (2.6) and (2.7) we obtain I,, = h,, = 2K3(log
Z, + log Z,) - 4K3
-(l-2,‘6){Q-*Q2-Q-IQ} +[Q+a(6)]P’[Q+26-3] where Q = aQ/M
log Q + 46-*Q-‘Q - [Q+a(6)]-2{Q--6-2}2
-Q,
and
Q = i3*Q/&S2 = -S-*Q
log Q + S-‘Q
- S-‘Q’-““B(6)[6P2
log Q + S-IQ + S-2Q’P”6B(S) log Q - { +)/Q}Q]
+ K3Q’-“‘C(6),
and C(6) = Zr’/“(log Z,)’ + Z:/‘(log Z2)2. Representing Z,, j = 1, 2, Q and its derivatives as functions of U and V, and referring to the expectations of functions of U and V given in Lemmas 2.2, A.l-A.3, we obtain i,, as stated in Theorem 2.1.
Appendix
B
n, be i.i.d. GBVE(l, 1, a), and V, as defined in (3.10). Then nl”Vn Lemma B.l. Let (Xl,, X2,), i= l,,.., is asymptotically normal with mean 0 and variance qf, where CT: is given in (3.3). Proof. If (X,, X2) - GBVE(l, 1, S) and T= min( XI, X,) we have Cov( X,, X2) = p and Cov( XI, T) = ulr 2T26 - 2-’ + T2(6 + l)/r(26 + l), so according to the multivariate CLT, the asymptotic distribution of J/*,(X, - l), (X2 - l), (T, - 2-6)] IS normal with mean 0 and covariance matrix Z = (a,,) and a,, = a,, for i, j = 1, 2, 3, u,, = u2* = 1, u33 = 2-*‘, ur2 = p and u,~ = u23 = ulr_ Using this, the stated result follows from the fact that 2,/x, + E(Z,) w.p. 1, and E(Z,)
=E[TZ(X,
=E(T)P(X,
< X2) =2-(‘+l).
Here, we have used the fact that T and Z( XI < X,) are independent. 48
0
Lemma B.2. Let W,(y) be as defined in (3.8). Then, n”2W,,(X) Proof. For arbitrary
0 < 7 < co, consider
N,(T) =
a sequence
3 0.
of neighborhoods
of 1 as
(y: Il(Y,‘>VT’) - (1~1)l(--“*J~
where I/. 11denotes the Euclidean norm (the choice of the reciprocal notational simplicity). For an arbitrary fixed E’ > 0, we define P,,(T)
Considering
= Pr
[
sup ln1’2W,(y)( VEUT)
the event [x E N,(T)]
(l/y,,
l/y,)
is only
for
1, Pin=
Pr[XeN,(r)].
and its complement,
> E’]
Pr[l n”%QX))
> E’
scale
we then have the inequality
+ P*,(T).
(B.1)
We are to show that for any E > 0, the right-hand side of (B.l) is smaller than E for all sufficiently large n. Since n’/‘(X - 1) is bounded in probability, we can choose 7 such that P2,,( r) < SE for all sufficiently large n. Henceforth 7 is fixed at this chosen value, and we proceed to show that pi,,(~) -+ 0 as n + cc. In order to bound n1/2W,(y) in the set y E N,(r) we refer to the definition of W(y) in (3.6). In the set S,,S’ we have (n/n)x,
GX*
T= X2,
G (Yl/Y2)X2?
which implies T(l/y, - l/y*) < X,/y1 - X,/y, < 0. Similarly, in the set S;S, we have T = X, and T(l/y, - l/y,) < X2/y2 - X,/y, d 0. We note that both these sets are empty if yi = y2. If y, # y2 one of the sets is empty. Since 1 - ~n-‘/~ < l/y, < 1 + rn-‘I*, i = 1, 2, we have 11/y, - l/y, 1 < 2rn-1/2, so IW(Y)I~2~n Hence,
-1/2T[Z(S,‘)
+Z(S,‘S)]
for YEN,(T)
for each i = 1,. . . , n, we have
(B.2) We now note that S,,S:: = 1x2,
< Xl, * X1,/y,
< X2,/y,]
c
[I
< &/X2,
< 1 + T’npl”],
for sufficiently large n and T’ = 37, and likewise, S&?, c [l - T’n-“‘* -c X1,/X2, < 11. Denoting A, = [l -r’np’/2, 1 + *‘n-1’2], and referring to (B.2) we then have,_for y E N,(T), 1w(y) I < T’n-“*T,Z( X1,/X2, E A,). This implies sup, E N?ol r?*W,, ( y) 1Q T ‘v,, where V, is the average of i.i.d. y = q.Z( Xu/X2, E n. Now, considering V = TZ( X,/X, E A,) and denoting p,, = Pr( X,/X, E A,), we have by A,), i= l,..., the Cauchy_Schwarz inequality, E(F,) = E(V) < { E(~)P,}“~. Since E(T*) < cc and p, -+ 0 as n + co, we have E( V,) ---f0 which implies that supyE Nn(7jln”‘W, ( y) 1% 0, and the proof is concluded. 0
Acknowledgement The research for this article AFOSR-87-0256.
was sponsored
by the Air Force
Office
of Scientific
Research
under
Grant
49
References Abramowitz, M. and I.E. Stegun, Handbook of Mathematical Funcrions (National Bureau of Standards, Washington, DC). Bain, L.J. (1978) Statistical Analysis of Reliability and LijeTesting Models (Dekker, New York). Block, H.W. and A.P. Basu (1974), A continuous bivariate exponential extension, J. Amer. Statist. Assoc. 69, 10311037. Freund, J.E. (1961) A bivariate extension of exponential distribution, J. Amer. Statist. Assoc. 56, 971-977. Gross, A.J. and C.F. Lam (1981) Paired observations from a survival distribution, Btometrics 37, 5044511. Gumbel, E.J. (1960) Bivariate exponential distribution, J. Amer. Statist. Assoc. 55, 698-707. Gupta, S’S,, A.S. Qureishi and B.K. Shah (1967) Best linear unbiased estimators of the parameters of the logistic distribution using order statistics, Technometrics 9, 43-56. Hougaard, P. (1986a), Survival models for heterogeneous populations derived from stable distributions, Biometrika 73, 387-396. Hougaard, P. (1986b), A class of multivariate failure time distributions, Biometrika 73, 671-678.
50
Johnson, N.L. and S. Kotz (1972) Distributions in Statrstics: Continuous Multivariate Distributions (Wiley, New York). Lee, L. (1979), Multivariate distributions having Weibull properties, J. Multivariate Anal. 9, 267-277. Lu, J.C. and G.K. Bhattacharyya (1988), Some bivariate extensions of the Weibull distribution, Tech. Rept. No. 821, Univ. of Wisconsin (Madison, WI). Mann, N.R., R.E. Schafer and N.D. Singpurwalla (1974) Methods for Statistical Analysis of Reliability and Lifetime Data (Wiley, New York). Marshall, A.W. and I. Olkin (1967) A multivariate exponential distribution, J. Amer. Statist. Assoc. 62, 30-44. Sarkar, S.K. (1987) A continuous bivariate exponential distribution, J. Amer. Statist. Assoc. 82, 667-675. Self, S.G. and K.Y. Liang (1987) Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions, J. Amer. Statist. Assoc. 82, 605-610.