- Email: [email protected]
A Bayesian analysis of record statistics from the Gompertz model
Applied Mathematics and Computation 145 (2003) 307–320 www.elsevier.com/locate/amc
A Bayesian analysis of record statistics from the Gompertz model Zeinhum F. Jaheen
*
Department of Mathematics, University of Assiut, Assiut, 71516, Egypt
Abstract The Gompertz distribution has been used as a growth model and it can be used to fit tumor growth. Record values can be viewed as order statistics from a sample whose size is determined by the values and the order of occurrence of observations. Based on record values from the two-parameter Gompertz distribution, Bayes estimators for the two unknown parameters are obtained by using Laplace approximation. These estimates are obtained based on the squared error and LINEX loss functions. Predictions for future upper record values from the Gompertz model are obtained from a Bayesian approach. The maximum likelihood and Bayes estimates are compared via Monte Carlo simulation study and a numerical example is given to illustrate the results of prediction. Ó 2002 Elsevier Inc. All rights reserved. Keywords: Growth model; Record values; Maximum likelihood; Bayes estimation; Laplace approximation; Squared error loss; LINEX loss; Bayes prediction; Numerical computations
1. Introduction The Gompertz model was formulated by Benjamin Gompertz [10] to fit mortality tables. It has been used as a growth model and also can be used to fit tumor growth. Ahuja and Nash [2] showed that the Gompertz distribution is related by a simple transformation to certain distributions in the family of distributions derived by Pearson. Garg et al. [9] obtained the maximum
* Address: Department of Mathematics, Umm Al-Qura University (TAIF), P.O. Box 3970, Taif, Saudi Arabia. E-mail address: [email protected] (Z.F. Jaheen).
0096-3003/$ - see front matter Ó 2002 Elsevier Inc. All rights reserved. doi:10.1016/S0096-3003(02)00489-7
308
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
likelihood estimates of the parameters of the Gompertz model. Osman [13] used the two-parameter Gompertz distribution to derive a compound Gompertz distribution assuming that one of the parameters is a random variable following the gamma distribution. He studied the properties of his distribution and suggested its use for modeling lifetime data and analyzing the survivals in heterogeneous populations. Applications and more recent survey for the Gompertz model are given by Al-Hussaini et al. [4]. The probability density function (pdf) and cumulative distribution function (cdf) of the Gompertz distribution (denoted by Gompða; bÞ) are given, respectively, by b f ðtÞ ¼ b exp at ½expðatÞ 1 ; t > 0; ða > 0; b > 0Þ; ð1:1Þ a b F ðtÞ ¼ 1 exp ½expðatÞ 1 : ð1:2Þ a The pdf of the Gompertz distribution is unimodal. It has positive skewness and an increasing hazard rate function. The Gompertz distribution can be shown to be a truncated form of the type I extreme value distribution [11]. One disadvantage when using squared error loss, is that it penalises overestimation or underestimation [14]. Overestimation of a parameter can lead to more severe or less severe consequences than underestimation, or vice versa. Subsequently, the use of an asymmetrical loss function, which associates greater importance to overestimation or underestimation, can be considered for the estimation of the parameters. The LINEX loss function is defined as LðDÞ ¼ eaD aD 1;
a 6¼ 0;
ð1:3Þ
where D ¼ /^ðhÞ /ðhÞ, the scalar estimation error, if /ðhÞ is estimated by /^ðhÞ [16]. The sign of a represents the direction and its magnitude represents the degree of symmetry. First, for a ¼ 1 the LINEX loss function is quite asymmetric about zero with overestimation being more costly than underestimation. Second, if a < 0, LðDÞ rises exponentially when D < 0 (underestimation) and almost linearly when D > 0 (overestimation). For a closed to zero, the LINEX is approximately squared error loss and therefore almost symmetric. The posterior expectation of the LINEX loss function in (1.3) is ^ E/ ð/^ðhÞ /ðhÞÞ ¼ ea/ðhÞ E/ ðea/ðhÞ Þ að/^ðhÞ E/ ð/ðhÞÞÞ 1;
ð1:4Þ
where E/ ðÞ denoting posterior expectation with respect to the posterior density of /. The Bayes estimator of /ðhÞ, denoted by /^BL , of the function /ðhÞ under the LINEX loss function is the value /^ðhÞ which minimizes (1.4), it is 1 /^BL ¼ ln E/ ðea/ðhÞ Þ ; ð1:5Þ a provided that E/ ðea/ðhÞ Þ exists and is finite (see [7]).
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
309
Chandler [8] introduced the study of record values and documented many of the basic properties of records. Record values can be viewed as order statistics from a sample whose size is determined by the values and the order of occurrence of observations. In a little over thirty years, a large number of publications devoted to records have appeared. This is possibly due to the fact that we encounter this notion frequently in daily life, especially in singling out record values from a set of others and in registering and recalling record values. Record values can be viewed as order statistics from a sample whose size is determined by the values and the order of occurrence of observations. Let X1 ; X2 ; . . . is a sequence of independent and identically distributed random variables with cdf F ðxÞ and pdf f ðxÞ. Set Yn ¼ maxðminÞX1 ; X2 ; . . . ; Xn , n P 1. We say Xj is an upper (lower) record of this sequence if Yj > ð<ÞYj1 ; j > 1. By definition, X1 is an upper as well as a lower record value. One can transform from upper record values to lower records by replacing the original sequence of random variables by Xj , j P 1 or (if P ðXi > 0 ¼ 1Þ for all i) by 1=Xi ; i P 1; the lower record values of this sequence will correspond to the upper record values of the original sequence. The notations XU ðnÞ and XLðnÞ are used for the nth upper and lower records, respectively. For more details on record values (see [1,5,6]). In this paper, Bayes estimates for the two unknown parameters a and b of the Gompertz distribution are obtained based on upper record values. The estimates are obtained based on the squared error and LINEX loss functions by using the approximation form of Tierney and Kadane [15]. Prediction bounds for future upper record values are obtained from a Bayesian approach. The maximum likelihood and Bayes estimates are compared via Monte Carlo simulation study. Numerical example is given to illustrate the results of prediction.
2. Estimation of the parameters In this section, we shall be concerned with estimation of the two unknown parameters a and b of the Gompertz model based on record values. Suppose we observe n upper record values XU ð1Þ ; XU ð2Þ ; . . . ; XUðnÞ from the Gompertz distribution with pdf given by (1.1). The likelihood function (LF) (see [1]) is given by Lða; bjxÞ ¼
n1 Y
H ðxi Þf ðxn Þ;
ð2:1Þ
i¼1
where x ¼ ðx1 ; x2 ; . . . ; xn Þ and H ðÞ is the hazard function corresponding to the pdf f ðÞ. It follows, from (1.1), (1.2) and (2.1), that
310
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320 n nax
Lða; bjxÞ ¼ b e where x ¼ 1n
Pn
i¼1
exp
b axn ðe 1Þ ; a
ð2:2Þ
xi .
2.1. Maximum likelihood estimation The natural logarithm of the LF (2.2) is b ‘ ln Lða; bjxÞ ¼ n lnðbÞ þ nax ðeaxn 1Þ: a
ð2:3Þ
Assuming that the parameters a and b are unknown, the maximum likelihood estimate (MLE) of the parameter b can be shown to be b^ ¼
n^ a ; 1
ð2:4Þ
ea^xn
where a^ is the MLE of the parameter a which can be obtained as a solution of the following non-linear equation ½^ aðx xn Þ þ 1ea^xn ð^ ax þ 1Þ ¼ 0:
ð2:5Þ
2.2. Bayes estimation Under the assumption that both of the parameters a and b are unknown, we may consider the joint prior density as a product of a conditional density of b for given a (which is taken to be the conjugate gamma prior when a is known) and a two parameter inverted gamma density for a. So that the joint prior density of a and b can be written as gða; bÞ ¼ g1 ðbjaÞg2 ðaÞ;
ð2:6Þ
where g1 ðajbÞ ¼
lm bm1 ebl=a ; CðmÞam
b > 0; ðl > 0; m > 0Þ;
ð2:7Þ
and g2 ðaÞ ¼
cd 1 c=a e ; CðdÞ adþ1
a > 0; ðc > 0; d > 0Þ:
ð2:8Þ
Multiplying g1 ðajbÞ by g2 ðaÞ, we obtain the joint prior density of a and b, given by (2.6), by lm cd bm1 1 ðbl þ cÞ ; a > 0; b > 0: ð2:9Þ gða; bÞ ¼ exp a CðmÞCðdÞ amþdþ1
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
311
The joint posterior density function of a; b given the data, denoted by qða; bjxÞ, can be written as Z Z qða; bjxÞ ¼ Lða; bjxÞgða; bÞ= Lða; bjxÞgða; bÞ da db; a
¼
A amþdþ1
b
nþm1
exp
b
b B1 ða; xn Þ þ B2 ða; xÞ ; a
ð2:10Þ
where B1 ða; xn Þ ¼ eaxn þ l 1; B2 ða; xÞ ¼ nax ac ;
ð2:11Þ
and A is the normalizing constant given by Z 1 ðnþmÞ and1 expfB2 ða; xÞg½B1 ða; xn Þ da: A1 ¼ Cðn þ mÞ
ð2:12Þ
0
2.2.1. Bayes estimators based on squared error loss function Under a squared error loss function, the Bayes estimator /~ of a function /ða; bÞ is given by Z Z /~ E½/ða; bÞjx ¼ /ða; bÞqða; bjxÞ da db ¼
a
Z Z
b
Z Z Lða; bjxÞgða; bÞ da db;
/ða; bÞLða; bjxÞgða; bÞ= a
b
a
ð2:13Þ
b
where Lða; bjxÞ and gða; bÞ are given, respectively, by (2.2) and (2.9). Generally, the ratio of two integrals given by (2.13) can not be obtained in a simple closed form. In this case, we can use numerical integration technique, which can be computationally intensive, especially in high dimensional parameter space. Instead, one can use Laplace approximation form due to Tierney and Kadane [15]. In the following, a review of Tierney and KadaneÕs approximation form is given. 2.2.2. The approximation form of Tierney and Kadane Tierney and Kadane [15] gave an approximate form for the evaluation of the ratio of integrals of the form (2.13) by writing the two expressions, LðkÞ ¼
1 ln qðkjxÞ n
and
So that (2.13) takes the form
L ðkÞ ¼ LðkÞ þ
1 ln /ðkÞ: n
ð2:14Þ
312
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
1=2 det X /^B ¼ exp½nL ðk^ Þ Lðk^Þ; det X 1=2 det X ¼ /ðk^ Þqðk^ jxÞ=qðk^jxÞ; det X
ð2:15Þ
where k^ and k^ maximize L ðkÞ and LðkÞ, respectively, and X and X are the negatives of the inverses of the matrices of second derivatives of L ðkÞ and LðkÞ, at k^ and k^, respectively. Here, we use the approximation form (2.15) to obtain Bayes estimators for the parameters a and b. The regularity condition required for using Tierney–KadaneÕs form (2.15) is that the posterior density should be unimodal. This condition is satisfied for our posterior density given by (2.10), see Appendix A. In our case, k ða; bÞ and the functions Lða; bÞ and L ða; bÞ are obtained from (2.14) and given by
Lða; bÞ ¼ 1n ln A ðm þ d þ 1Þ ln a þ ðn þ m 1Þ ln b ba B1 ða; xn Þ þ B2 ða; xn Þ ; L ða; bÞ ¼ Lða; bÞ þ 1n ln /ða; bÞ:
ð2:16Þ The mode of the posterior density (2.10), denoted by ðaM ; bM Þ, is the solution of the following non-linear equations 9 oL 1 m þ d þ 1 bB1 ða; xn Þ bxn eaxn > ¼ þ þ nxn ; > L1 ða; bÞ ¼ = oa n a a2 a oL 1 n þ m 1 B1 ða; xn Þ > > ; ¼ : L2 ða; bÞ ¼ ob n b a
ð2:17Þ
It follows, From (2.17), that aM ðn þ m 1Þ b^M ¼ ; B1 ðaM ; xn Þ
ð2:18Þ
and aM is a solution of the following non-linear equation hðaÞ ¼
nd2 xn eaxn ðn þ m 1Þ axn þ nxn ¼ 0: a e þl1
ð2:19Þ
It has been theoretically shown, see the Appendix A, that the equation hðaÞ has a unique root and then the posterior density (2.10) is a unimodal. To obtain det X, we first derive the second derivatives of Lða; bÞ with respect to a and b as
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
o2 L 1 m þ d þ 1 2bB1 ða; xn Þ 2beaxn ¼ þ oa2 n a2 a3 a2 o2 L 1 B1 ða; xn Þ xn eaxn ¼ L12 ¼ ; oaob n a2 a o2 L nþm1 ; L22 ¼ 2 ¼ ob nb2
L11 ¼
313
9 b 2 axn > > xn e ; > > > a > > = > > > > > > > ;
ð2:20Þ
where B1 ða; xn Þ is given by (2.11). Therefore, 2
det X ¼ 1=½L11 L22 ðL12 Þ ;
ð2:21Þ
evaluated at the posterior mode ð^ aM ; b^M Þ. Now, the Bayes estimates of a and b are computed as follows: (i) If /ða; bÞ ¼ a, we have from (2.16), a L ða; bÞ
¼ Lða; bÞ þ
1 ln a: n
ð2:22Þ
Thus, a L11
¼ L11
1 ; na2
a L12
¼ L12
and
a L22
¼ L22 :
ð2:23Þ
Hence, 2
det Xa ¼ 1=½a L11 a L22 ða L12 Þ ;
ð2:24Þ
evaluated at the mode of a L ða; bÞ, denoted by solving the two equations L1 ða; bÞ þ
1 ¼0 na
and
ð^ a1 ; b^1 Þ,
L2 ða; bÞ ¼ 0:
which is obtained by
ð2:25Þ
Substituting (2.21) and (2.24) in (2.15), the Bayes estimator of a takes the form !1=2 2 L11 L22 ðL12 Þ qð^ a1 ; b^1 jxÞ=qð^ aM ; b^M jxÞ; ð2:26Þ a^BS ¼ a^1 2 a L11 a L22 ða L12 Þ where qða; bjxÞ is the posterior density given by (2.10) evaluated at the modes ð^ a1 ; b^1 Þ and ð^ aM ; b^M Þ of the functions a L ða; bÞ and Lða; bÞ. (ii) If /ða; bÞ ¼ b, we have from (2.16), b L ða; bÞ
¼ Lða; bÞ þ
1 ln b: n
ð2:27Þ
Thus, b L11
¼ L11
1 ; nb2
b L12
¼ L12
and
b L22
¼ L22 :
ð2:28Þ
314
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
Hence, 2
det Xb ¼ 1=½b L11 b L22 ðb L12 Þ ;
ð2:29Þ
evaluated at the mode of b L ða; bÞ, denoted by ð^ a2 ; b^2 Þ, which is obtained by solving the two equations L1 ða; bÞ þ
1 ¼ 0; nb
L2 ða; bÞ ¼ 0:
ð2:30Þ
Substituting (2.21) and (2.29) in (2.15), the Bayes estimator of b takes the form !1=2 2 L11 L22 ðL12 Þ ^ ^ bBS ¼ b2 qð^ a2 ; b^2 jxÞ=qð^ aM ; b^M jxÞ; ð2:31Þ 2 L L ð L Þ b 11 b 22 b 12 where qða; bjxÞ is the posterior density given by (2.10) evaluated at the modes ð^ a2 ; b^2 Þ and ð^ aM ; b^M Þ of the functions a L ða; bÞ and Lða; bÞ. 2.2.3. Bayes estimators based on LINEX loss function Under the LINEX loss function (1.3), the Bayes estimator /~ of a function /ða; bÞ is given by (1.5), which can be rewritten as 1 /^BL ¼ ln Eðea/ jxÞ a Z Z Z Z 1 a/ e Lða; bjxÞgða; bÞ= Lða; bjxÞgða; bÞ da db : ¼ ln a a b a b ð2:32Þ The ratio of integrals (2.32) can also be computed by using Tierney and KadaneÕs approximate form (2.15). In the following, the Bayes estimators of a and b are computed as follows: (i) If /ða; bÞ ¼ a, we obtain from (2.16), a L ða; bÞ
1 ¼ Lða; bÞ aa: n
ð2:33Þ
Thus, a L11
¼ L11 ;
a L12
¼ L12
and
a L22
¼ L22 :
ð2:34Þ
Hence, 2
det X a ¼ 1=½a L11 a L22 ða L12 Þ ;
ð2:35Þ
^ evaluated at the mode of a L ða; bÞ, denoted by ð^ a 3 ; b3 Þ, which is obtained by solving the two equations a L1 ða; bÞ ¼ 0 and L2 ða; bÞ ¼ 0: ð2:36Þ n
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
315
Substituting (2.21) and (2.35) in (2.32), the Bayes estimator of a takes the form 2 31 , !1=2 2 1 L11 L22 ðL12 Þ ^ qð^ a qð^ aM ; b^M jxÞ5 ; a^BL ¼ ln 4a^ 3 3 ; b3 jxÞ 2 a L L ð L Þ a 12 a 11 a 22 ð2:37Þ where qða; bjxÞ is the posterior density given by (2.10) evaluated at the modes ^ ð^ a aM ; b^M Þ of the functions a L ða; bÞ and Lða; bÞ. 3 ; b3 Þ and ð^ (ii) If /ða; bÞ ¼ b, we have from (2.16), b L ða; bÞ
1 ¼ Lða; bÞ ab: n
ð2:38Þ
Thus, b L11
¼ L11 ;
b L12
¼ L12
and
b L22
¼ L22 :
ð2:39Þ
Hence, 2 det X b ¼ 1=½b L11 b L22 ðb L12 Þ ;
ð2:40Þ
evaluated at the mode of b L ða; bÞ, denoted by ða 4 ; b4 Þ, which is obtained by solving the two equations
L1 ða; bÞ ¼ 0;
L2 ða; bÞ
a ¼ 0: n
ð2:41Þ
Substituting (2.21) and (2.40) in Eq. (2.32), the Bayes estimator of b takes the form 2 31 , !1=2 2 1 L11 L22 ðL12 Þ ^ b^BL ¼ ln 4b^ qð^ a qð^ aM ; b^M jxÞ5 ; 3 4 ; b4 jxÞ 2 a L L ð L Þ b 11 b 22 b 22 ð2:42Þ where qða; bjxÞ is the posterior density given by (2.10) evaluated at the modes ^ ð^ a aM ; b^M Þ of the functions a L ða; bÞ and Lða; bÞ. 4 ; b4 Þ and ð^
3. Prediction of future record values Suppose that we have n upper records XU ð1Þ ¼ x1 ; XU ð2Þ ¼ x2 ; . . . ; XU ðnÞ ¼ xn from the Gomp(a; b) distribution. Based on such a record sample, Bayesian prediction is needed for the sth upper record, 1 < n < s. Let Y XU ðsÞ be the sth upper record value, 1 < n < s. The conditional pdf of Y for given xn is given, (see [1]), by
316
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
½vðyÞ vðxn Þsn1 f ðy; hÞ ; 1 F ðxn ; hÞ Cðs nÞ
f ðyjxn ; hÞ ¼
ð3:1Þ
where vðÞ ¼ lnð1 F ðÞÞ. For the Gomp(a; b) distribution, with pdf given by (1.1), the function f ðyjxn ; hÞ is then given by sn1 beay b b B3 ða; y; xn Þ f ðyjxn ; a; bÞ ¼ exp B3 ða; y; xn Þ ; ð3:2Þ a Cðs nÞ a where B3 ða; y; xn Þ ¼ eay eaxn :
ð3:3Þ
Following Aitchison and Dunsmore [3], the Bayes predictive density function of y given x is given by Z f ðyjxÞ ¼ f ðyjxn ; hÞqðhjxÞ dh: ð3:4Þ H
Substituting (2.10) and (3.2) in (3.4), we get f ðyjxÞ ¼ A
Cðs þ mÞ Iðy; xÞ; Cðs nÞ
ð3:5Þ
where Iðy; xÞ ¼
Z
1
and eay expfB2 ða; xÞg½B3 ða; y; xn Þsn1 ½B1 ða; yÞðsþmÞ da:
0
ð3:6Þ Bayesian prediction bounds for Y ¼ XU ðsÞ , given the previous data x, are obtained by evaluating PrðY P kjxÞ, for some positive k. It follows, from (3.5), that Z 1 PrðY P kjxÞ ¼ f ðyjxÞ dy; k
¼A
Ij ðk; xÞ Cðs þ mÞ sn1 j sn1 Xj¼0 ð1Þ ; j Cðs nÞ nþmþj
ð3:7Þ
where Ij ðk; xÞ ¼
Z 0
1
"
# j fB ða; x Þg 1 n da: and1 expfB2 ða; xÞg fB1 ða; kÞgnþmþj
ð3:8Þ
It can be easily shown that f ðyjxÞ is a density function on the positive half of the real line by proving that PrðY P xn jxÞ ¼ 1.
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
317
A 100s% Bayesian prediction interval for Y XU ðsÞ is such that P ½LLðxÞ < Y < ULðxÞ ¼ s;
ð3:9Þ
where LLðxÞ and ULðxÞ are the lower and upper limits satisfying 9 P ½Y > LLðxÞjx ¼ ð1 þ sÞ=2; = and ; P ½Y > ULðxÞjx ¼ ð1 sÞ=2:
ð3:10Þ
Special case: For the special case, when s ¼ n þ 1, which is practically of special interest, Bayesian prediction bounds for the future upper record value Ynþ1 ¼ XU ðnþ1Þ are obtained from PrðYnþ1 P k1 jxÞ ¼ gðk1 Þ=gðxn Þ where gðk1 Þ ¼
Z
1
ð3:11Þ
ðnþmÞ
and1 expfB2 ða; xÞg½B1 ða; k1 Þ
da;
ð3:12Þ
0
and gðxn Þ is obtained from (3.12) with k1 being replaced by xn . 4. Numerical computations In the following, the maximum likelihood and Bayes (squared error and LINEX) estimates are compared based on a Monte Carlo simulation study and a numerical example is given to illustrate the results of prediction. 4.1. Monte Carlo simulation The estimates obtained in Section 2 are computed and compared according to the following steps: 1. For a given vector of prior parameters ðl; m; c; dÞ, we generate a and b from the joint prior density (2.9). The IMSL [17] is used in the generation of the gamma random variates. 2. For given a; b obtained in step (1), we generate n ¼ ð5; 10; 15Þ upper record values from the Gompða; bÞ with pdf (1.1). 3. The ML estimate of a is computed by solving the nonlinear equation (2.5) by using ZSPOW routine from the IMSL [17] library. Substituting the MLE of a in (2.4) gives the MLE of b. 4. The Bayes estimates of a and b are computed from (2.26), (2.31), (2.37) and (2.42). 5. The squared deviations ðw wÞ2 are computed for different sizes n where () stands for an estimate (ML or Bayes) and w stands for the parameter (a or b).
318
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
Table 1 ER of the estimates of a for different values of l, m, c, d, n, and a (l, m, c, d)
n
ERð^ aML Þ
ERð^ aBS Þ
ERð^ aBL Þ a¼4
a¼8
a ¼ 16
(1.5, 3.0, 1.0, 2.0)
5 10 15
0.1425 0.1138 0.1057
0.1337 0.1124 0.1027
0.1352 0.1131 0.1036
0.1327 0.1119 0.1018
0.1286 0.1108 0.1009
(1.3, 2.5, 2.4, 1.6)
5 10 15
0.0532 0.0460 0.0418
0.0517 0.0399 0.0376
0.0636 0.0415 0.0402
0.0493 0.0365 0.0325
0.0458 0.0330 0.0307
a¼4
a¼8
a ¼ 16
Table 2 ER of the Estimates of b for different values of l, m, c, d, n, and a (l, m, c, d)
n
ERðb^ML Þ
ERðb^BS Þ
ERðb^BL Þ
(1.5, 3.0, 1.0, 2.0)
5 10 15
0.2047 0.1889 0.1637
0.1913 0.1854 0.1612
0.2007 0.1876 0.1629
0.1751 0.1738 0.1564
0.1729 0.1715 0.1513
(1.3, 2.5, 2.4, 1.6)
5 10 15
0.0746 0.0673 0.0607
0.0687 0.0658 0.0546
0.0723 0.0662 0.0580
0.0664 0.0643 0.0531
0.0653 0.0620 0.0502
6. The above steps are repeated 1000 times and the estimated risk (ER) is computed by averaging the squared deviations over the 1000 repetitions. The computational results are displayed in Tables 1 and 2. 4.2. Numerical example In this subsection we are interested in finding the lower and upper 95% prediction bounds for the next upper record value in the sample of size n ¼ 10 when both of the parameters a and b are unknown and prior information about a and b suggested that l ¼ 1:5, m ¼ 3:0, c ¼ 1:0 and d ¼ 2:0 for the joint prior density given in (2.9). An upper record sample of size 10 is generated from the Gompða ¼ 1:483; b ¼ 5:572Þ model given by (1.1) and written in order form as: 0:12528; 0:21211; 0:22784; 0:26063; 0:65258; 0:66056; 0:68255; 0:79385; 0:83778; 0:92206 Using our results in equation (3.11), the lower and upper 95% prediction bounds for xU ðnþ1Þ , the next record value, are 0.86413 and 2.41785, respectively.
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
319
5. Concluding remarks In this paper, the maximum likelihood and Bayes methods of estimation are used for estimation of the parameters of the Gompertz distribution based on record values. Bayesian prediction bounds are obtained for the future record values. It has been noticed, from Tables 1 and 2, that the estimated risks of the estimates decrease as n increases and the Bayes estimates have the smallest estimated risks as compared with their corresponding maximum likelihood estimates. The estimated risks of the Bayes estimates that are obtained based on LINEX loss function are smaller than the corresponding estimated risks of the estimates, which are obtained, based on squared error loss function. The estimated risks decrease as a increases. Different values of the prior parameters l, m, c and d rather than those appearing in the above tables have been considered but did not change the previous conclusion. If the prior parameters are unknown, the empirical Bayes approach may be used to estimate such parameters, (see, for example, [12]). Appendix A. Unimodality of the posterior density To show that the posterior density function defined by (2.10) is unimodal, it suffices to show that the equation hðaÞ ¼ 0, where hðaÞ is given by (2.19), has a unique root, say a . It then follows, from (2.18), that the posterior density qða; bjxÞ has only one mode ða ; b Þ. If the functional relationship y ¼ hðaÞ is monotone decreasing and crosses the horizontal a-axis, then it will cross the a-axis in only one point a satisfying hða Þ ¼ 0. The curve of y ¼ hðaÞ is monotone decreasing if dy=da < 0, for all positive values of a. It follows, from (2.19), that dhðaÞ ¼ ðD1 þ D2 Þ; da
ðA:1Þ
where D1 ¼
nd2 >0 a2
D2 ¼ ðn þ m 1Þðl 1Þ
if
9 > > =
n > d þ 2; x2n eaxn
ðeaxn þ l 1Þ2
>0
for all
> n þ m > 1; l > 1: > ; ðA:2Þ
Therefore, dy=da < 0 for all a > 0, n > d þ 2 and l > 1. On the other hand, it can be shown that hðaÞ tends to þ1 as a ! 0þ and to a negative value as a ! 1, then the function y ¼ hðaÞ is continuous on (0; 1).
320
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
Therefore, y ¼ hðaÞ is a continuous function which decreases monotonically from positive to negative values. So that the curve of y ¼ hðaÞ must cross the horizontal a-axis at exactly one point a .
References [1] M. Ahsanullah, Record Statistics, Nova Science Publishers, Commack, NY, 1995. [2] J.C. Ahuja, S.W. Nash, The generalized Gompertz–Verhulst family of distributions, Sankhya, Part A 29 (1967) 141–156. [3] J. Aitchison, I.R. Dunsmore, Statistical Prediction Analysis, Cambridge University Press, 1975. [4] E.K. Al-Hussaini, G.R. Al-Dayian, S.A. Adham, On finite mixture of two-component Gompertz lifetime model, J. Statist. Comput. Simul. 67 (1) (2000) 1–15. [5] B.C. Arnold, N. Balakrishnan, H.N. Nagaraja, A First Course In Order Statistics, John Wiley and Sons, New York, 1992. [6] B.C. Arnold, N. Balakrishnan, H.N. Nagaraja, Records, John Wiley and Sons, New York, 1998. [7] R. Calabria, G. Pulcini, Point estimation under asymmetric loss functions for left truncated exponential samples, Commun. Statist. Theor. Meth. 25 (3) (1996) 285–600. [8] K.N. Chandler, The distribution and frequency of record values, J. Roy. Stat. Soc. B 14 (1952) 220–228. [9] M.L. Garg, B.R. Rao, K. Redmond, Maximum likelihood estimation of the parameters of the Gompertz survival function, J. Roy. Stat. Soc. Ser. Appl. Stat. 19 (1970) 152–159. [10] B. Gompertz, On the nature of the function expressive of the law of human mortality, and on a new model of determining the value of life contingencies, Philos. Trans. Roy. Soc. Lond. 115 (1825) 513–585. [11] N.L. Johnson, S. Kotz, N. Balakrishnan, Continuous Univariate Distributions, vol. 2, second ed., John Wiley and Sons, New York, 1995. [12] J.L. Maritz, T. Lwin, Empirical Bayes Methods, second ed., Chapman & Hall, London, 1989. [13] M.I. Osman, A new model for analyzing the survival of heterogenous data, Ph.D. thesis, Case Western Reserve University, Cleveland, OH, 1987. [14] A.L. Rukhin, M.M.N. Ananda, Risk behavior of variance estimators in multivariate normal distribution, Stat. Prob. Lett. 13 (1992) 159–166. [15] L. Tierney, J.B. Kadane, Accurate approximations for posterior moments and marginal densities, J. Amer. Stat. Assoc. 81 (1986) 82–86. [16] H.R. Varian, A Bayesian approach to real estate assessment, in: Studies in Bayesian Econometrics and Statistics in Honer of Leonard J. Savage, North Holland, Amsterdam, 1975, pp. 195–208. [17] IMSL Reference Manual, IMSL, Houston, TX, 1984.
A Bayesian analysis of record statistics from the Gompertz model Zeinhum F. Jaheen
*
Department of Mathematics, University of Assiut, Assiut, 71516, Egypt
Abstract The Gompertz distribution has been used as a growth model and it can be used to fit tumor growth. Record values can be viewed as order statistics from a sample whose size is determined by the values and the order of occurrence of observations. Based on record values from the two-parameter Gompertz distribution, Bayes estimators for the two unknown parameters are obtained by using Laplace approximation. These estimates are obtained based on the squared error and LINEX loss functions. Predictions for future upper record values from the Gompertz model are obtained from a Bayesian approach. The maximum likelihood and Bayes estimates are compared via Monte Carlo simulation study and a numerical example is given to illustrate the results of prediction. Ó 2002 Elsevier Inc. All rights reserved. Keywords: Growth model; Record values; Maximum likelihood; Bayes estimation; Laplace approximation; Squared error loss; LINEX loss; Bayes prediction; Numerical computations
1. Introduction The Gompertz model was formulated by Benjamin Gompertz [10] to fit mortality tables. It has been used as a growth model and also can be used to fit tumor growth. Ahuja and Nash [2] showed that the Gompertz distribution is related by a simple transformation to certain distributions in the family of distributions derived by Pearson. Garg et al. [9] obtained the maximum
* Address: Department of Mathematics, Umm Al-Qura University (TAIF), P.O. Box 3970, Taif, Saudi Arabia. E-mail address: [email protected] (Z.F. Jaheen).
0096-3003/$ - see front matter Ó 2002 Elsevier Inc. All rights reserved. doi:10.1016/S0096-3003(02)00489-7
308
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
likelihood estimates of the parameters of the Gompertz model. Osman [13] used the two-parameter Gompertz distribution to derive a compound Gompertz distribution assuming that one of the parameters is a random variable following the gamma distribution. He studied the properties of his distribution and suggested its use for modeling lifetime data and analyzing the survivals in heterogeneous populations. Applications and more recent survey for the Gompertz model are given by Al-Hussaini et al. [4]. The probability density function (pdf) and cumulative distribution function (cdf) of the Gompertz distribution (denoted by Gompða; bÞ) are given, respectively, by b f ðtÞ ¼ b exp at ½expðatÞ 1 ; t > 0; ða > 0; b > 0Þ; ð1:1Þ a b F ðtÞ ¼ 1 exp ½expðatÞ 1 : ð1:2Þ a The pdf of the Gompertz distribution is unimodal. It has positive skewness and an increasing hazard rate function. The Gompertz distribution can be shown to be a truncated form of the type I extreme value distribution [11]. One disadvantage when using squared error loss, is that it penalises overestimation or underestimation [14]. Overestimation of a parameter can lead to more severe or less severe consequences than underestimation, or vice versa. Subsequently, the use of an asymmetrical loss function, which associates greater importance to overestimation or underestimation, can be considered for the estimation of the parameters. The LINEX loss function is defined as LðDÞ ¼ eaD aD 1;
a 6¼ 0;
ð1:3Þ
where D ¼ /^ðhÞ /ðhÞ, the scalar estimation error, if /ðhÞ is estimated by /^ðhÞ [16]. The sign of a represents the direction and its magnitude represents the degree of symmetry. First, for a ¼ 1 the LINEX loss function is quite asymmetric about zero with overestimation being more costly than underestimation. Second, if a < 0, LðDÞ rises exponentially when D < 0 (underestimation) and almost linearly when D > 0 (overestimation). For a closed to zero, the LINEX is approximately squared error loss and therefore almost symmetric. The posterior expectation of the LINEX loss function in (1.3) is ^ E/ ð/^ðhÞ /ðhÞÞ ¼ ea/ðhÞ E/ ðea/ðhÞ Þ að/^ðhÞ E/ ð/ðhÞÞÞ 1;
ð1:4Þ
where E/ ðÞ denoting posterior expectation with respect to the posterior density of /. The Bayes estimator of /ðhÞ, denoted by /^BL , of the function /ðhÞ under the LINEX loss function is the value /^ðhÞ which minimizes (1.4), it is 1 /^BL ¼ ln E/ ðea/ðhÞ Þ ; ð1:5Þ a provided that E/ ðea/ðhÞ Þ exists and is finite (see [7]).
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
309
Chandler [8] introduced the study of record values and documented many of the basic properties of records. Record values can be viewed as order statistics from a sample whose size is determined by the values and the order of occurrence of observations. In a little over thirty years, a large number of publications devoted to records have appeared. This is possibly due to the fact that we encounter this notion frequently in daily life, especially in singling out record values from a set of others and in registering and recalling record values. Record values can be viewed as order statistics from a sample whose size is determined by the values and the order of occurrence of observations. Let X1 ; X2 ; . . . is a sequence of independent and identically distributed random variables with cdf F ðxÞ and pdf f ðxÞ. Set Yn ¼ maxðminÞX1 ; X2 ; . . . ; Xn , n P 1. We say Xj is an upper (lower) record of this sequence if Yj > ð<ÞYj1 ; j > 1. By definition, X1 is an upper as well as a lower record value. One can transform from upper record values to lower records by replacing the original sequence of random variables by Xj , j P 1 or (if P ðXi > 0 ¼ 1Þ for all i) by 1=Xi ; i P 1; the lower record values of this sequence will correspond to the upper record values of the original sequence. The notations XU ðnÞ and XLðnÞ are used for the nth upper and lower records, respectively. For more details on record values (see [1,5,6]). In this paper, Bayes estimates for the two unknown parameters a and b of the Gompertz distribution are obtained based on upper record values. The estimates are obtained based on the squared error and LINEX loss functions by using the approximation form of Tierney and Kadane [15]. Prediction bounds for future upper record values are obtained from a Bayesian approach. The maximum likelihood and Bayes estimates are compared via Monte Carlo simulation study. Numerical example is given to illustrate the results of prediction.
2. Estimation of the parameters In this section, we shall be concerned with estimation of the two unknown parameters a and b of the Gompertz model based on record values. Suppose we observe n upper record values XU ð1Þ ; XU ð2Þ ; . . . ; XUðnÞ from the Gompertz distribution with pdf given by (1.1). The likelihood function (LF) (see [1]) is given by Lða; bjxÞ ¼
n1 Y
H ðxi Þf ðxn Þ;
ð2:1Þ
i¼1
where x ¼ ðx1 ; x2 ; . . . ; xn Þ and H ðÞ is the hazard function corresponding to the pdf f ðÞ. It follows, from (1.1), (1.2) and (2.1), that
310
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320 n nax
Lða; bjxÞ ¼ b e where x ¼ 1n
Pn
i¼1
exp
b axn ðe 1Þ ; a
ð2:2Þ
xi .
2.1. Maximum likelihood estimation The natural logarithm of the LF (2.2) is b ‘ ln Lða; bjxÞ ¼ n lnðbÞ þ nax ðeaxn 1Þ: a
ð2:3Þ
Assuming that the parameters a and b are unknown, the maximum likelihood estimate (MLE) of the parameter b can be shown to be b^ ¼
n^ a ; 1
ð2:4Þ
ea^xn
where a^ is the MLE of the parameter a which can be obtained as a solution of the following non-linear equation ½^ aðx xn Þ þ 1ea^xn ð^ ax þ 1Þ ¼ 0:
ð2:5Þ
2.2. Bayes estimation Under the assumption that both of the parameters a and b are unknown, we may consider the joint prior density as a product of a conditional density of b for given a (which is taken to be the conjugate gamma prior when a is known) and a two parameter inverted gamma density for a. So that the joint prior density of a and b can be written as gða; bÞ ¼ g1 ðbjaÞg2 ðaÞ;
ð2:6Þ
where g1 ðajbÞ ¼
lm bm1 ebl=a ; CðmÞam
b > 0; ðl > 0; m > 0Þ;
ð2:7Þ
and g2 ðaÞ ¼
cd 1 c=a e ; CðdÞ adþ1
a > 0; ðc > 0; d > 0Þ:
ð2:8Þ
Multiplying g1 ðajbÞ by g2 ðaÞ, we obtain the joint prior density of a and b, given by (2.6), by lm cd bm1 1 ðbl þ cÞ ; a > 0; b > 0: ð2:9Þ gða; bÞ ¼ exp a CðmÞCðdÞ amþdþ1
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
311
The joint posterior density function of a; b given the data, denoted by qða; bjxÞ, can be written as Z Z qða; bjxÞ ¼ Lða; bjxÞgða; bÞ= Lða; bjxÞgða; bÞ da db; a
¼
A amþdþ1
b
nþm1
exp
b
b B1 ða; xn Þ þ B2 ða; xÞ ; a
ð2:10Þ
where B1 ða; xn Þ ¼ eaxn þ l 1; B2 ða; xÞ ¼ nax ac ;
ð2:11Þ
and A is the normalizing constant given by Z 1 ðnþmÞ and1 expfB2 ða; xÞg½B1 ða; xn Þ da: A1 ¼ Cðn þ mÞ
ð2:12Þ
0
2.2.1. Bayes estimators based on squared error loss function Under a squared error loss function, the Bayes estimator /~ of a function /ða; bÞ is given by Z Z /~ E½/ða; bÞjx ¼ /ða; bÞqða; bjxÞ da db ¼
a
Z Z
b
Z Z Lða; bjxÞgða; bÞ da db;
/ða; bÞLða; bjxÞgða; bÞ= a
b
a
ð2:13Þ
b
where Lða; bjxÞ and gða; bÞ are given, respectively, by (2.2) and (2.9). Generally, the ratio of two integrals given by (2.13) can not be obtained in a simple closed form. In this case, we can use numerical integration technique, which can be computationally intensive, especially in high dimensional parameter space. Instead, one can use Laplace approximation form due to Tierney and Kadane [15]. In the following, a review of Tierney and KadaneÕs approximation form is given. 2.2.2. The approximation form of Tierney and Kadane Tierney and Kadane [15] gave an approximate form for the evaluation of the ratio of integrals of the form (2.13) by writing the two expressions, LðkÞ ¼
1 ln qðkjxÞ n
and
So that (2.13) takes the form
L ðkÞ ¼ LðkÞ þ
1 ln /ðkÞ: n
ð2:14Þ
312
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
1=2 det X /^B ¼ exp½nL ðk^ Þ Lðk^Þ; det X 1=2 det X ¼ /ðk^ Þqðk^ jxÞ=qðk^jxÞ; det X
ð2:15Þ
where k^ and k^ maximize L ðkÞ and LðkÞ, respectively, and X and X are the negatives of the inverses of the matrices of second derivatives of L ðkÞ and LðkÞ, at k^ and k^, respectively. Here, we use the approximation form (2.15) to obtain Bayes estimators for the parameters a and b. The regularity condition required for using Tierney–KadaneÕs form (2.15) is that the posterior density should be unimodal. This condition is satisfied for our posterior density given by (2.10), see Appendix A. In our case, k ða; bÞ and the functions Lða; bÞ and L ða; bÞ are obtained from (2.14) and given by
Lða; bÞ ¼ 1n ln A ðm þ d þ 1Þ ln a þ ðn þ m 1Þ ln b ba B1 ða; xn Þ þ B2 ða; xn Þ ; L ða; bÞ ¼ Lða; bÞ þ 1n ln /ða; bÞ:
ð2:16Þ The mode of the posterior density (2.10), denoted by ðaM ; bM Þ, is the solution of the following non-linear equations 9 oL 1 m þ d þ 1 bB1 ða; xn Þ bxn eaxn > ¼ þ þ nxn ; > L1 ða; bÞ ¼ = oa n a a2 a oL 1 n þ m 1 B1 ða; xn Þ > > ; ¼ : L2 ða; bÞ ¼ ob n b a
ð2:17Þ
It follows, From (2.17), that aM ðn þ m 1Þ b^M ¼ ; B1 ðaM ; xn Þ
ð2:18Þ
and aM is a solution of the following non-linear equation hðaÞ ¼
nd2 xn eaxn ðn þ m 1Þ axn þ nxn ¼ 0: a e þl1
ð2:19Þ
It has been theoretically shown, see the Appendix A, that the equation hðaÞ has a unique root and then the posterior density (2.10) is a unimodal. To obtain det X, we first derive the second derivatives of Lða; bÞ with respect to a and b as
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
o2 L 1 m þ d þ 1 2bB1 ða; xn Þ 2beaxn ¼ þ oa2 n a2 a3 a2 o2 L 1 B1 ða; xn Þ xn eaxn ¼ L12 ¼ ; oaob n a2 a o2 L nþm1 ; L22 ¼ 2 ¼ ob nb2
L11 ¼
313
9 b 2 axn > > xn e ; > > > a > > = > > > > > > > ;
ð2:20Þ
where B1 ða; xn Þ is given by (2.11). Therefore, 2
det X ¼ 1=½L11 L22 ðL12 Þ ;
ð2:21Þ
evaluated at the posterior mode ð^ aM ; b^M Þ. Now, the Bayes estimates of a and b are computed as follows: (i) If /ða; bÞ ¼ a, we have from (2.16), a L ða; bÞ
¼ Lða; bÞ þ
1 ln a: n
ð2:22Þ
Thus, a L11
¼ L11
1 ; na2
a L12
¼ L12
and
a L22
¼ L22 :
ð2:23Þ
Hence, 2
det Xa ¼ 1=½a L11 a L22 ða L12 Þ ;
ð2:24Þ
evaluated at the mode of a L ða; bÞ, denoted by solving the two equations L1 ða; bÞ þ
1 ¼0 na
and
ð^ a1 ; b^1 Þ,
L2 ða; bÞ ¼ 0:
which is obtained by
ð2:25Þ
Substituting (2.21) and (2.24) in (2.15), the Bayes estimator of a takes the form !1=2 2 L11 L22 ðL12 Þ qð^ a1 ; b^1 jxÞ=qð^ aM ; b^M jxÞ; ð2:26Þ a^BS ¼ a^1 2 a L11 a L22 ða L12 Þ where qða; bjxÞ is the posterior density given by (2.10) evaluated at the modes ð^ a1 ; b^1 Þ and ð^ aM ; b^M Þ of the functions a L ða; bÞ and Lða; bÞ. (ii) If /ða; bÞ ¼ b, we have from (2.16), b L ða; bÞ
¼ Lða; bÞ þ
1 ln b: n
ð2:27Þ
Thus, b L11
¼ L11
1 ; nb2
b L12
¼ L12
and
b L22
¼ L22 :
ð2:28Þ
314
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
Hence, 2
det Xb ¼ 1=½b L11 b L22 ðb L12 Þ ;
ð2:29Þ
evaluated at the mode of b L ða; bÞ, denoted by ð^ a2 ; b^2 Þ, which is obtained by solving the two equations L1 ða; bÞ þ
1 ¼ 0; nb
L2 ða; bÞ ¼ 0:
ð2:30Þ
Substituting (2.21) and (2.29) in (2.15), the Bayes estimator of b takes the form !1=2 2 L11 L22 ðL12 Þ ^ ^ bBS ¼ b2 qð^ a2 ; b^2 jxÞ=qð^ aM ; b^M jxÞ; ð2:31Þ 2 L L ð L Þ b 11 b 22 b 12 where qða; bjxÞ is the posterior density given by (2.10) evaluated at the modes ð^ a2 ; b^2 Þ and ð^ aM ; b^M Þ of the functions a L ða; bÞ and Lða; bÞ. 2.2.3. Bayes estimators based on LINEX loss function Under the LINEX loss function (1.3), the Bayes estimator /~ of a function /ða; bÞ is given by (1.5), which can be rewritten as 1 /^BL ¼ ln Eðea/ jxÞ a Z Z Z Z 1 a/ e Lða; bjxÞgða; bÞ= Lða; bjxÞgða; bÞ da db : ¼ ln a a b a b ð2:32Þ The ratio of integrals (2.32) can also be computed by using Tierney and KadaneÕs approximate form (2.15). In the following, the Bayes estimators of a and b are computed as follows: (i) If /ða; bÞ ¼ a, we obtain from (2.16), a L ða; bÞ
1 ¼ Lða; bÞ aa: n
ð2:33Þ
Thus, a L11
¼ L11 ;
a L12
¼ L12
and
a L22
¼ L22 :
ð2:34Þ
Hence, 2
det X a ¼ 1=½a L11 a L22 ða L12 Þ ;
ð2:35Þ
^ evaluated at the mode of a L ða; bÞ, denoted by ð^ a 3 ; b3 Þ, which is obtained by solving the two equations a L1 ða; bÞ ¼ 0 and L2 ða; bÞ ¼ 0: ð2:36Þ n
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
315
Substituting (2.21) and (2.35) in (2.32), the Bayes estimator of a takes the form 2 31 , !1=2 2 1 L11 L22 ðL12 Þ ^ qð^ a qð^ aM ; b^M jxÞ5 ; a^BL ¼ ln 4a^ 3 3 ; b3 jxÞ 2 a L L ð L Þ a 12 a 11 a 22 ð2:37Þ where qða; bjxÞ is the posterior density given by (2.10) evaluated at the modes ^ ð^ a aM ; b^M Þ of the functions a L ða; bÞ and Lða; bÞ. 3 ; b3 Þ and ð^ (ii) If /ða; bÞ ¼ b, we have from (2.16), b L ða; bÞ
1 ¼ Lða; bÞ ab: n
ð2:38Þ
Thus, b L11
¼ L11 ;
b L12
¼ L12
and
b L22
¼ L22 :
ð2:39Þ
Hence, 2 det X b ¼ 1=½b L11 b L22 ðb L12 Þ ;
ð2:40Þ
evaluated at the mode of b L ða; bÞ, denoted by ða 4 ; b4 Þ, which is obtained by solving the two equations
L1 ða; bÞ ¼ 0;
L2 ða; bÞ
a ¼ 0: n
ð2:41Þ
Substituting (2.21) and (2.40) in Eq. (2.32), the Bayes estimator of b takes the form 2 31 , !1=2 2 1 L11 L22 ðL12 Þ ^ b^BL ¼ ln 4b^ qð^ a qð^ aM ; b^M jxÞ5 ; 3 4 ; b4 jxÞ 2 a L L ð L Þ b 11 b 22 b 22 ð2:42Þ where qða; bjxÞ is the posterior density given by (2.10) evaluated at the modes ^ ð^ a aM ; b^M Þ of the functions a L ða; bÞ and Lða; bÞ. 4 ; b4 Þ and ð^
3. Prediction of future record values Suppose that we have n upper records XU ð1Þ ¼ x1 ; XU ð2Þ ¼ x2 ; . . . ; XU ðnÞ ¼ xn from the Gomp(a; b) distribution. Based on such a record sample, Bayesian prediction is needed for the sth upper record, 1 < n < s. Let Y XU ðsÞ be the sth upper record value, 1 < n < s. The conditional pdf of Y for given xn is given, (see [1]), by
316
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
½vðyÞ vðxn Þsn1 f ðy; hÞ ; 1 F ðxn ; hÞ Cðs nÞ
f ðyjxn ; hÞ ¼
ð3:1Þ
where vðÞ ¼ lnð1 F ðÞÞ. For the Gomp(a; b) distribution, with pdf given by (1.1), the function f ðyjxn ; hÞ is then given by sn1 beay b b B3 ða; y; xn Þ f ðyjxn ; a; bÞ ¼ exp B3 ða; y; xn Þ ; ð3:2Þ a Cðs nÞ a where B3 ða; y; xn Þ ¼ eay eaxn :
ð3:3Þ
Following Aitchison and Dunsmore [3], the Bayes predictive density function of y given x is given by Z f ðyjxÞ ¼ f ðyjxn ; hÞqðhjxÞ dh: ð3:4Þ H
Substituting (2.10) and (3.2) in (3.4), we get f ðyjxÞ ¼ A
Cðs þ mÞ Iðy; xÞ; Cðs nÞ
ð3:5Þ
where Iðy; xÞ ¼
Z
1
and eay expfB2 ða; xÞg½B3 ða; y; xn Þsn1 ½B1 ða; yÞðsþmÞ da:
0
ð3:6Þ Bayesian prediction bounds for Y ¼ XU ðsÞ , given the previous data x, are obtained by evaluating PrðY P kjxÞ, for some positive k. It follows, from (3.5), that Z 1 PrðY P kjxÞ ¼ f ðyjxÞ dy; k
¼A
Ij ðk; xÞ Cðs þ mÞ sn1 j sn1 Xj¼0 ð1Þ ; j Cðs nÞ nþmþj
ð3:7Þ
where Ij ðk; xÞ ¼
Z 0
1
"
# j fB ða; x Þg 1 n da: and1 expfB2 ða; xÞg fB1 ða; kÞgnþmþj
ð3:8Þ
It can be easily shown that f ðyjxÞ is a density function on the positive half of the real line by proving that PrðY P xn jxÞ ¼ 1.
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
317
A 100s% Bayesian prediction interval for Y XU ðsÞ is such that P ½LLðxÞ < Y < ULðxÞ ¼ s;
ð3:9Þ
where LLðxÞ and ULðxÞ are the lower and upper limits satisfying 9 P ½Y > LLðxÞjx ¼ ð1 þ sÞ=2; = and ; P ½Y > ULðxÞjx ¼ ð1 sÞ=2:
ð3:10Þ
Special case: For the special case, when s ¼ n þ 1, which is practically of special interest, Bayesian prediction bounds for the future upper record value Ynþ1 ¼ XU ðnþ1Þ are obtained from PrðYnþ1 P k1 jxÞ ¼ gðk1 Þ=gðxn Þ where gðk1 Þ ¼
Z
1
ð3:11Þ
ðnþmÞ
and1 expfB2 ða; xÞg½B1 ða; k1 Þ
da;
ð3:12Þ
0
and gðxn Þ is obtained from (3.12) with k1 being replaced by xn . 4. Numerical computations In the following, the maximum likelihood and Bayes (squared error and LINEX) estimates are compared based on a Monte Carlo simulation study and a numerical example is given to illustrate the results of prediction. 4.1. Monte Carlo simulation The estimates obtained in Section 2 are computed and compared according to the following steps: 1. For a given vector of prior parameters ðl; m; c; dÞ, we generate a and b from the joint prior density (2.9). The IMSL [17] is used in the generation of the gamma random variates. 2. For given a; b obtained in step (1), we generate n ¼ ð5; 10; 15Þ upper record values from the Gompða; bÞ with pdf (1.1). 3. The ML estimate of a is computed by solving the nonlinear equation (2.5) by using ZSPOW routine from the IMSL [17] library. Substituting the MLE of a in (2.4) gives the MLE of b. 4. The Bayes estimates of a and b are computed from (2.26), (2.31), (2.37) and (2.42). 5. The squared deviations ðw wÞ2 are computed for different sizes n where () stands for an estimate (ML or Bayes) and w stands for the parameter (a or b).
318
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
Table 1 ER of the estimates of a for different values of l, m, c, d, n, and a (l, m, c, d)
n
ERð^ aML Þ
ERð^ aBS Þ
ERð^ aBL Þ a¼4
a¼8
a ¼ 16
(1.5, 3.0, 1.0, 2.0)
5 10 15
0.1425 0.1138 0.1057
0.1337 0.1124 0.1027
0.1352 0.1131 0.1036
0.1327 0.1119 0.1018
0.1286 0.1108 0.1009
(1.3, 2.5, 2.4, 1.6)
5 10 15
0.0532 0.0460 0.0418
0.0517 0.0399 0.0376
0.0636 0.0415 0.0402
0.0493 0.0365 0.0325
0.0458 0.0330 0.0307
a¼4
a¼8
a ¼ 16
Table 2 ER of the Estimates of b for different values of l, m, c, d, n, and a (l, m, c, d)
n
ERðb^ML Þ
ERðb^BS Þ
ERðb^BL Þ
(1.5, 3.0, 1.0, 2.0)
5 10 15
0.2047 0.1889 0.1637
0.1913 0.1854 0.1612
0.2007 0.1876 0.1629
0.1751 0.1738 0.1564
0.1729 0.1715 0.1513
(1.3, 2.5, 2.4, 1.6)
5 10 15
0.0746 0.0673 0.0607
0.0687 0.0658 0.0546
0.0723 0.0662 0.0580
0.0664 0.0643 0.0531
0.0653 0.0620 0.0502
6. The above steps are repeated 1000 times and the estimated risk (ER) is computed by averaging the squared deviations over the 1000 repetitions. The computational results are displayed in Tables 1 and 2. 4.2. Numerical example In this subsection we are interested in finding the lower and upper 95% prediction bounds for the next upper record value in the sample of size n ¼ 10 when both of the parameters a and b are unknown and prior information about a and b suggested that l ¼ 1:5, m ¼ 3:0, c ¼ 1:0 and d ¼ 2:0 for the joint prior density given in (2.9). An upper record sample of size 10 is generated from the Gompða ¼ 1:483; b ¼ 5:572Þ model given by (1.1) and written in order form as: 0:12528; 0:21211; 0:22784; 0:26063; 0:65258; 0:66056; 0:68255; 0:79385; 0:83778; 0:92206 Using our results in equation (3.11), the lower and upper 95% prediction bounds for xU ðnþ1Þ , the next record value, are 0.86413 and 2.41785, respectively.
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
319
5. Concluding remarks In this paper, the maximum likelihood and Bayes methods of estimation are used for estimation of the parameters of the Gompertz distribution based on record values. Bayesian prediction bounds are obtained for the future record values. It has been noticed, from Tables 1 and 2, that the estimated risks of the estimates decrease as n increases and the Bayes estimates have the smallest estimated risks as compared with their corresponding maximum likelihood estimates. The estimated risks of the Bayes estimates that are obtained based on LINEX loss function are smaller than the corresponding estimated risks of the estimates, which are obtained, based on squared error loss function. The estimated risks decrease as a increases. Different values of the prior parameters l, m, c and d rather than those appearing in the above tables have been considered but did not change the previous conclusion. If the prior parameters are unknown, the empirical Bayes approach may be used to estimate such parameters, (see, for example, [12]). Appendix A. Unimodality of the posterior density To show that the posterior density function defined by (2.10) is unimodal, it suffices to show that the equation hðaÞ ¼ 0, where hðaÞ is given by (2.19), has a unique root, say a . It then follows, from (2.18), that the posterior density qða; bjxÞ has only one mode ða ; b Þ. If the functional relationship y ¼ hðaÞ is monotone decreasing and crosses the horizontal a-axis, then it will cross the a-axis in only one point a satisfying hða Þ ¼ 0. The curve of y ¼ hðaÞ is monotone decreasing if dy=da < 0, for all positive values of a. It follows, from (2.19), that dhðaÞ ¼ ðD1 þ D2 Þ; da
ðA:1Þ
where D1 ¼
nd2 >0 a2
D2 ¼ ðn þ m 1Þðl 1Þ
if
9 > > =
n > d þ 2; x2n eaxn
ðeaxn þ l 1Þ2
>0
for all
> n þ m > 1; l > 1: > ; ðA:2Þ
Therefore, dy=da < 0 for all a > 0, n > d þ 2 and l > 1. On the other hand, it can be shown that hðaÞ tends to þ1 as a ! 0þ and to a negative value as a ! 1, then the function y ¼ hðaÞ is continuous on (0; 1).
320
Z.F. Jaheen / Appl. Math. Comput. 145 (2003) 307–320
Therefore, y ¼ hðaÞ is a continuous function which decreases monotonically from positive to negative values. So that the curve of y ¼ hðaÞ must cross the horizontal a-axis at exactly one point a .
References [1] M. Ahsanullah, Record Statistics, Nova Science Publishers, Commack, NY, 1995. [2] J.C. Ahuja, S.W. Nash, The generalized Gompertz–Verhulst family of distributions, Sankhya, Part A 29 (1967) 141–156. [3] J. Aitchison, I.R. Dunsmore, Statistical Prediction Analysis, Cambridge University Press, 1975. [4] E.K. Al-Hussaini, G.R. Al-Dayian, S.A. Adham, On finite mixture of two-component Gompertz lifetime model, J. Statist. Comput. Simul. 67 (1) (2000) 1–15. [5] B.C. Arnold, N. Balakrishnan, H.N. Nagaraja, A First Course In Order Statistics, John Wiley and Sons, New York, 1992. [6] B.C. Arnold, N. Balakrishnan, H.N. Nagaraja, Records, John Wiley and Sons, New York, 1998. [7] R. Calabria, G. Pulcini, Point estimation under asymmetric loss functions for left truncated exponential samples, Commun. Statist. Theor. Meth. 25 (3) (1996) 285–600. [8] K.N. Chandler, The distribution and frequency of record values, J. Roy. Stat. Soc. B 14 (1952) 220–228. [9] M.L. Garg, B.R. Rao, K. Redmond, Maximum likelihood estimation of the parameters of the Gompertz survival function, J. Roy. Stat. Soc. Ser. Appl. Stat. 19 (1970) 152–159. [10] B. Gompertz, On the nature of the function expressive of the law of human mortality, and on a new model of determining the value of life contingencies, Philos. Trans. Roy. Soc. Lond. 115 (1825) 513–585. [11] N.L. Johnson, S. Kotz, N. Balakrishnan, Continuous Univariate Distributions, vol. 2, second ed., John Wiley and Sons, New York, 1995. [12] J.L. Maritz, T. Lwin, Empirical Bayes Methods, second ed., Chapman & Hall, London, 1989. [13] M.I. Osman, A new model for analyzing the survival of heterogenous data, Ph.D. thesis, Case Western Reserve University, Cleveland, OH, 1987. [14] A.L. Rukhin, M.M.N. Ananda, Risk behavior of variance estimators in multivariate normal distribution, Stat. Prob. Lett. 13 (1992) 159–166. [15] L. Tierney, J.B. Kadane, Accurate approximations for posterior moments and marginal densities, J. Amer. Stat. Assoc. 81 (1986) 82–86. [16] H.R. Varian, A Bayesian approach to real estate assessment, in: Studies in Bayesian Econometrics and Statistics in Honer of Leonard J. Savage, North Holland, Amsterdam, 1975, pp. 195–208. [17] IMSL Reference Manual, IMSL, Houston, TX, 1984.