Inference for the generalized Rayleigh distribution based on progressively censored data

Journal of Statistical Planning and Inference 141 (2011) 3313–3322

Contents lists available at ScienceDirect

Journal of Statistical Planning and Inference journal homepage: www.elsevier.com/locate/jspi

Inference for the generalized Rayleigh distribution based on progressively censored data Mohammad Z. Raqab a, Mohamed T. Madi b, a b

Department of Statistics and Operations Research, King Saud University, Riyadh 11451, Saudi Arabia Department of Statistics, UAE University, Al-Ain, United Arab Emirates

a r t i c l e i n f o

abstract

Article history: Received 17 May 2010 Received in revised form 11 March 2011 Accepted 15 April 2011 Available online 28 April 2011

In this paper, and based on a progressive type-II censored sample from the generalized Rayleigh (GR) distribution, we consider the problem of estimating the model parameters and predicting the unobserved removed data. Maximum likelihood and Bayesian approaches are used to estimate the scale and shape parameters. The Gibbs and Metropolis samplers are used to predict the life lengths of the removed units in multiple stages of the progressively censored sample. Artificial and real data analyses have been performed for illustrative purposes. & 2011 Elsevier B.V. All rights reserved.

Keywords: Generalized Rayleigh distribution Maximum likelihood estimation Bayesian estimation Bayesian prediction Gibbs and Metropolis sampling Importance sampling

1. Introduction Recently, the type-II progressively censoring scheme has received considerable interest in the literature. It is a generalization of type-II right censoring scheme. It can be described as follows. n units are placed on a life-testing experiment and only mð o nÞ are completely observed until failure. The censoring occurs progressively in m stages. These m stages offer failure times of the m completely observed units. At the time of the first failure (the first stage), R1 of the n 1 surviving units are randomly removed (censored) from the experiment. Continuing on, at the time of the second failure (the second stage), R2 of the n2R1 surviving units are randomly withdrawn (censored). Finally, at the time of the mth failure (the mth stage), all the remaining Rm ¼ nmR1     Rm1 surviving units are withdrawn. We will refer to this as progressive type-II right censoring scheme ðR1 ,R2 , . . . ,Rm Þ. It is clear that this scheme includes the conventional type-II right censoring scheme ðR1 ¼ R2 ¼    ¼ Rm1 ¼ 0, Rm ¼ nmÞ and the complete sampling scheme ðR1 ¼ R2 ¼    ¼ Rm1 ¼ 0,n ¼ mÞ. This scheme is considered an important method of obtaining data in lifetime studies. There are many situations where the removal of units prior to failure is preplanned in order to save time and money associated with testing. Cohen (1963, 1966) was one of the earliest to study the more general censoring scheme when removal of the units is also allowed in between. Although several inference techniques for estimating parameters from different lifetime distributions based on progressive type-II censored data have appeared in the literature, the GR model has not been considered yet. It is observed that the two parameter generalized Rayleigh (GR) distribution is always right skewed and can be used quite effectively to analyze any skewed data. Shapes of the different pdfs of the GR distribution can be found in Raqab and Kundu (2006).

 Corresponding author.

E-mail addresses: [email protected], [email protected] (M.Z. Raqab), [email protected] (M.T. Madi). 0378-3758/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2011.04.016

3314

M.Z. Raqab, M.T. Madi / Journal of Statistical Planning and Inference 141 (2011) 3313–3322

The GR distribution with parameters a and l 40 has the cumulative distribution function (cdf) 2

Fðx; a, lÞ ¼ ð1eðlxÞ Þa ,

a, l 4 0,

ð1:1Þ

and probability density function (pdf) 2

2

2

f ðx; a, lÞ ¼ 2al xeðlxÞ ð1eðlxÞ Þa1 ,

a, l 4 0,

ð1:2Þ

respectively. Here a and l 4 0 are the shape and inverse scale parameters, respectively. From now on the generalized Rayleigh distribution with parameters a and l will be denoted by GR(a, lÞ. Several aspects of the GR distribution have been studied by Surles and Padgett (2001), Raqab and Kundu (2006) and Kundu and Raqab (2005, 2007). For some general references on Burr type X distribution, the readers are referred to Sartawi and Abu-Salih (1991), Ahmad et al. (1997), and the references cited there. Let X1 ,X2 , . . . ,Xn denote the failure times of n-independent units placed on a life-testing experiment. Assume these Xi ’s come from a common continuous cumulative distribution function Fðx; yÞ and probability density function f ðx; yÞ. Let X1:m:n rX2:m:n r    rXm:m:n denote the above-mentioned m progressively type-II right censored sample. Recently, several articles, on estimating the unknown parameters for different distributions under this censoring scheme, have appeared. See for example Balakrishnan and Aggrawala (2000), Balakrishnan et al. (2004), Ng et al. (2002), Ng (2005), Mousa and Al-Sagheer (2006), Balakrishnan (2007) and Basak et al. (2009). Kim et al. (2011) studied the two-parameter exponentiated Weibull distribution under type-II progressive censoring. They derived the MLE and Bayes estimators of the two shape parameters, as well as the reliability function. The MLE are derived by solving the likelihood equations via a Newton–Raphson method, while the Bayes estimators are obtained using the Lindley (1980) approximation. In this paper, we estimate the scale and shape parameters involved in the GR model and discuss the prediction of the lifetime lengths of all censored units in all m stages of censoring. In Section 2, we describe the EM algorithm for determining the maximum likelihood estimators of the parameters. Another useful method for determining the MLE is the Newton–Raphson method. The method is illustrated in Section 5. The Bayes estimates for the shape and scale parameters a and l, respectively, are derived in Section 3 using importance sampling. In Section 4, we implement Gibbs and Metropolis sampling to provide sample-based estimates for the predictive density functions of the parameters as well as the times to failure of the Rj ðj ¼ 1,2, . . . ,mÞ units still surviving at the time of observation Xj:m:n ðj ¼ 1,2, . . . ,mÞ. Section 5 presents a numerical simulation and data analysis to illustrate the so-obtained results. 2. Maximum likelihood estimators Suppose n independent items are placed on test. The ordered m failures are observed under the type-II progressively P censoring plan R ¼ ðR1 , . . . ,Rm Þ, where Ri Z 0,n ¼ m þ m j ¼ 1 Rj . Let Y ¼ ðX1:m:n ,X2:m:n , . . . ,Xm:m:n Þ with X1:m:n r X2:m:n r    r Xm:m:n denotes the progressive type-II right censored data from a population with pdf and cdf given in (1.1) and (1.2), respectively. The likelihood function based on a progressive type-II censored sample is given by Lða, l; x1:m:n , . . . ,xm:m:n Þ ¼ c

m Y

f ðxi:m:n ; a, lÞ½1Fðxi:m:n ; a, lÞRi ,

ð2:1Þ

i¼1

where c ¼ nðn1R1 Þðn2R1 R2 Þ    ðnm þ 1R1     Rm1 Þ: From (1.1), (1.2) and (2.1), we write the likelihood function of a and l based on progressive type-II censored sample Y as follows: Lða, ljyÞpam l

2m

2

expf½ml y2m þða1ÞDl ðyÞ þ Ta, l ðyÞg,

ð2:2Þ

where y2m ¼

m m X 2 1X x2 ,D ðyÞ ¼ lnð1eðlxj:m:n Þ Þ m j ¼ 1 j:m:n l j¼1

and Ta, l ðyÞ ¼

m X

2 Rj ln½1ð1eðlxj:m:n Þ Þa :

j¼1

Therefore the log-likelihood function can be written as 2

LogLpmlna þ 2mlnlml y2m ða1ÞDl ðyÞTa, l ðyÞ:

ð2:3Þ

Differentiating (2.3) with respect to a and l, we obtain the likelihood equations: m

a

Dl ðyÞ

2 2 m X Ri ð1eðlxi:m:n Þ Þa lnð1eðlxi:m:n Þ Þ

1ð1eðlxi:m:n Þ Þa 2

i¼1

¼0

ð2:4Þ

M.Z. Raqab, M.T. Madi / Journal of Statistical Planning and Inference 141 (2011) 3313–3322

3315

and ða1Þl

2

ðlxi:m:n Þ2 i:m:n e 2 1eðlxi:m:n Þ

m X x2 i¼1

2

al

2 ðlxi:m:n Þ2 ð1eðlxi:m:n Þ Þa1 i:m:n e 2 a 1ð1eðlxi:m:n Þ Þ

m X Ri x2 i¼1

2

ml y2m þ m ¼ 0:

ð2:5Þ

For complete sample, the last term of Eq. (2.4) has to be cancelled and the MLE of a is obtained as a function of l and Eq. (2.5) can be solved numerically by setting Ri ¼ 0, for all i. In our progressive censored data, the MLEs of a and l have to be obtained by solving the two dimensional equations in (2.4) and (2.5). For simplicity, we propose the EM algorithm, suggested by Dempster et al. (1977) and used by Ng et al. (2002), to compute the MLEs of a and l. Let Z ¼ ðZ1 , . . . ,Zm Þ, with Zj ¼ ðZj1 , . . . ,ZjRj Þ, j ¼ 1,2, . . . ,m, be the censored data. We consider the censored data as missing data. The combination ðY,ZÞ ¼ W forms the complete data set. The log-likelihood function based on W is m X

2

Hðw; a, lÞp2nlnl þnlnal

m X

x2j:m:n þða1Þ

j¼1

2

lnð1eðlxj:m:n Þ Þl

2

j¼1

Rj m X X

z2jk þ ða1Þ

j¼1k¼1

Rj m X X

2

lnð1eðlzjk Þ Þ:

ð2:6Þ

j¼1k¼1

In the E-step, one requires to compute the pseudo-likelihood function. This can be obtained from Hðw; a, lÞ by replacing any function of Zjk (say gðzjk Þ) by EðgðZjk ÞjZjk 4xj:m:n Þ. Therefore the pseudo-log-likelihood function is 2

H ðw; a, lÞp2nlnl þnlnal

m X

x2j:m:n þða1Þ

j¼1

þ ða1Þ

Rj m X X

m X

2

lnð1eðlxj:m:n Þ Þl

j¼1

2

Rj m X X

2 EðZjk jZjk 4 xj:m:n Þ

j¼1k¼1

2

E½lnð1eðlzjk Þ ÞjZjk 4 xj:m:n :

ð2:7Þ

j¼1k¼1

Given Xj:m:n ¼ xj:m:n , the conditional distribution of Zjk follows a truncated GR distribution with left truncation at xj:m:n . That is, fZjY ðzj jYÞ ¼ fW ðzj Þ=½1FW ðxj:m:n Þ, zj 4 xj:m:n . The last two terms in the pseudo-log-likelihood function are evaluated as follows: Z

a

2 E1 ðn; a, lÞ ¼ EðZjk jZjk 4 nÞ ¼

l2 ½1Fðn; a, lÞ

eðlnÞ

2

ðlnuÞð1uÞa1 du

0

and 2

E2 ðn; a, lÞ ¼ E½lnð1eðlZjk Þ ÞjZjk 4 n ¼

    2 2 1 1 1 ð1eðlnÞ Þa lnð1eðlnÞ Þ  : ½1Fðn; a, lÞ a a

Secondly, the M-step involves the maximization of the pseudo-likelihood function by replacing E1 ðn; a, lÞ and E2 ðn; a, lÞ ðkÞ ðk þ 1Þ in (2.7). Assume that at the kth stage, the estimates of ða, lÞ are ðaðkÞ , l Þ, then ðaðk þ 1Þ , l Þ can be obtained by maximizing m X

2

H ðw; a, lÞp2nlnl þ nlnal

x2j:m:n þ ða1Þ

j¼1

m X

2

2

lnð1eðlxj:m:n Þ Þl

j¼1

m X

ðkÞ

Rj E1 ðxj:m:n ; aðkÞ , l Þþ ða1Þ

j¼1

m X

ðkÞ

Rj E2 ðxj:m:n ; aðkÞ , l Þ,

j¼1

ð2:8Þ ðkÞ

with respect to a and l. From (2.8), the estimate of a can be obtained as a function of aðkÞ and l :

a^ M ðlÞ ¼ Pm

j¼1

ðlxj:m:n Þ2

lnð1e

n : P ðkÞ ðkÞ Þþ m j ¼ 1 Rj E2 ðxj:m:n ; a , l Þ

ð2:9Þ

Plugging (2.9) into (2.8), we have H ðw; a, lÞH ðw; a^ M ðlÞ, lÞ ¼ nln



a

a^ M ðlÞ



 n

a a^ M ðlÞ

 1 r 0:

The last inequality holds by the fact that ln x rx1, with equality iff x ¼ 1. Therefore, the maximization of H ðw; a^ M ðlÞ, lÞ can be obtained easily by solving 2

hðlÞ ¼ l ,

ð2:10Þ

where 2

31   m 2 2 m m X X xj:m:n eðlxj:m:n Þ ^ M ðlÞ1 X a 1 1 ðkÞ 2 ðkÞ x  þ R E1 ðxj:m:n ; a , l Þ5 : ðhðlÞ ¼ 4 ðlxj:m:n Þ2 n j ¼ 1 j:m:n nj¼1 j n j ¼ 1 1e ðk þ 1Þ

ðk þ 1Þ

is obtained, aðk þ 1Þ is obtained as aðk þ 1Þ ¼ a^ ðl Þ. Once l The maximum likelihood estimators can also be obtained using the Newton Raphson method, by solving the equations in (2.4) and (2.5) using the second order partial derivatives of the log-likelihood function. Details can be found in Ng et al. (2002). The method is illustrated in Examples 1 and 2 of Section 5.

3316

M.Z. Raqab, M.T. Madi / Journal of Statistical Planning and Inference 141 (2011) 3313–3322

3. Posterior distributions and Bayes estimation In this section we present the posterior densities of the parameters a and l based on progressive censored data and then obtain the corresponding Bayes estimates of these parameters. A natural choice for the priors of a and l (Raqab and Madi, 2009) would be to assume that the two quantities are independent gamma Gða0 ,b0 Þ and generalized exponential power GEPða1 ,b1 Þ distributions, respectively, with the following densities: ba00 a0 1 b0 a a e Gða0 Þ

ga0 ,b0 ðaÞ ¼

and

ha1 ,b1 ðlÞ ¼

2ba11 2a1 1 b1 l2 l e , Gða1 Þ

ð3:1Þ

where a, l 4 0 and a0 ,b0 ,a1 ,b1 are chosen to reflect prior knowledge about a and l. By combining (3.1) and (2.2), we obtain the joint posterior density of a and l

pða, ljyÞpga0 þ m,b0 þ Dl ðyÞ ðaÞha

1

þ m,b1 þ my2m

ðlÞðb0 þ Dl ðyÞÞða0 þ mÞ expfDl ðyÞTa, l ðyÞg:

ð3:2Þ

In light of the progressive sample, we update the prior information about the shape and scale parameters via the posterior model. The marginal posterior density of l is given by

pðljyÞpha

1

þ m,b1 þ my2m

ðlÞ,Ql ðyÞ,

ð3:3Þ

where Ql ðyÞ ¼ ðb0 þDl ðyÞÞða0 þ mÞ expfDl ðyÞgEa ðexpfTa, l ðyÞgÞ and Ea denotes the expectation with respect to Gða0 þm,b0 þ Dl ðyÞÞ: It follows from (3.3) and using importance sampling, that Bayes estimate of l is R1 0 lQl ðyÞha1 þ m,b1 þ my2m ðlÞ dl E ½lQl ðyÞ , ð3:4Þ l^ B ¼ EðljyÞ ¼ R 1 ¼ l El ½Ql ðyÞ Q ðyÞh ð l Þ d l l 2 0 a þ m,b þ my 1

1

m

where El denotes the expectation with respect to GEPða1 þ m,b1 þ my2m Þ: The marginal posterior density of a given l and y is

pðajy, lÞpga0 þ m,b0 þ Dl ðyÞ ðaÞexpfTa, l ðyÞg:

ð3:5Þ

From (3.5), the marginal posterior of a is Z 1 Z 1 pðajyÞ ¼ pða, ljyÞ dl ¼ pðajl,yÞpðljyÞ dl ¼ Eljy ½ga0 þ m,b0 þ Dl ðyÞ ðaÞexpfTa, l ðyÞg 0

¼

0

El ½Ql ðyÞexpfTa, l ðyÞgga0 þ m,b0 þ Dl ðyÞ ðaÞ : El ½Ql ðyÞ

Using the fact that Eðajl,yÞ ¼ Ea ½aexpfTa, l ðyÞg=Ea ½expfTa, l ðyÞg, we obtain the Bayes estimate of a as

a^ B ¼ EðajyÞ ¼ EEðajl,yÞ ¼

El ½Ql ðyÞEðajl,yÞ : El ½Ql ðyÞ

ð3:6Þ

4. Sample-based prediction Prediction of unobserved or censored observations is of interest, especially in actuarial, medical, and engineering sciences. See for example Balakrishnan and Rao (1997), Kaminsky and Nelson (1998), and Basak et al. (2006). In this section we consider the problem of predicting the censored data in the progressively censored sample from the GR distribution. Mainly, we are interested in the posterior density of the jth order statistic from a sample of size Rj removed items at stage j, Zk:Rj ðk ¼ 1,2, . . . ,Rj ; j ¼ 1,2, . . . ,mÞ based on observing the progressively type-II right censored sample, Y ¼ ðx1:m:n ,x2:m:n , . . . ,xm:m:n Þ: The posterior predictive density of Zk:Rj given the observed censored data is Z 1Z 1 fZk:R jY ðzk:Rj , a, lÞpða, ljyÞ da dl, zk:Rj 4xj:m:n : pðzk:Rj jdataÞ ¼ 0

0

j

Here fZk:R jY ðzk:Rj , a, lÞ is the conditional density of Zk:Rj and a, l given Y ¼ y. Using the Markovian property of progressively j type-II right censored order statistics, we have fZk:R jY ðzk:Rj , a, lÞ ¼ fZk:R jXj:m:n ðzk:Rj , a, lÞ ¼ c½Fðzk:Rj ÞFðxj:m:n Þk1 ½1Fðzk:Rj ÞRj k j

j

f ðzk:Rj Þ ½1Fðxj:m:n ÞRj

,

zk:Rj 4 xj:m:n ,

ð4:1Þ

M.Z. Raqab, M.T. Madi / Journal of Statistical Planning and Inference 141 (2011) 3313–3322

3317

where c is defined in (2.1). Therefore the conditional density of Zjk and a, l given the progressively censored sample Y ¼ y from the GR distribution, after algebraic manipulation, can be written as 2

fZk:R jY ðzk:Rj , a, lÞpal zk:Rj expfðlzk:Rj Þ2 g j

j k k1 R X X

ð1Þl þ i



k1

i¼0l¼0

i



rj k



l

 expfaðiZl ðxj:m:n Þ þ ðk þ liÞZl ðzk:Rj ÞÞ

þ Rj ba, l ðzk:Rj Þ þ Zl ðzk:Rj Þg,

ð4:2Þ

where

Zl ðzÞ ¼ lnð1expfðlzÞ2 gÞ and ba, l ðzÞ ¼ ln½1ð1expfðlzÞ2 gÞa : From (3.2), (4.1) and (4.2), the posterior predictive density of Zk:Rj at any point zk:Rj 4xj:m:n is obtained to be

pða, l,zjyÞpam þ a0 l2ðm þ a1 þ 1Þ1 zk:Rj expfl2 ðb1 þ z2k:Rj þ my2m Þ þ Dl ðyÞTa, l ðyÞ þ Zl ðzk:Rj Þ þ Rj ba, l ðzk:Rj Þg j 1 k1 R X X

ð1Þi þ l



k1

i¼0l¼0



Rj k

i



l

 expfaðb0 þ Dl ðyÞ þiZl ðxj:m:n Þ þðk þliÞZl ðzk:Rj ÞÞg:

ð4:3Þ

From (4.3), it is obvious that the form of the predictive density function pðzk:Rj jyÞ will not be tractable and the computation of the predictive Bayes estimate EðZk:Rj jYÞ will be impossible. Consequently, we opt for stochastic simulation procedures, namely, the Gibbs and Metropolis samplers (Gilks et al., 1995) to generate samples from the predictive distributions. For this, let us consider the problem of predicting Z ¼ ðZ1 ,Z2 , . . . ,Zm Þ, where Zj ¼ ðZ1:Rj , . . . ,ZRj :Rj Þ,j ¼ 1,2, . . . ,m, where Zk:Rj is the kth order statistic from a sample of size Rj removed observations at stage j. By forming the product of the extended predictive likelihood and the joint prior of a and l, the full Bayesian model is expressed as

pða, l,ZjyÞpan þ a0 1 l2ðn þ a1 Þ1 expf½ðb1 þ ny2n Þl2 þ aðb0 þ Dl ðy,zÞÞDl ðy,zÞg, where m X

n ¼ mþ

Rj

and

j¼1

0 !1 Rj m X X 1 2 2 y2n ¼ @ x þ z A, n j ¼ 1 j:m:n k ¼ 1 jk

Dl ðy,zÞ ¼ Dl ðyÞ þ Dl ðzÞ,

and Dl ðzÞ ¼

Rj m X X

lnð1expfðlzk:Rj Þ2 gÞ:

j¼1k¼1

Setting

wðkÞ j

¼ ðZ1 , . . . ,Zj1 ,ZðkÞ ,Zj þ 1 , . . . ,Zm Þ, with ZðkÞ ¼ ðZ1:Rj , . . . ,Zk1:Rj ,Zk þ 1:Rj , . . . ,ZRj :Rj Þ, the full conditional distribution j j

of Zk:Rj ð1r k rRj Þ, is found to be 8 2 > al zk:Rj expfðlzk:Rj Þ2 gð1expfðlzk:Rj Þ2 gÞa1 Iðzk1:Rj o zk:Rj o zk þ 1:Rj Þ > > > , > > < ½1expfðlzk þ 1:Rj Þ2 ga ½1expfðlzk1:Rj Þ2 ga ðkÞ pðzk:Rj jy,wj , a, lÞ ¼ 2 2 2 a1 > > > al zk:Rj expfðlzk:Rj Þ gð1expfðlzk:Rj Þ gÞ Iðzk:Rj 4 zk1:Rj Þ > > , > : 1½1expfðlzk1:Rj Þ2 ga

kaRj ,

k ¼ Rj ,

where zj0 ¼ xj:m:n The full conditional distribution of l is

pðljy, aÞpl2ðn þ a1 Þ1 expfðb1 þ ny2n Þl2 gexpfða1ÞDl ðy,zÞg, and the conditional distribution of ajy, l is Gða0 þ n,b0 þ Dl ðy,zÞÞ: 5. Data analysis and discussion To illustrate the above procedures, we present the analysis of two data sets. The first data set is artificial and the second one is a real life one. The computations are performed using Visual Fortran 5.0 and Minitab 15. The graphs are drawn using S-PLUS 6.0. The procedures can be applied easily for any data set. Example 1 (Simulated data). The following progressively type-II right censored data set is generated from a GR distribution with a ¼ 2:0 and l ¼ 1:5: i

1

2

3

4

5

6

7

8

9

10

ri yi

0 0.351

0 0.362

2 0.409

0 0.411

0 0.624

1 0.633

0 0.669

0 0.687

0 0.785

3 0.788

3318

M.Z. Raqab, M.T. Madi / Journal of Statistical Planning and Inference 141 (2011) 3313–3322

Using the EM algorithm outlined in Section 2, the maximum likelihood estimates turned out to be

a^ ¼ 1:969, l^ ¼ 1:584: After seven iterations of the Newton Raphson iterative procedure, the MLEs of a and l are obtained to be

a^ ¼ 2:688, l^ ¼ 1:844 With the same starting values and level of accuracy, it took the EM algorithm 18 iterations to converge. For the Bayesian approach, we assume that a1 ¼ a2 ¼ 1:5 and b1 ¼ b2 ¼ 1:0: The Bayes point estimates of a and l are determined as follows. 1. Ten thousand l values are generated from GEPða1 þm,my 2m þ b1 Þ: 2. For each li , generate 10,000 a values from Gða0 þ m,b0 þ Dli ðyÞÞ. 3. Compute Ea ðexpfTa, l ðyÞgÞ and Ea ðaexpfTa, l ðyÞgÞ by averaging expfTa, l ðyÞg and aexpfTa, l ðyÞÞg with respect to the simulated a values in step 2. 4. Compute Ql ðyÞ and Eðajli ,yÞ: 5. Average the numerators and the denominators of (3.4) and (3.6) with respect to the l values simulated in step 1. The resulting Bayes estimates for a and l are found to be a^ B ¼ 1:966, l^ B ¼ 1:575. After setting the initial values for a, l and Z, a sampler single chain, with a pre-determined number of iterations, is run and used as input in Raftery and Lewis (1992) Fortran program to determine the required number of iterations that should be run to attain convergence. Subsequent to convergence, 1000 draws of equally spaced variates were collected for the parameters a and l as well as for the order statistics Z. Table 1 presents the 100pth (p ¼0.005, 0.025, 0.05,0.5, 0.95, 0.975, 0.995) simulated percentiles of a, l, and the order statistics zk:Rj . We can then establish simulated confidence intervals for these quantities based on different reliability coefficients. For example, 95% confidence intervals for a, l, were, respectively, (0.864,3.905) and (1.020,2.062). From the iterative method, the Bayesian predictive estimates of the order statistics of censored observations and prediction intervals for Zk:Rj are presented in Table 2. For example, the prediction estimates of the last R10 ¼ 3 order statistics still surviving at the observation x10:10:16 are obtained to be z^ 1:R10 ¼ 0:886,

z^ 2:R10 ¼ 1:005,

z^ 3:R10 ¼ 1:199:

Note that the intervals for a and l contain the true values of the parameters used to generate the data. Example 2 (Real data). In this example, we analyze the ball bearing data, which was reported by Lawless (1982, p. 288) and represents the number of million revolutions before failing for each of 23 ball bearings in a life test. The progressively censored sample we use is as follows: i

1

2

3

4

5

6

7

8

9

10

11

12

ri yi

0 17.88

0 28.92

0 33.00

3 41.52

0 51.84

0 51.96

0 54.12

3 55.56

0 68.88

0 84.12

0 93.12

5 98.64

The above data were studied by several authors. Gupta and Kundu (2001) used this data set to compare gamma, Weibull and generalized exponential models. Meintanis (2008) showed that the GR distribution works quite well for these ball bearing data. The maximum likelihood estimates based on the above data were determined to be l^ ¼ 0:011, a^ ¼ 1:100, using the EM algorithm, and l^ ¼ 0:011, a^ ¼ 1:167, using the NR method. The Bayes estimates using the same prior parameter values as in the previous example turned out to be l^ B ¼ 0:012, a^ B ¼ 1:37. Further, 95% prediction intervals for a and l were determined to be, respectively, (0.706,2.169) and (0.008,0.016). The Bayes predictive point estimates and 95% prediction intervals for zk:Rj are presented in Table 3 with R4 ¼ 3, R8 ¼ 3, R12 ¼ 5. To get an idea about the distributional Table 1 Simulated percentiles of the estimated distributions of a, l, and Zk:Rj . p

0.005

0.025

0.05

0.5

0.95

0.975

0.995

a l

0.648 0.839

0.864 1.020

0.960 1.094

1.760 1.550

3.450 1.993

3.905 2.062

4.630 2.172

Z1:R3 Z2:R3

0.411 0.459

0.420 0.524

0.435 0.568

0.613 0.880

1.016 1.487

1.117 1.639

1.315 1.919

Z1:R6

0.634

0.641

0.649

0.848

1.352

1.557

1.771

Z1:R10 Z2:R10 Z3:R10

0.789 0.796 0.847

0.791 0.816 0.875

0.793 0.828 0.902

0.853 0.964 1.155

1.065 1.305 1.707

1.121 1.356 1.823

1.247 1.569 2.102

M.Z. Raqab, M.T. Madi / Journal of Statistical Planning and Inference 141 (2011) 3313–3322

3319

Table 2 Prediction values and intervals for Zk:Rj ,k ¼ 1, . . . ,Rj ,j ¼ 1, . . . ,m. Predicted values

95% Prediction intervals

Z1:R3 Z2:R3

0.659 0.941

(0.422,1.126) (0.511,1.624)

Z1:R6

0.940

(0.647,1.611)

Z1:R10 Z2:R10 Z3:R10

0.886 1.005 1.199

(0.791,1.153) (0.819,1.363) (0.878,1.784)

Table 3 Prediction values and intervals for Zk:Rj ,k ¼ 1, . . . ,Rj ,j ¼ 1, . . . ,m. Predicted values

95% Prediction intervals

Z1:R4 Z2:R4 Z3:R4

63.07 85.50 118.56

(42.12,107.27) (50.40,140.66) (65.93,198.41)

Z1:R8 Z2:R8 Z3:R8

73.37 94.81 124.84

(56.09,114.96) (60.64,155.65) (71.96,209.72)

Z1:R12 Z2:R12 Z3:R12 Z4:R12 Z5:R12

105.78 113.74 123.58 135.79 157.93

(98.85,126.85) (100.42,144.99) (103.22,166.14) (107.60,186.59) (115.22,230.98)

1.0

density

0.8

0.6

0.4

0.2

0.0 0.0

0.5

1.0

1.5

α

2.0

2.5

3.0

Fig. 1. Estimate of the density function of a.

behavior of a, l and Zk:R4 with R4 ¼ 3,k ¼ 1,2, the kernel method for estimating the predictive densities is used and the resulting graphs are displayed in Figs. 1–4, respectively. To see how the MLEs and the Bayes estimators compare, we carried out a Monte Carlo simulation. We compared the Bayes estimators to the MLEs in terms of bias and mean squared error (MSE), for different sample sizes and censoring schemes, including two schemes used in Kim et al. (2011). For a particular n, m and a censoring scheme, we generate a progressively censored sample from the GR distribution with a ¼ 2:0 and l ¼ 1:5: In each case, we compute the MLEs and the Bayes estimators of a and l: We replicate the process 1000 times and compute the average bias and MSE. The results, up to three decimal places, are reported in Tables 4 and 5. The MLE estimates are computed using both the EM algorithm and the NR methods. The bias and MSE values for the NR method are displayed between parentheses.

M.Z. Raqab, M.T. Madi / Journal of Statistical Planning and Inference 141 (2011) 3313–3322

150

density

100

50

0 0.005

0.010

0.015

0.020

λ Fig. 2. Estimate of the density function of l.

0.025

density

0.020

0.015

0.010

0.005

0.000 50

100

150

200

z1:R4 Fig. 3. Estimate of the predictive density function for Z1:R4 .

0.015

density

3320

0.010

0.005

0.000 50

100

150 z2:R4

200

Fig. 4. Estimate of the predictive density function for Z2:R4 .

250

M.Z. Raqab, M.T. Madi / Journal of Statistical Planning and Inference 141 (2011) 3313–3322

3321

Table 4 The average biases and MSEs for the Bayes and MLEs of a. n

m

20 20 20 20 20 25 25 25 30 30 30 30 30

Censoring scheme

10 10 10 10 10 10 10 10 10 10 15 15 15

(10, 9*0) (2,0,2,0,2,0,2,0,2,0) (9*0,10) 10*1 (5*0,5*2) (15,9*0) (5,5,5,7*0) (9*0,15) (0,0,5,5,5,5,4*0) (4*0,10,10,4*0) (0,0,5,5,5,10*0) (14*0,15) (15,14*0)

Bayes

MLE

Bias

MSE

Bias

MSE

 0.048  0.056 0.173 0.222 0.013  0.042  0.048 0.190  0.045  0.049 0.023 0.330 0.003

0.271 0.259 0.924 1.156 0.984 0.260 0.271 0.799 0.230 0.228 0.240 0.795 0.243

 0.451(0.889)  0.419(0.173)  0.096(0.053) 0.867(0.089) 1.033(  0.464)  0.588(0.695)  0.609(0.677)  0.126(0.258)  0.719(  0.384)  0.727(  0.728)  0.476(0.450) 0.158(0.467)  0.495(1.199)

0.284(2.163) 0.266(0.717) 0.864(0.589) 1.252(0.554) 1.542(0.439) 0.376(1.283) 0.397(1.206) 0.716(0.527) 0.527(0.282) 0.538(0.601) 0.274(0.599) 0.741(0.642) 0.295(2.492)

Table 5 The average biases and MSEs for the Bayes and MLEs of l. n

m

20 20 20 20 20 25 25 25 30 30 30 30 30

Censoring scheme

10 10 10 10 10 10 10 10 10 10 15 15 15

(10, 9*0) (2,0,2,0,2,0,2,0,2,0) (9*0,10) 10*1 (5*0,5*2) (15,9*0) (5,5,5,7*0) (9*0,15) (0,0,5,5,5,5,4*0) (4*0,10,10,4*0) (0,0,5,5,5,10*0) (14*0,15) (15,14*0)

Bayes

MLE

Bias

MSE

Bias

MSE

 0.063  0.075  0.063  0.081  0.161  0.060  0.063  0.032  0.072  0.074  0.035 0.071  0.038

0.036 0.037 0.147 0.215 0.198 0.035 0.035 0.123 0.038 0.039 0.028 0.104 0.027

 0.235(  0.161)  0.210(  0.257)  0.192(  0.285) 0.254(  0.266) 0.266(  0.371)  0.324(  0.237)  0.342(  0.286)  0.179(  0.170)  0.409(  0.524)  0.415(  0.770)  0.248(  0.307)  0.020(0.065)  0.250(  0.138)

0.080(0.059) 0.074(0.098) 0.176(0.113) 0.170(0.102) 0.166(0.169) 0.123(0.080) 0.134(0.104) 0.149(0.170) 0.181(0.289) 0.186(0.606) 0.077(0.112) 0.106(0.109) 0.079(0.040)

Table 6 The average biases and MSEs for the Bayes estimator using non-informative prior. n

m

20 20 25 25 30 30 30

10 10 10 10 10 10 15

Censoring scheme

(2,0,2,0,2,0,2,0,2,0) (9*0,10) (5,5,5,7*0) (9*0,15) (0,0,5,5,5,5,4*0) (4*0,10,10,4*0) (14*0,15)

a

l Bias

MSE

Bias

MSE

 0.362  0.303 1.626  0.250  0.387  0.504  0.032

0.368 0.237 0.551 0.172 0.346 0.469 0.101

 0.141  0.120 0.437 0.008 0.043  0.275 0.169

1.349 1.012 1.775 0.825 1.208 1.280 0.758

It can be seen, from these numerical results, that, overall, the Bayes estimators perform better than the MLEs, both in terms of bias and MSE, except for right censoring, where the MLE of a compares quite well with its Bayes estimator. It can also be seen that, overall, our MLE and Bayes estimators for the shape parameter a resulted in better MSE values than those of Kim et al. (2011) for their shape parameter y. We also used the non-informative prior pða, lÞp1=ðalÞ which corresponds to a0 ¼ b0 ¼ a1 ¼ b1 ¼ 0: The partial results in Table 6 show that the corresponding Bayes estimators yielded higher MSE values than the maximum likelihood estimators.

3322

M.Z. Raqab, M.T. Madi / Journal of Statistical Planning and Inference 141 (2011) 3313–3322

6. Conclusions In this article, we have considered the problem of estimation and prediction for the two-parameter GR distribution when the data are progressively type-II censored. The maximum likelihood estimators of the model parameters have been derived. We have also proposed a Bayesian approach to estimating the model parameters and predicting the times to failure of units censored in multiple stages in a progressively censored sample. Importance sampling is used to estimate the shape and scale parameters involved in the model and the Gibbs and Metropolis samplers were simple and useful in predicting the times to failure of the Rj ðj ¼ 1,2, . . . ,mÞ units still surviving at the observation xj:m:n ðj ¼ 1,2, . . . ,mÞ: We have compared the MLEs and Bayes estimates obtained by numerical simulation in terms of the average biases and MSEs for different censoring schemes. It is observed that overall, the Bayes estimators perform better, when compared with the MLEs, except when a non-informative prior is used.

Acknowledgments The authors are thankful to two referees for their insightful comments and suggestions that led to substantive improvements in the article. The authors are grateful to the deanship of scientific research at King Saud University and the Faculty of Business and Economics at UAE University for supporting this research project. References Ahmad, K.E., Fakhry, M.E., Jaheen, Z.F., 1997. Empirical Bayes estimation of PðX o YÞ and characterization of Burr type X model. Journal of Statistical Planning and Inference 64, 297–308. Balakrishnan, N., 2007. Progressive methodology: an appraisal (with discussions). Test 16 (2), 211–259. Balakrishnan, N., Aggrawala, R., 2000. Progressive Censoring, Theory, Methods and Applications. Birkhauser, Boston. Balakrishnan, N., Kannan, N., Lin, C.T., Wu, S.J.S., 2004. Inference for the extreme value distribution under progressive type-II censoring. Journal of Statistical Computation and Simulation 25, 25–45. Balakrishnan, N., Rao, C.R., 1997. Large sample approximations to best linear unbiased estimation and best linear unbiased prediction based on progressively censored samples and some applications. In: Panchapakesan, S., Balakrishnan, N. (Eds.), Advances in Statistical Decision Theory and Applications. Birkhauser, Boston, pp. 431–444. Basak, I., Basak, P., Balakrishnan, N., 2006. On some predictors of times to failure of censored items in progressively censored samples. Computational Statistics and Data Analysis 50, 1313–1337. Basak, P., Basak, I., Balakrishnan, N., 2009. Estimation for the three-parameter lognormal distribution based on progressively censored data. Computational Statistics and Data Analysis 53, 3580–3592. Cohen, A.C., 1963. Progressive censored samples in life testing. Technometrics 5, 327–329. Cohen, A.C., 1966. Life testing and early failure. Technometrics 8, 539–549. Dempster, A.P., Laird, N.M., Rubin, D.D., 1977. Maximum likelihood from incomplete data via EM algorithm. Journal of the Royal Statistical Society, Series B 39, 1–38. Gilks, W.R., Richardson, S., Spiegelhalter, D.J., 1995. Markov Chain Monte Carlo in Practise. Chapman & Hall, London. Gupta, R.D., Kundu, D., 2001. Exponentiated exponential distribution, an alternative to Gamma and Weibull distributions. Biometrical Journal 43 (1), 117–130. Kaminsky, K.S., Nelson, P.I., 1998. Prediction of order statistics. In: Balakrishnan, N., Rao, C.R. (Eds.), Handbook of Statistics, Order Statistics: Applications, vol. 17. North-Holland, Amsterdam, pp. 431–450. Kim, C., Jung, J., Chung, Y., 2011. Bayesian estimation for the exponentiated Weibull model under type II progressive censoring. Statistical Papers 52 (1), 53–70. Kundu, D., Raqab, M.Z., 2005. Generalized Rayleigh distribution; different methods of estimation. Computational Statistics and Data Analysis 49, 187–200 . Kundu, D., Raqab, M.Z., 2007. Discriminating between the log-normal and generalized Rayleigh distribution. Statistics 41 (6), 505–515. Lawless, J.F., 1982. Statistical Models and Methods for Lifetime Data. Wiley, New York. Lindley, D.V., 1980. Approximate Bayesian methods. Trabajos de Stadistica 21, 223–237. Meintanis, S.G., 2008. A new approach of goodness-of-fit testing for exponentiated laws applied to the generalized Rayleigh distribution. Computational Statistics and Data Analysis 52 (5), 2496–2503. Mousa, M.A., Al-Sagheer, S.A., 2006. Statistical inference for the Rayleigh model based on progressively type-II censored data. Statistics 40 (2), 149–157. Ng, H.K.T., 2005. Parameter estimation for a modified Weibull distribution for progressively type-II censored samples. IEEE Transactions on Reliability 54 (3), 374–380. Ng, H.K.T., Chan, C.S., Balakrishnan, N., 2002. Estimation of parameters from progressively ordered censored data using EM algorithm. Computational Statistics and Data Analysis 39, 371–386. Raftery, A.E., Lewis, S., 1992. How many iterations in the Gibbs sampler? In: Bernardo, J.M., Berger, J., Dawid, A.P., Smith, A.F.M. (Eds.), Bayesian Statistics, vol. 4, Oxford, UK, pp. 763–773. Raqab, M.Z., Kundu, D., 2006. Burr type X distribution: revisited. Journal of Probability and Statistical Sciences 4, 179–193. Raqab, M.Z., Madi, M.T., 2009. Bayesian analysis for the exponentiated Rayleigh distribution. Metron, International Journal of Statistics LXVII (3), 269–288. Sartawi, H.A., Abu-Salih, M.S., 1991. Bayes prediction bounds for the Burr type X model. Communications in Statistics—Theory and Methods 20, 2307–2330. 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.