- Email: [email protected]
Inference for the Weibull distribution with progressive hybrid censoring
Computational Statistics and Data Analysis 56 (2012) 451–467
Contents lists available at SciVerse ScienceDirect
Computational Statistics and Data Analysis journal homepage: www.elsevier.com/locate/csda
Inference for the Weibull distribution with progressive hybrid censoring Chien-Tai Lin ∗ , Cheng-Chieh Chou, Yen-Lung Huang Department of Mathematics, Tamkang University, Tamsui 25137, Taiwan
article
info
Article history: Received 4 April 2010 Received in revised form 25 July 2011 Accepted 3 September 2011 Available online 16 September 2011 Keywords: Gibbs sampling Lindley’s approximation Markov Chain Monte Carlo method Metropolis–Hastings algorithm Tierney–Kadane’s approximation
abstract Recently, progressive hybrid censoring schemes have become quite popular in life-testing and reliability studies. In this paper, we investigate the maximum likelihood estimation and Bayesian estimation for a two-parameter Weibull distribution based on adaptive Type-I progressively hybrid censored data. The Bayes estimates of the unknown parameters are obtained by using the approximation forms of Lindley (1980) and Tierney and Kadane (1986) as well as two Markov Chain Monte Carlo methods under the assumption of gamma priors. Computational formulae for the expected number of failures is provided and it can be used to determine the optimal adaptive Type-I progressive hybrid censoring schemes under a pre-determined budget of experiment. © 2011 Elsevier B.V. All rights reserved.
1. Introduction A new progressive hybrid censoring scheme—adaptive Type-I progressive hybrid censoring scheme (adaptive Type-I PHCS), introduced by Lin and Huang (2011), is a mixture of Type-I and Type-II progressive censoring schemes, and it can be described as follows. Suppose n identical units are placed on test with progressive censoring scheme (R1 , R2 , . . . , Rm ), 1 ≤ m ≤ n, and the experiment is terminated at time T , where T ∈ (0, ∞) and integers Ri ’s are fixed in advance. At the time of the first failure Y1:m:n , R1 of the remaining units are randomly removed. Similarly, at the time of the second failure Y2:m:n , R2 of the remaining units are randomly removed and so on. Let D denote the number of failures that occur before time T . If the mth failure Ym:m:n occurs before time T , the process will not stop at the time Ym:m:n , but continue to observe ∑D failures (without any further withdrawals) up to time T . Then, at time T all the remaining R⋆D = n − D − i=1 Ri units are removed and the experiment is terminated. Note that the progressive censoring scheme in this case will become (R1 , R2 , . . . , Rm , Rm+1 , . . . , RD ), where Rm = Rm+1 = · · · = RD = 0. On the other hand, the process when Ym:m:n > T will have a progressive censoring scheme as (R1 , R2 , . . . , RD ). The work on progressive hybrid censoring has become quite popular in life-testing and reliability studies. Kundu and Joarder (2006) and Childs et al. (2008) have considered a Type-I PHCS in the context of life-testing experiment is terminated at time min{Ym:m:n , T }. Kundu and Joarder (2006) discussed the maximum likelihood and Bayesian estimation for the exponential distribution, whereas Childs et al. (2008) dealt with the derivation of the exact distribution of the maximum likelihood estimator (MLE) of the unknown parameter. In the same paper, Childs et al. (2008) has also proposed another Type-II PHCS that would terminate the experiment at the random time max{Ym:m:n , T }, and derived the corresponding exact distribution of MLE. Note that on the case Ym:m:n < T , the process under this censoring scheme will be performed exactly the same as our proposed scheme. Ng et al. (2009) and Lin et al. (2009) have investigated point and interval estimation for exponential and Weibull lifetimes under an adaptive Type-II PHCS. When Ym:m:n < T , adaptive Type-II PHCS is equivalent to Type-I PHCS. In contrast, no more units will be removed whenever the experimental time passes time T , and all the remaining
∗
Corresponding author. Tel.: +886 2 26263976; fax: +886 2 26209916. E-mail address: [email protected] (C.-T. Lin).
0167-9473/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.csda.2011.09.002
452
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
Rm = n − R1 −· · ·− RD − m units are removed from the experiment immediately after the mth failure has occurred. Recently, Lin and Huang (2011) studied point and interval estimation for exponential lifetimes under an adaptive Type-I PHCS, and investigated Bayesian sampling plans under different progressive hybrid censoring schemes when a general loss function is used. Upon applying the Lindley’s 1980 approximation and Gibbs sampling procedure, Kundu (2007, 2008) and Banerjee and Kundu (2008) have studied Bayesian inference of parameters for a Weibull distribution based on Type-I hybrid, Type-II hybrid, or Type-II progressive censored sample. The Lindley’s approximation is generally accurate enough in these works, but, as Lindley has pointed out, the required evaluation of third derivatives of the posterior can be rather tedious—in particular, in problems with several parameters. Thus, in this paper, an easily computable approximation of Tierney and Kadane (1986) is considered to approximate the Bayes estimates of the parameters for a Weibull distribution based on an adaptive Type-I progressive hybrid censored sample. Note that the marginal posterior density function of the shape parameter cannot be obtained analytically. For this reason, Kundu (2007) has also used a Gibbs sampling technique along with the approach for generating random variates from a log-concave density (see Devroye, 1984) to construct the Bayes estimates and credible intervals for the unknown parameters. In order to make a comparison to this procedure, the Metropolis–Hastings algorithm with normal proposal distribution is then employed. The rest of this article is organized as follows. In Section 2, we derive the MLEs of two unknown parameters of Weibull distribution and the corresponding asymptotic variances and covariance. In Section 3, we cover the Bayesian analysis based on the approximation forms of Lindley (1980) and Tierney and Kadane (1986) as well as two Markov Chain Monte Carlo (MCMC) methods. In Section 4, we give the computational formulae for the expected number of failures, which will be used to find the optimal adaptive Type-I PHCS in Section 6. Numerical comparisons of the MLEs and the Bayes estimates are presented in Section 5. In Section 6, we propose a procedure to determine the optimality adaptive Type-I PHCS when a cost constraint given in Eq. (14), which includes the sampling cost, the time-consuming cost, and the salvage, is used. 2. Maximum likelihood estimation In this section we discuss the MLEs of the parameters µ and σ based on the data observed under the proposed censoring scheme and derive their asymptotic variance–covariance matrix by using the results of Ng et al. (2002) and Lin et al. (2006). 2.1. Maximum likelihood estimators Consider lifetime of units in the life test Yi , i = 1, . . . , n, to be independent and identically distributed as Weibull with probability density function (pdf) δ
fY (y; δ, λ) = λδ yδ−1 e−λy ,
y > 0,
where δ > 0 is the shape parameter and λ > 0 is the scale parameter. Note that in the derivation of MLEs and expected number of failures in Section 4, instead of working with the parametric model for Yi , it is often more convenient to work with the equivalent model for the log-lifetimes Xi = ln Yi for i = 1, 2, . . . , n, which follows the extreme value distribution with parameters µ = − ln λ/δ and σ = 1/δ . Now, suppose that an adaptive Type-I progressive hybrid censored extreme value sample with a given adaptive Type-I PHCS (R1 , . . . , Rm ) is observed. We denote the resulting failure times by X1:m:n , X2:m:n , . . . , Xm:m:n , Xm+1:n , . . . , XD:n when Xm:m:n < TL , and X1:m:n , X2:m:n , . . . , XD:m:n when Xm:m:n > TL , where TL ≡ ln T . Then the likelihood function of the given sample is given by L(µ, σ ) = cD
D ∏
⋆
fX (xi )[1 − FX (xi )]Ri [1 − FX (TL )]RD ,
(1)
i=1
where fX (x; µ, σ ) =
1
σ
[ exp
x−µ
FX (x; µ, σ ) = 1 − exp − exp D ∏
− exp
σ [
cD =
x−µ
m −
,
]
,
−∞ < x < ∞,
−∞ < x < ∞,
Rj ,
(2)
j =i
i=1
and R⋆D = n − D −
σ ]
σ
γi with γi = m − i + 1 +
x−µ
∑D
i=1
Ri with Rm = Rm+1 = · · · = RD = 0 if D ≥ m. Here we use the usual conventions that
∏0
j =1
αj ≡ 1
and j=i αj ≡ 0. Also, for simplicity of notation, we use xi instead of xi:m:n or xi:n , and yi instead of yi:m:n or yi:n in the remaining discussion.
∑i−1
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
453
From (1), the log-likelihood function given D ≥ 1 is ln L(µ, σ ) = ln(cD ) − D ln σ +
D −
ξi −
D −
i=1
(Ri + 1) exp(ξi ) − R⋆D exp(tL ),
i=1
where
ξi = (xi − µ)/σ and tL = (TL − µ)/σ .
(3)
Thus, the MLE of σ , denoted as σˆ , can be obtained by D ∑
D ∑
xi
i=1
D
=
(Ri + 1)xi exp
xi
+ R⋆D TL exp
xi
+ R⋆D exp
σ
i =1 D ∑
(Ri + 1) exp
i=1
σ
TL
σ
− σ.
(4)
TL
σ
Since Eq. (4) cannot be solved analytically, some numerical methods such as a bisection algorithm or other iterative procedures are needed. Alternatively, we can follow the work of Lin et al. (2009) to obtain the approximate maximum likelihood estimates of the parameters, and use them as the initial values in the iterative procedure for solving Eq. (4). After obtaining σˆ , the MLE of µ, denoted as µ ˆ , can be easily computed as
D ∑
(Ri + 1) exp
i =1 µ ˆ = σˆ ln
xi σˆ
+ R⋆D exp
D
TL σˆ .
Note that the MLEs do not exist when D = 0. 2.2. Asymptotic variances and covariance Using the missing-information principle, we can compute the observed Fisher information matrix. For θ = (µ, σ )′ , we define X , Z and W to be the observed, missing and complete data, and IX (θ), IW |X (θ) and IW (θ) to be the corresponding Fisher information matrix, respectively. Then IX (θ) = IW (θ) − IW |X (θ).
(5)
The complete information matrix for the extreme value distribution (see Ng et al., 2002) is IW (θ) =
n
1−γ c2
1 1−γ
σ2
,
where γ = 0.577215665 is the Euler’s constant and c 2 is π 2 /6 + (1 − γ 2 ) = 1.823680661. Now IW |X (θ) =
D −
(j)
∗ Rj IW |X (θ) + R∗D IW |X (θ),
(6)
j=1
(j)
∗ where IW |X (θ) and IW |X (θ) are the information matrix of a single observation for the truncated extreme value distribution with left truncation at xj and TL , respectively. The logarithm of the pdf of the truncated extreme value distribution with left truncation at xj is
ln fZj (zj |Zj > xj , µ, σ ) = − log σ + exp(ξj ) +
zj − µ
− exp
σ
zj − µ
σ
,
where ξj is as defined in Eq. (3). Thus, the three second partial derivatives with respect to µ and σ are
∂ 2 ln fZj ∂µ2 ∂ 2 ln fZj ∂µ∂σ ∂ ln fZj
= =
2
∂σ 2
=
1
[
exp(ξj ) − exp
σ2 1
zj − µ
]
σ
[
exp(ξj ) + ξj exp(ξj ) − exp
σ2 1
,
σ2
ξj2 exp(ξj ) + 2ξj exp(ξj ) −
−2
zj − µ
σ
exp
zj − µ
σ
zj − µ
σ
zj − µ
+2
−
zj − µ
σ
2
σ
exp
zj − µ
σ
exp
zj − µ
σ
+1 .
zj − µ
σ
] +1 ,
454
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
Taking the negative of the expected value of these three second partial derivatives with the results (see Ng et al., 2002 and Lin et al., 2006)
[ E
Zj − µ
σ
[ E
exp
] ∞ − exp[(p + 1)ξj ] Zj > xj , µ, σ = ψ(1) exp(eξj ) − [ψ(1) + ξj − ψ(p + 2)] , Γ (p + 2) p=0
Zj − µ
] ∞ − exp[(p + 2)ξj ] Zj > xj , µ, σ = exp(eξj ) − , Γ (p + 3)
σ p=0 ] ∞ − Zj − µ Zj − µ exp[(p + 2)ξj ] Zj > xj , µ, σ = ψ(2) exp(eξj ) − [ψ(2) + ξj − ψ(p + 3)] E exp , σ σ Γ (p + 3) p=0 2 Zj − µ Zj − µ E exp Zj > xj , µ, σ = [ψ ′ (2) + ψ 2 (2)] exp(eξj ) − {[ξj + ψ(2) − ψ(p + 3)]2 σ σ ∞ − exp[(p + 2)ξj ] , + ψ ′ (2) − ψ ′ (p + 3)} Γ (p + 3) p=0 [
where ψ(t ) =
d dt
ln Γ (t ) and ψ ′ (t ) =
d dt
ψ(t ) are the digamma and trigamma functions (see Abramowitz and Stegun,
(j) 1964), the missing information matrix IW |X (θ) can be easily computed. Replacing ξj with tL defined in Eq. (3), the information ∗ IW |X (θ) can also be computed with the same arguments. Then, from (5) and (6), the asymptotic variance–covariance matrix
of the maximum likelihood estimate θ can be obtained by inverting the observed information matrix IX ( θ). 3. Bayesian estimation In this section, we consider the Bayesian estimation under the assumption that the prior distribution of λ is Gamma(a, b) with pdf given by
π1 (λ|a, b) =
ba
λa−1 e−λb , λ > 0, Γ (a) where a > 0 and b > 0; and the prior of δ, π2 (δ), has the support on (0, ∞), and it is independent of the prior of λ. Thus, the joint posterior density function of δ and λ is given by π (δ, λ)ℓ(y |δ, λ) , π (δ, λ|y ) ∝ ∞ ∞ π (δ, λ)ℓ(y |δ, λ)dδ dλ 0 0 where π (δ, λ) = π1 (λ|a, b)π2 (δ) is the joint prior density function and
ℓ(y |δ, λ) = cD δ λ
D D
D ∏
δ−1
yi
D − exp −λ (Ri + 1)yδi + R⋆D T δ
i =1
i=1
is the likelihood function based on the adaptive Type-I progressively hybrid censored sample y = (y1 , . . . , yD ). Then, under a squared error loss, the Bayes estimate of any function of δ and λ, say g (δ, λ), is g (δ, λ)eQ (δ,λ) dδ dλ 0∞ ∞ , (7) π (δ, λ)ℓ(y |δ, λ)dδ dλ eQ (δ,λ) dδ dλ 0 0 0 0 where Q (δ, λ) = ln π (δ, λ) + ln ℓ(y |δ, λ) ≡ ρ(δ, λ) + L(δ, λ). Note that Eq. (7) cannot be obtained analytically even when π2 (δ) is known. Therefore, we adopt two approximations—Lindley’s approximation and Tierney and Kadane’s
∞∞
E (g (δ, λ)|y ) =
0
0
g (δ, λ)π (δ, λ)ℓ(y |δ, λ)dδ dλ
∞∞
∞∞
=
0
approximation, as well as two MCMC methods for approximating (7). 3.1. Lindley’s approximation Lindley (1980) has proposed approximations for the ratio in (7). The basic idea is to expand L(δ, λ) and ρ(δ, λ) of (7) into a Taylor series expansion about the MLE of (δ, λ), or Q (δ, λ) of (7) about the posterior mode of (δ, λ). The integrations all involve the moments that capture the first-order error terms of the normal approximation. This approximation can be viewed as an application of the Laplace method for integrals, and it will produce reasonable results as long as the posterior is unimodal or at least dominated by a single mode. For the two-parameter case, (δ, λ) = (λ1 , λ2 ), Lindley’s approximation of (7) is the form g (λ1 , λ2 ) +
1− 2
gij σij +
−
Uj ρj +
1
1
2
2
1 2
L30 σ11 U1 +
+ L12 (σ22 U1 + 2σ12 U2 ) + L03 σ22 U2 ,
1 2
L21 (2σ12 U1 + σ11 U2 ) (8)
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
455
where Uk = j gj σkj with σij being the (i, j)th elements of the inverse of the Fisher information matrix, and each suffix of variable g , ρ , or L denotes differentiation once with respect to the variable having that suffix. Evaluating all the expressions in (8) at the MLEs of λ1 and λ2 produces the approximation gˆBL to (7). In our case, with the specification of δ ∼ Gamma(α, β), we have
∑
D −
L(δ, λ) = ln(cD ) + D ln δ + D ln λ + (δ − 1)
D − δ ⋆ δ ln(yi ) − λ (Ri + 1)yi + RD T .
i=1
i=1
The MLEs of δ and λ are the solutions of the equations L10 =
L01 =
D
D −
ln(yi ) − λ
D −
δ
+
D
D − δ ⋆ δ − (Ri + 1)yi + RD T = 0.
λ
δ
(Ri + 1)yi ln(yi ) + RD T ln T
δ
i=1
⋆
= 0,
i=1
i =1
Now, to apply Lindley’s form in (8), we first obtain L20 = −
L11
D − ˆ ˆ δ 2 ⋆ δ 2 − λˆ (Ri + 1)yi (ln yi ) + RD T (ln T ) ,
D
δˆ 2 i=1 D − δˆ ⋆ δˆ =− (Ri + 1)yi ln(yi ) + RD T ln T ,
L02 = −
i=1
L30 =
L21
2D
δˆ 3
D − ˆδ ˆ (Ri + 1)yi (ln yi )3 + R⋆D T δ (ln T )3 , − λˆ
D
λˆ 2
,
L03 =
i=1
D − δˆ 2 ⋆ δˆ 2 =− (Ri + 1)yi (ln yi ) + RD T (ln T ) ,
2D
λˆ 3
,
L12 = 0;
i=1
and W
σ11 =
SW − V 2
,
S
σ22 =
SW − V 2
,
σ12 = σ21 =
−V SW − V 2
,
where S=
V =
D
δˆ 2
+ λˆ
D −
D − δˆ 2 ⋆ δˆ 2 (Ri + 1)yi (ln yi ) + RD T (ln T ) , i=1
(Ri + 1)yδi ln yi + R⋆D T δ ln T , ˆ
ˆ
W =
i=1
D
λˆ 2
.
For the posterior mean of δ , we have g (δ, λ) = δ and hence g1 = 1, g2 = 0, gij = 0 for i, j = 1, 2. Thus, 1
3
1
2
2
2
2 δˆBL = δˆ + ρ1 σ11 + ρ2 σ12 + L30 σ11 + L21 σ11 σ12 + L03 σ12 σ22 ,
where ρ1 =
α−1 δˆ
− β and ρ2 =
a−1
λˆ
− b. Similarly, for the posterior mean of λ, we have g (δ, λ) = λ and
1
1
1
2
2
2
2 2 + σ11 σ22 ) + L03 σ22 . λˆ BL = λˆ + ρ1 σ12 + ρ2 σ22 + L30 σ11 σ12 + L21 (2σ12
3.2. Tierney and Kadane’s approximation Setting H (δ, λ) = Q (δ, λ)/n and H ∗ (δ, λ) = [ln g (δ, λ) + Q (δ, λ)]/n, Tierney and Kadane (1986) reexpressed the expression in (7) for the posterior mean of g as
∞∞ E (g (δ, λ)|y ) = 0∞ 0∞ 0
0
∗ enH (δ,λ) dδ dλ
enH (δ,λ) dδ dλ
,
and then applied Laplace’s method to produce the approximation gˆBT (δ, λ) =
[
det Σ ∗ det Σ
]1/2
˜ ∗ ) − H (δ˜ , λ) ˜ exp n H ∗ (δ˜ ∗ , λ
(9)
456
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
˜ ∗ ) and (δ, ˜ λ) ˜ maximize H ∗ (δ, λ) and H (δ, λ), respectively, and Σ ∗ and Σ are the negatives of to E (g (δ, λ)|y ), where (δ˜ ∗ , λ ∗ ˜ ∗ ) and (δ, ˜ λ) ˜ , respectively. the inverse Hessians of H (δ, λ) and H (δ, λ) at (δ˜ ∗ , λ In our case, we have H (δ, λ) =
1
k + (α + D − 1) ln δ + (a + D − 1) ln λ − δβ + (δ − 1)
n
D −
ln(yi ) − λ
i=1
D −
δ
⋆
δ
(Ri + 1)yi + RD T + b
,
i=1
˜ can be obtained by solving the following two equations where k is a constant; therefore, (δ˜ , λ) D D − − α+D−1 −β + ln(yi ) − λ (Ri + 1)yδi ln(yi ) + R⋆D T δ ln T = 0, δ i =1 i =1 D − ∂ H (δ, λ) 1 a+D−1 δ ⋆ δ = − (Ri + 1)yi + RD T + b = 0, ∂λ n λ i=1 ∂ H (δ, λ) 1 = ∂δ n
and, from the second derivatives of H (δ, λ),
D − ∂ 2 H (δ, λ) 1 α+D−1 δ 2 ⋆ δ 2 = −λ − (Ri + 1)yi (ln yi ) + RD T (ln T ) , ∂δ 2 n δ2 i =1 ∂ 2 H (δ, λ) 1 a+D−1 = − , ∂λ2 n λ2 D − ∂ 2 H (δ, λ) 1 = − (Ri + 1)yδi ln(yi ) + R⋆D T δ ln T , ∂λ∂δ n i=1 ˜ λ) ˜ is given by the determinant of the negative of the inverse Hessian of H (δ, λ) at (δ, 2 −1 det Σ = (H11 H22 − H12 ) ,
where H11
H22
H12
D − 1 α+D−1 ∂ 2 H (δ, λ) ˜δ ˜ = (Ri + 1)yi (ln yi )2 + R⋆D T δ (ln T )2 , + λ˜ =− ∂δ 2 δ, n δ˜ 2 ˜ λ˜ i=1 1 a+D−1 ∂ 2 H (δ, λ) = =− , ∂λ2 δ, n λ˜ 2 ˜ λ˜ D 1 − ∂ 2 H (δ, λ) δ˜ ⋆ δ˜ = = − (Ri + 1)yi ln(yi ) + RD T ln T . ∂λ∂δ δ, n i=1 ˜ λ˜
ˆ BT in Eq. (9) can then be Now, following the same arguments with g (δ, λ) = δ and λ, respectively, in H ∗ (δ, λ), δˆ BT and λ obtained straightforwardly. 3.3. MCMC methods Based on the independent gamma priors λ ∼ Gamma(a, b) and δ ∼ Gamma(α, β), the posterior pdfs of λ and δ are given by
λ|δ, y ∼ Gamma a + D,
D −
δ
⋆
δ
(Ri + 1)yi + RD T + b
(10)
i=1
and
π (δ|λ, y ) ∝ δ
α+D−1
D ∏
δ−1
yi
D − exp −δβ − λ (Ri + 1)yδi + R⋆D T δ
i=1
.
i=1
Noting that the density function π (δ|λ, y ) is not known, but it follows from Lemma 1 of Kundu (2007) that
π (δ|y ) ∝ δ
α+D−1 −δβ
e
D ∏ i =1
δ−1
yi
−(a+D) D − δ ⋆ δ (Ri + 1)yi + RD T + b
(11)
i =1
can be easily identified to be a log-concave function. Therefore, we may adopt the Gibbs sampling technique of Kundu (2007) to generate samples from π (δ|y ) in (11) using the approach of Devroye (1984) and in turn to obtain the Bayes estimates of δ and λ. The detailed procedure can be described as follows:
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
457
Step 1. Set t = 1. Step 2. Generate δ (t ) from π (δ|y ) using the simple algorithm proposed by Devroye (1984): (a) Compute c = π(m|y ), where m is the mode of π (δ|y ). (b) Generate u from Uniform(0, 2), and v from Uniform(0, 1). (c) If u ≤ 1 then set w = u and d = v ; otherwise, set w = 1 − ln(u − 1) and d = v(u − 1). (d) Set w = m + w/c. (e) Accept w if d ≤ π(w|y )/c; otherwise, go to (d). Step 3. Generate λ(t ) from π(λ|δ (t ) , y ) as given in (10). Step 4. Set t = t + 1. Step 5. Repeat Steps 2–4, K times, and obtain (δ (i) , λ(i) ) for i = 1, . . . , K . Alternatively, we can use a Metropolis–Hastings step to generate random samples from the posterior density π (δ|y ). By experimentation, we observed that π (δ|y ) is similar to a normal distribution. The steps of this Metropolis–Hastings algorithm are carried out as follows: Step 1. Set t = 1. Step 2. Generate δ (t ) from π (δ|y ) using the Metropolis–Hastings algorithm with the proposal distribution q(δ) ≡ I (δ > ˆ Var(δ)) ˆ : 0)N (δ, (a) Let v = δ (t −1) . Here we set δ (0) ≡ δˆ . (b) Generate w from the proposal distribution q.
π(w|y )q(v)
(c) Let p(v, w) = min 1, π(v|y )q(w) . (d) Accept w with probability p(v, w) or accept v with probability 1 − p(v, w). Step 3. Generate λ(t ) from π (λ|δ (t ) , y ) as given in (10). Step 4. Set t = t + 1. Step 5. Repeat Steps 2–4, K times, and obtain (δ (i) , λ(i) ) for i = 1, . . . , K . the Bayes estimates of δ and λ with respect to squared error loss function, respectively, are δˆ BM = Eˆ (δ|y ) = ∑KHence, ∑ (i) δ / K , and λˆ BM = Eˆ (λ|y ) = Ki=1 λ(i) /K . i =1 4. Expected number of failures E (D) In some applications, it is of interest to have the average number of failures (or effective sample size) in the experiment for a particular life testing plan. For a Weibull distribution, we can show that the probability mass function of D for a prefixed stopping time T = exp(TL ) under adaptive Type-I PHCS is given by (see Appendix)
] [ TL − µ for d = 0, exp − n exp σ [ ] − d−1 TL − µ c exp −γ exp ci,d−1 (R1 + 1, . . . , Rd−1 + 1) d d+1 σ i =0 d ∑ T −µ P (D = d) = (Rj +1) exp Lσ 1−exp − j=d−i × for d = 1, . . . , m − 1, d ∑ (Rj +1) [ j=d−i ] − [ ] m TL − µ ai,m TL − µ ∗ ∗ 1 − exp −(γ − γ ) exp cm exp −γm exp i m σ γi − γm∗ σ i =1 where cd and γi are given in Eq. (2), γm∗ = n − m − ai , d =
d ∏
1
k=1
γk − γi
,
∑m−1 i=1
for d = m,
Ri ,
1 ≤ i ≤ d ≤ m,
(12)
k̸=i
and ci,r (α1 , . . . , αr ) =
i ∏
r −i+j
(−1)i
∑
j=1 k=r −i+1
αk
r −i r −i
∏∑ j=1 k=j
, αk
i = 1, 2, . . . , r ,
(13)
458
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
and P (D = m + d)
T −µ d−1 − m − exp −γj exp TLσ−µ exp −(γm∗ − d) exp L σ − c ( 1 , . . . , 1 ) a i,d−1 j,m cd′ cm i + 1 γj − (γm∗ − d) i=0 j=1 T −µ exp −(γm∗ − d + i + 1) exp L σ − exp −γj exp TLσ−µ − for d = 1, . . . , γm∗ − 1, γj − (γm∗ − d) − (i + 1) = T −µ ∗ γ− m m −1 − 1 − exp −γj exp L σ ∗ −1 (1, . . . , 1)aj,m c i ,γ m cγ′ ∗ cm m i+1 γj i = 0 j = 1 TL −µ TL −µ − exp −γ exp exp −( i + 1 ) exp j σ σ − for d = γm∗ , γj − (i + 1) where cd′ = γm∗ (γm∗ − 1) · · · (γm∗ − d + 1). Then the expected number of failures E (D) = without any difficulty.
∑m+γm∗ d=1
d P (D = d) can be computed
5. Simulation study The simulation study has been carried out to examine the performance of the MLEs and Bayes estimators — Lindley’s approximation (Lindley), Tierney and Kadane’s approximation (T&K), and two MCMC methods (which we abbreviate as MCMC1 and MCMC2), for n = 30 and m = 5, TL = ln T = −0.5 and 0.5, and parameters (δ, λ) = (1, 1) and (1.25, 0.5353) under three different sampling schemes: Scheme 1: R1 = · · · = Rm−1 = 0 and Rm = n − m; Scheme 2: R1 = n − m and R2 = · · · = Rm = 0; and Scheme 3: R1 = · · · = Rm−1 = 1 and Rm = n − 2m + 1. The Fortran and R programs and the details of our simulated results are available upon request. For each setting, the average biases and mean squared errors (MSEs) of MLEs and Bayes estimates based on 10 000 simulations when D > 0 are plotted in Figs. 1–6. Here the Bayes estimates are all computed by assuming that δ and λ have Gamma(α, β) and Gamma(a, b) priors, respectively, including the noninformative gamma priors, prior 0 : α = β = a = b = 0, and informative priors, prior k : α = β = a = b = 2k, where k = 1, . . . , 10, and a squared error loss function is used. The associated expected number of failures E (D) and the average number of observed failures based on 10 000 replications, details are not reported, are also computed. They agree closely to each other and this suggests that the Monte Carlo simulation study provides an alternative approach to access the number of failures occurring for the experiment under a particular setting. In terms of minimum biases and MSEs, the performance of the MCMC2 estimators for the parameter λ is generally best followed by the MCMC1 and T&K estimators and then the Lindley estimators and MLEs; in contrast, the ordering of the performance of the estimators (from best to worst) for the parameter δ is T&K, MCMC2, MCMC1, Lindley estimators and MLEs. It can be further observed that the MSEs of the Lindley estimators are usually larger than other Bayes estimators when the prior k gets larger. Compared to the MLEs, the MCMC2, MCMC1 and T&K estimators in most of the cases do greatly improve the efficiency of the estimates of δ and λ. On the whole, the MCMC2 and T&K estimators perform quite well for both parameters, whereas the Lindley estimators and MLEs perform rather poorly. Thus, it would seem reasonable to recommend the use of the MCMC2 and T&K estimators, especially the MCMC2 estimators, for estimating the unknown parameters δ and λ. Similar findings can be also observed for larger m, but the corresponding figures are not presented here for brevity. 6. Planning life test Due to limited available funding, it is very important to find the optimal censoring scheme in practice. It is often to choose the particular censoring scheme which provides more information of the unknown parameters among the class of possible schemes. Our criterion used to compare the possible adaptive progressive hybrid censoring schemes to determine the optimum one under the limited budget is based on the estimation precision of the logarithm of the pth (0 < p < 1) quantile of Weibull distribution, ln Tp = ln[− ln(1 − p)/λ]/δ , as was used by Kundu (2008) as well. Let C1 be the cost for inspecting an item, C2 be the salvage incurred by an unfailed item in the inspection, C3 be the cost per unit time used for life test. Then, the average total cost involved for a given adaptive progressive hybrid censoring scheme is Cb = C1 n − C2 (n − E (D)) + C3 T = (C1 − C2 )n + C2 E (D) + C3 T ,
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
0.1
0.2
MLE Lindley T&K MCMC1 MCMC2
-0.1
0.0
Bias(lambda)
0.1 0.0 -0.1
Bias(delta)
0.2
MLE Lindley T&K MCMC1 MCMC2
0
5
10
15
20
0
5
5
10
15
MSE(lambda)
MLE Lindley T&K MCMC1 MCMC2
0
10
15
20
Prior
0.00 0.05 0.10 0.15 0.20 0.25
0.00 0.05 0.10 0.15 0.20 0.25
Prior
MSE(delta)
459
20
MLE Lindley T&K MCMC1 MCMC2
0
5
Prior
10
15
20
15
20
15
20
Prior
0
Bias(lambda)
5
10
15
MLE Lindley T&K MCMC1 MCMC2
-0.2 0.0 0.2 0.4 0.6
MLE Lindley T&K MCMC1 MCMC2
-0.2 0.0 0.2 0.4 0.6
Bias(delta)
(a) δ = 1 and λ = 1.
20
0
5
MLE Lindley T&K MCMC1 MCMC2
0.4 0.2 0.0
0.2
0.4
MSE(lambda)
0.6
MLE Lindley T&K MCMC1 MCMC2
0.6
10 Prior
0.0
MSE(delta)
Prior
0
5
10
15
20
0
5
Prior
10 Prior
(b) δ = 1.25 and λ = 0.5353. Fig. 1. The biases and MSEs of the estimators for different priors under Scheme 1 when n = 30, m = 5, TL = −0.5.
where n is positive integer, and E (D) ≥ 1 and T > 0. Thus, the problem boils down as follows: for a given budget Cb , among the class of all possible (n, m, (R1 , . . . , Rm ), T )’s, which satisfy
(C1 − C2 )n + C2 E (D) + C3 T ≤ Cb ,
(14)
we want to find the best scheme (n , m , (R1 , . . . , Rm ) , T ), which maximizes the information measure ∗
I (P ) =
∗
1
[∫
Vposterior(P ) (ln Tp )dp
Edata
∗
]−1
∗
,
(15)
0
where Vposterior(P ) (ln Tp ) denotes the posterior variance of ln Tp for a given censoring scheme P = (R1 , . . . , Rm ) and the information measure in (15) is independent of p. Because of the restriction of the integer variables n, m, and (R1 , . . . , Rm ), the nonlinear mixed integer programming (see Grossmann, 2002) can be used to find the optimal solution. From (14) and E (D) ≥ 1, we have (C1 − C2 )n + C2 + C3 T ≤ Cb . C −C Since T > 0, we further have (C1 − C2 )n + C2 ≤ Cb , which implies n ≤ Cb −C2 . Thus, the upper bound of n is obtained. For a 1
2
given value of n, it is clear that C2 E (D) + C3 T ≤ Cb − (C1 − C2 )n. Since E (D) ≥ 1, we therefore obtain T ≤
Cb −(C1 −C2 )n−C2 , C3
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
1.0 0.5
20
5
10 Prior
15
20
2.0
MLE Lindley T&K MCMC1 MCMC2
0.0
0.0
0.5
1.0
1.5
2.0
MLE Lindley T&K MCMC1 MCMC2
0
1.5
15
1.0
10 Prior
MSE(lambda)
5
0.5
0
MSE(delta)
MLE Lindley T&K MCMC1 MCMC2
0.0
0.5 0.0
Bias(delta)
1.0
MLE Lindley T&K MCMC1 MCMC2
Bias(lambda)
460
0
5
10 Prior
15
20
0
5
10 Prior
15
20
1.5 1.0 0.5
Bias(lambda)
-0.5 20
5
10 Prior
15
20
MLE Lindley T&K MCMC1 MCMC2
0.0
0.0
0.5
1.0
1.5
MLE Lindley T&K MCMC1 MCMC2
0
1.5
15
1.0
10 Prior
MSE(lambda)
5
0.5
0
MSE(delta)
MLE Lindley T&K MCMC1 MCMC2
0.0
0.5
1.0
MLE Lindley T&K MCMC1 MCMC2
-0.5 0.0
Bias(delta)
1.5
(a) δ = 1 and λ = 1.
0
5
10 Prior
15
20
0
5
10 Prior
15
20
(b) δ = 1.25 and λ = 0.5353. Fig. 2. The biases and MSEs of the estimators for different priors under Scheme 2 when n = 30, m = 5, TL = −0.5.
the upper bound of T (which is denoted as Tn ). Since larger T yields larger information measure, we shall find the largest value of T , Tn,m,(R1 ,...,Rm ) , which is less than the upper bound Tn and also satisfies the constraint (14) for each given n, m and (R1 , . . . , Rm ). Therefore, the optimization can be determined in a finite number of steps as follows: Step 1. Compute the upper bound of the number of test units n˜ =
C b −C 2 C 1 −C 2
, where [y] is the greatest integer that is less than
or equal to y. Step 2. Set n = 2. Step 3. Compute the upper bound of the length of inspection interval Tn =
Cb −(C1 −C2 )n−C2 . C3
Step 4. Set m = 1. Step 5. Apply the bisection method to obtain Tn,m,(R1 ,...,Rm ) and compute the corresponding information measure, I (P )n,m,(R1 ,...,Rm ) , for all possible choices (n, m, (R1 , . . . , Rm )). Step 6. Set m = m + 1. If m ≤ n − 1 go to Step 5, else go to Step 7. Step 7. Record the particular scheme (n, mn , (R1 , . . . , Rm )n , Tn ), which has the largest information measure, I (P )n , among all possible choices of (n, m, (R1 , . . . , Rm ), Tn,m,(R1 ,...,Rm ) ) for 1 ≤ m ≤ n − 1 in Table 1.
0.3 0.2 0.1
Bias(lambda)
0.10
MSE(lambda)
MLE Lindley T&K MCMC1 MCMC2
0.20
0
5
10 Prior
15
20
MLE Lindley T&K MCMC1 MCMC2
0.00
0.00
MSE(delta)
20
0.30
15
0.20
10 Prior
MLE Lindley T&K MCMC1 MCMC2
0.10
5
0.30
0
461
-0.1 0.0
0.1
0.2
MLE Lindley T&K MCMC1 MCMC2
-0.1 0.0
Bias(delta)
0.3
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
0
5
10 Prior
15
20
0
5
10 Prior
15
20
10 Prior
15
20
10 Prior
15
20
Bias(lambda) 10 Prior
15
20
0
0.6
MLE Lindley T&K MCMC1 MCMC2
0.0
0.0
5
0.4
0.4
MSE(lambda)
MLE Lindley T&K MCMC1 MCMC2
0.6
MLE Lindley T&K MCMC1 MCMC2
0.8
5
0.2
MSE(delta)
0.8
0
-0.2 0.0 0.2 0.4 0.6 0.8
MLE Lindley T&K MCMC1 MCMC2
0.2
Bias(delta)
-0.2 0.0 0.2 0.4 0.6 0.8
(a) δ = 1 and λ = 1.
0
5
10 Prior
15
20
0
5
(b) δ = 1.25 and λ = 0.5353. Fig. 3. The biases and MSEs of the estimators for different priors under Scheme 3 when n = 30, m = 5, TL = −0.5.
Step 8. Set n = n + 1. If n ≤ n˜ go to Step 3, else go to Step 9. Step 9. Obtain the optimal scheme (n∗ , m∗ , (R1 , . . . , Rm )∗ , T ∗ ), which has the largest information measure, I (P )∗ , among all possible choices of (n, mn , (R1 , . . . , Rm )n , Tn ) for 1 ≤ n ≤ n˜ in Table 1. Note that it is not possible to compute Edata
1 0
Vposterior(P ) (ln Tp )dp analytically, so we conduct a Monte Carlo simulation
study along with the MCMC2 method discussed in Section 3.3 to approximate I (P ). More precisely, for fixed p, we first approximate Vposterior(P ) (ln Tp ) = Eposterior(P ) (ln Tp )2 − (Eposterior(P ) (ln Tp ))2 by using the MCMC2 method with K = 1000 to obtain the Bayes estimates of g1 (δ, λ) = ln Tp (δ, λ) = (ψ − ln λ)/δ and g2 (δ, λ) = (ln Tp (δ, λ))2 = [ψ 2 − 2ψ ln λ + (ln λ)2 ]/δ 2 , where ψ = ln(− ln(1 − p)), and then compute
1 0
Vposterior(P ) (ln Tp )dp. Upon repeating the same procedure 10 000 times, the average of
an approximation for I (P ) = Edata
1 0
Vposterior(P ) (ln Tp )dp .
1 0
Vposterior(P ) (ln Tp )dp’s shall be
15
0.06 0.02
20
MLE Lindley T&K MCMC1 MCMC2
0
5
10 Prior
15
0
MSE(lambda)
10 Prior
5
10 Prior
15
0.00 0.02 0.04 0.06 0.08 0.10
5
0.00 0.02 0.04 0.06 0.08 0.10
0
MSE(delta)
MLE Lindley T&K MCMC1 MCMC2
-0.02
Bias(lambda)
0.06 0.02 -0.02
Bias(delta)
MLE Lindley T&K MCMC1 MCMC2
0.10
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467 0.10
462
20
20
MLE Lindley T&K MCMC1 MCMC2
0
5
10 Prior
15
20
10 Prior
15
20
10 Prior
15
20
5
0.0 0.2 0.4 0.6
Bias(lambda) 10 Prior
15
20
0
5
MSE(lambda)
MLE Lindley T&K MCMC1 MCMC2
0
10 Prior
15
20
0.00 0.05 0.10 0.15 0.20
0.00 0.05 0.10 0.15 0.20
0
MSE(delta)
MLE Lindley T&K MCMC1 MCMC2
-0.4
0.0 0.2 0.4 0.6
MLE Lindley T&K MCMC1 MCMC2
-0.4
Bias(delta)
(a) δ = 1 and λ = 1.
5
MLE Lindley T&K MCMC1 MCMC2
0
5
(b) δ = 1.25 and λ = 0.5353. Fig. 4. The biases and MSEs of the estimators for different priors under Scheme 1 when n = 30, m = 5, TL = 0.5.
For illustrative purposes, let us consider the case Cb = 100, C1 = 12, C2 = 3, and C3 = 10 when (δ, λ) = (1, 1) and (1.25, 0.5353), respectively. Assuming that both δ and λ have independent Gamma(20, 20) prior distributions, we report the largest information measure I (P )n and the corresponding censoring scheme (n, mn , (R1 , . . . , Rm )n , Tn ) for 2 ≤ n ≤ n˜ = 10 in Table 1. Therefore, the optimal progressive censoring scheme (n∗ , m∗ , (R1 , . . . , Rm )∗ , T ∗ ) is (7,3, (0,0,4), 1.9107) for (δ, λ) = (1, 1), and (6,4, (0,0,0,2), 3.0147) for (δ, λ) = (1.25, 0.5353). Finally, we have also checked the effect of efficiency if we depart from the optimal scheme slightly. For instance, we depart from the optimal scheme (7, 3, (0, 0, 4), 1.9107) for (δ, λ) = (1, 1) slightly, say (0,1,3), which has an information measure 8.8763, then the relative efficiency is 8.8763/9.1181 = 0.9735; and at (1,0,3), which has an information measure 8.8445, the relative efficiency becomes 8.8445/9.1181 = 0.9700; they are near optimal. This shows that the efficiency does not change much if we depart from the optimal scheme only slightly. The same conclusion can also be made for (δ, λ) = (1.25, 0.5353). Acknowledgments We express our sincere thanks to the editor, associate editor, and the reviewers for their useful suggestions which led to an improvement over an earlier version of this manuscript. The first author thanks the National Science of Council of Taiwan (NSC Grant Number 99-2118-M-032-011-MY3) for funding this research.
Bias(lambda)
0
5
10 Prior
15
20
MLE Lindley T&K MCMC1 MCMC2
0.0
0.5
1.0
MLE Lindley T&K MCMC1 MCMC2
0.0
MSE(delta)
20
1.5
15
MLE Lindley T&K MCMC1 MCMC2
1.0
10 Prior
463
0.5
5
1.5
0
-0.2 0.0 0.2 0.4 0.6 0.8
MLE Lindley T&K MCMC1 MCMC2
MSE(lambda)
Bias(delta)
-0.2 0.0 0.2 0.4 0.6 0.8
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
0
5
10 Prior
15
20
0
5
10 Prior
15
20
10 Prior
15
20
10 Prior
15
20
10 Prior
15
0.2 0.4 0.6 0.8 1.0
20
0
0
MSE(lambda)
MLE Lindley T&K MCMC1 MCMC2
5
10 Prior
15
0.0 0.1 0.2 0.3 0.4 0.5 0.6
5
0.0 0.1 0.2 0.3 0.4 0.5 0.6
0
MSE(delta)
MLE Lindley T&K MCMC1 MCMC2
-0.2
Bias(lambda)
MLE Lindley T&K MCMC1 MCMC2
-0.2
Bias(delta)
0.2 0.4 0.6 0.8 1.0
(a) δ = 1 and λ = 1.
5
MLE Lindley T&K MCMC1 MCMC2
0
20
5
(b) δ = 1.25 and λ = 0.5353. Fig. 5. The biases and MSEs of the estimators for different priors under Scheme 2 when n = 30, m = 5, TL = 0.5.
Appendix. Probability Mass Function of D For extreme value distribution, the pdf of Xd:m:n , d = 1, . . . , m, is given by (see Balakrishnan and Aggarwala, 2000) fXd:m:n (x) = cd
d −
ai,d
i=1
1
σ
[ exp
x−µ
σ
− γi exp
x−µ
σ
]
,
−∞ < x < ∞,
where cd , γi , and ai,d are given in Eqs. (2) and (12). Let D be the number of failures up to time T and TL = ln T . Clearly, TL −µ Pr(D = 0) = Pr(X1:m:n > TL ) = exp −n exp . It then follows from the joint density of Xr :m:n and Xs:m:n for σ 1 ≤ r < s ≤ m (see Balakrishnan, 2007), fXr :m:n ,Xs:m:n (xr , xs ) = cs
r −1 s− r −1 − − i=0
k=0
ci,r −1 (R1 + 1, . . . , Rr −1 + 1)ck,s−r −1 (Rr +1 + 1, . . . , Rs−1 + 1)
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
10 Prior
15
0.10 0.05
Bias(lambda)
0.00
5
10 Prior
15
20
MLE Lindley T&K MCMC1 MCMC2
0.00
0.04
MSE(lambda)
MLE Lindley T&K MCMC1 MCMC2
0.08
0
0.00
MSE(delta)
20
0.12
5
0.12
0
MLE Lindley T&K MCMC1 MCMC2
0.08
0.05 0.00
Bias(delta)
0.10
MLE Lindley T&K MCMC1 MCMC2
0.04
464
0
5
10 Prior
15
20
0
5
10 Prior
15
20
10 Prior
15
20
10 Prior
15
20
10 Prior
15
0.0 0.2 0.4 0.6 0.8
20
0
5
MSE(lambda)
MLE Lindley T&K MCMC1 MCMC2
0
10 Prior
15
20
5
0.00 0.05 0.10 0.15 0.20
5
0.00 0.05 0.10 0.15 0.20
0
MSE(delta)
MLE Lindley T&K MCMC1 MCMC2
-0.4
Bias(lambda)
MLE Lindley T&K MCMC1 MCMC2
-0.4
Bias(delta)
0.0 0.2 0.4 0.6 0.8
(a) δ = 1 and λ = 1.
MLE Lindley T&K MCMC1 MCMC2
0
5
(b) δ = 1.25 and λ = 0.5353. Fig. 6. The biases and MSEs of the estimators for different priors under Scheme 3 when n = 30, m = 5, TL = 0.5. s−∑ 1−k
× f (xr )[1 − F (xr )]
(Rj +1)−1
j=r −i
f (xs )[1 − F (xs )]
n−
s−∑ 1−k
(Rj +1)−1
j=1
,
xr < xs ,
where ci,j is given in Eq. (13), to have P (D = d) = P (Xd:m:n < TL < Xd+1:m:n )
= cd+1
d−1 −
ci,d−1 (R1 + 1, . . . , Rd−1 + 1)
i =0
∫
∞
∫
TL
× TL
= cd+1
d−1 − i =0
d ∑
f (xd )[1 − F (xd )]
(Rj +1)−1
j=d−i
−∞
ci,d−1 (R1 + 1, . . . , Rd−1 + 1)
f (xd+1 )[1 − F (xd+1 )]γd+1 −1 dxd dxd+1
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
465
Table 1 The information measure I (P )n for (n, mn , (R1 , . . . , Rm )n , Tn ) when 2 ≤ n ≤ n˜ = 10. n
mn
(R1 , . . . , Rm )n
Tn
Tn
E ( D)
I (P )n
1 2 1 3 5 3 4 2 4
(1) (1,0) (3) (1,0,1) (0,0,0,0,1) (0,0,4) (1,0,0,3) (5,2) (0,6,0,0)
7.9000 7.0000 6.1000 5.2000 4.3000 3.4000 2.5000 1.6000 0.7000
7.6003 6.7006 5.2066 4.3150 2.8991 1.9107 1.2749 1.0532 0.3612
1.9990 1.9982 3.9781 3.9499 5.6696 5.9642 5.0837 2.8227 2.1293
4.2435 4.8165 7.3065 7.6797 9.0157 9.1181 8.3129 7.0320 5.3858
(1) (1,0) (0,0,1) (1,0,0,0) (0,0,0,2) (0,1,0,0,1) (0,2,0,0,1) (1,5,0) (0,6,1)
7.9000 7.0000 6.1000 5.2000 4.3000 3.4000 2.5000 1.6000 0.7000
7.6007 6.7014 5.2176 4.3394 3.0147 2.2762 1.6078 1.2083 0.5034
1.9977 1.9953 3.9413 3.8687 5.2843 4.7460 3.9741 2.3058 1.6552
3.6199 4.1332 6.0566 6.4624 7.5273 7.4703 6.7446 6.3858 4.4670
δ = 1 and λ = 1 2 3 4 5 6 7 8 9 10
δ = 1.25 and λ = 0.5353 2 3 4 5 6 7 8 9 10
1 2 3 4 4 5 5 3 3
×
d ∑ TL −µ 1 − exp − (Rj + 1) exp σ TL −µ exp −γ exp d + 1 j=d−i σ d ∑
(Rj + 1)
γd+1
j=d−i
] − [ d−1 TL − µ ci,d−1 (R1 + 1, . . . , Rd−1 + 1) = cd exp −γd+1 exp σ i=0 d ∑ TL −µ 1 − exp − (Rj + 1) exp σ j=d−i × , for d = 1, . . . , m − 1. d ∑ (Rj + 1) j=d−i
Given Xm:m:n = xm , Xm+1:n is now distributed as the first order statistic of a random sample of size γm∗ = n − m − from a truncated extreme value distribution with cumulative distribution function ] [ xm − µ x−µ ∗ ∗ + γm exp , xm < x < ∞. FXm+1:n (x|Xm:m:n = xm ) = 1 − exp −γm exp
σ
∑m−1 i=1
Ri
σ
Therefore, we have P (D = m) = P (Xm:m:n < TL ≤ Xm+1:n )
∫
∞
P (x < TL ≤ Xm+1:n |Xm:m:n = x)fXm:m:n (x)dx
= −∞
∫
TL
=
[1 − FXm+1:n (TL |Xm:m:n = x)]fXm:m:n (x)dx
−∞
[ ] − [ ] m TL − µ TL − µ ai , m ∗ = cm exp −γm∗ exp 1 − exp −(γ − γ ) exp . i m σ γi − γm∗ σ i=1 In a similar fashion, Xm+r :n and Xm+s:n (for 1 ≤ r < s ≤ γm∗ ) are distributed as the r th and sth order statistics of a random sample of size γm∗ from a truncated exponential distribution with joint density function (see Balakrishnan et al., 2002) fXr :n ,Xs:n |Xm:m:n (xr , xs |xm ) = cs′
r −1 s− r −1 − − i=0
[ ×
f (xr ) 1 − F (xm )
][
f (xs ) 1 − F ( xm )
ci,r −1 (1, . . . , 1)ck,s−r −1 (1, . . . , 1)
k=0
][
1 − F (xr ) 1 − F (xm )
]s−r −1−k+i [
1 − F (xs ) 1 − F (xm )
]γm∗ −s+k
,
xm < xr < xs < ∞,
466
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
where cs′ = γm∗ (γm∗ − 1) · · · (γm∗ − s + 1). Then, for d = 1, . . . , γm∗ − 1, we have P (Xm+d:n < x < Xm+d+1:n |Xm:m:n = xm ) ′
= cd+1
d−1 −
] ∗ xm − µ ci,d−1 (1, . . . , 1) exp γm − d + i + 1 exp σ [
i =0 ∞∫ x
∫ ×
f (y)f (z )[1 − F (y)]i [1 − F (z )]γm −d−1 dydz ∗
xm
x
] ∗ xm − µ = cd+1 ci,d−1 (1, . . . , 1) exp γm − d + i + 1 exp σ i =0 xm −µ x−µ x−µ exp −(i + 1) exp − exp −(i + 1) exp σ exp −(γm∗ − d) exp σ σ × i+1 γm∗ − d [ ] xm − µ x−µ = cd′ exp (γm∗ − d) exp − exp σ σ x −µ d−1 − 1 − exp (i + 1) exp mσ − exp x−µ σ . × ci,d−1 (1, . . . , 1) i+1 i=0 d−1 −
′
[
It leads to P (D = m + d) = P (Xm+d:n < TL ≤ Xm+d+1:n ) ∞
∫
P (x < Xm+d:n < TL ≤ Xm+d+1:n |Xm:m:n = x)fXm:m:n (x)dx
= −∞
= cd′ cm
TL −µ TL −µ ∗ d−1 − m − exp −γ exp exp −(γ − d ) exp − j m σ σ ci,d−1 (1, . . . , 1)aj,m
γj − (γm∗ − d) T −µ exp −(γm∗ − d + i + 1) exp L σ − exp −γj exp TLσ−µ − , for d = 1, . . . , γm∗ − 1. γj − (γm∗ − d) − (i + 1) i+1
i=0 j=1
Finally, we use the representation of Balakrishnan et al. (2002) for the pdf of Xγm∗ :n to have P (D = m + γm∗ ) = P (Xγm∗ :n < TL ) ∞
∫ =
−∞ TL
∫
P (x < Xγm∗ :n < TL |Xm:m:n = x)fXm:m:n (x)dx TL
∫
= −∞
= cγ′ m∗ cm
×
x
cγ′ ∗
∗ γ− m −1
m
x−µ
]
σ
i =0
∗ γ− m m −1 −
f (y)[1 − F (y)]i dyfXm:m:n (x)dx
ci,γm∗ −1 (1, . . . , 1)aj,m
i+1 i=0 j=1 1 − exp −γj exp TLσ−µ
[
ci,γm∗ −1 (1, . . . , 1) exp (i + 1) exp
γj
exp −(i + 1) exp
−
TL −µ
σ
− exp −γj exp TLσ−µ
γj − (i + 1)
.
References Abramowitz, M., Stegun, I.A. (Eds.), 1964. Handbook of Mathematical Functions, with Formulas, Graphs and Mathematical Tables. Dover, New York. Balakrishnan, N., 2007. Progressive censoring methodology: an appraisal (with discussions). Test 16, 211–296. Balakrishnan, N., Aggarwala, R., 2000. Progressive Censoring: Theory, Methods and Applications. Birkhäuser, Boston. Balakrishnan, N., Childs, A., Chandrasekar, B., 2002. An efficient computational method for moments of order statistics under progressive censoring. Statistics and Probability Letters 60, 359–365. Banerjee, A., Kundu, D., 2008. Inference based on type-II hybrid censored data from a Weibull distribution. IEEE Transactions on Reliability 57, 369–378. Childs, A., Chandrasekar, B., Balakrishnan, N., 2008. Exact likelihood inference for an exponential parameter under progressive hybrid censoring schemes. In: Vonta, F., Nikulin, M., Limnios, N., Huber-Carol, C. (Eds.), Statistical Models and Methods for Biomedical and Technical Systems. Birkhäuser, pp. 323–334. Devroye, L., 1984. A simple algorithm for generating random variates with a log-concave density. Computing 33, 247–257.
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
467
Grossmann, I.E., 2002. Review of nonlinear mixed-integer and disjunctive programming techniques. Optimization and Engineering 3, 227–252. Kundu, D., 2007. On hybrid censoring Weibull distribution. Journal of Statistical Planning and Inference 137, 2127–2142. Kundu, D., 2008. Bayesian inference and life testing plan for the Weibull distribution in presence of progressive censoring. Technometrics 50, 144–154. Kundu, D., Joarder, A., 2006. Analysis of type-II progressively hybrid censored data. Computational Statistics & Data Analysis 50, 2509–2528. Lin, C.T., Huang, Y.L., 2011. On progressive hybrid censored exponential distribution. Journal of Statistical Computation and Simulation, First published on: 21 June 2011 (iFirst). Lin, C.T., Ng, H.K.T., Chan, P.S., 2009. Statistical inference of type-II progressively hybrid censored data with Weibull lifetimes. Communications in Statistics— Theory and Methods 38, 1710–1729. Lin, C.T., Wu, S.J.S., Balakrishnan, N., 2006. Inference for log-gamma distribution based on progressively type-II censored data. Communications in Statistics— Theory and Methods 35, 1271–1292. Lindley, D.V., 1980. Approximate Bayesian method. Trabajos de Estadistica 31, 223–237. Ng, H.K.T., Chan, P.S., Balakrishnan, N., 2002. Estimation of parameters from progressively censored data using EM algorithm. Computational Statistics & Data Analysis 39, 371–386. Ng, H.K.T., Kundu, D., Chan, P.S., 2009. Statistical analysis of exponential lifetimes under an adaptive type-II progressive censoring scheme. Naval Research Logistics 56, 687–698. Tierney, L., Kadane, J.B., 1986. Accurate approximations for posterior moments and marginal densities. Journal of the American Statistical Association 81, 82–86.
Contents lists available at SciVerse ScienceDirect
Computational Statistics and Data Analysis journal homepage: www.elsevier.com/locate/csda
Inference for the Weibull distribution with progressive hybrid censoring Chien-Tai Lin ∗ , Cheng-Chieh Chou, Yen-Lung Huang Department of Mathematics, Tamkang University, Tamsui 25137, Taiwan
article
info
Article history: Received 4 April 2010 Received in revised form 25 July 2011 Accepted 3 September 2011 Available online 16 September 2011 Keywords: Gibbs sampling Lindley’s approximation Markov Chain Monte Carlo method Metropolis–Hastings algorithm Tierney–Kadane’s approximation
abstract Recently, progressive hybrid censoring schemes have become quite popular in life-testing and reliability studies. In this paper, we investigate the maximum likelihood estimation and Bayesian estimation for a two-parameter Weibull distribution based on adaptive Type-I progressively hybrid censored data. The Bayes estimates of the unknown parameters are obtained by using the approximation forms of Lindley (1980) and Tierney and Kadane (1986) as well as two Markov Chain Monte Carlo methods under the assumption of gamma priors. Computational formulae for the expected number of failures is provided and it can be used to determine the optimal adaptive Type-I progressive hybrid censoring schemes under a pre-determined budget of experiment. © 2011 Elsevier B.V. All rights reserved.
1. Introduction A new progressive hybrid censoring scheme—adaptive Type-I progressive hybrid censoring scheme (adaptive Type-I PHCS), introduced by Lin and Huang (2011), is a mixture of Type-I and Type-II progressive censoring schemes, and it can be described as follows. Suppose n identical units are placed on test with progressive censoring scheme (R1 , R2 , . . . , Rm ), 1 ≤ m ≤ n, and the experiment is terminated at time T , where T ∈ (0, ∞) and integers Ri ’s are fixed in advance. At the time of the first failure Y1:m:n , R1 of the remaining units are randomly removed. Similarly, at the time of the second failure Y2:m:n , R2 of the remaining units are randomly removed and so on. Let D denote the number of failures that occur before time T . If the mth failure Ym:m:n occurs before time T , the process will not stop at the time Ym:m:n , but continue to observe ∑D failures (without any further withdrawals) up to time T . Then, at time T all the remaining R⋆D = n − D − i=1 Ri units are removed and the experiment is terminated. Note that the progressive censoring scheme in this case will become (R1 , R2 , . . . , Rm , Rm+1 , . . . , RD ), where Rm = Rm+1 = · · · = RD = 0. On the other hand, the process when Ym:m:n > T will have a progressive censoring scheme as (R1 , R2 , . . . , RD ). The work on progressive hybrid censoring has become quite popular in life-testing and reliability studies. Kundu and Joarder (2006) and Childs et al. (2008) have considered a Type-I PHCS in the context of life-testing experiment is terminated at time min{Ym:m:n , T }. Kundu and Joarder (2006) discussed the maximum likelihood and Bayesian estimation for the exponential distribution, whereas Childs et al. (2008) dealt with the derivation of the exact distribution of the maximum likelihood estimator (MLE) of the unknown parameter. In the same paper, Childs et al. (2008) has also proposed another Type-II PHCS that would terminate the experiment at the random time max{Ym:m:n , T }, and derived the corresponding exact distribution of MLE. Note that on the case Ym:m:n < T , the process under this censoring scheme will be performed exactly the same as our proposed scheme. Ng et al. (2009) and Lin et al. (2009) have investigated point and interval estimation for exponential and Weibull lifetimes under an adaptive Type-II PHCS. When Ym:m:n < T , adaptive Type-II PHCS is equivalent to Type-I PHCS. In contrast, no more units will be removed whenever the experimental time passes time T , and all the remaining
∗
Corresponding author. Tel.: +886 2 26263976; fax: +886 2 26209916. E-mail address: [email protected] (C.-T. Lin).
0167-9473/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.csda.2011.09.002
452
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
Rm = n − R1 −· · ·− RD − m units are removed from the experiment immediately after the mth failure has occurred. Recently, Lin and Huang (2011) studied point and interval estimation for exponential lifetimes under an adaptive Type-I PHCS, and investigated Bayesian sampling plans under different progressive hybrid censoring schemes when a general loss function is used. Upon applying the Lindley’s 1980 approximation and Gibbs sampling procedure, Kundu (2007, 2008) and Banerjee and Kundu (2008) have studied Bayesian inference of parameters for a Weibull distribution based on Type-I hybrid, Type-II hybrid, or Type-II progressive censored sample. The Lindley’s approximation is generally accurate enough in these works, but, as Lindley has pointed out, the required evaluation of third derivatives of the posterior can be rather tedious—in particular, in problems with several parameters. Thus, in this paper, an easily computable approximation of Tierney and Kadane (1986) is considered to approximate the Bayes estimates of the parameters for a Weibull distribution based on an adaptive Type-I progressive hybrid censored sample. Note that the marginal posterior density function of the shape parameter cannot be obtained analytically. For this reason, Kundu (2007) has also used a Gibbs sampling technique along with the approach for generating random variates from a log-concave density (see Devroye, 1984) to construct the Bayes estimates and credible intervals for the unknown parameters. In order to make a comparison to this procedure, the Metropolis–Hastings algorithm with normal proposal distribution is then employed. The rest of this article is organized as follows. In Section 2, we derive the MLEs of two unknown parameters of Weibull distribution and the corresponding asymptotic variances and covariance. In Section 3, we cover the Bayesian analysis based on the approximation forms of Lindley (1980) and Tierney and Kadane (1986) as well as two Markov Chain Monte Carlo (MCMC) methods. In Section 4, we give the computational formulae for the expected number of failures, which will be used to find the optimal adaptive Type-I PHCS in Section 6. Numerical comparisons of the MLEs and the Bayes estimates are presented in Section 5. In Section 6, we propose a procedure to determine the optimality adaptive Type-I PHCS when a cost constraint given in Eq. (14), which includes the sampling cost, the time-consuming cost, and the salvage, is used. 2. Maximum likelihood estimation In this section we discuss the MLEs of the parameters µ and σ based on the data observed under the proposed censoring scheme and derive their asymptotic variance–covariance matrix by using the results of Ng et al. (2002) and Lin et al. (2006). 2.1. Maximum likelihood estimators Consider lifetime of units in the life test Yi , i = 1, . . . , n, to be independent and identically distributed as Weibull with probability density function (pdf) δ
fY (y; δ, λ) = λδ yδ−1 e−λy ,
y > 0,
where δ > 0 is the shape parameter and λ > 0 is the scale parameter. Note that in the derivation of MLEs and expected number of failures in Section 4, instead of working with the parametric model for Yi , it is often more convenient to work with the equivalent model for the log-lifetimes Xi = ln Yi for i = 1, 2, . . . , n, which follows the extreme value distribution with parameters µ = − ln λ/δ and σ = 1/δ . Now, suppose that an adaptive Type-I progressive hybrid censored extreme value sample with a given adaptive Type-I PHCS (R1 , . . . , Rm ) is observed. We denote the resulting failure times by X1:m:n , X2:m:n , . . . , Xm:m:n , Xm+1:n , . . . , XD:n when Xm:m:n < TL , and X1:m:n , X2:m:n , . . . , XD:m:n when Xm:m:n > TL , where TL ≡ ln T . Then the likelihood function of the given sample is given by L(µ, σ ) = cD
D ∏
⋆
fX (xi )[1 − FX (xi )]Ri [1 − FX (TL )]RD ,
(1)
i=1
where fX (x; µ, σ ) =
1
σ
[ exp
x−µ
FX (x; µ, σ ) = 1 − exp − exp D ∏
− exp
σ [
cD =
x−µ
m −
,
]
,
−∞ < x < ∞,
−∞ < x < ∞,
Rj ,
(2)
j =i
i=1
and R⋆D = n − D −
σ ]
σ
γi with γi = m − i + 1 +
x−µ
∑D
i=1
Ri with Rm = Rm+1 = · · · = RD = 0 if D ≥ m. Here we use the usual conventions that
∏0
j =1
αj ≡ 1
and j=i αj ≡ 0. Also, for simplicity of notation, we use xi instead of xi:m:n or xi:n , and yi instead of yi:m:n or yi:n in the remaining discussion.
∑i−1
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
453
From (1), the log-likelihood function given D ≥ 1 is ln L(µ, σ ) = ln(cD ) − D ln σ +
D −
ξi −
D −
i=1
(Ri + 1) exp(ξi ) − R⋆D exp(tL ),
i=1
where
ξi = (xi − µ)/σ and tL = (TL − µ)/σ .
(3)
Thus, the MLE of σ , denoted as σˆ , can be obtained by D ∑
D ∑
xi
i=1
D
=
(Ri + 1)xi exp
xi
+ R⋆D TL exp
xi
+ R⋆D exp
σ
i =1 D ∑
(Ri + 1) exp
i=1
σ
TL
σ
− σ.
(4)
TL
σ
Since Eq. (4) cannot be solved analytically, some numerical methods such as a bisection algorithm or other iterative procedures are needed. Alternatively, we can follow the work of Lin et al. (2009) to obtain the approximate maximum likelihood estimates of the parameters, and use them as the initial values in the iterative procedure for solving Eq. (4). After obtaining σˆ , the MLE of µ, denoted as µ ˆ , can be easily computed as
D ∑
(Ri + 1) exp
i =1 µ ˆ = σˆ ln
xi σˆ
+ R⋆D exp
D
TL σˆ .
Note that the MLEs do not exist when D = 0. 2.2. Asymptotic variances and covariance Using the missing-information principle, we can compute the observed Fisher information matrix. For θ = (µ, σ )′ , we define X , Z and W to be the observed, missing and complete data, and IX (θ), IW |X (θ) and IW (θ) to be the corresponding Fisher information matrix, respectively. Then IX (θ) = IW (θ) − IW |X (θ).
(5)
The complete information matrix for the extreme value distribution (see Ng et al., 2002) is IW (θ) =
n
1−γ c2
1 1−γ
σ2
,
where γ = 0.577215665 is the Euler’s constant and c 2 is π 2 /6 + (1 − γ 2 ) = 1.823680661. Now IW |X (θ) =
D −
(j)
∗ Rj IW |X (θ) + R∗D IW |X (θ),
(6)
j=1
(j)
∗ where IW |X (θ) and IW |X (θ) are the information matrix of a single observation for the truncated extreme value distribution with left truncation at xj and TL , respectively. The logarithm of the pdf of the truncated extreme value distribution with left truncation at xj is
ln fZj (zj |Zj > xj , µ, σ ) = − log σ + exp(ξj ) +
zj − µ
− exp
σ
zj − µ
σ
,
where ξj is as defined in Eq. (3). Thus, the three second partial derivatives with respect to µ and σ are
∂ 2 ln fZj ∂µ2 ∂ 2 ln fZj ∂µ∂σ ∂ ln fZj
= =
2
∂σ 2
=
1
[
exp(ξj ) − exp
σ2 1
zj − µ
]
σ
[
exp(ξj ) + ξj exp(ξj ) − exp
σ2 1
,
σ2
ξj2 exp(ξj ) + 2ξj exp(ξj ) −
−2
zj − µ
σ
exp
zj − µ
σ
zj − µ
σ
zj − µ
+2
−
zj − µ
σ
2
σ
exp
zj − µ
σ
exp
zj − µ
σ
+1 .
zj − µ
σ
] +1 ,
454
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
Taking the negative of the expected value of these three second partial derivatives with the results (see Ng et al., 2002 and Lin et al., 2006)
[ E
Zj − µ
σ
[ E
exp
] ∞ − exp[(p + 1)ξj ] Zj > xj , µ, σ = ψ(1) exp(eξj ) − [ψ(1) + ξj − ψ(p + 2)] , Γ (p + 2) p=0
Zj − µ
] ∞ − exp[(p + 2)ξj ] Zj > xj , µ, σ = exp(eξj ) − , Γ (p + 3)
σ p=0 ] ∞ − Zj − µ Zj − µ exp[(p + 2)ξj ] Zj > xj , µ, σ = ψ(2) exp(eξj ) − [ψ(2) + ξj − ψ(p + 3)] E exp , σ σ Γ (p + 3) p=0 2 Zj − µ Zj − µ E exp Zj > xj , µ, σ = [ψ ′ (2) + ψ 2 (2)] exp(eξj ) − {[ξj + ψ(2) − ψ(p + 3)]2 σ σ ∞ − exp[(p + 2)ξj ] , + ψ ′ (2) − ψ ′ (p + 3)} Γ (p + 3) p=0 [
where ψ(t ) =
d dt
ln Γ (t ) and ψ ′ (t ) =
d dt
ψ(t ) are the digamma and trigamma functions (see Abramowitz and Stegun,
(j) 1964), the missing information matrix IW |X (θ) can be easily computed. Replacing ξj with tL defined in Eq. (3), the information ∗ IW |X (θ) can also be computed with the same arguments. Then, from (5) and (6), the asymptotic variance–covariance matrix
of the maximum likelihood estimate θ can be obtained by inverting the observed information matrix IX ( θ). 3. Bayesian estimation In this section, we consider the Bayesian estimation under the assumption that the prior distribution of λ is Gamma(a, b) with pdf given by
π1 (λ|a, b) =
ba
λa−1 e−λb , λ > 0, Γ (a) where a > 0 and b > 0; and the prior of δ, π2 (δ), has the support on (0, ∞), and it is independent of the prior of λ. Thus, the joint posterior density function of δ and λ is given by π (δ, λ)ℓ(y |δ, λ) , π (δ, λ|y ) ∝ ∞ ∞ π (δ, λ)ℓ(y |δ, λ)dδ dλ 0 0 where π (δ, λ) = π1 (λ|a, b)π2 (δ) is the joint prior density function and
ℓ(y |δ, λ) = cD δ λ
D D
D ∏
δ−1
yi
D − exp −λ (Ri + 1)yδi + R⋆D T δ
i =1
i=1
is the likelihood function based on the adaptive Type-I progressively hybrid censored sample y = (y1 , . . . , yD ). Then, under a squared error loss, the Bayes estimate of any function of δ and λ, say g (δ, λ), is g (δ, λ)eQ (δ,λ) dδ dλ 0∞ ∞ , (7) π (δ, λ)ℓ(y |δ, λ)dδ dλ eQ (δ,λ) dδ dλ 0 0 0 0 where Q (δ, λ) = ln π (δ, λ) + ln ℓ(y |δ, λ) ≡ ρ(δ, λ) + L(δ, λ). Note that Eq. (7) cannot be obtained analytically even when π2 (δ) is known. Therefore, we adopt two approximations—Lindley’s approximation and Tierney and Kadane’s
∞∞
E (g (δ, λ)|y ) =
0
0
g (δ, λ)π (δ, λ)ℓ(y |δ, λ)dδ dλ
∞∞
∞∞
=
0
approximation, as well as two MCMC methods for approximating (7). 3.1. Lindley’s approximation Lindley (1980) has proposed approximations for the ratio in (7). The basic idea is to expand L(δ, λ) and ρ(δ, λ) of (7) into a Taylor series expansion about the MLE of (δ, λ), or Q (δ, λ) of (7) about the posterior mode of (δ, λ). The integrations all involve the moments that capture the first-order error terms of the normal approximation. This approximation can be viewed as an application of the Laplace method for integrals, and it will produce reasonable results as long as the posterior is unimodal or at least dominated by a single mode. For the two-parameter case, (δ, λ) = (λ1 , λ2 ), Lindley’s approximation of (7) is the form g (λ1 , λ2 ) +
1− 2
gij σij +
−
Uj ρj +
1
1
2
2
1 2
L30 σ11 U1 +
+ L12 (σ22 U1 + 2σ12 U2 ) + L03 σ22 U2 ,
1 2
L21 (2σ12 U1 + σ11 U2 ) (8)
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
455
where Uk = j gj σkj with σij being the (i, j)th elements of the inverse of the Fisher information matrix, and each suffix of variable g , ρ , or L denotes differentiation once with respect to the variable having that suffix. Evaluating all the expressions in (8) at the MLEs of λ1 and λ2 produces the approximation gˆBL to (7). In our case, with the specification of δ ∼ Gamma(α, β), we have
∑
D −
L(δ, λ) = ln(cD ) + D ln δ + D ln λ + (δ − 1)
D − δ ⋆ δ ln(yi ) − λ (Ri + 1)yi + RD T .
i=1
i=1
The MLEs of δ and λ are the solutions of the equations L10 =
L01 =
D
D −
ln(yi ) − λ
D −
δ
+
D
D − δ ⋆ δ − (Ri + 1)yi + RD T = 0.
λ
δ
(Ri + 1)yi ln(yi ) + RD T ln T
δ
i=1
⋆
= 0,
i=1
i =1
Now, to apply Lindley’s form in (8), we first obtain L20 = −
L11
D − ˆ ˆ δ 2 ⋆ δ 2 − λˆ (Ri + 1)yi (ln yi ) + RD T (ln T ) ,
D
δˆ 2 i=1 D − δˆ ⋆ δˆ =− (Ri + 1)yi ln(yi ) + RD T ln T ,
L02 = −
i=1
L30 =
L21
2D
δˆ 3
D − ˆδ ˆ (Ri + 1)yi (ln yi )3 + R⋆D T δ (ln T )3 , − λˆ
D
λˆ 2
,
L03 =
i=1
D − δˆ 2 ⋆ δˆ 2 =− (Ri + 1)yi (ln yi ) + RD T (ln T ) ,
2D
λˆ 3
,
L12 = 0;
i=1
and W
σ11 =
SW − V 2
,
S
σ22 =
SW − V 2
,
σ12 = σ21 =
−V SW − V 2
,
where S=
V =
D
δˆ 2
+ λˆ
D −
D − δˆ 2 ⋆ δˆ 2 (Ri + 1)yi (ln yi ) + RD T (ln T ) , i=1
(Ri + 1)yδi ln yi + R⋆D T δ ln T , ˆ
ˆ
W =
i=1
D
λˆ 2
.
For the posterior mean of δ , we have g (δ, λ) = δ and hence g1 = 1, g2 = 0, gij = 0 for i, j = 1, 2. Thus, 1
3
1
2
2
2
2 δˆBL = δˆ + ρ1 σ11 + ρ2 σ12 + L30 σ11 + L21 σ11 σ12 + L03 σ12 σ22 ,
where ρ1 =
α−1 δˆ
− β and ρ2 =
a−1
λˆ
− b. Similarly, for the posterior mean of λ, we have g (δ, λ) = λ and
1
1
1
2
2
2
2 2 + σ11 σ22 ) + L03 σ22 . λˆ BL = λˆ + ρ1 σ12 + ρ2 σ22 + L30 σ11 σ12 + L21 (2σ12
3.2. Tierney and Kadane’s approximation Setting H (δ, λ) = Q (δ, λ)/n and H ∗ (δ, λ) = [ln g (δ, λ) + Q (δ, λ)]/n, Tierney and Kadane (1986) reexpressed the expression in (7) for the posterior mean of g as
∞∞ E (g (δ, λ)|y ) = 0∞ 0∞ 0
0
∗ enH (δ,λ) dδ dλ
enH (δ,λ) dδ dλ
,
and then applied Laplace’s method to produce the approximation gˆBT (δ, λ) =
[
det Σ ∗ det Σ
]1/2
˜ ∗ ) − H (δ˜ , λ) ˜ exp n H ∗ (δ˜ ∗ , λ
(9)
456
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
˜ ∗ ) and (δ, ˜ λ) ˜ maximize H ∗ (δ, λ) and H (δ, λ), respectively, and Σ ∗ and Σ are the negatives of to E (g (δ, λ)|y ), where (δ˜ ∗ , λ ∗ ˜ ∗ ) and (δ, ˜ λ) ˜ , respectively. the inverse Hessians of H (δ, λ) and H (δ, λ) at (δ˜ ∗ , λ In our case, we have H (δ, λ) =
1
k + (α + D − 1) ln δ + (a + D − 1) ln λ − δβ + (δ − 1)
n
D −
ln(yi ) − λ
i=1
D −
δ
⋆
δ
(Ri + 1)yi + RD T + b
,
i=1
˜ can be obtained by solving the following two equations where k is a constant; therefore, (δ˜ , λ) D D − − α+D−1 −β + ln(yi ) − λ (Ri + 1)yδi ln(yi ) + R⋆D T δ ln T = 0, δ i =1 i =1 D − ∂ H (δ, λ) 1 a+D−1 δ ⋆ δ = − (Ri + 1)yi + RD T + b = 0, ∂λ n λ i=1 ∂ H (δ, λ) 1 = ∂δ n
and, from the second derivatives of H (δ, λ),
D − ∂ 2 H (δ, λ) 1 α+D−1 δ 2 ⋆ δ 2 = −λ − (Ri + 1)yi (ln yi ) + RD T (ln T ) , ∂δ 2 n δ2 i =1 ∂ 2 H (δ, λ) 1 a+D−1 = − , ∂λ2 n λ2 D − ∂ 2 H (δ, λ) 1 = − (Ri + 1)yδi ln(yi ) + R⋆D T δ ln T , ∂λ∂δ n i=1 ˜ λ) ˜ is given by the determinant of the negative of the inverse Hessian of H (δ, λ) at (δ, 2 −1 det Σ = (H11 H22 − H12 ) ,
where H11
H22
H12
D − 1 α+D−1 ∂ 2 H (δ, λ) ˜δ ˜ = (Ri + 1)yi (ln yi )2 + R⋆D T δ (ln T )2 , + λ˜ =− ∂δ 2 δ, n δ˜ 2 ˜ λ˜ i=1 1 a+D−1 ∂ 2 H (δ, λ) = =− , ∂λ2 δ, n λ˜ 2 ˜ λ˜ D 1 − ∂ 2 H (δ, λ) δ˜ ⋆ δ˜ = = − (Ri + 1)yi ln(yi ) + RD T ln T . ∂λ∂δ δ, n i=1 ˜ λ˜
ˆ BT in Eq. (9) can then be Now, following the same arguments with g (δ, λ) = δ and λ, respectively, in H ∗ (δ, λ), δˆ BT and λ obtained straightforwardly. 3.3. MCMC methods Based on the independent gamma priors λ ∼ Gamma(a, b) and δ ∼ Gamma(α, β), the posterior pdfs of λ and δ are given by
λ|δ, y ∼ Gamma a + D,
D −
δ
⋆
δ
(Ri + 1)yi + RD T + b
(10)
i=1
and
π (δ|λ, y ) ∝ δ
α+D−1
D ∏
δ−1
yi
D − exp −δβ − λ (Ri + 1)yδi + R⋆D T δ
i=1
.
i=1
Noting that the density function π (δ|λ, y ) is not known, but it follows from Lemma 1 of Kundu (2007) that
π (δ|y ) ∝ δ
α+D−1 −δβ
e
D ∏ i =1
δ−1
yi
−(a+D) D − δ ⋆ δ (Ri + 1)yi + RD T + b
(11)
i =1
can be easily identified to be a log-concave function. Therefore, we may adopt the Gibbs sampling technique of Kundu (2007) to generate samples from π (δ|y ) in (11) using the approach of Devroye (1984) and in turn to obtain the Bayes estimates of δ and λ. The detailed procedure can be described as follows:
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
457
Step 1. Set t = 1. Step 2. Generate δ (t ) from π (δ|y ) using the simple algorithm proposed by Devroye (1984): (a) Compute c = π(m|y ), where m is the mode of π (δ|y ). (b) Generate u from Uniform(0, 2), and v from Uniform(0, 1). (c) If u ≤ 1 then set w = u and d = v ; otherwise, set w = 1 − ln(u − 1) and d = v(u − 1). (d) Set w = m + w/c. (e) Accept w if d ≤ π(w|y )/c; otherwise, go to (d). Step 3. Generate λ(t ) from π(λ|δ (t ) , y ) as given in (10). Step 4. Set t = t + 1. Step 5. Repeat Steps 2–4, K times, and obtain (δ (i) , λ(i) ) for i = 1, . . . , K . Alternatively, we can use a Metropolis–Hastings step to generate random samples from the posterior density π (δ|y ). By experimentation, we observed that π (δ|y ) is similar to a normal distribution. The steps of this Metropolis–Hastings algorithm are carried out as follows: Step 1. Set t = 1. Step 2. Generate δ (t ) from π (δ|y ) using the Metropolis–Hastings algorithm with the proposal distribution q(δ) ≡ I (δ > ˆ Var(δ)) ˆ : 0)N (δ, (a) Let v = δ (t −1) . Here we set δ (0) ≡ δˆ . (b) Generate w from the proposal distribution q.
π(w|y )q(v)
(c) Let p(v, w) = min 1, π(v|y )q(w) . (d) Accept w with probability p(v, w) or accept v with probability 1 − p(v, w). Step 3. Generate λ(t ) from π (λ|δ (t ) , y ) as given in (10). Step 4. Set t = t + 1. Step 5. Repeat Steps 2–4, K times, and obtain (δ (i) , λ(i) ) for i = 1, . . . , K . the Bayes estimates of δ and λ with respect to squared error loss function, respectively, are δˆ BM = Eˆ (δ|y ) = ∑KHence, ∑ (i) δ / K , and λˆ BM = Eˆ (λ|y ) = Ki=1 λ(i) /K . i =1 4. Expected number of failures E (D) In some applications, it is of interest to have the average number of failures (or effective sample size) in the experiment for a particular life testing plan. For a Weibull distribution, we can show that the probability mass function of D for a prefixed stopping time T = exp(TL ) under adaptive Type-I PHCS is given by (see Appendix)
] [ TL − µ for d = 0, exp − n exp σ [ ] − d−1 TL − µ c exp −γ exp ci,d−1 (R1 + 1, . . . , Rd−1 + 1) d d+1 σ i =0 d ∑ T −µ P (D = d) = (Rj +1) exp Lσ 1−exp − j=d−i × for d = 1, . . . , m − 1, d ∑ (Rj +1) [ j=d−i ] − [ ] m TL − µ ai,m TL − µ ∗ ∗ 1 − exp −(γ − γ ) exp cm exp −γm exp i m σ γi − γm∗ σ i =1 where cd and γi are given in Eq. (2), γm∗ = n − m − ai , d =
d ∏
1
k=1
γk − γi
,
∑m−1 i=1
for d = m,
Ri ,
1 ≤ i ≤ d ≤ m,
(12)
k̸=i
and ci,r (α1 , . . . , αr ) =
i ∏
r −i+j
(−1)i
∑
j=1 k=r −i+1
αk
r −i r −i
∏∑ j=1 k=j
, αk
i = 1, 2, . . . , r ,
(13)
458
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
and P (D = m + d)
T −µ d−1 − m − exp −γj exp TLσ−µ exp −(γm∗ − d) exp L σ − c ( 1 , . . . , 1 ) a i,d−1 j,m cd′ cm i + 1 γj − (γm∗ − d) i=0 j=1 T −µ exp −(γm∗ − d + i + 1) exp L σ − exp −γj exp TLσ−µ − for d = 1, . . . , γm∗ − 1, γj − (γm∗ − d) − (i + 1) = T −µ ∗ γ− m m −1 − 1 − exp −γj exp L σ ∗ −1 (1, . . . , 1)aj,m c i ,γ m cγ′ ∗ cm m i+1 γj i = 0 j = 1 TL −µ TL −µ − exp −γ exp exp −( i + 1 ) exp j σ σ − for d = γm∗ , γj − (i + 1) where cd′ = γm∗ (γm∗ − 1) · · · (γm∗ − d + 1). Then the expected number of failures E (D) = without any difficulty.
∑m+γm∗ d=1
d P (D = d) can be computed
5. Simulation study The simulation study has been carried out to examine the performance of the MLEs and Bayes estimators — Lindley’s approximation (Lindley), Tierney and Kadane’s approximation (T&K), and two MCMC methods (which we abbreviate as MCMC1 and MCMC2), for n = 30 and m = 5, TL = ln T = −0.5 and 0.5, and parameters (δ, λ) = (1, 1) and (1.25, 0.5353) under three different sampling schemes: Scheme 1: R1 = · · · = Rm−1 = 0 and Rm = n − m; Scheme 2: R1 = n − m and R2 = · · · = Rm = 0; and Scheme 3: R1 = · · · = Rm−1 = 1 and Rm = n − 2m + 1. The Fortran and R programs and the details of our simulated results are available upon request. For each setting, the average biases and mean squared errors (MSEs) of MLEs and Bayes estimates based on 10 000 simulations when D > 0 are plotted in Figs. 1–6. Here the Bayes estimates are all computed by assuming that δ and λ have Gamma(α, β) and Gamma(a, b) priors, respectively, including the noninformative gamma priors, prior 0 : α = β = a = b = 0, and informative priors, prior k : α = β = a = b = 2k, where k = 1, . . . , 10, and a squared error loss function is used. The associated expected number of failures E (D) and the average number of observed failures based on 10 000 replications, details are not reported, are also computed. They agree closely to each other and this suggests that the Monte Carlo simulation study provides an alternative approach to access the number of failures occurring for the experiment under a particular setting. In terms of minimum biases and MSEs, the performance of the MCMC2 estimators for the parameter λ is generally best followed by the MCMC1 and T&K estimators and then the Lindley estimators and MLEs; in contrast, the ordering of the performance of the estimators (from best to worst) for the parameter δ is T&K, MCMC2, MCMC1, Lindley estimators and MLEs. It can be further observed that the MSEs of the Lindley estimators are usually larger than other Bayes estimators when the prior k gets larger. Compared to the MLEs, the MCMC2, MCMC1 and T&K estimators in most of the cases do greatly improve the efficiency of the estimates of δ and λ. On the whole, the MCMC2 and T&K estimators perform quite well for both parameters, whereas the Lindley estimators and MLEs perform rather poorly. Thus, it would seem reasonable to recommend the use of the MCMC2 and T&K estimators, especially the MCMC2 estimators, for estimating the unknown parameters δ and λ. Similar findings can be also observed for larger m, but the corresponding figures are not presented here for brevity. 6. Planning life test Due to limited available funding, it is very important to find the optimal censoring scheme in practice. It is often to choose the particular censoring scheme which provides more information of the unknown parameters among the class of possible schemes. Our criterion used to compare the possible adaptive progressive hybrid censoring schemes to determine the optimum one under the limited budget is based on the estimation precision of the logarithm of the pth (0 < p < 1) quantile of Weibull distribution, ln Tp = ln[− ln(1 − p)/λ]/δ , as was used by Kundu (2008) as well. Let C1 be the cost for inspecting an item, C2 be the salvage incurred by an unfailed item in the inspection, C3 be the cost per unit time used for life test. Then, the average total cost involved for a given adaptive progressive hybrid censoring scheme is Cb = C1 n − C2 (n − E (D)) + C3 T = (C1 − C2 )n + C2 E (D) + C3 T ,
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
0.1
0.2
MLE Lindley T&K MCMC1 MCMC2
-0.1
0.0
Bias(lambda)
0.1 0.0 -0.1
Bias(delta)
0.2
MLE Lindley T&K MCMC1 MCMC2
0
5
10
15
20
0
5
5
10
15
MSE(lambda)
MLE Lindley T&K MCMC1 MCMC2
0
10
15
20
Prior
0.00 0.05 0.10 0.15 0.20 0.25
0.00 0.05 0.10 0.15 0.20 0.25
Prior
MSE(delta)
459
20
MLE Lindley T&K MCMC1 MCMC2
0
5
Prior
10
15
20
15
20
15
20
Prior
0
Bias(lambda)
5
10
15
MLE Lindley T&K MCMC1 MCMC2
-0.2 0.0 0.2 0.4 0.6
MLE Lindley T&K MCMC1 MCMC2
-0.2 0.0 0.2 0.4 0.6
Bias(delta)
(a) δ = 1 and λ = 1.
20
0
5
MLE Lindley T&K MCMC1 MCMC2
0.4 0.2 0.0
0.2
0.4
MSE(lambda)
0.6
MLE Lindley T&K MCMC1 MCMC2
0.6
10 Prior
0.0
MSE(delta)
Prior
0
5
10
15
20
0
5
Prior
10 Prior
(b) δ = 1.25 and λ = 0.5353. Fig. 1. The biases and MSEs of the estimators for different priors under Scheme 1 when n = 30, m = 5, TL = −0.5.
where n is positive integer, and E (D) ≥ 1 and T > 0. Thus, the problem boils down as follows: for a given budget Cb , among the class of all possible (n, m, (R1 , . . . , Rm ), T )’s, which satisfy
(C1 − C2 )n + C2 E (D) + C3 T ≤ Cb ,
(14)
we want to find the best scheme (n , m , (R1 , . . . , Rm ) , T ), which maximizes the information measure ∗
I (P ) =
∗
1
[∫
Vposterior(P ) (ln Tp )dp
Edata
∗
]−1
∗
,
(15)
0
where Vposterior(P ) (ln Tp ) denotes the posterior variance of ln Tp for a given censoring scheme P = (R1 , . . . , Rm ) and the information measure in (15) is independent of p. Because of the restriction of the integer variables n, m, and (R1 , . . . , Rm ), the nonlinear mixed integer programming (see Grossmann, 2002) can be used to find the optimal solution. From (14) and E (D) ≥ 1, we have (C1 − C2 )n + C2 + C3 T ≤ Cb . C −C Since T > 0, we further have (C1 − C2 )n + C2 ≤ Cb , which implies n ≤ Cb −C2 . Thus, the upper bound of n is obtained. For a 1
2
given value of n, it is clear that C2 E (D) + C3 T ≤ Cb − (C1 − C2 )n. Since E (D) ≥ 1, we therefore obtain T ≤
Cb −(C1 −C2 )n−C2 , C3
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
1.0 0.5
20
5
10 Prior
15
20
2.0
MLE Lindley T&K MCMC1 MCMC2
0.0
0.0
0.5
1.0
1.5
2.0
MLE Lindley T&K MCMC1 MCMC2
0
1.5
15
1.0
10 Prior
MSE(lambda)
5
0.5
0
MSE(delta)
MLE Lindley T&K MCMC1 MCMC2
0.0
0.5 0.0
Bias(delta)
1.0
MLE Lindley T&K MCMC1 MCMC2
Bias(lambda)
460
0
5
10 Prior
15
20
0
5
10 Prior
15
20
1.5 1.0 0.5
Bias(lambda)
-0.5 20
5
10 Prior
15
20
MLE Lindley T&K MCMC1 MCMC2
0.0
0.0
0.5
1.0
1.5
MLE Lindley T&K MCMC1 MCMC2
0
1.5
15
1.0
10 Prior
MSE(lambda)
5
0.5
0
MSE(delta)
MLE Lindley T&K MCMC1 MCMC2
0.0
0.5
1.0
MLE Lindley T&K MCMC1 MCMC2
-0.5 0.0
Bias(delta)
1.5
(a) δ = 1 and λ = 1.
0
5
10 Prior
15
20
0
5
10 Prior
15
20
(b) δ = 1.25 and λ = 0.5353. Fig. 2. The biases and MSEs of the estimators for different priors under Scheme 2 when n = 30, m = 5, TL = −0.5.
the upper bound of T (which is denoted as Tn ). Since larger T yields larger information measure, we shall find the largest value of T , Tn,m,(R1 ,...,Rm ) , which is less than the upper bound Tn and also satisfies the constraint (14) for each given n, m and (R1 , . . . , Rm ). Therefore, the optimization can be determined in a finite number of steps as follows: Step 1. Compute the upper bound of the number of test units n˜ =
C b −C 2 C 1 −C 2
, where [y] is the greatest integer that is less than
or equal to y. Step 2. Set n = 2. Step 3. Compute the upper bound of the length of inspection interval Tn =
Cb −(C1 −C2 )n−C2 . C3
Step 4. Set m = 1. Step 5. Apply the bisection method to obtain Tn,m,(R1 ,...,Rm ) and compute the corresponding information measure, I (P )n,m,(R1 ,...,Rm ) , for all possible choices (n, m, (R1 , . . . , Rm )). Step 6. Set m = m + 1. If m ≤ n − 1 go to Step 5, else go to Step 7. Step 7. Record the particular scheme (n, mn , (R1 , . . . , Rm )n , Tn ), which has the largest information measure, I (P )n , among all possible choices of (n, m, (R1 , . . . , Rm ), Tn,m,(R1 ,...,Rm ) ) for 1 ≤ m ≤ n − 1 in Table 1.
0.3 0.2 0.1
Bias(lambda)
0.10
MSE(lambda)
MLE Lindley T&K MCMC1 MCMC2
0.20
0
5
10 Prior
15
20
MLE Lindley T&K MCMC1 MCMC2
0.00
0.00
MSE(delta)
20
0.30
15
0.20
10 Prior
MLE Lindley T&K MCMC1 MCMC2
0.10
5
0.30
0
461
-0.1 0.0
0.1
0.2
MLE Lindley T&K MCMC1 MCMC2
-0.1 0.0
Bias(delta)
0.3
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
0
5
10 Prior
15
20
0
5
10 Prior
15
20
10 Prior
15
20
10 Prior
15
20
Bias(lambda) 10 Prior
15
20
0
0.6
MLE Lindley T&K MCMC1 MCMC2
0.0
0.0
5
0.4
0.4
MSE(lambda)
MLE Lindley T&K MCMC1 MCMC2
0.6
MLE Lindley T&K MCMC1 MCMC2
0.8
5
0.2
MSE(delta)
0.8
0
-0.2 0.0 0.2 0.4 0.6 0.8
MLE Lindley T&K MCMC1 MCMC2
0.2
Bias(delta)
-0.2 0.0 0.2 0.4 0.6 0.8
(a) δ = 1 and λ = 1.
0
5
10 Prior
15
20
0
5
(b) δ = 1.25 and λ = 0.5353. Fig. 3. The biases and MSEs of the estimators for different priors under Scheme 3 when n = 30, m = 5, TL = −0.5.
Step 8. Set n = n + 1. If n ≤ n˜ go to Step 3, else go to Step 9. Step 9. Obtain the optimal scheme (n∗ , m∗ , (R1 , . . . , Rm )∗ , T ∗ ), which has the largest information measure, I (P )∗ , among all possible choices of (n, mn , (R1 , . . . , Rm )n , Tn ) for 1 ≤ n ≤ n˜ in Table 1. Note that it is not possible to compute Edata
1 0
Vposterior(P ) (ln Tp )dp analytically, so we conduct a Monte Carlo simulation
study along with the MCMC2 method discussed in Section 3.3 to approximate I (P ). More precisely, for fixed p, we first approximate Vposterior(P ) (ln Tp ) = Eposterior(P ) (ln Tp )2 − (Eposterior(P ) (ln Tp ))2 by using the MCMC2 method with K = 1000 to obtain the Bayes estimates of g1 (δ, λ) = ln Tp (δ, λ) = (ψ − ln λ)/δ and g2 (δ, λ) = (ln Tp (δ, λ))2 = [ψ 2 − 2ψ ln λ + (ln λ)2 ]/δ 2 , where ψ = ln(− ln(1 − p)), and then compute
1 0
Vposterior(P ) (ln Tp )dp. Upon repeating the same procedure 10 000 times, the average of
an approximation for I (P ) = Edata
1 0
Vposterior(P ) (ln Tp )dp .
1 0
Vposterior(P ) (ln Tp )dp’s shall be
15
0.06 0.02
20
MLE Lindley T&K MCMC1 MCMC2
0
5
10 Prior
15
0
MSE(lambda)
10 Prior
5
10 Prior
15
0.00 0.02 0.04 0.06 0.08 0.10
5
0.00 0.02 0.04 0.06 0.08 0.10
0
MSE(delta)
MLE Lindley T&K MCMC1 MCMC2
-0.02
Bias(lambda)
0.06 0.02 -0.02
Bias(delta)
MLE Lindley T&K MCMC1 MCMC2
0.10
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467 0.10
462
20
20
MLE Lindley T&K MCMC1 MCMC2
0
5
10 Prior
15
20
10 Prior
15
20
10 Prior
15
20
5
0.0 0.2 0.4 0.6
Bias(lambda) 10 Prior
15
20
0
5
MSE(lambda)
MLE Lindley T&K MCMC1 MCMC2
0
10 Prior
15
20
0.00 0.05 0.10 0.15 0.20
0.00 0.05 0.10 0.15 0.20
0
MSE(delta)
MLE Lindley T&K MCMC1 MCMC2
-0.4
0.0 0.2 0.4 0.6
MLE Lindley T&K MCMC1 MCMC2
-0.4
Bias(delta)
(a) δ = 1 and λ = 1.
5
MLE Lindley T&K MCMC1 MCMC2
0
5
(b) δ = 1.25 and λ = 0.5353. Fig. 4. The biases and MSEs of the estimators for different priors under Scheme 1 when n = 30, m = 5, TL = 0.5.
For illustrative purposes, let us consider the case Cb = 100, C1 = 12, C2 = 3, and C3 = 10 when (δ, λ) = (1, 1) and (1.25, 0.5353), respectively. Assuming that both δ and λ have independent Gamma(20, 20) prior distributions, we report the largest information measure I (P )n and the corresponding censoring scheme (n, mn , (R1 , . . . , Rm )n , Tn ) for 2 ≤ n ≤ n˜ = 10 in Table 1. Therefore, the optimal progressive censoring scheme (n∗ , m∗ , (R1 , . . . , Rm )∗ , T ∗ ) is (7,3, (0,0,4), 1.9107) for (δ, λ) = (1, 1), and (6,4, (0,0,0,2), 3.0147) for (δ, λ) = (1.25, 0.5353). Finally, we have also checked the effect of efficiency if we depart from the optimal scheme slightly. For instance, we depart from the optimal scheme (7, 3, (0, 0, 4), 1.9107) for (δ, λ) = (1, 1) slightly, say (0,1,3), which has an information measure 8.8763, then the relative efficiency is 8.8763/9.1181 = 0.9735; and at (1,0,3), which has an information measure 8.8445, the relative efficiency becomes 8.8445/9.1181 = 0.9700; they are near optimal. This shows that the efficiency does not change much if we depart from the optimal scheme only slightly. The same conclusion can also be made for (δ, λ) = (1.25, 0.5353). Acknowledgments We express our sincere thanks to the editor, associate editor, and the reviewers for their useful suggestions which led to an improvement over an earlier version of this manuscript. The first author thanks the National Science of Council of Taiwan (NSC Grant Number 99-2118-M-032-011-MY3) for funding this research.
Bias(lambda)
0
5
10 Prior
15
20
MLE Lindley T&K MCMC1 MCMC2
0.0
0.5
1.0
MLE Lindley T&K MCMC1 MCMC2
0.0
MSE(delta)
20
1.5
15
MLE Lindley T&K MCMC1 MCMC2
1.0
10 Prior
463
0.5
5
1.5
0
-0.2 0.0 0.2 0.4 0.6 0.8
MLE Lindley T&K MCMC1 MCMC2
MSE(lambda)
Bias(delta)
-0.2 0.0 0.2 0.4 0.6 0.8
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
0
5
10 Prior
15
20
0
5
10 Prior
15
20
10 Prior
15
20
10 Prior
15
20
10 Prior
15
0.2 0.4 0.6 0.8 1.0
20
0
0
MSE(lambda)
MLE Lindley T&K MCMC1 MCMC2
5
10 Prior
15
0.0 0.1 0.2 0.3 0.4 0.5 0.6
5
0.0 0.1 0.2 0.3 0.4 0.5 0.6
0
MSE(delta)
MLE Lindley T&K MCMC1 MCMC2
-0.2
Bias(lambda)
MLE Lindley T&K MCMC1 MCMC2
-0.2
Bias(delta)
0.2 0.4 0.6 0.8 1.0
(a) δ = 1 and λ = 1.
5
MLE Lindley T&K MCMC1 MCMC2
0
20
5
(b) δ = 1.25 and λ = 0.5353. Fig. 5. The biases and MSEs of the estimators for different priors under Scheme 2 when n = 30, m = 5, TL = 0.5.
Appendix. Probability Mass Function of D For extreme value distribution, the pdf of Xd:m:n , d = 1, . . . , m, is given by (see Balakrishnan and Aggarwala, 2000) fXd:m:n (x) = cd
d −
ai,d
i=1
1
σ
[ exp
x−µ
σ
− γi exp
x−µ
σ
]
,
−∞ < x < ∞,
where cd , γi , and ai,d are given in Eqs. (2) and (12). Let D be the number of failures up to time T and TL = ln T . Clearly, TL −µ Pr(D = 0) = Pr(X1:m:n > TL ) = exp −n exp . It then follows from the joint density of Xr :m:n and Xs:m:n for σ 1 ≤ r < s ≤ m (see Balakrishnan, 2007), fXr :m:n ,Xs:m:n (xr , xs ) = cs
r −1 s− r −1 − − i=0
k=0
ci,r −1 (R1 + 1, . . . , Rr −1 + 1)ck,s−r −1 (Rr +1 + 1, . . . , Rs−1 + 1)
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
10 Prior
15
0.10 0.05
Bias(lambda)
0.00
5
10 Prior
15
20
MLE Lindley T&K MCMC1 MCMC2
0.00
0.04
MSE(lambda)
MLE Lindley T&K MCMC1 MCMC2
0.08
0
0.00
MSE(delta)
20
0.12
5
0.12
0
MLE Lindley T&K MCMC1 MCMC2
0.08
0.05 0.00
Bias(delta)
0.10
MLE Lindley T&K MCMC1 MCMC2
0.04
464
0
5
10 Prior
15
20
0
5
10 Prior
15
20
10 Prior
15
20
10 Prior
15
20
10 Prior
15
0.0 0.2 0.4 0.6 0.8
20
0
5
MSE(lambda)
MLE Lindley T&K MCMC1 MCMC2
0
10 Prior
15
20
5
0.00 0.05 0.10 0.15 0.20
5
0.00 0.05 0.10 0.15 0.20
0
MSE(delta)
MLE Lindley T&K MCMC1 MCMC2
-0.4
Bias(lambda)
MLE Lindley T&K MCMC1 MCMC2
-0.4
Bias(delta)
0.0 0.2 0.4 0.6 0.8
(a) δ = 1 and λ = 1.
MLE Lindley T&K MCMC1 MCMC2
0
5
(b) δ = 1.25 and λ = 0.5353. Fig. 6. The biases and MSEs of the estimators for different priors under Scheme 3 when n = 30, m = 5, TL = 0.5. s−∑ 1−k
× f (xr )[1 − F (xr )]
(Rj +1)−1
j=r −i
f (xs )[1 − F (xs )]
n−
s−∑ 1−k
(Rj +1)−1
j=1
,
xr < xs ,
where ci,j is given in Eq. (13), to have P (D = d) = P (Xd:m:n < TL < Xd+1:m:n )
= cd+1
d−1 −
ci,d−1 (R1 + 1, . . . , Rd−1 + 1)
i =0
∫
∞
∫
TL
× TL
= cd+1
d−1 − i =0
d ∑
f (xd )[1 − F (xd )]
(Rj +1)−1
j=d−i
−∞
ci,d−1 (R1 + 1, . . . , Rd−1 + 1)
f (xd+1 )[1 − F (xd+1 )]γd+1 −1 dxd dxd+1
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
465
Table 1 The information measure I (P )n for (n, mn , (R1 , . . . , Rm )n , Tn ) when 2 ≤ n ≤ n˜ = 10. n
mn
(R1 , . . . , Rm )n
Tn
Tn
E ( D)
I (P )n
1 2 1 3 5 3 4 2 4
(1) (1,0) (3) (1,0,1) (0,0,0,0,1) (0,0,4) (1,0,0,3) (5,2) (0,6,0,0)
7.9000 7.0000 6.1000 5.2000 4.3000 3.4000 2.5000 1.6000 0.7000
7.6003 6.7006 5.2066 4.3150 2.8991 1.9107 1.2749 1.0532 0.3612
1.9990 1.9982 3.9781 3.9499 5.6696 5.9642 5.0837 2.8227 2.1293
4.2435 4.8165 7.3065 7.6797 9.0157 9.1181 8.3129 7.0320 5.3858
(1) (1,0) (0,0,1) (1,0,0,0) (0,0,0,2) (0,1,0,0,1) (0,2,0,0,1) (1,5,0) (0,6,1)
7.9000 7.0000 6.1000 5.2000 4.3000 3.4000 2.5000 1.6000 0.7000
7.6007 6.7014 5.2176 4.3394 3.0147 2.2762 1.6078 1.2083 0.5034
1.9977 1.9953 3.9413 3.8687 5.2843 4.7460 3.9741 2.3058 1.6552
3.6199 4.1332 6.0566 6.4624 7.5273 7.4703 6.7446 6.3858 4.4670
δ = 1 and λ = 1 2 3 4 5 6 7 8 9 10
δ = 1.25 and λ = 0.5353 2 3 4 5 6 7 8 9 10
1 2 3 4 4 5 5 3 3
×
d ∑ TL −µ 1 − exp − (Rj + 1) exp σ TL −µ exp −γ exp d + 1 j=d−i σ d ∑
(Rj + 1)
γd+1
j=d−i
] − [ d−1 TL − µ ci,d−1 (R1 + 1, . . . , Rd−1 + 1) = cd exp −γd+1 exp σ i=0 d ∑ TL −µ 1 − exp − (Rj + 1) exp σ j=d−i × , for d = 1, . . . , m − 1. d ∑ (Rj + 1) j=d−i
Given Xm:m:n = xm , Xm+1:n is now distributed as the first order statistic of a random sample of size γm∗ = n − m − from a truncated extreme value distribution with cumulative distribution function ] [ xm − µ x−µ ∗ ∗ + γm exp , xm < x < ∞. FXm+1:n (x|Xm:m:n = xm ) = 1 − exp −γm exp
σ
∑m−1 i=1
Ri
σ
Therefore, we have P (D = m) = P (Xm:m:n < TL ≤ Xm+1:n )
∫
∞
P (x < TL ≤ Xm+1:n |Xm:m:n = x)fXm:m:n (x)dx
= −∞
∫
TL
=
[1 − FXm+1:n (TL |Xm:m:n = x)]fXm:m:n (x)dx
−∞
[ ] − [ ] m TL − µ TL − µ ai , m ∗ = cm exp −γm∗ exp 1 − exp −(γ − γ ) exp . i m σ γi − γm∗ σ i=1 In a similar fashion, Xm+r :n and Xm+s:n (for 1 ≤ r < s ≤ γm∗ ) are distributed as the r th and sth order statistics of a random sample of size γm∗ from a truncated exponential distribution with joint density function (see Balakrishnan et al., 2002) fXr :n ,Xs:n |Xm:m:n (xr , xs |xm ) = cs′
r −1 s− r −1 − − i=0
[ ×
f (xr ) 1 − F (xm )
][
f (xs ) 1 − F ( xm )
ci,r −1 (1, . . . , 1)ck,s−r −1 (1, . . . , 1)
k=0
][
1 − F (xr ) 1 − F (xm )
]s−r −1−k+i [
1 − F (xs ) 1 − F (xm )
]γm∗ −s+k
,
xm < xr < xs < ∞,
466
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
where cs′ = γm∗ (γm∗ − 1) · · · (γm∗ − s + 1). Then, for d = 1, . . . , γm∗ − 1, we have P (Xm+d:n < x < Xm+d+1:n |Xm:m:n = xm ) ′
= cd+1
d−1 −
] ∗ xm − µ ci,d−1 (1, . . . , 1) exp γm − d + i + 1 exp σ [
i =0 ∞∫ x
∫ ×
f (y)f (z )[1 − F (y)]i [1 − F (z )]γm −d−1 dydz ∗
xm
x
] ∗ xm − µ = cd+1 ci,d−1 (1, . . . , 1) exp γm − d + i + 1 exp σ i =0 xm −µ x−µ x−µ exp −(i + 1) exp − exp −(i + 1) exp σ exp −(γm∗ − d) exp σ σ × i+1 γm∗ − d [ ] xm − µ x−µ = cd′ exp (γm∗ − d) exp − exp σ σ x −µ d−1 − 1 − exp (i + 1) exp mσ − exp x−µ σ . × ci,d−1 (1, . . . , 1) i+1 i=0 d−1 −
′
[
It leads to P (D = m + d) = P (Xm+d:n < TL ≤ Xm+d+1:n ) ∞
∫
P (x < Xm+d:n < TL ≤ Xm+d+1:n |Xm:m:n = x)fXm:m:n (x)dx
= −∞
= cd′ cm
TL −µ TL −µ ∗ d−1 − m − exp −γ exp exp −(γ − d ) exp − j m σ σ ci,d−1 (1, . . . , 1)aj,m
γj − (γm∗ − d) T −µ exp −(γm∗ − d + i + 1) exp L σ − exp −γj exp TLσ−µ − , for d = 1, . . . , γm∗ − 1. γj − (γm∗ − d) − (i + 1) i+1
i=0 j=1
Finally, we use the representation of Balakrishnan et al. (2002) for the pdf of Xγm∗ :n to have P (D = m + γm∗ ) = P (Xγm∗ :n < TL ) ∞
∫ =
−∞ TL
∫
P (x < Xγm∗ :n < TL |Xm:m:n = x)fXm:m:n (x)dx TL
∫
= −∞
= cγ′ m∗ cm
×
x
cγ′ ∗
∗ γ− m −1
m
x−µ
]
σ
i =0
∗ γ− m m −1 −
f (y)[1 − F (y)]i dyfXm:m:n (x)dx
ci,γm∗ −1 (1, . . . , 1)aj,m
i+1 i=0 j=1 1 − exp −γj exp TLσ−µ
[
ci,γm∗ −1 (1, . . . , 1) exp (i + 1) exp
γj
exp −(i + 1) exp
−
TL −µ
σ
− exp −γj exp TLσ−µ
γj − (i + 1)
.
References Abramowitz, M., Stegun, I.A. (Eds.), 1964. Handbook of Mathematical Functions, with Formulas, Graphs and Mathematical Tables. Dover, New York. Balakrishnan, N., 2007. Progressive censoring methodology: an appraisal (with discussions). Test 16, 211–296. Balakrishnan, N., Aggarwala, R., 2000. Progressive Censoring: Theory, Methods and Applications. Birkhäuser, Boston. Balakrishnan, N., Childs, A., Chandrasekar, B., 2002. An efficient computational method for moments of order statistics under progressive censoring. Statistics and Probability Letters 60, 359–365. Banerjee, A., Kundu, D., 2008. Inference based on type-II hybrid censored data from a Weibull distribution. IEEE Transactions on Reliability 57, 369–378. Childs, A., Chandrasekar, B., Balakrishnan, N., 2008. Exact likelihood inference for an exponential parameter under progressive hybrid censoring schemes. In: Vonta, F., Nikulin, M., Limnios, N., Huber-Carol, C. (Eds.), Statistical Models and Methods for Biomedical and Technical Systems. Birkhäuser, pp. 323–334. Devroye, L., 1984. A simple algorithm for generating random variates with a log-concave density. Computing 33, 247–257.
C.-T. Lin et al. / Computational Statistics and Data Analysis 56 (2012) 451–467
467
Grossmann, I.E., 2002. Review of nonlinear mixed-integer and disjunctive programming techniques. Optimization and Engineering 3, 227–252. Kundu, D., 2007. On hybrid censoring Weibull distribution. Journal of Statistical Planning and Inference 137, 2127–2142. Kundu, D., 2008. Bayesian inference and life testing plan for the Weibull distribution in presence of progressive censoring. Technometrics 50, 144–154. Kundu, D., Joarder, A., 2006. Analysis of type-II progressively hybrid censored data. Computational Statistics & Data Analysis 50, 2509–2528. Lin, C.T., Huang, Y.L., 2011. On progressive hybrid censored exponential distribution. Journal of Statistical Computation and Simulation, First published on: 21 June 2011 (iFirst). Lin, C.T., Ng, H.K.T., Chan, P.S., 2009. Statistical inference of type-II progressively hybrid censored data with Weibull lifetimes. Communications in Statistics— Theory and Methods 38, 1710–1729. Lin, C.T., Wu, S.J.S., Balakrishnan, N., 2006. Inference for log-gamma distribution based on progressively type-II censored data. Communications in Statistics— Theory and Methods 35, 1271–1292. Lindley, D.V., 1980. Approximate Bayesian method. Trabajos de Estadistica 31, 223–237. Ng, H.K.T., Chan, P.S., Balakrishnan, N., 2002. Estimation of parameters from progressively censored data using EM algorithm. Computational Statistics & Data Analysis 39, 371–386. Ng, H.K.T., Kundu, D., Chan, P.S., 2009. Statistical analysis of exponential lifetimes under an adaptive type-II progressive censoring scheme. Naval Research Logistics 56, 687–698. Tierney, L., Kadane, J.B., 1986. Accurate approximations for posterior moments and marginal densities. Journal of the American Statistical Association 81, 82–86.