Inference about reliability parameter with gamma strength and stress

Journal of Statistical Planning and Inference 142 (2012) 848–854

Contents lists available at SciVerse ScienceDirect

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

Inference about reliability parameter with gamma strength and stress Kai Huang, Jie Mi n, Zeyi Wang Department of Mathematics and Statistics, Florida International University, Miami, FL 33199, USA

a r t i c l e in f o

abstract

Article history: Received 23 December 2010 Received in revised form 3 October 2011 Accepted 11 October 2011 Available online 20 October 2011

The statistical inference about the reliability parameter R involving independent gamma stress and strength is considered. Assuming the two shape parameters are known arbitrary real numbers, the UMVUE of R is obtained. The performances of the UMVUE and the MLE are compared numerically based on extensive Monte Carlo simulation. Simulation studies indicate that the performance of the two estimators are about the same. The MLE is preferred due to its computational simplicity. & 2011 Elsevier B.V. All rights reserved.

Keywords: Strength Stress Gamma distribution UMVUE MLE MSE Bias

1. Introduction Extensive research has been conducted on the stress–strength model. This model involves two independent random variables X and Y, and the parameter of interest is the probability R ¼ PðX 4YÞ. This model has wide applicability. Naturally, from the engineering point of view X and Y can represent the strength of a structure and the stress imposed on it. With this interpretation the probability R is often called the reliability parameter which will be used in this article. In addition to its applications in engineering it is also applied in many other fields. For example, in a medical application let X represent the response for a control group and Y represent the response for a treatment group (see Hauck et al., 2000; Simonoff et al., 1986). Then in this case the probability PðX 4 YÞ measures the effect of the treatment. This functional can also be used for bioequivalence assessment (see Wellek, 1993). It is frequently used to assess the effectiveness of diagnostic markers in distinguishing between diseased and healthy individuals (Reiser, 2000). In biology, this probability is useful in estimating heritability of a generic trait (Schwarz and Wearden, 1959). An extensive review of this topic is presented in Kotz et al. (2003). For most recent work on the topic it is referred to Baklizi (2008a,b), Eryilmaz (2008a, 2008b, 2010, 2011), Krishnamoorthy et al. (2007, 2009), Krishnamoorthy and Lin (2010), Kundu and Raqab (2009), Rezaei et al. (2010), Saracoglu and Kaya (2007) and the references therein. In the literature various distributions of X and Y are considered which include exponential, normal, Weibull, etc. Constantine et al. (1986) consider the case when both X and Y have gamma distributions. It was assumed there that the two scale parameters are unknown but the two shape parameters are known integers. Under these assumptions they

n

Corresponding author. Tel.: þ 1 305 348 2602; fax: þ1 305 348 6895. E-mail addresses: khuang@fiu.edu (K. Huang), mi@fiu.edu (J. Mi), [email protected] (Z. Wang).

0378-3758/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2011.10.005

K. Huang et al. / Journal of Statistical Planning and Inference 142 (2012) 848–854

849

derived the uniformly minimum variance unbiased estimator (UMVUE), the maximum likelihood estimator (MLE) of R, and obtained an exact confidence interval for R. Under the same assumptions, Constantine et al. (1989, 1990) conducted bootstrap study and compared 10 methods of constructing confidence intervals for R, mainly by using the Mann–Whitney statistic estimator and bootstrap method. Their results are summarized in Kotz et al. (2003). Using normal approximation to gamma distribution, Krishnamoorthy et al. (2008) proposed normal-based approaches for assessing R. In this article we will study the same model but will relax the assumption requiring the two known shape parameters to be integers, i.e., we allow the two shape parameters to be any positive real numbers. The article is organized as follows. The UMVUE of R is derived in Section 2. In Section 3 we present results of numerical studies based on Monte Carlo simulations. It is observed numerically that MLE of R is biased, it is superior to the UMVUE since its MSE is nearly the same as, or somewhat smaller, compared to that of the UMVUE, and additionally its computation is much easier. 2. Point estimation of reliability parameter Let the strength X and stress Y have gamma distributions with density function f X ðxÞ ¼

lk k1 lx x e GðkÞ

and

f Y ðyÞ ¼

tg GðgÞ

yg1 ety

respectively, where both the shape parameters k and g are known but the scale parameters l and t are unknown. Our quantity of interest is the reliability parameter R ¼ PðX 4 YÞ. Suppose X ¼ fX 1 , . . . ,X m g is a sample from the strength population and Y ¼ fY 1 , . . . ,Y n g is a sample from the stress population. It is shown in Constantine et al. (1986) that R can be given by the incomplete beta function or the hypergeometric series. In particular, if at least one of the shape parameters k and g is an integer, then R can be expressed as a finite sum. Here we use different but simpler approach to derive the similar but modified results as follows. Theorem 1. The reliability parameter R is given by 1Fðg=ðkrÞ; 2k,2gÞ, where Fð; 2k,2gÞ denotes the F-distribution with numerator degrees of freedom 2k and denominator degrees of freedom 2g, and is strictly increasing in r ¼ t=l. Moreover, if k is an integer, then it holds that     k1 X Gðj þ gÞ r g 1 j R¼ 1þr Gðj þ 1ÞGðgÞ 1 þ r j¼0 and if g is an integer then R can also be expressed as a finite sum by noting R ¼ 1PðY oXÞ. Proof. By definition it follows that   X 41 : R ¼ PðX 4 YÞ ¼ P Y Note that 2lX  w2 ð2kÞ, 2tY  w2 ð2gÞ and they are independent of each other. Thus         2lX l 2lX=ð2kÞ g1 g g 4 4 R¼P ; 2k,2g , ¼P ¼ P x4 ¼ 1F 2tY 2tY=ð2gÞ t kr kr kr where random variable x  Fð2k,2gÞ. From this it is obvious that R ¼ RðrÞ strictly increases in r 40. Now suppose k is an integer. Note that for any y 40 Z 1 k k1 X l k1 lx ðlyÞj ly x e dx ¼ e : PðX 4 yÞ ¼ GðkÞ j! y j¼0 Therefore, R ¼ PðX 4 YÞ ¼

Z

1

PðX 4 yÞ 0

¼

k1 X

j g

l t Gðj þ gÞ jþg j ¼ 0 j!GðgÞðl þ tÞ

Z

1 0

tg g1 ty y e dy ¼ GðgÞ

Z

1 0

0 @

k1 X ðlyÞj j¼0

j!

1 ely A

tg g1 ty y e dy GðgÞ

    k1 X ðl þ tÞj þ g j þ g1 ðl þ tÞy Gðj þ gÞ r g 1 j y e dy ¼ : j!GðgÞ 1 þ r 1þr Gðj þ gÞ j¼0

&

b can be easily obtained Due to the invariance property of the maximum likelihood estimator the MLE of R, denoted as R, Pm Pn b b b ¼t b=l where l ¼ k=X , t b ¼ g=Y , X ¼ X =m and Y ¼ Y =n. from replacing r by r i¼1

i

j¼1

j

Assuming both shape parameters k and g are known integers, Constantine et al. (1986) obtained the UMVUE of R ~ for any, not expressed in terms of the hypergeometric series. Here we will derive the UMVUE of R, denoted as R, necessarily integers, shape parameters k and g.

850

K. Huang et al. / Journal of Statistical Planning and Inference 142 (2012) 848–854

The following result about conditional density is needed for deriving the UMVUE of R. Tong (1977) gave a general method of obtaining conditional density for exponential families. But it is not easy to apply that method to the gamma case, so we present a short way for this purpose as below. P Lemma 1. The conditional density function of X1, given U  m i ¼ 1 X i ¼ u is GðmkÞ xk1  x ðm1Þk1  f X9U ðx9uÞ ¼ 1 , 0 ox ou: u GðkÞGððm1ÞkÞ uk Pm P Proof. Let U ¼ m i ¼ 1 X i ,W ¼ i ¼ 2 X i . Then U ¼ X 1 þ W and the joint cumulative distribution function of ðX 1 ,UÞ is Gðx,uÞ  PðX 1 rx,U r uÞ ¼ PðX 1 rx,X 1 þ W r uÞ. Note that U and W follow gamma distributions with density functions f U ðuÞ ¼

lmk mk1 lu u e , u40 GðmkÞ

and f W ðwÞ ¼

lðm1Þk Gððm1ÞkÞ

wðm1Þk1 elw ,

w4 0:

Thus Gðx,uÞ ¼

Z

x

PðW r usÞ

0

lk k1 ls s e ds: GðkÞ

From this we have mk

@2 G l ¼ xk1 ðuxÞðm1Þk1 elu , @u@x GðkÞGððm1ÞkÞ

0 ox o u o 1:

Therefore, the conditional density function of X1 given U ¼

Pm

i¼1

X i ¼ u is

mk

f X 1 9U ðx9uÞ ¼ ¼

l GðmkÞ GðmkÞ xk1 ðuxÞðm1Þk1 xk1 ðuxÞðm1Þk1 elu mk  ¼ mk1  l u GðkÞGððm1ÞkÞ G ðkÞ G ððm1ÞkÞ umk1 l u e 1 xk1  x ðm1Þk1  1 , Bðk,ðm1ÞkÞ uk u

0 ox o u o 1,

where Bða,bÞ is the beta function, i.e., Bða,bÞ ¼ GðaÞGðbÞ=Gða þ bÞ.

&

Theorem 2. The UMVUE R~ of R based on samples fX 1 , . . . ,X m g and fY 1 , . . . ,Y n g is given by 8 n m X X R 1 sk1 ð1sÞðm1Þk1 > > > 0 Ius=v ðg,ðn1ÞgÞ ds if v ¼ yj Z xi ¼ u, > > Bðk,ðm1ÞkÞ < j¼1 i¼1 ~ R¼ n m X X > R v=u sk1 ð1sÞðm1Þk1 > > Ius=v ðg,ðn1ÞgÞ ds þ 1I1v=u ðk,ðm1ÞkÞ if v ¼ yj o xi ¼ u, > 0 > : Bðk,ðm1ÞkÞ j¼1 i¼1 where Ix ða,bÞ is the regularized incomplete beta function: Z x a1 t ð1tÞb1 dt: Ix ða,bÞ ¼ Bða,bÞ 0 Proof. Obviously the indicator function ( 1 if X 1 4 Y 1 , IðX 1 4Y 1 Þ ¼ 0 if X 1 r Y 1 P is an unbiased estimator of R ¼ PðX 4 YÞ. Furthermore, ðU,VÞ is the sufficient and complete statistic, where U ¼ m i ¼ 1 Xi Pn and V ¼ j ¼ 1 Y j . Hence, according to the Blackwell–Rao and Lehmann–Scheffe theorems EðIðX 1 4Y 1 Þ9U,VÞ is the unique UMVUE of R. We have R~ ¼ EðIðX 1 4 Y 1 Þ9U ¼ u,V ¼ vÞ ¼ PðX 1 4 Y 1 9U ¼ u,V ¼ vÞ: P According to Lemma 1 the conditional density function of Y1 given V ¼ nj¼ 1 Y j ¼ v is GðngÞ yg1  yðn1Þg1 1 yg1  yðn1Þg1   1 ¼ 1 , 0 oy o v o 1: f Y 1 9V ¼ v Bðg,ðn1ÞgÞ vg v GðgÞGððn1ÞgÞ vg

K. Huang et al. / Journal of Statistical Planning and Inference 142 (2012) 848–854

851

So, denoting x4v ¼ minfx,vg we obtain Z 1Z 1 Iðy ovÞ yg1  yðn1Þg1 Iðx ouÞIðy oxÞ xk1  x ðm1Þk1 1  1 dx dy R~ ¼ g k Bðg,ðn1ÞgÞ v v Bðk,ðm1ÞkÞ u u 0 0   Z Z u yg1  x4v 1 xk1  x ðm1Þk1 1 yðn1Þg1 dy dx: 1  1 ¼ k u Bðg,ðn1ÞgÞ v v v 0 Bðk,ðm1ÞkÞ u 0

ð1Þ

Let us consider two cases: Case 1: u rv. In this case x r u implies x rv and x4v ¼ x. Hence Z x  Z u 1 xk1  x ðm1Þk1 1 yg1  yðn1Þg1 R~ ¼ 1  1 dy dx g k u v 0 Bðk,ðm1ÞkÞ u 0 Bðg,ðn1ÞgÞ v ! Z us=v Z 1 k1 s ð1sÞðm1Þk1 1  t g1 ð1tÞðn1Þg1 dt ds: ¼ Bðg,ðn1ÞgÞ Bðk,ðm1ÞkÞ 0 0

ð2Þ

Case 2: u 4v. Now the integral in (1) can be expressed as the sum of integrals on ð0,vÞ and (v,u). On the interval x 2 ð0,vÞ certainly x4v ¼ x, and on the interval x 2 ðv,uÞ, we have x4v ¼ v. Thus it follows that Z x  Z v  yg1  1 xk1  x ðm1Þk1 1 yðn1Þg1 dy R~ ¼ dx 1  1 k u v v 0 Bðk,ðm1ÞkÞ u 0 Bðg,ðn1ÞgÞ v   Z u Z yg1  v 1 xk1  x ðm1Þk1 1 yðn1Þg1 dy dx þ 1  1 k u v v v Bðk,ðm1ÞkÞ u 0 Bðg,ðn1ÞgÞ v ! Z Z Z v=u k1 su=v g1 1 k1 s ð1sÞðm1Þk1 t ð1tÞðn1Þg1 s ð1sÞðm1Þk1 dt ds þ ds: ð3Þ ¼ Bðk,ðm1ÞkÞ Bð g ,ðn1Þ g Þ Bðk,ðm1ÞkÞ 0 0 v=u From (2) and (3) we obtain the desired result.

&

3. Numerical studies and conclusion We will conduct simulation studies in which the known shape parameters are k ¼ 2:5, g ¼ 4:5. The scale parameters l and t are chosen in such a way that the reliability parameter R ¼ PðX 4YÞ equals a given value based on Theorem 1. Using N ¼10,000 replications, the following figures show the MSEs of the UMVUE R~ and the counterparts of MLE R^ corresponding to various sample sizes m¼n ¼5,10,20,40,80 and 160. Also, as an estimation of bias, the differences R~ R and R^ R between the average R~ , R^ and the true value of R 2 ð0; 1Þ are displayed. ^ while the right one shows the estimated In each of Figs. 1–5, the left plot shows the MSEs, of the estimator R~ and R, ~ ^ ~ ^ biases of R R and R R of R and R. From these figures, it is obvious that the average of the 10,000 values of the UMVUE R~ is ~ As far as the MLE, it is very close to the true value of R. It is certainly consistent with the unbiasedness property of R.

MSE

0.035

UMVUE

UMVUE

0.03

MLE

UMVUE MLE

0.015 0.01

0.025

UMVUE

0.005

MLE

0

0.02

−0.005

0.015

−0.01 −0.015

0.01

MLE

−0.02

0.005 0

Difference

0.02

−0.025 0

0.2

0.4

0.6

0.8

1

−0.03

0

0.2

0.4

0.6

0.8

Fig. 1. A plot of the MSE (left panel) and differences R~ R and R^ R (right panel) for k ¼ 2:5, g ¼ 4:5, N ¼ 100; 000, m ¼ n ¼ 5.

1

852

K. Huang et al. / Journal of Statistical Planning and Inference 142 (2012) 848–854

MSE

0.016

UMVUE

0.014

UMVUE

MLE

MLE 0.005

MLE

0.012

Difference

0.01

UMVUE

UMVUE

0.01

0

0.008 −0.005

0.006 0.004

MLE

−0.01

0.002 0

0

0.2

0.4

0.6

0.8

−0.015

1

0

0.2

0.4

0.6

0.8

1

Fig. 2. A plot of the MSE (left panel) and differences R~ R and R^ R (right panel) for k ¼ 2:5, g ¼ 4:5, N ¼ 100; 000, m ¼ n ¼ 10.

8

x 10−3

MSE

6

UMVUE MLE

7

x 10−3

Difference UMVUE MLE

4

6

2

5

UMVUE

0

4 −2

3

−4

2

MLE −6

1 0

0

0.2

0.4

0.6

0.8

1

−8

0

0.2

0.4

0.6

0.8

1

Fig. 3. A plot of the MSE (left panel) and differences R~ R and R^ R (right panel) for k ¼ 2:5, g ¼ 4:5, N ¼ 100; 000, m ¼ n ¼ 20.

4

x 10−3

MSE

3

UMVUE MLE

3.5

x 10−3

Difference UMVUE MLE

2

3

UMVUE

1

2.5

0

2 −1

1.5

−2

1

−3

0.5 0

0

0.2

0.4

0.6

0.8

1

−4

MLE

0

0.2

0.4

0.6

0.8

Fig. 4. A plot of the MSE (left panel) and differences R~ R and R^ R (right panel) for k ¼ 2:5, g ¼ 4:5, N ¼ 100; 000, m ¼ n ¼ 40.

1

K. Huang et al. / Journal of Statistical Planning and Inference 142 (2012) 848–854

x 10−3

x 10−3

MSE

2

MLE

UMVUE MLE

1.5

1.6

1

1.4

0.5

1.2

UMVUE

0

1

−0.5

0.8

−1

0.6

−1.5

0.4

−2

0.2 0

Difference

2

UMVUE

1.8

853

0

0.2

0.4

0.6

0.8

1

−2.5

MLE 0

0.2

0.4

0.6

0.8

1

Fig. 5. A plot of the MSE (left panel) and differences R~ R and R^ R (right panel) for k ¼ 2:5, g ¼ 4:5, N ¼ 100; 000, m ¼ n ¼ 80.

MSE−−−UMVUE

0.4

1 0.9

0.35

0.8

0.3

0.7

0.25

0.6

0.2

0.5

0.15

0.4 0.3

0.1

0.2

0.05 0

MSE−−−MLE

x 10−3

0.1 0

0.2

0.4

0.6

0.8

1

UMVUE

0.1

0

1

0

0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

MLE

x 10−3

0.5

−0.1 0

−0.2 −0.3

−0.5

−0.4 −1

−0.5 −0.6

0

0.2

0.4

0.6

0.8

1

−1.5

0

0.2

0.4

^ Fig. 6. A plot of the MSE of UMVUE (top, left panel), and MSE of MLE (top, right panel), and differences RR (bottom, left panel) and RR (bottom, right panel) for k ¼ 2:5, g ¼ 4:5, N ¼ 100; 000, m ¼ n ¼ 160.

854

K. Huang et al. / Journal of Statistical Planning and Inference 142 (2012) 848–854

notable that R^ overestimates when R is small and underestimates when R is large, but its bias is always small and actually is only of order 10  3 if sample sizes m ¼ n Z 20. The rate of convergence of the mean squared error is an important characteristics of an estimator. Our extensive numerical results showed that the MSE of the MLE is reduced by a factor of two when the sample size is doubled from m ¼n¼ 5 until m ¼n¼ 5  214 ¼81,920; here we have reported the results only for sample sizes up to m ¼n¼ 160. This indicates that the estimator converge to the true value of R at the rate of Oðn1 Þ, where n is the sample size. Usually, the MSEs of estimators such as UMVUE or MLE for a certain parameter should decrease as the sample sizes ^ However, this pattern of decreasing MSE holds only when the sample increase. This pattern has been observed for MLE R. ~ The reason for that is because we have to evaluate beta function and incomplete sizes are not too large for the UMVUE R. beta function with large parameters ðn1Þg and ðm1Þk according to Theorem 2, these numerical computation will accumulate big errors. Moreover, the integration in the formula of R~ makes the error even bigger because its integrand involves incomplete beta function. This phenomena is displayed in Fig. 6 in which we see that the bias of R~ is negative all ~ the time. This is one of the key factors showing that the MLE R^ is superior to the UMVUE R. From the above example, it is observed that even though R^ is a biased estimator of R its MSE is about the same as that of ~ Considering this R~ and actually is smaller than the MSE of R~ most of the time. This indicates that R^ is more efficient than R. and the fact that the computation of R~ based on Theorem 2 is more complicated than that of R^ based on Theorem 1, the MLE R^ is recommended for estimating the reliability parameter R.

Acknowledgments The authors wish to thank the reviewers and associate editor for their helpful comments and suggestions to improve the manuscript. References Baklizi, A., 2008a. Estimation of PrðX o YÞ using record values in the one and two parameter exponential distributions. Communications in Statistics—Theory and Methods 37, 692–698. Baklizi, A., 2008b. Likelihood and Bayesian estimation of PrðX o YÞ using lower record values from the generalized exponential distribution. Computational Statistics & Data Analysis 52, 3468–3473. Constantine, K., Karson, M., Tse, S., 1986. Estimation of PðY o XÞ in the gamma case. Communication in Statistics—Simulation and Computation 15, 365–388. Constantine, K., Karson, M., Tse, S., 1989. Bootstrapping estimators of PðY o XÞ in the gamma case. Journal of Statistical Computation and Simulation 19, 217–231. Constantine, K., Karson, M., Tse, S., 1990. Confidence interval estimation of PðY o XÞ in the gamma case. Communication in Statistics—Simulation and Computation 19, 225–244. Eryilmaz, S., 2008a. Consecutive k-out-of-n: G system in stress–strength setup. Communication in Statistics—Simulation and Computation 37, 579–589. Eryilmaz, S., 2008b. Multivariate stress–strength reliability model and its evaluation for coherent structures. Journal of Multivariate Analysis 99, 1878–1887. Eryilmaz, S., 2010. On system reliability in stress–strength setup. Statistics and Probability Letters 80, 834–839. Eryilmaz, S., 2011. A new perspective to stress–strength models. Annals of the Institute of Statistical Mathematics 63 (1), 101–115. Hauck, W.W., Hyslop, T., Anderson, S., 2000. Generalized treatment elects for clinical trials. Statistics in Medicine 19, 887–899. Kotz, S., Lumelsdii, Y., Pensky, M., 2003. The Stress–Strength Model and its Generalization. World Scientific, Singapore. Krishnamoorthy, K., Mathew, T., Mukherjee, S., 2008. Normal-based methods for a gamma distribution: prediction and tolerance intervals and stress– strength reliability. Technometrics 50, 69–78. Krishnamoorthy, K., Mukherjee, S., Guo, H., 2007. Inference on reliability in two-parameter exponential stress–strength model. Metrika 65, 261–273. Krishnamoorthy, K., Lin, Y., Xia, Y., 2009. Confidence limits and prediction limits for a Weibull distribution based on the generalized variable approach. Journal of Statistical Planning and Inference 139, 2675–2684. Krishnamoorthy, K., Lin, Y., 2010. Confidence limits for stress–strength reliability involving Weibull models. Journal of Statistical Planning and Inference 140, 1754–17644. Kundu, D., Raqab, M.Z., 2009. Estimation of R ¼ PðY o XÞ for three-parameter Weibull distribution. Statistics & Probability Letters 79, 1839–1846. Reiser, B., 2000. Measuring the effectiveness of diagnostic markers in the presence of measurement error through the use of ROC curves. Statistics in Medicine 19, 2115–2129. Rezaei, S., Tahmasbi, R., Mahmoodi, M., 2010. Estimation of PðY o XÞ for generalized Pareto distribution. Journal of Statistical Planning and Inference 140, 480–494. Saracoglu, B., Kaya, M.F., 2007. Maximum likelihood estimation and confidence intervals of system reliability for Gomperts distribution in stress–strength models. Selcuk Journal of Applied Mathematics 8, 25–36. Schwarz, L., Wearden, S., 1959. A distribution-free asymptotic method of estimating, testing and setting confidence limits to heritability. Biometrics 15, 227–235. Simonoff, J.S., Hochberg, Y., Reiser, B., 1986. Alternative estimation procedures for PrðX o YÞ in categorized data. Biometrics 42, 895–907. Tong, H., 1977. On the estimation of PrðY r XÞ for exponential families. IEEE Transactions on Reliability 26, 54–56. Wellek, S., 1993. Basing the analysis of comparative bioavailability trials on an individualized statistical definition of equivalence. Biometrical Journal 35, 47–55.