Estimation of P[Y<X] for generalized Pareto distribution

Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

Contents lists available at ScienceDirect

Journal of Statistical Planning and Inference journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / j s p i

Estimation of P[Y < X] for generalized Pareto distribution Sadegh Rezaeia , Rasool Tahmasbib, ∗ , Manijeh Mahmoodib a b

Department of Statistics, Shahid Chamran University of Ahvaz, Ahvaz, Iran Department of Statistics, Amirkabir University of Technology (Polytechnic), Tehran, Iran

A R T I C L E

I N F O

Article history: Received 28 December 2008 Received in revised form 23 July 2009 Accepted 24 July 2009 Available online 4 August 2009 Keywords: Bayes estimator Generalized Pareto distribution Maximum likelihood estimator Bootstrap confidence intervals Asymptotic distributions

A B S T R A C T

This paper deals with the estimation of P[Y < X] when X and Y are two independent generalized Pareto distributions with different parameters. The maximum likelihood estimator and its asymptotic distribution are obtained. An asymptotic confidence interval of P[Y < X] is constructed using the asymptotic distribution. Assuming that the common scale parameter is known, MLE, UMVUE, Bayes estimation of R and confidence interval are obtained. The ML estimator of R, asymptotic distribution and Bayes estimation of R in general case is also studied. Monte Carlo simulations are performed to compare the different proposed methods. © 2009 Elsevier B.V. All rights reserved.

1. Introduction In reliability contexts, inferences about R = P[Y < X], where X and Y independent distributions, are a subject of interest. For example, in mechanical reliability of a system, if X is the strength of a component which is subject to stress Y, then R is a measure of system performance. The system fails, if at any time the applied stress is greater than its strength. The problem of estimating of R have been widely used in the statistical literature. The maximum likelihood estimator (MLE) of P[Y < X], when X and Y are normally distributed, has been considered by Downtown (1973), Govidarajulu (1967), Woodward and Kelley (1977) and Owen et al. (1977). Tong (1974, 1977) considered the estimation of P[Y < X] when X and Y are independent exponential random variables. Awad et al. (1981) considered the MLE of R, when X and Y have bivariate exponential distribution. Ahmad et al. (1997) and Surles and Padgett (1998, 2001) considered the estimation of P[Y < X], where X and Y are Burr type X random variable. The gamma case is studied in Constantine and Karson (1986). The theoretical and practical results on the theory and applications of the stress–strength relationships in industrial and economic systems during the last decades are collected and digested in Kotz et al. (2003). The class of life-time distributions (in particular, exponential and gamma) is considered by Nadarajah (2003). Estimation of P[Y < X] from logistic (Nadarajah, 2004a), Laplace (Nadarajah, 2004b), exponential case with common location parameter (Baklizi and El-Masri, 2004), Burr type III (Mokhlis, 2005), beta (Nadarajah, 2005a), gamma (Nadarajah, 2005b), bivariate exponential (Nadarajah and Kotz, 2006) and Weibull (Kundu and Gupta, 2006) distributions are also studied. Kundu and Gupta (2005) considered the estimation of P[Y < X], when X and Y are independent generalized exponential distribution. Recently, inferences on reliability in two-parameter exponential stress–strength model (Krishnamoorthy et al., 2007) and ML estimation of systemreliability for Gompertz distribution (Saraçoglu and Kaya, 2007) are considered. Kakade et al. (2008) studied exponentiated Gumbel case and Abd Elfattah and Marwa (to appear) studied exponential case based on censored samples.

∗ Corresponding author. E-mail addresses: [email protected] (S. Rezaei), [email protected], [email protected] (R. Tahmasbi), [email protected] (M. Mahmoodi). 0378-3758/$ - see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2009.07.024

S. Rezaei et al. / Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

481

In this paper, we focus on estimation of R = P[Y < X], where X and Y follow the generalized Pareto (GP) distribution with different parameters. We obtain the MLE of R and its asymptotic distribution. The asymptotic distribution is used to construct an asymptotic confidence interval. Two bootstrap confidence intervals of R are also proposed. Assuming that the common scale parameter is known, the maximum likelihood estimator of P[Y < X], confidence intervals, UMVUE and Bayes estimation of R are obtained. This paper is organized as follows: In the next Section, the GP distribution is introduced. Estimation of R with common scale parameter is given in Section 3. In this section, the ML estimator of R, asymptotic distribution and bootstrap confidence intervals are presented. Estimation of R if the common scale parameters are known is discussed in Section 4. In this section MLE, UMVUE and Bayes estimation of R are discussed. In Section 5, the estimation of R in general case is studied. The ML estimator of R, asymptotic distribution and Bayes estimation of R are presented in Section 5. The different proposed methods have been compared using Monte Carlo simulations and their results have been reported in Section 6. In Section 7, a numerical example is illustrated and the results of different methods are compared. 2. Generalized Pareto distribution A random variable X is said to have generalized Pareto distribution, if its probability density function (pdf) is given by f(,,) (x) =

1





1+

x−

−(1/ +1) ,



where ,  ∈ R and  ∈ (0, +∞). For convenience, we reparametrized this distribution by defining /  ≡ , 1/  ≡  and  ≡ 0. Therefore, f (x) = (1 + x)−(+1) ;

x > 0.

The cumulative distribution function is defined by F(x) = 1 − (1 + x)− for  > 0 and  > 0. Here  and  are the shape and scale parameters, respectively. It is also well known that this distribution has decreasing failure rate (DFR) property. This distribution is also known as Pareto distribution of the second type or Lomax distribution. Shi et al. (1999) used generalized Pareto distribution to estimate the size of the maximum inclusion in clean steels and application of this distribution to reinsurance is discussed by Kremer (1997). 3. Estimation of R with common scale parameter In this section, we investigate the properties of R, when the common scale parameter , is the same, and the general case is studied in Section 5. 3.1. Maximum likelihood estimator of R To investigate the properties of R, denote by GP(, ) the distribution of reparametrized GP. Let X ∼ GP(, ) and Y ∼ GP(, ), where X and Y are independent random variables. Therefore,  R = P[Y < X] =

0

 =

=

∞

0

∞

P(Y < X|X = x)P(X = x) dx

(1 + x)−(+1) (1 − (1 + x)− ) dx

 . +

(1)

To compute the MLE of R, suppose X1 , X2 , . . . , Xn is a random sample from GP(, ) and Y1 , Y2 , . . . , Ym is a random sample from GP(, ). Therefore, the log-likelihood function of the observed sample is L(, , ) = n ln  + m ln  + (n + m) ln  − ( + 1)

n  i=1

ln(1 + xi ) − ( + 1)

m  j=1

ln(1 + yj ).

(2)

482

S. Rezaei et al. / Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

The MLEs of ,  and  say ˆ , ˆ and ˆ , respectively, can be obtained as the solutions of n *L n  = − ln(1 + xi ), *  i=1

(3)

m *L m  = − ln(1 + yj ), *  j=1

(4)

n m   yj *L n + m xi = − ( + 1) − ( + 1) .  1 + xi 1 +  yj * i=1 j=1

(5)

From (3)–(5), we obtain n , ln(1 + xi ) m ˆ = m ln(1 +  yj ) j=1

ˆ = n

(6)

i=1

(7)

and ˆ can be obtained as the solution of the non-linear equation f () =

n+m



n m n m     yj yj xi xi n m × × − m − − . 1 +  x 1 +  y 1 +  x 1 + yj ln(1 +  x ) ln(1 +  y ) i j i i j i=1 j=1

− n

i=1

j=1

i=1

(8)

j=1

Therefore, ˆ can be obtained as a solution of the non-linear equation of the form g() = , where ⎡ g() = (n + m)⎣ n i=1

⎤−1 n m n m     yj yj n xi xi m ⎦ . × × + m + + 1 + xi 1 + yj 1 + xi 1 +  yj ln(1 + xi ) j=1 ln(1 + yj ) i=1

j=1

i=1

(9)

j=1

Since, ˆ is a fixed point solution of the non-linear equation (9), therefore, it can be obtained using an iterative scheme as follows: g((j) ) = (j+1) ,

(10)

where (j) is the jth iterate of ˆ . The iteration procedure should be stopped when |(j) − (j+1) | is sufficiently small. Once we obtain ˆ , then ˆ and ˆ can be obtained from (6) and (7), respectively. Since ML estimators are invariant, so the MLE of R becomes Rˆ =

ˆ ˆ + ˆ

.

(11)

3.2. Asymptotic distribution In this section, the asymptotic distribution of hˆ = (ˆ , ˆ , ˆ ) and the asymptotic distribution of Rˆ are obtained. Based on the ˆ the asymptotic confidence interval of R is derived. We denote the expected Fisher information asymptotic distribution of R, matrix of h = (, , ) as J(h) = E(I; ), where I = [Iij ]i,j=1,2,3 is the observed information matrix. I.e., ⎡

2

* L ⎢ *2 ⎢ ⎢ 2 ⎢ * L I() = − ⎢ ⎢ ⎢ ** ⎢ 2 ⎣ * L ** It is easy to see that I11 =

n

2

,

I12 = I21 = 0,

2

* L ** 2 * L *2 2 * L **

2

⎤

* L ** ⎥ ⎥ ⎥ 2 * L ⎥ ⎥. ⎥ ** ⎥ ⎥ 2 * L ⎦ *2

(12)

S. Rezaei et al. / Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

I13 = I31 =

n 

xi , 1 + xi

i=1

I22 =

m

2

,

I23 = I32 =

I33 =

m 

yj

j=1

1 + yj

n+m

2

− ( + 1)

,

n 

x2i

i=1

(1 + xi )2

− ( + 1)

m 

y2j

j=1

(1 + yj )2

.

Using the integrals of the form 

∞ 0

xr−1 (1 + x)− dx =  B(r, − r) −r

for 0 < r < , where B(x, y) is beta function, we have J11 = J22 =

n

2 m

2

, ,

J12 = J21 = 0, J13 = J31 = J23 = J32 = J33 =

1

2

n



B(2, ),

m



B(2, ),

[(n + m) − n( + 1)B(3, ) − m( + 1)B(3, )].

Theorem 1. As n → ∞ and m → ∞ and n/m → p then √ √ √ [ n(ˆ − ), m(ˆ − ), n(ˆ − )] → N3 (0, U−1 (, , )), where ⎡

u11 U(, , ) = ⎣ 0 u31

0 u22 u32

⎤ u13 u23 ⎦ u33

and u11 =

1 1 J11 = 2 , n 

u13 = u31 = u22 =

1 1 J22 = 2 , m 

u23 = u32 = u33 =

 1 J13 = B(2, ), n 

√ p  1 J23 = √ × B(2, ), n p 

( + 1) ( + 1) 1 p+1 − B(3, ) − B(3, ). J33 = 2 2 2 n p  p

Proof. The proof follows from the asymptotic normality of MLE. 

483

484

S. Rezaei et al. / Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

Theorem 2. As n → ∞ and m → ∞ so that n/m → p then √ n(Rˆ − R) → N(0, B),

(13)

where B=

1 k( + )

4

√ 2 [ (u22 u33 − u223 ) − 2 pu23 u31 + 2 p(u11 u33 − u213 )]

and k = u11 u22 u33 − u11 u23 u32 − u13 u22 u31 . Proof. It is clear that √ √ Var[ n(Rˆ − R)] = E[ n(Rˆ − R)]2 2

√ √  n(ˆ − ) −  n(ˆ − ) =E ( + )(ˆ + ˆ )  ⎤ ⎡ √ √ n 2 √ ˆ n 2 √ 2 ˆ − )]2 − ˆ − ) m(ˆ − )]  [ m(  −  )] +  [ n(   [ n(  2 ⎥ ⎢m m ⎥. =E⎢ ⎦ ⎣ 2 2 ˆ ( + ) (ˆ + ) Using Theorem 1, the consistency and asymptotic normality of MLE, the proof is complete.



Note that Theorem 2 can be used to construct asymptotic confidence intervals. To compute the confidence interval of R, the variance B needs to be estimated. To estimate it, the empirical Fisher information matrix and the MLE estimates of ,  and  are used, as follows: uˆ 11 =

1

ˆ 2

,

uˆ 13 = uˆ 31 = uˆ 22 =

1 2 ˆ

ˆ B(3, ˆ ), ˆ

,

ˆ 1 uˆ 23 = uˆ 32 = √ × B(2, ˆ ), p ˆ uˆ 33 =

p+1 2 pˆ

−

ˆ (ˆ + 1) 1 B(3, ˆ ) − B(3, ˆ ). 2 ˆ pˆ

3.3. Bootstrap confidence intervals It is clear that the confidence intervals based on the asymptotic results do not perform very well for small sample size. So, two confidence intervals based on the parametric bootstrap methods are proposed: (i) percentile bootstrap method (Efron, 1982) and (ii) bootstrap-t method (Hall, 1988). The algorithms for estimating the confidence intervals of R using both methods are illustrated below. Algorithm of the percentile bootstrap method: Step 1: From the sample {x1 , x2 , . . . , xn } and {y1 , y2 , . . . , ym }, compute ˆ , ˆ and ˆ . Step 2: Use ˆ and ˆ to generate a bootstrap sample {x∗1 , x∗2 , . . . , x∗n } and similarly use ˆ and ˆ to generate a sample {y∗1 , y∗2 , . . . , y∗m }. Based on {x∗1 , x∗2 , . . . , x∗n } and {y∗1 , y∗2 , . . . , y∗m } compute the bootstrap sample estimate of R using (11), say R∗ . Step 3: Repeat step 2, N boot times. Step 4: Let G(x)=P(Rˆ ∗ ⱕ x), be the cumulative distribution of Rˆ ∗ . Define Rˆ boot =G−1 (x) for a given x. The approximate 100(1− )% confidence interval of R is given by      1−

Rˆ boot , Rˆ boot . 2 2 Algorithm of the bootstrap-t method: Step 1: From the sample {x1 , x2 , . . . , xn } and {y1 , y2 , . . . , ym }, compute ˆ , ˆ , ˆ .

S. Rezaei et al. / Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

485

Step 2: Using ˆ and ˆ generate a bootstrap sample {x∗1 , x∗2 , . . . , x∗n } and similarly using ˆ and ˆ generate a sample {y∗1 , y∗2 , . . . , y∗m }. Based on {x∗1 , x∗2 , . . . , x∗n } and {y∗1 , y∗2 , . . . , y∗m } compute the bootstrap sample estimate of R using (11), say R∗ and following statistic √ ˆ∗ ˆ n(R − R) T∗ =  , Var(Rˆ ∗ ) where Var(Rˆ ∗ ) is obtained using the observed or expected Fisher information matrix. Step 3: Repeat step 2, N boot times. Step 4: For the T ∗ values obtained in step 2, determine the upper and lower bounds of the 100(1 − )% confidence interval of R as follows: let H(x) = P(T ∗ ⱕ x) be the cumulative distribution function of T ∗ . For a given x, define  ˆ −1 (x). Rˆ boot-t (x) = Rˆ + n−1/2 Var(R)H ˆ can be computed as same as computing the Var(Rˆ ∗ ). The approximate 100(1 − )% confidence interval of R is Here also, Var(R) given by    

 , Rˆ boot-t 1 − . Rˆ boot-t 2 2 4. Estimation of R if  is known In this section, we consider the estimation of R when  is known. Without loss of generality, we can assume that  = 1. 4.1. MLE of R Let X1 , X2 , . . . , Xn be a random sample from GP(, 1) and Y1 , Y2 . . . . , Ym be a random sample from GP(, 1). Based on Section 2, it is clear that the MLE of R will be Rˆ =

ˆ ˆ + ˆ

,

(14)

where n , ln(1 + xi ) m . ˆ = m ln(1 + yj ) j=1

ˆ = n

(15)

i=1

(16)

Therefore,  m ni=1 ln(1 + xi ) . Rˆ = m  n j=1 ln(1 + yj ) + m ni=1 ln(1 + xi ) It is easy to see that ln(1 + Xi ) follows an exponential distribution with mean −1 . Therefore, 2  2 2 m j=1 ln(1 + Yj ) ∼  (2m). So, Rˆ ∼

1 1+

(17) n i=1

ln(1 + Xi ) ∼ 2 (2n) and

 F 

or 1 − Rˆ R × ∼ F, 1−R Rˆ where the random variable F follows a F2m,2n distribution with 2m and 2n degrees of freedom. So, the probability density function (pdf) of Rˆ is as follows:  1 − x m−1 x × n+m ,  m 1 − x 1+ n x 



fRˆ (x) =

1 m x2 B(m, n) n

m

where 0 < x < 1 and ,  ⱖ 0.

486

S. Rezaei et al. / Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

The 100(1 − )% confidence interval of R can be obtained as

 1 1 , , 1+F × (1/ Rˆ − 1) 1 + F × (1/ Rˆ − 1) (1− /2;2n,2m)

(18)

( /2;2n,2m)

where F( /2;2n,2m) and F(1− /2;2n,2m) are the lower and upper /2th percentile points of a F distribution. 4.2. UMVUE of R Using the results of Tong (1974, 1977), the UMVUE of R is obtained. This method is also used by Kundu and Gupta (2005)   and Kakade et al. (2008). When the common scale parameter  is known ( ni=1 ln(1 + Xi ), m i=1 ln(1 + Yi )) is a jointly sufficient statistics for (, ). Let us define  1 if ln(1 + X1 ) < ln(1 + Y1 ), (X1 , Y1 ) = 0 o.w. ˜ can be Clearly, 1 − (X1 , Y1 ) is an unbiased estimator of R, since E(1 − (X1 , Y1 )) = P(Y < X) = R. Therefore, the UMVUE of R, say R, obtained as ⎛ ⎞⎤ ⎡ n m   ˜R = E ⎣1 − (X1 , Y1 )| ⎝ (19) ln(1 + Xi ), ln(1 + Yi )⎠⎦ . i=1

i=1

Since ln(1 + Xi ) and ln(1 + Yi ) are exponential random variables, so using the results of Tong (1974,1977) it follows that ⎧ i  m−1  (n − 1)!(m − 1)! T1 ⎪ ⎪ ⎪1 − (−1)i = if T1 ⱕ T2 , ⎨ (m − i − 1)!(n + i − 1)! T2 i=0 R˜ =   i n−1 ⎪  (n − 1)!(m − 1)! T2 ⎪ ⎪ ⎩ (−1)i if T2 ⱕ T1 , (m + i − 1)!(n − i − 1)! T 1 i=0   where T1 = ni=1 ln(1 + Xi ) and T2 = m i=1 ln(1 + Yi ).

(20)

4.3. Bayes estimation of R In this section, we obtain the Bayes estimation of R under assumption that the shape parameters  and  are random variables. It is assumed that  and  have independent gamma priors with the pdfs a

() =

b11 a1 −1 −b1   e ; (a1 )

() =

b22 a2 −1 −b2   e ; (a2 )

>0

(21)

 > 0,

(22)

and a

with the parameters  ∼ Gamma(a1 , b1 ) and  ∼ Gamma(a2 , b2 ). The posterior pdfs of  and  are as follows:

|data ∼ Gamma(a1 + n, b1 + T1 ),

(23)

|data ∼ Gamma(a2 + m, b2 + T2 ),

(24)

n

m

where T1 = i=1 ln(1 + xi ) and T2 = i=1 ln(1 + yi ). Since a priori  and  are independent, using (23) and (24), the posterior pdf of R becomes fR (r) = C

ra2 +m−1 (1 − r)a1 +n−1 ((b1 + T1 )(1 − r) + (b2 + T2 )r)(n+m+a1 +a2 )

for 0 < r < 1

(25)

and 0 otherwise, where C=

(n + m + a1 + a2 ) (b + T1 )a1 +n (b2 + T2 )a2 +m . (n + a1 )(m + a2 ) 1

There is no explicit expression for the posterior mean or median. On the other hand the posterior mode can be easily obtained: d rA2 −1 (1 − r)A1 −1 [2r2 (B2 − B1 ) + r(2B1 − 2B2 − A1 B2 − A2 B1 ) + A2 B1 ] fR (r) = , dr [B1 (1 − r) + B2 r](A1 +A2 +3) where B1 = b1 + T1 , B2 = b2 + T2 , A1 = a1 + n − 1 and A2 = a2 + m − 1.

S. Rezaei et al. / Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

487

Note that, for r ∈ (0, 1), (d/dr)fR (r) = 0 has only two roots. Using the fact that limr→0+ (d/dr)fR (r) > 0 and limr→1− (d/dr)fR (r) < 0, it easily follows that the density function fR (r) has a unique mode. The posterior mode can be obtained as the unique root of which lies in between 0 and 1 of the following quadratic: 2r2 (B1 − B2 ) + r(2B2 − 2B1 + A1 B2 + A2 B1 ) − A2 B1 = 0. Now, consider the following loss function:  L(a, b) =

1 if |a − b| > c, 0 if |a − b| ⱕ c.

(26)

It is known that Bayes estimates with respect to the above equation is the midpoint of the `modal interval' of length 2c of the posterior distribution (see Ferguson, 1967). Therefore, the posterior mode is an approximate Bayes estimator of R with respect to the above loss function when the constant c is small. As mentioned before, the Bayes estimate of R under squared error loss cannot be computed analytically. Alternatively, using the approximate method of Lindley (1980) and Ahmad et al. (1997), it can be easily seen that the approximate Bayes estimate of R, say Rˆ Bayes , under squared error loss function is ⎡

⎤

Rˆ Bayes = R˜ ⎣1 +

˜ R˜ 2 (˜ (n + a1 − 1) − ˜ (m + a2 − 2)) ⎦ 2

˜ (n + b1 − 1)(m + a2 − 1)

,

(27)

where R˜ =

˜

˜ + ˜

˜ =

,

n + a1 − 1 b1 + T1

˜ =

and

m + a2 − 1 . b2 + T2

5. Estimation of R in general case Computing the R when the scale parameter is different is considered, in this section. Surles and Padgett (1998, 2001) considered this case, also. In Surles and Padgett (2001) there is no exact expression for R, but they presented a bound for it. 5.1. Maximum likelihood estimator of R Let X ∼ GP(, 1 ) and Y ∼ GP(, 2 ), where X and Y are independent random variables. Therefore,  R= =

∞

0∞ 0

=1−

P(Y < X|X = x)P(X = x) dx

1 (1 + 1 x)−(+1) (1 − (1 + 2 x)− ) dx



∞ 0

=1−



1 (1 + 1 x)−(+1) (1 + 2 x)− dx

1 2

− 

∞

0

(1 + t)−(+1)



2 +t 1

−  dt.

Considering integral of the form, 

∞ 0

x −1 ( + x)− (x + )−Q dx = 

− −Q



 

B( ,  − + Q)2 F1 , ;  + Q; 1 − ,



where > 0,  > − Q and 2 F1 (.) is Gauss' hypergeometric function [ ]. Then, R=1−





2 2 F  + 1, 1;  +  + 1; 1 − .  +  1 2 1 1 

(28)

To compute the MLE of R, suppose X1 , X2 , . . . , Xn is a random sample from GP(, 1 ) and Y1 , Y2 , . . . , Ym is a random sample from GP(, 2 ). Therefore, the log-likelihood function of the observed sample is L(, , 1 , 2 ) = n ln  + n ln 1 − ( + 1)

n  i=1

ln(1 + 1 xi ) + m ln  + m ln 2 − ( + 1)

m  j=1

ln(1 + 2 yj ).

(29)

488

S. Rezaei et al. / Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

The MLEs of , , 1 and 2 say ˆ , ˆ , ˆ1 and ˆ2 , respectively, can be obtained as the solutions of n *L n  = − ln(1 + 1 xi ), *  i=1

(30)

m *L m  = − ln(1 + 2 yj ), *  j=1

(31)

n  *L xi n = − ( + 1) , 1 + 1 xi *1 1 i=1

(32)

m  *L yi m = − ( + 1) .  1 + 2 yi *2 2 i=1

(33)

From the above equations, we obtain n , ln(1 + 1 xi ) i=1 m ˆ = m j=1 ln(1 + 2 yj )

ˆ = n

(34) (35)

and ˆ1 and ˆ2 can be obtained as the solution of the non-linear equations f (1 ) =

n

1

  xi xi n × − 1 + 1 xi 1 + 1 xi i=1 ln(1 + 1 xi )

(36)

  yi yi m × − . 1 +  y 1 + 2 yi ln(1 +  y ) 2 i 2 i i=1

(37)

n

n

i=1

i=1

− n

and f (2 ) =

m

2

m

m

i=1

i=1

− m

By invariance property of the ML estimators, the MLE of R becomes   ˆ ˆ ˆR = 1 − ˆ 2 F ˆ + 1, 1; ˆ + ˆ + 1; 1 − 2 . 2 1 ˆ + ˆ ˆ1 ˆ1

(38)

5.2. Asymptotic distribution The asymptotic distribution of hˆ = (ˆ , ˆ , ˆ1 , ˆ2 ) is obtained and similar to Theorems 1 and 2, the asymptotic distribution of Rˆ could be obtained. We denote the expected Fisher information matrix of h = (, , 1 , 2 ) as J(h) = E(I; ), where I = [Iij ]i,j=1,2,3,4 is 2

the observed information matrix. I.e., Iij = −* L/ *i *j . It is easy to see that I11 =

n

2

,

I12 = 0, I13 =

n 

xi , 1 + 1 xi

i=1

I14 = 0, I22 =

m

2

,

I23 = 0, I24 =

I33 =

m 

yj

j=1

1 + 2 yj

n

21

,

− ( + 1)

n 

x2i

i=1

(1 + 1 xi )2

,

S. Rezaei et al. / Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

489

I34 = 0, I44 =

m

− ( + 1)

22

m 

y2j

j=1

(1 + 2 yj )2

.

Therefore, J11 =

n

2

,

J12 = 0, J13 =

n

1

B(2, ),

J14 = 0, J22 =

m

2

,

J23 = 0, J24 = J33 =

m

2 n

21

B(2, ),

−

n( + 1)

21

B(3, ),

J34 = 0, J44 =

m

22

−

n( + 1)

22

B(3, ).

Based on the above Fisher information matrix, it is possible to present confidence intervals of R based on the percentile bootstrap and bootstrap-t methods. But, since they are similar to those of mentioned in Section 3.3 and for saving space, we omit them. 5.3. Bayes estimation of R In this section, we obtain the Bayes estimation of R under assumption that the shape parameters  and  and scale parameters 1 and 2 are random variables. We mainly obtain the Bayes estimate of R under the squared error loss, and the corresponding credible interval by the Gibbs sampling technique. It is assumed that , , 1 and 2 have independent gamma priors with the parameters  ∼ Gamma(a1 , b1 ),  ∼ Gamma(a2 , b2 ) 1 ∼ Gamma(a3 , b3 ) and 2 ∼ Gamma(a4 , b4 ). Based on the above assumptions, we have the likelihood function of the observed data as n

L(data; , , 1 , 2 ) = n 1

n m   m m (1 + 1 xi )−(+1)  2 (1 + 2 yi )−(+1) . i=1

i=1

Therefore, the joint density of the data, , , 1 and 2 can be obtained as L(data, , , 1 , 2 ) = L(data; , , 1 , 2 ) () () (1 ) (2 ),

(39)

where (.) is the prior distribution. Therefore the joint posterior density of , , 1 and 2 given the data is L(, , 1 , 2 |data) =  ∞  ∞  ∞  ∞ 0

0

0

0

L(data, , , 1 , 2 ) . L(data, , , 1 , 2 ) d d d1 d2

(40)

Since (40) cannot be obtained analytically, we adopt the Gibbs sampling technique to compute the Bayes estimate of R, and the corresponding credible interval of R. The posterior pdfs of  and  are as follows: ⎛ ⎞ n  ⎝ (41) |, 1 , 2 , data ∼ Gamma a1 + n, b1 + ln(1 + 1 xi )⎠ , i=1

⎛

|, 1 , 2 , data ∼ Gamma ⎝a2 + m, b2 +

m  i=1

⎞ ln(1 + 2 yi )⎠

(42)

490

S. Rezaei et al. / Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

and n+a3 −1 −b3 1

f1 (1 |, , 2 , data) ∝ 1

e

n 

(1 + 1 xi )−(+1) ,

(43)

i=1 m+a4 −1 −b4 2

f2 (2 |, , 1 , data) ∝ 2

e

m  (1 + 2 yi )−(+1) .

(44)

i=1

The posterior pdfs of 1 and 2 are not known, but the plots of them show that they are similar to normal distribution. So to generate random numbers from these distributions, we use the Metropolis–Hastings method with normal proposal distribution. Therefore the algorithm of Gibbs sampling is as follows: (0) (0) (0) Step 1: Start with an initial guess ((0) ,  , 1 , 2 ). Step 2: Set t = 1.  (t−1) Step 3: Generate (t) from Gamma(a1 + n, b1 + ni=1 ln(1 + 1 xi )). m (t) (t−1) Step 4: Generate  from Gamma(a2 + m, b2 + i=1 ln(1 + 2 yi )). (t)

(t−1)

Step 5: Using Metropolis–Hastings, generate 1 from f1 with the N(1 , 1) proposal distribution. (t) (t−1) Step 6: Using Metropolis–Hastings, generate 2 from f2 with the N(2 , 1) proposal distribution. Step 7: Compute R(t) from (38). Step 8: Set t = t + 1. Step 9: Repeat steps 3–8, T times. (t−1) 2 ,  ) proposal distribution as follows: Note that in steps 5 and 6, we use the Metropolis–Hastings algorithm with q ∼ N(

a. b. c. d.

(t−1)

Let x =  . Generate y from the proposal distribution q. Let p(x, y) = min{1, f (y)/f (x)q(x)/q(y)}. Accept y with probability p(x, y) or accept x with probability 1 − p(x, y).

Now the approximate posterior mean, and posterior variance of R become T 1  (t) ˆ R E(R|data) = T

(45)

T 1  (t) ˆ ˆ (R − E(R|data))2 , V(R|data) = T

(46)

t=1

and

t=1

respectively. Based on T, and R values, using the method proposed by Chen and Shao (1999), the approximate highest posterior density (HPD) credible interval of R can be easily constructed. This method is also used by Kundu and Gupta (2006), but for generating random numbers from an unknown distribution, they used the method proposed by Devroye (1984) to generate a sample from a log-concave density function. 6. Simulation results In this section, we present some results based on Monte Carlo simulations to compare the performance of the different methods mainly for small sample sizes. We consider tree cases separately to draw inference on R, namely when (i) the common scale parameter  is unknown, (ii) the common scale parameter  is known and (iii) the scale parameter 1 and 2 are unknown. In the first two cases, we consider the following small sample sizes: m, n = 15, 25 and 50. We take different values for ,  and , also. For the first case, we assumed that the common scale parameter , is unknown. From the sample, we compute the estimate of , using the iterative algorithm (10). We have used the initial estimate to be 1 and the iterative process stops when the difference between the two consecutive iterates are less than 10−6 . Once we estimate , we estimate  and  using (6) and (7), respectively. Finally, we obtain the MLE of R using (11). We report the average biases and mean square error (MSE) of Rˆ over 1000 replications. We also compute the 95% confidence intervals based on the asymptotic distribution of Rˆ through the estimation of B using two different methods: using the expected information matrix (Theorem 2) and replacing the parameters by their estimates (denoted by CIe ), and estimating B from the observed information matrix (denoted by CIo ). We report the average confidence intervals and the coverage percentages, cp, based on 1000 replications. All the results are reported in Table 1. Some of the points are clear from this experiment. Even for small sample sizes, the performance of the MLEs are quite ˆ It is observed that when m = n and m, n increase, then MSE(R) ˆ and biases decrease. It satisfactory in terms of biases and MSE(R).

S. Rezaei et al. / Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

491

Table 1 Simulation results and estimation of the parameters when  is unknown from 1000 samples. (n, m)

(, , )

(ˆ , ˆ , ˆ )

Bias(R)

ˆ MSE(R)

CIe

CIo

cp

(15, 15)

(1,1,1) (1.5,2,1) (2,1.5,0.5) (2,1.5,1.5) (0.5,0.5,1) (0.5,1.5,3)

(1.18,1.13,0.95) (1.95,2.75,0.88) (2.70,1.98,0.43) (2.60,1.92,1.31) (0.55,0.54,1.02) (0.57,1.82,2.82)

0.0001 0.0016 0.0008 0.0025 0.0016 0.0082

0.0091 0.0092 0.0094 0.0089 0.0096 0.0064

(0.321,0.679) (0.397,0.749) (0.253,0.606) (0.250,0.603) (0.320,0.677) (0.612,0.904)

(0.321,0.679) (0.399,0.747) (0.255,0.604) (0.251,0.601) (0.320,0.677) (0.634,0.882)

0.941 0.933 0.930 0.934 0.934 0.880

(25, 25)

(1,1,1) (1.5,2,1) (2,1.5,0.5) (2,1.5,1.5) (0.5,0.5,1) (0.5,1.5,3)

(1.06,1.08,0.97) (1.74,2.33,0.94) (2.40,1.81,0.46) (2.37,1.75,1.41) (0.52,0.52,1.00) (0.53,1.68,2.98)

0.0010 0.0006 0.0000 0.0029 0.0013 0.0045

0.0056 0.0052 0.0052 0.0056 0.0051 0.0036

(0.362,0.639) (0.435,0.709) (0.292,0.565) (0.289,0.562) (0.360,0.637) (0.641,0.868)

(0.362,0.639) (0.436,0.707) (0.293,0.564) (0.290,0.561) (0.360,0.637) (0.710,0.799)

0.937 0.937 0.938 0.931 0.951 0.890

(25, 50)

(1,1,1) (1.5,2,1) (2,1.5,0.5) (2,1.5,1.5) (0.5,0.5,1) (0.5,1.5,3)

(1.06,1.07,0.97) (1.70,2.23,0.94) (2.28,1.65,0.47) (2.17,1.60,1.45) (0.51,0.51,0.97) (0.53,1.58,2.92)

0.0001 0.0020 0.0026 0.0057 0.0021 0.0007

0.0040 0.0037 0.0039 0.0037 0.0038 0.0028

(0.380,0.620) (0.451,0.688) (0.308,0.544) (0.305,0.541) (0.378,0.618) (0.651,0.850)

(0.381,0.618) (0.448,0.691) (0.299,0.554) (0.295,0.551) (0.434,0.562) (0.678,0.824)

0.942 0.952 0.955 0.960 0.968 0.916

(50, 50)

(1,1,1) (1.5,2,1) (2,1.5,0.5) (2,1.5,1.5) (0.5,0.5,1) (0.5,1.5,3)

(1.05,1.04,0.97) (1.62,2.18,0.93) (2.12,1.58,0.48) (2.16,1.63,1.43) (0.51,0.51,1.00) (0.52,1.59,2.97)

0.0002 0.0033 0.0042 0.0001 0.0005 0.0044

0.0026 0.0028 0.0026 0.0026 0.0025 0.0019

(0.402,0.598) (0.478,0.672) (0.328,0.521) (0.332,0.526) (0.403,0.599) (0.674,0.835)

(0.402,0.598) (0.480,0.670) (0.329,0.520) (0.333,0.524) (0.403,0.599) (0.659,0.849)

0.941 0.924 0.935 0.938 0.957 0.977

Table 2 Simulation results and estimation of the parameters when  is known from 1000 samples. (n, m)

R

Rˆ

CIMLE

R˜

Rˆ Bayes

cp

(15,15)

0.500 0.571 0.429 0.667 0.333

0.497 0.571 0.431 0.666 0.336

(0.328,0.667) (0.396,0.730) (0.271,0.605) (0.496,0.802) (0.199,0.507)

0.497 0.573 0.429 0.671 0.331

0.497 0.569 0.433 0.661 0.341

0.948 0.952 0.947 0.947 0.961

(25,25)

0.500 0.571 0.429 0.667 0.333

0.499 0.569 0.429 0.665 0.336

(0.365,0.633) (0.432,0.696) (0.302,0.565) (0.534,0.775) (0.226,0.467)

0.499 0.570 0.427 0.668 0.333

0.499 0.567 0.430 0.662 0.339

0.947 0.950 0.939 0.952 0.949

(25,50)

0.500 0.571 0.429 0.667 0.333

0.498 0.566 0.425 0.661 0.331

(0.386,0.620) (0.452,0.682) (0.319,0.549) (0.553,0.763) (0.238,0.448)

0.501 0.570 0.427 0.666 0.331

0.501 0.568 0.429 0.661 0.335

0.950 0.953 0.953 0.948 0.954

verifies the consistency property of the MLE estimators of R. Comparing the average confidence lengths and coverage percentages, it is observed that estimating B either by substituting the parameters by their estimates in Theorem 2 or using observed Fisher information matrix, leads to very similar estimates. The confidence intervals based on the MLEs work quite well unless the sample size is very small, say (15, 15). It is observed that when the sample size is increase, then the coverage percentages of the confidence intervals based on the asymptotic results reach the nominal level, 95%. Considering the case when the common scale parameter is known, we obtain the estimates of R using the MLE and UMVUE methods. The ML estimation of R, Rˆ and its UMVUE estimation, R˜ are reported in Table 2. Since for the MLE, the exact distribution is also known therefore it can be used to construct confidence intervals, also. We denote this confidence interval by CIMLE . To compute different Bayes estimates, we prefer to use the non-informative prior, because we do not have any prior information on R. On the other hand, the non-informative prior, i.e., a1 = a2 = b1 = b2 = 0 provides prior distributions which are not proper, we adopt the suggestion of Congdon (2001, p. 20) and Kundu and Gupta (2005), i.e., choose a1 = a2 = b1 = b2 = 0.0001. The average estimation of Bayes estimates, Rˆ Bayes , is presented in Table 2, under the same prior distributions and based on 1000 replications. In this case as expected, for all the methods when m and n increase, then the average biases decrease. In the third case when the scale parameters 1 and 2 are different and unknown, we obtain the estimates of R using the MLE method. The ML estimation of R, Rˆ are reported in Table 3. To compute different Bayes estimates, we prefer to use the

492

S. Rezaei et al. / Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

Table 3 Simulation results and estimation of the parameters when 1 and 2 are different from 1000 samples. (1 = 1). (n, m)





2

R

Rˆ

ˆ E(R)

ˆ V(R)

(25,25)

1 1 1 2 2 2

1.5 1.5 1.5 1.5 1.5 1.5

0.5 2 3 0.5 2 3

0.7124 0.4787 0.4086 0.5619 0.3032 0.2393

0.7082 0.4787 0.4067 0.5599 0.3000 0.2356

0.7119 0.4784 0.4091 0.5607 0.3034 0.2393

0.0213 0.0251 0.0315 0.0205 0.0272 0.0269

(25,50)

1 1 1 2 2 2

1.5 1.5 1.5 1.5 1.5 1.5

0.5 2 3 0.5 2 3

0.7124 0.4787 0.4086 0.5619 0.3032 0.2393

0.7119 0.4805 0.4076 0.5603 0.3011 0.2363

0.7116 0.4798 0.4096 0.5622 0.3029 0.2398

0.0253 0.0246 0.0267 0.0230 0.0221 0.0225

(50,50)

1 1 1 2 2 2

1.5 1.5 1.5 1.5 1.5 1.5

0.5 2 3 0.5 2 3

0.7124 0.4787 0.4086 0.5619 0.3032 0.2393

0.7129 0.4780 0.4089 0.5614 0.3036 0.2440

0.7122 0.4784 0.4083 0.5613 0.3030 0.2390

0.0181 0.0191 0.0217 0.0180 0.0223 0.0220

Table 4 The data generated with n = m = 20,  = 1,  = 1.4 and  = 1. x

y

x

y

x

y

x

y

0.0024 1.0863 19.3212 6.5356 6.1724

0.0359 0.7108 2.0499 4.7506 2.4195

4.9743 0.0648 1.7140 0.4599 1.8520

5.0290 1.4026 0.2223 4.9126 0.0538

0.1250 0.0185 10.4438 0.3367 0.0088

0.5266 0.7632 0.7383 1.3550 0.0206

0.9618 1.8178 0.7337 0.0296 0.1922

3.3423 0.5894 1.1491 0.6478 0.0824

Table 5 Parameters estimation and bootstrap confidence intervals with N = 500 boot times for the data presented in Table 3.



(ˆ , ˆ , ˆ )

Rˆ

Rˆ ∗

Bˆ e

Bˆ o

CIboot-p

CIboot-t

Unknown Known

(1.13,1.31,1.02) (1.15,1.32,1.00)

0.5346 0.5353

0.5379 0.5361

0.1242 0.1245

0.1243 0.1247

(0.368,0.718) (0.380,0.691)

(0.349,0.686) (0.359,0.702)

non-informative prior, because we do not have any prior information on R. So, we adopt the suggestion of Congdon (2001, p. 20) ˆ and Kundu and Gupta (2005), i.e., choose a1 = a2 = b1 = b2 = 0.0001. The average and variance estimations of Bayes estimates, E(R) ˆ and V(R), are presented in Table 3, under the same prior distributions and based on 1000 replications. In this case as expected, for all the methods when m and n increase, then the average biases decrease. 7. Numerical example Since the confidence intervals based on the asymptotic results for small sample sizes do not perform very well results, so we present an analysis based on two bootstrap methods. To do this, we simulate 20 numbers from GE(1, 1) and 20 numbers from GE(1.4, 1), reported in Table 4. A complete analysis of these data are presented in this section. In the first row of Table 5, it is assumed that  is unknown. The MLE of (, , ) is obtained using (6)–(8). The iterative procedure is stops whenever two consecutive values are less than 10−6 . Using the percentile bootstrap method, we present the mean of 500 bootstrap samples of R by Rˆ ∗ and its 95% confidence interval by CIboot-p . We also present the confidence interval of bootstrap-t method. To compute CIboot-t , the Bˆ e and Bˆ o are computed using the expected and observed information matrix, respectively. In the second row of Table 5, it is assumed that  is known and is equal to 1. The parameters are estimated using (15)–(17). As same as above, the bootstrap confidence intervals are obtained. The performance of the bootstrap confidence intervals are quite well. But construction of bootstrap confidence intervals are computationally more demanding than the asymptotic confidence intervals. Using (20) the UMVUE of R becomes R˜ = 0.5362 and using (27) the Bayes estimation becomes Rˆ Bayes = 0.5344. The profile likelihood function of  is plotted in Fig. 1. It is an upside down function and has a unique maximum. The posterior probability density function of R for the given data set, is plotted in Fig. 2. It is almost symmetric, as expected.

S. Rezaei et al. / Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

493

−86 −88

log L (λ; α, β)

−90 −92 −94 −96 −98 −100 0

0.5

1

1.5

2

2.5

λ Fig. 1. The profile likelihood of  for the given data set presented in Table 3.

x 10−48 2.5

2

f (r)

1.5

1

0.5

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

r Fig. 2. Posterior probability density function of R for the given data set presented in Table 3.

8. Conclusion In this paper, we have addressed the problem of estimating P[Y < X] for the generalized Pareto distributions. We consider the cases when the common scale parameter is known or unknown and when the scale parameters are different. When the common scale parameter is unknown, it is observed that the maximum likelihood estimator works quite well. We can use the asymptotic distribution of the maximum likelihood estimator to construct confidence intervals which work well even for small sample sizes. Based on the simulation results, we recommend to use the parametric Bootstrap percentile method, when the sample size is very small. When the common scale parameter is known we propose maximum likelihood estimator and uniformly minimum variance unbiased estimator. The data analysis indicates that the confidence interval based on the exact distribution of the MLE and the corresponding credible intervals based on the non-informative priors are almost identical. References Abd Elfattah, A.M., Marwa, O.M. Estimating of P(Y < X) in the exponential case based on censored samples. Electronic Journal of Statistics, to appear, ISSN:1935–7524.

494

S. Rezaei et al. / Journal of Statistical Planning and Inference 140 (2010) 480 -- 494

Ahmad, K.E., Fakhry, M.E., Jaheen, Z.F., 1997. Empirical Bayes estimation of P(Y < X) and characterizations of the Burr-type X model. Journal of Statistical Planning and Inference 64, 297–308. Awad, A.M., Azzam, M.M., Hamadan, M.A., 1981. Some inference result in P(Y < X) in the bivariate exponential model. Communications in Statistics—Theory and Methods 10, 2515–2524. Baklizi, A., El-Masri, A.Q., 2004. Shrinkage estimation of P(X < Y) in the exponential case with common location parameter. Metrika 59 (2), 163–171. Chen, M.H., Shao, Q.M., 1999. Monte Carlo estimation of Bayesian credible and HPD intervals. Journal of Computational and Graphical Statistics 8, 69–92. Congdon, P., 2001. Bayesian Statistical Modeling. Wiley, New York. Constantine, K., Karson, M., 1986. The estimation of P(Y < X) in gamma case. Communications in Statistics—Computations and Simulations 15, 365–388. Devroye, L., 1984. A simple algorithm for generating random variates with a log-concave density. Computing 33, 247–257. Downtown, F., 1973. The estimation of P(Y < X) in the normal case. Technometrics 15, 551–558. Efron, B., 1982. The jackknife, the bootstrap and other resampling plans. In: CBMS-NSF Regional Conference Series in Applied Mathematics, vol. 38, SIAM, Phiadelphia, PA. Ferguson, T.S., 1967. Mathematical Statistics: A Decision Theoretic Approach. Academic Press, New York. Govidarajulu, Z., 1967. Tow sided confidence limits for P(X > Y) based on normal samples of X and Y. Sankhya B 29, 35–40. Hall, P., 1988. Theoretical comparison of bootstrap confidence intervals. Annals of Statistics 16, 927–953. Kakade, C.S., Shirke, D.T., Kundu, D., 2008. Inference for P(Y < X) in exponentiated Gumbel distribution. Journal of Statistics and Applications 3 (1–2), 121–133. Kotz, S., Lumelskii, Y., Pensky, M., 2003. The Stress–Strength Model and its Generalizations: Theory and Applications. World Scientific, Singapore. Kremer, E., 1997. A characterization of the generalized Pareto-distribution with an application to reinsurance. Bl¨atter der DGVFM 23 (1), 17–19. Krishnamoorthy, K., Mukherjee, S., Guo, H., 2007. Inference on reliability in two-parameter exponential stress–strength model. Metrika 65 (3), 261–273. Kundu, D., Gupta, R.D., 2005. Estimation of P[Y < X] for generalized exponential distribution. Metrika 61, 291–308. Kundu, D., Gupta, R.D., 2006. Estimation of P[Y < X] for Weibull distributions. IEEE Transactions on Reliability 55 (2), 270–280. Lindley, D.V., 1980. Approximate Bayes method. Trabajos de Estadistica 3, 281–288. Mokhlis, N.A., 2005. Reliability of a stress–strength model with Burr type III distributions. Communications in Statistics—Theory and Methods 34 (7), 1643–1657. Nadarajah, S., 2003. Reliability for lifetime distributions. Mathematical and Computer Modelling 37 (7–8), 683–688. Nadarajah, S., 2004a. Reliability for logistic distributions. Elektronnoe Modelirovanie 26 (3), 65–82. Nadarajah, S., 2004b. Reliability for Laplace distributions. Mathematical Problems in Engineering 2004 (2), 169–183. Nadarajah, S., 2005a. Reliability for some bivariate beta distributions. Mathematical Problems in Engineering (1), 101–111. Nadarajah, S., 2005b. Reliability for some bivariate gamma distributions. Mathematical Problems in Engineering (2), 151–163. Nadarajah, S., Kotz, .S., 2006. Reliability for some bivariate exponential distributions. Mathematical Problems in Engineering 2006, 1–14. Owen, D.B., Craswell, K.J., Hanson, D.L., 1977. Non-parametric upper confidence bounds for P(Y < X) and confidence limits for P(Y < X) when X and Y are normal. Journal of the American Statistical Association 59, 906–924. Saraçoglu, B., Kaya, M.F., 2007. Maximum likelihood estimation and confidence intervals of system reliability for Gompertz distribution in stress–strength models. Selçuk Journal of Applied Mathematics 8 (2), 25–36. Shi, G., Atkinson, H.V., Sellars, C.M., Anderson, C.W., 1999. Application of the generalized Pareto distribution to the estimation of the size of the maximum inclusion in clean steels. Acta Materialia 47 (5), 1455–1468. Surles, J.G., Padgett, W.J., 1998. Inference for P(Y < X) in the Burr type X model. Journal of Applied Statistical Sciences 7, 225–238. Surles, J.G., Padgett, W.J., 2001. Inference for reliability and stress–strength for a scaled Burr-type X distribution. Lifetime Data Analysis 7, 187–200. Tong, H., 1974. A note on the estimation of P(Y < X) in the exponential case. Technometrics 16, 625. Tong, H., 1977. On the estimation of P(Y < X) for exponential families. IEEE Transactions on Reliability 26, 54–56. Woodward, W.A., Kelley, G.D., 1977. Minimum variance unbiased estimation of P(Y < X) in the normal case. Technometrics 19, 95–98.