- Email: [email protected]
Estimation for the three-parameter lognormal distribution based on progressively censored data
Computational Statistics and Data Analysis 53 (2009) 3580–3592
Contents lists available at ScienceDirect
Computational Statistics and Data Analysis journal homepage: www.elsevier.com/locate/csda
Estimation for the three-parameter lognormal distribution based on progressively censored data Prasanta Basak a , Indrani Basak a,∗ , N. Balakrishnan b a b
Penn State Altoona, Altoona, PA, United States McMaster University, Hamilton, Canada
article
info
Article history: Received 27 February 2009 Accepted 27 March 2009 Available online 1 April 2009
a b s t r a c t Some work has been done in the past on the estimation of parameters of the threeparameter lognormal distribution based on complete and censored samples. In this article, we develop inferential methods based on progressively Type-II censored samples from a three-parameter lognormal distribution. In particular, we use the EM algorithm as well as some other numerical methods to determine maximum likelihood estimates (MLEs) of parameters. The asymptotic variances and covariances of the MLEs from the EM algorithm are computed by using the missing information principle. An alternative estimator, which is a modification of the MLE, is also proposed. The methodology developed here is then illustrated with some numerical examples. Finally, we also discuss the interval estimation based on large-sample theory and examine the actual coverage probabilities of these confidence intervals in case of small samples by means of a Monte Carlo simulation study. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Lognormal distribution is one of the distributions commonly used for modeling lifetimes or reaction-times, and is particularly useful for modeling data which are long-tailed and positively skewed. It has been discussed extensively by many authors including Cohen (1951, 1988), Hill (1963), Harter and Moore (1966), Crow and Shimizu (1988), Johnson et al. (1994), Munro and Wixley (1970) and Rukhin (1984). It is well-known that there is a close relationship between normal and lognormal distributions. If X = log(Y − γ ) is normally distributed with mean µ and standard deviation σ , then the distribution of Y becomes a three-parameter lognormal distribution with parameter θ = (γ , µ, σ ). The probability density function of such a three-parameter lognormal distribution is
[log(y − γ ) − µ]2 √ 1 exp − , f (y; θ) = σ 2π (y − γ ) 2σ 2 0, otherwise.
γ < y < ∞, σ > 0, −∞ < µ < ∞
(1.1)
In (1.1), σ 2 and µ are the variance and mean of the underlying normal variable X , but become the shape and scale parameters of the lognormal variable Y . It is more convenient to use β = exp(µ) and w = exp(σ 2 ) as the scale and shape parameters of the lognormal variable Y , respectively. Also, γ is the threshold (location) parameter of the lognormal random variable.
∗
Corresponding author. Tel.: +1 814 9495263; fax: +1 814 949 5456. E-mail address: [email protected] (I. Basak).
0167-9473/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.csda.2009.03.015
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
3581
When the threshold parameter γ is known, the parameter estimation can be done using the well-known results for normal distribution by making the transformation from Y to X . Estimation methods become more complex when γ is unknown. Adding to that complexity is the fact that quite often the data on lifetimes or reaction-times come with ‘‘censoring.’’ Censoring occurs when exact survival times are known only for a portion of the individuals or items under study. The complete survival times may not have been observed by the experimenter either intentionally or unintentionally, and there are numerous examples for each kind; see, for example, Nelson (1982) and Balakrishnan and Cohen (1991). In this article, we consider a general scheme of progressive Type-II right censoring. Under this scheme, n units are placed on a life-testing experiment and only m (< n) are completely observed until failure. The censoring occurs progressively in m stages. These m stages offer failure times of the m completely observed units. At the time of the first failure (the first stage), R1 of the n − 1 surviving units are randomly withdrawn (censored intentionally) from the experiment, R2 of the n − 2 − R1 surviving units are withdrawn at the time of the second failure (the second stage), and so on. Finally, at the time of the m-th failure (the m-th stage), all the remaining Rm = n − m − R1 − · · · − Rm−1 surviving units are withdrawn. This scheme (R1 , R2 , . . . , Rm ) is referred to as progressive Type-II right censoring scheme. It is clear that this scheme includes the conventional Type-II right censoring scheme (when R1 = R2 = · · · = Rm−1 = 0 and Rm = n − m) and complete sampling scheme (when n = m and R1 = R2 = · · · = Rm = 0). The ordered lifetime data which arise from such a progressive Type-II right censoring scheme are called progressively Type-II right censored order statistics. For theory, methods and applications of progressive censoring, readers are referred to the book by Balakrishnan and Aggarwala (2000) and the recent discussion paper by Balakrishnan (2007). Suppose n independent lognormally distributed units are placed on a life-testing experiment. Let Y1:m:n ≤ · · · ≤ Ym:m:n denote the above mentioned m progressively Type-II right censored order statistics. For ease in notation, let us use Yj (j = 1, . . . , m) to denote these Yj:m:n ’s. Note that we observe only Y = (Y1 , . . . , Ym ). The purpose of this article is to discuss different estimation procedures for the parameter θ based on the progressively Type-II right censored order statistics Y. Recently, Ng et al. (2002) considered two-parameter lognormal distribution (which does not include the threshold parameter) and discussed Newton–Raphson algorithm as well as EM algorithm for finding the MLEs. Inclusion of the third parameter γ introduces an unusual feature in the likelihood function. Hill (1963) has shown that there exist paths along which the likelihood function of any ordered sample y1 , . . . , yn from the three-parameter lognormal distribution tends to ∞ as (γ , µ, σ ) approaches (y1 , −∞, ∞). Thus, the global maximization of the likelihood leads to the unreasonable estimate (y1 , −∞, ∞) although, in fact, the likelihood at the point is zero. But Hill (1963) also retained the idea that reasonable estimates could be obtained by solving the likelihood equations. To get these reasonable estimates, Cohen (1951), Cohen and Whitten (1988) and Harter and Moore (1966) equated partial derivatives of the likelihood function to zero and solved the resulting equations. These estimates are called local maximum likelihood estimates (LMLEs). Harter and Moore (1966) and Calitz (1973) noted that although these LMLEs are not true MLEs according to the usual definition, they are reasonable estimates and appear to possess most of the desirable properties associated with MLEs. However, it is noted in Harter and Moore (1966) and Cohen and Whitten (1988) that sometimes the likelihood function may have no clearly defined local maximum for small samples and so LMLEs fail to produce estimates in that case. In the past, some work has been done on estimation methods for the three-parameter lognormal distribution based on complete and censored samples; see, for example, the books by Cohen and Whitten (1988) and Balakrishnan and Cohen (1991). In Section 2, we describe the Newton–Raphson algorithm for determining the MLEs of the parameter θ based on a progressively censored sample. The second derivatives of the log-likelihood are required in order to use the algorithm. These computations are complicated when data are progressively censored. Another viable alternative to Newton–Raphson algorithm is the well-known EM algorithm and in Section 3 we discuss how that can be used to determine the MLEs in this case. Asymptotic variances and covariances of the maximum likelihood estimates generated through the EM algorithm are given in Section 4. Section 5 describes some simplified estimation methods which yields simple alternative estimators. All these methods discussed in this article are then illustrated with some numerical examples in Section 6. With these illustrative examples, we also discuss the interval estimation based on large-sample theory and then examine the actual coverage probabilities in case of small samples through Monte Carlo simulations. 2. Newton–Raphson algorithm One of the standard methods of determining the maximum likelihood estimates is the Newton–Raphson algorithm. In this section, we describe the Newton–Raphson algorithm for finding the MLEs numerically when life-times are distributed as a three-parameter lognormal distribution with parameter θ . These MLEs are local MLEs, as mentioned earlier, corresponding to the Newton–Raphson algorithm, and would be denoted here by LMLE1. The log-likelihood function log L(θ) = log L based on the progressively Type-II right censored order statistics Y, is log L(θ) = const . − m log σ −
m X j =1
log(y −γ )−µ
log(yj − γ ) −
m 1X
2 j =1
Ψ0j2 +
m X
¯ Ψ0j Rj log Φ
,
(2.1)
j =1
j ¯ = 1−Φ where Ψ0j = ; see Balakrishnan and Aggarwala (2000). In (2.1) and throughout this article, φ and Φ σ denote the probability density and survival function of the standard normal distribution, respectively. Three likelihood
3582
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
equations which need to be solved simultaneously for the required estimate θˆ = (µ, ˆ σˆ , γˆ ) are as follows: m m X 1 ∂ log L 1 X 1 = Ψ0j + Rj L0j + =0 ∂γ σ j=1 yj − γ y − γ j j =1 " # m m X 1 X ∂ log L = Ψ0j + Rj L0j = 0 ∂µ σ j =1 j =1 " # m m X X ∂ log L 1 2 = −m + Ψ0j + Rj Ψ0j L0j = 0 ∂σ σ j =1 j =1
(2.2)
¯ (Ψ0j ) is the hazard function of the standard normal distribution at Ψ0j . in which L0j = φ(Ψ0j )/Φ In the Newton–Raphson algorithm, the simultaneous solution is obtained through an iterative procedure. In each iterative step, the corrections a, b, c to the previous estimates γ0 , µ0 , σ0 produce new estimates γˆ , µ, ˆ σˆ as γˆ = γ0 + a,
µ ˆ = µ0 + b and σˆ = σ0 + c .
The iteration method is based on Taylor series expansions of the estimating equations in (2.2) in the neighborhood of the previous simultaneous estimates. Neglecting powers of a, b and c above the first order and using Taylor’s theorem, we get the following equations which need to be solved for a, b, c:
∂ 2 log L ∂ 2 log L + b +c ∂γ 2 0 ∂γ ∂µ 0 ∂ 2 log L ∂ 2 log L +b +c a ∂µ∂γ 0 ∂µ2 0 ∂ 2 log L ∂ 2 log L a + b +c ∂σ ∂γ ∂σ ∂µ
a
0
0
∂ 2 log L ∂ log L = − ∂γ ∂σ 0 ∂γ 0 ∂ 2 log L ∂ log L =− ∂µ∂σ 0 ∂µ 0 ∂ 2 log L ∂ log L =− , ∂σ 2 0 ∂σ 0
(2.3)
where the notation A|0 , for any partial derivative A, means the partial derivative evaluated at (γ0 , µ0 , σ0 ). The second derivatives needed in (2.3) are as follows: m m m ∂ 2 log L X 1 1 X 1 1 X 1 = + ( Ψ + R L ) − (1 − Rj Ψ0j L0j + Rj L20j ) 0j j 0j 2 2 2 ∂γ 2 ( y − γ ) σ ( y − γ ) σ ( y − γ )2 j j j j=1 j =1 j =1 " # m X ∂ 2 log L 1 = − m + Rj L0j (L0j − Ψ0j ) ∂µ2 σ2 j =1 " # m X 2 ∂ 2 log L 1 3 2 2 = 2 m− 3Ψ0j + Rj (2Ψ0j L0j − Ψ0j L0j + Ψ0j L0j ) ∂σ 2 σ j =1 m 1 1 X ∂ 2 log L =− 2 1 − Rj Ψ0j L0j + Rj L20j ∂γ ∂µ σ j =1 yj − γ m 1 ∂ 2 log L 1 X =− 2 2Ψ0j + Rj (L0j − Ψ0j2 L0j + Ψ0j L20j ) ∂γ ∂σ σ j =1 yj − γ m ∂ 2 log L 1 X =− 2 2Ψ0j + Rj (L0j − Ψ0j2 L0j + Ψ0j L20j ) . ∂µ∂σ σ j =1
(2.4)
Using the Newton–Raphson method, Ng et al. (2002) discussed methods of finding MLEs while Balakrishnan et al. (2003) discussed the construction of confidence intervals of µ and σ for the two-parameter lognormal distribution. Mi and Balakrishnan (2003) have made use of the fact that the lognormal density is log-concave in order to establish that the MLEs of µ and σ do exist and are unique. Hence, in that case, the EM algorithm and the Newton–Raphson algorithm will converge to the same values. In order to obtain LMLE1 for the three-parameter lognormal distribution, we suggest a variation of the above method, which was used by Cohen (1951). Calitz (1973) found this method of Cohen (1951) produced better convergence in the estimation process. Instead of solving all three equations in (2.2) simultaneously, the suggested method is to start with a γ(0) which is less than y1 and use the second and third equations of (2.2) to get µ(0) = µ(γ(0) ) and σ(0) = σ (γ(0) ). The first equation is used as the test equation. If the left hand side of this equation equals zero when (γ(0) , µ(0) , σ(0) ) substituted there, then no further iteration is required. Otherwise, another γ(1) is chosen which is less than y1 yielding µ(1) = µ(γ(1) )
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
3583
and σ(1) = σ (γ(1) ) and look for the sign change in the left hand side of the first equation of (2.2). Finally, interpolating on γ values, first find the final estimate γˆ and then the final estimates µ ˆ = µ(γˆ ) and σˆ = σ (γˆ ) are obtained using the second and third equations of (2.2). The final estimates (γˆ , µ, ˆ σˆ ) when substituted in the first equation of (2.2) should produce zero or close to zero. As noted by Cohen and Whitten (1988), usually a single root γˆ will be found. In the event that multiple roots of γˆ occur, Cohen and Whitten (1988) suggested to use the root which results in the closest agreement between y¯ and 2
ˆ σˆ2 E[ [Y ] = γˆ + eµ+ .
3. EM Algorithm In the case of progressively censored samples, an alternative to the Newton–Raphson algorithm is the use of EM algorithm for numerically finding the MLEs. One advantage of the EM algorithm is that asymptotic variances and covariances of the EM algorithm estimates can be computed, which is discussed in Section 4. EM algorithm, introduced by Dempster et al. (1977), is a very popular tool to handle any missing or incomplete data situation; readers are referred to the book by McLachlan and Krishnan (1997) for a detailed discussion on EM algorithm and its applications. This algorithm is an iterative method which has two steps. In the E-step, it replaces any missing data by its expected value and in the M-step the log-likelihood function is maximized with the observed data and expected value of the incomplete data, producing an update of the parameter estimates. The MLEs of the parameters are obtained by repeating the E- and M-steps until convergence occurs. Since progressive censoring model can be viewed as a missing data problem, EM algorithm can be applied to obtain the MLEs of the parameters in this case. Let us denote the censored data vector as Z = (Z1 , Z2 , . . . , Zm ), where the j-th stage censored data vector Zj is a 1 ×Rj vector, Zj = (Zj1 , Zj2 , . . . , ZjRj ), for j = 1, 2, . . . , m. The complete data set is then obtained by combining the observed data Y and the censored data Z. E-step of the algorithm requires the computation of the conditional expectation of functions of censored data vector Z, conditional on the observed data vector Y and the current value of the parameters. In particular, one computes the conditional expectation of the log-likelihood E [log L(Y, Z, θ )|Y = y] as E [log L(Y, Z, θ )|Y = y] = const . − n log σ −
m X
log(yj − γ ) −
m 1X
log(yj − γ ) − µ
2
2 j =1 σ " # 2 Rj Rj m X m X X 1X log(Zjk − γ ) − µ E log(Zjk − γ ) | Zjk > yj − E − | Zjk > yj . 2 j=1 k=1 σ j=1 k=1 j =1
(3.1)
The above conditional expectations are obtained using the result that given Yj = yj , Zj ’s have a left-truncated distribution F , truncated at yj . More specifically, the conditional probability density of Z, given Y, is given by (see Balakrishnan and Aggarwala (2000))
fZ|Y (z|y; θ) =
Rj m Y Y
fZjk |Yj (zjk |yj ; θ),
(3.2)
j=1 k=1
where fZjk |Yj (zjk |yj ; θ) =
f (zjk ; θ)
(3.3)
1 − F (yj ; θ)
and f (zjk ; θ) is given by (1.1) and F denotes the corresponding cumulative distribution function. In the M-step of the (h + 1)-th iteration, we will denote the updated estimates of the parameter θ as θ (h+1) . This θ (h+1) maximizes the loglikelihood function involving the observed data Y, conditional expectation of the log-likelihood function of censored data vector Z given the observed data vector Y, and the h-th iteration value of the parameter θ (h) . As a starting value θ (0) , one can use a γ(0) < y1 and µ(0) and σ(0) are computed on the basis of the so-called ‘‘pseudo-complete’’ sample which involves observed data Y and the censored observations at the j-th step Zj all taken to be yj . Thus, µ(0) = µ(γ(0) ) and σ(0) = σ (γ(0) ) are then given by
µ(0) = σ(0)
m 1X
n
(Rj + 1) log yj − γ(0)
1 v j=" # . u m u1 X 2 2 t = (Rj + 1) log yj − γ(0) − µ(0) n j=1
(3.4)
3584
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
Starting with the initial estimates in (3.4), the (h + 1)-th iteration value θ (h+1) is obtained using the h-th iteration value θ h as follows:
µ(h+1) = σ(h+1)
1
( m X
n
log yj − γ(h) +
m X
Rj E log(Z − γ(h) )|Z > yj ; θ (h)
)
j =1
j=1
v ( ) u m m u1 X X 2 2 t = Rj E log (Z − γ(h) )|Z > yj ; µ(h+1) , σ(h) , γ(h) − µ2(h+1) log yj − γ(h) + n
j =1
j =1
.
(3.5)
The conditional expectations in the above expression (3.5) can be obtained as follows:
)
E log(Z − γ(h) )|Z > yj ; θ (h) = σ(h) L0j(h) + µ(h)
E log2 (Z − γ(h) )|Z > yj ; µ(h+1) , σ(h) , γ(h) = σ(2h) 1 + Ψ0j∗ L∗0j + 2σ(h) µ(h+1) L∗0j + µ2(h+1)
where log yj − γ(h) − µ(h)
Ψ0j(h) = Ψ0j (θ (h) ) =
σ(h)
,
φ(Ψ0j(h) ) , ¯ (Ψ0j(h) ) Φ log yj − γ(h) − µ(h+1) Ψ0j∗ = , σ(h) φ(Ψ0j∗ ) L∗0j = . ¯ (Ψ0j∗ ) Φ L0j(h) = L0j (θ (h) ) =
γ(h+1) is obtained by solving the following equation for γ : m X µ(h+1) − σ(2h+1)
−
m X
Rj E
j =1 ∗ θ (h+1) =
m X + µ(h+1) − σ(2h+1) Rj E
− yj − γ j =1 log(Z − γ ) Z > yj ; γ , θ ∗ (h+1) = 0, Z −γ j =1
m X log(yj − γ )
1
yj − γ
j =1
Z > yj ; γ , θ ∗ (h+1) Z −γ 1
(3.6)
(µ(h+1) , σ(h+1) ). The conditional expectations in the above expression (3.6) can be obtained as follows: σ2 (h+1) 1 ∗ −µ(h+1) 2 P0j(h+1) (γ ) E Z > yj ; γ , θ (h+1) = e Z −γ 2 σ (h+1) log(Z − γ ) −µ(h+1) Z > yj ; γ , θ ∗ 2 E σ(h+1) P0j(h+1) (γ ) + (µ(h+1) − σ(2h+1) )P0j(h+1) (γ ) , (h+1) = e Z −γ
where
where log yj − γ − µ(h+1)
Ψ0j(h+1) (γ ) =
Ψ0j (γ , θ ∗(h+1) )
=
P0j(h+1) (γ ) = P0j (γ , θ ∗(h+1) ) =
σ(h+1)
,
¯ (Ψ0j(h+1) (γ ) + σˆ (h+1) ) Φ . ¯ (Ψ0j(h+1) (γ )) Φ
The final MLEs of the parameters are local MLEs corresponding to the EM algorithm, and would be denoted here by LMLE2. In order to obtain LMLE2, we suggest a variation of the above method. As mentioned before, one can start with a γ(0) which is less than y1 and get µ(0) = µ(γ(0) ) and σ(0) = σ (γ(0) ) by using (3.4) and then one gets µ(1) and σ(1) by using (3.5). The suggestion is to treat (3.6) as the test equation. If the left hand side of that equation equals zero when (γ(0) , µ(1) , σ(1) ) substituted there, then no further iteration is needed. Otherwise, the suggested method involves choosing another γ(1) which is less than y1 and repeating the procedure looking for the sign change in (3.6). Finally, interpolating on γ values, one gets the final estimate (γˆ , µ, ˆ σˆ ) which when substituted in (3.6) should produce zero or close to zero. 4. Asymptotic variances and covariances of the EM algorithm estimates Asymptotic variances and covariances of the MLEs when the EM algorithm is used can be obtained by using the missing information principle of Louis (1982) and Tanner (1993). This principle is basically Observed information = Complete information − Missing information. Based on this principle, Louis (1982) developed the procedure of finding the observed information matrix when EM algorithm is used to find MLEs in an incomplete data situation. We adopt this principle in the situation of progressive
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
3585
Type II right censoring. We will denote the complete, observed and missing (expected) information by I (θ), IY (θ) and IZ |Y (θ), respectively. The complete information I (θ) is given by I (θ) = −E
∂ 2 E [log L(Y, Z, θ)|Y = y]
(4.1)
∂θ 2
in which E [log L(Y, Z, θ)|Y = y] is given by (3.1). In (4.1), the expectation is taken with respect to both Y and Z. The Fisher information matrix for a single observation which is censored at the time of the j-th failure is given by (j) IZ |Y (θ)
" = −E
∂ 2 log fZjk |Yj (zjk |yj ; θ)
# (4.2)
∂θ 2 (j)
in which fZjk |Yj (zjk |yj ; θ) is given by (3.3). In (4.2), the expectation is taken with respect to zjk so that IZ |Y (θ) is a function of yj and θ . Then, the missing (expected) information is simply IZ|Y (θ) =
m X
(j)
Rj IZ |Y (θ),
(4.3)
j =1
(j)
where IZ |Y (θ) is given by (4.2). The observed information IY (θ) is then obtained as follows: IY (θ) = I (θ) − IZ|Y (θ),
(4.4)
where I (θ) and IZ|Y (θ) are given by (4.1) and (4.3), respectively, and are derived below. Finally, upon inverting the observed information matrix IY (θ) in (4.4), one gets the asymptotic variances and covariances of the MLEs when EM algorithm is used. 4.1. Complete information matrix I (θ) The log-likelihood function log L∗ (θ) based on n uncensored observations yi , i = 1, . . . , n, is given by log L∗ (θ) =
n X
log f (yi ; θ),
(4.5)
i=1
where f (yi ; θ) is as given in (1.1). On differentiating the log-likelihood in (4.5) and equating to zero, one obtains the estimating equations given by (2.2) with m = n and Rj = 0; j = 1, 2, . . . , m. Negative of the second derivatives of the log-likelihood function log L∗ (θ) are obtained by appropriately differentiating the first derivatives in (2.2) with m = n and Rj = 0, j = 1, 2, . . . , m, and are given by (2.4) with m = n and Rj = 0, j = 1, 2, . . . , m. For the three-parameter lognormal distribution, it can be shown that
1
σ2
= e 2 −µ Yi − γ 1 2 E = e2(σ −µ) (Yi − γ )2 log(Yi − γ ) σ2 E = (µ − σ 2 )e 2 −µ Yi − γ log(Yi − γ ) 2 E = (µ − 2σ 2 )e2(σ −µ) . (Yi − γ )2
E
(4.6)
One gets the complete information matrix I (θ) in (4.1), by using (2.4) with m = n, Rj = 0, j = 1, 2, . . . , m, and (4.6), as follows:
1 + σ2
2(σ 2 −µ)
σ2 e 1 σ 2 −µ I (θ) = n e 2 σ2 2 2 σ − e 2 −µ σ
1
σ2
e
σ 2 −µ 2
1
2
− e σ
σ 2 −µ
2
σ2
0
0
σ2
2
.
(4.7)
Expressions for the asymptotic variances and covariances of the MLE of θ in the uncensored case are obtained by inverting
−1 2 the matrix I (θ) in (4.7). Denoting β = eµ , w = eσ and H = w(1 + σ 2 ) − (1 + 2σ 2 ) , the asymptotic variances and
3586
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
covariances are:
σ 2 β2 H n w σ2 V (µ) ˆ = (1 + H )
n 2 σ 2 V (σˆ ) = (1 + 2σ H ) 2n 2 σ β Cov(γˆ , µ) ˆ = − √ H n w σ3 β Cov(γˆ , σˆ ) = √ H n w 3 σ Cov(µ, ˆ σˆ ) = − H . V (γˆ ) =
n It is worthwhile to mention here that, for the uncensored case, Cohen (1951) has obtained the asymptotic variances and covariances of the MLEs of the parameters (γ , β, σ ), while those corresponding to the parameters (γ , µ, σ 2 ) were obtained by Hill (1963). 4.2. Missing information matrix IZ|Y (θ) The logarithm of the density function of an observation zjk = z censored at yj , the time of the j-th failure, is given by [see Eq. (3.2)]
log fz |yj (z |yj ; θ) = const . − log σ − log(z − γ ) − log 1 − Φ
log(yj − γ ) − µ
σ
−
1 2σ 2
[log(z − γ ) − µ]2 . (4.8)
Differentiating (4.8) with respect to γ , µ and σ , one gets
∂ log fz |yj ∂γ ∂ log fz |yj ∂µ ∂ log fz |yj ∂σ
= = =
1 log(z − γ ) − µ
σ 2 1
z−γ log(z − γ ) − µ
σ
σ
1
σ
"
+
1 z −γ
−
1
1
σ yj − γ
L0j ,
− L0j ,
log(z − γ ) − µ
2
σ
− 1 + ψ0j L0j
#
.
(4.9)
By using the properties of the left-truncated log-normal distribution, it can then be shown that E log(Z − γ ) − µ|Z > yj ; θ = σ L0j ,
E {log(Z − γ ) − µ}2 |Z > yj ; θ = σ 2 1 + Ψ0j L0j , E {log(Z − γ ) − µ}3 |Z > yj ; θ = σ 3 2 + Ψ0j2 , E {log(Z − γ ) − µ}4 |Z > yj ; θ = σ 4 3 1 + Ψ0j L0j + Ψ0j3 L0j , 1 σ2 E |Z > yj ; θ = e 2 −µ LL1j , Z −γ " # 2 1 2 E Z > yj ; θ = e2(σ −µ) LL2j , Z −γ 2 log(Z − γ ) Z > yj ; θ = e σ2 −µ σ L1j + (µ − σ 2 )LL1j , E Z −γ " # 2 log(Z − γ ) 2 E Z > yj ; θ = e2(σ −µ) σ 2 Ψ2j + 2σ (µ − 2σ 2 ) L2j + σ 2 + (µ − 2σ 2 )2 LL2j , Z −γ log(Z − γ ) Z > yj ; θ = e2(σ 2 −µ) σ L2j + (µ − 2σ 2 )LL2j , E 2 (Z − γ ) σ2 (log(Z − γ ))2 E Z > y ; θ = e 2 −µ σ 2 Ψ1j + 2σ (µ − σ 2 ) L2j + σ 2 + (µ − σ 2 )2 LL2j j Z −γ σ2 (log(Z − γ ))3 E Z > y ; θ = e 2 −µ σ 2 Ψ1j2 + 2σ 3 − 3σ (µ − σ 2 )2 L1j j Z −γ + 3(µ − σ 2 )σ 2 Ψ1j + 6σ (µ − σ 2 )2 L2j − 2(µ − σ 2 )3 LL1j + 3(µ − σ 2 ) σ 2 + (µ − σ 2 )2 LL2j ,
(4.10)
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
3587
in which Ψ0j , Ψ1j , Ψ2j , L0j , L1j , L2j , LL1j and LL2j are given by log yj − γ − µ
Ψ0j =
σ
log(yj − γ ) − (µ − σ 2 )
, σ 2 log(yj − γ ) − (µ − 2σ ) Ψ2j = , σ φ(Ψ1j ) φ(Ψ0j ) , L1j = , L0j = ¯ (Ψ0j ) ¯ (Ψ0j ) Φ Φ ¯ (Ψ1j ) ¯ (Ψ2j ) Φ Φ LL1j = , LL2j = . ¯ (Ψ0j ) ¯ (Ψ0j ) Φ Φ Ψ1j =
,
φ(Ψ2j ) L2j = , ¯ (Ψ0j ) Φ
(4.11)
(j)
The Fisher information matrix based on one observation zjk censored at yj , IZ |Y (θ), in (4.2), can be obtained using (4.9)–(4.11) as follows:
∂ log fZ |yj
2 Ψ e2(σ −µ) 2 + σ2j2
2
1
2(σ 2 −µ)
1
2e
σ 2 −µ 2
− L2j + e 1 + 2 LL2j − 2 L0j L1j σ 2 (yj − γ ) σ σ σ (yj − γ ) ∂ log fZ |yj 2 1 E = 2 1 + Ψ0j L0j − L20j ∂µ σ ∂ log fZ |yj 2 1 E = 2 2 + Ψ0j L0j 1 − Ψ0j L0j + Ψ0j2 ∂σ σ 2 ∂ log fZ |yj ∂ log fZ |yj 1 σ = 3 e 2 −µ σ (σ + 1)L1j − σ 2 (σ + 1)LL1j E ∂γ ∂µ σ σ L0j + 1 − σ 2 −σ L0j LL1j − L0j yj − γ " ∂ log fZ |yj ∂ log fZ |yj 1 2µ (σ − 1)L0j + µ = L0j E ∂γ ∂σ σ σ 2 (yj − γ ) ) ( 2µ2 + Ψ1j2 1 − Ψ0j L0j + 4µ 2µ2 σ2 + e 2 −µ − 3 + + − 3σ L1j 2 σ σ σ n µ µ o µ µ σ2 − − 1 − 4 − 1 σ + µ − + e 2 −µ 2Ψ1j L2j σ 2 σ 2 σ2 σ 2 2 2 µ − σ 2 − µσ µ µ µ σ σ −µ 2 2 2 2 2 + 2e 2 −µ − 1 − − σ + µ LL + 2e σ + (µ − σ ) LL 1j 2j σ2 σ σ2 σ4 ∂ log fZ |yj ∂ log fZ |yj 1 = 2 L0j + Ψ0j L0j (Ψ0j − L0j ) . E ∂µ ∂σ σ E
∂γ
=
L2 2 0j
(4.12)
Using IZ|Y (θ) from (4.12) and I (θ) from (4.7), the observed information can be obtained from (4.4). Finally, one can get the asymptotic variances and covariances of the MLEs, when the EM algorithm is used, by inverting this observed information matrix IY (θ). 5. Alternative estimators As an alternative estimator, one can use modified MLE (MMLE). MMLEs use first order statistic for estimating γ and are therefore easier to compute than the MLEs. For small samples when LMLEs fail to converge, these MMLEs (which always exist) produce reasonable estimates. For progressively censored data, one can have two versions of MMLEs. One version corresponds to the Newton–Raphson algorithm and the other version corresponds to the EM algorithm, and they will be ∂ log L denoted by MMLE1 and MMLE2, respectively. In both versions, the likelihood equation ∂γ = 0 is replaced by
γ + eµ+σ Φ
−1 [ 1 ] n+1
Pm
= y1
(5.1)
in which n = m + j=1 Rj is the total sample size. Eq. (5.1) is used as a test equation and it replaces the first equation in (2.2) for the Newton–Raphson algorithm version and Eq. (3.6) for the EM algorithm version. One starts with a first approximation γ1 and make the transformation log(yi − γ1 ). One then calculates conditional estimates µ1 = µ(γ1 ) and σ1 = σ (γ1 ) by using the second and third equations of (2.2) for the Newton–Raphson algorithm and (3.5) for the EM algorithm. The values (γ1 , µ1 , σ1 ) are substituted in Eq. (5.1). If the test equation is satisfied, then no further iteration is required. Otherwise, a
3588
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
Table 1 Progressively censored data for Example 1. i
1
2
3
4
5
6
7
8
9
10
11
12
Yi Ri
0.265 0
0.269 0
0.297 2
0.315 0
0.338 0
0.379 2
0.392 0
0.402 0
0.412 0
0.416 0
0.418 0
0.484 4
Table 2 Estimates and standard deviations of estimates for γ , µ and σ based on data in Table 1.
γ Estimator LMLE2
µ
γˆ
σγˆ
0.16
0.09
σ
µ ˆ −1.21
σµˆ
σˆ
σσˆ
0.19
0.29
0.05
Table 3 Progressively censored data for Example 2. i
1
2
3
4
5
6
7
Yi Ri
152.7 0
172.0 0
172.5 0
173.3 1
193.0 0
204.7 0
234.9 2
Table 4 Estimates and standard deviations of estimates for γ , µ and σ based on data in Table 3.
γ Estimator LMLE2
µ
σ
γˆ
σγˆ
µ ˆ
σµˆ
σˆ
σσˆ
110.30
10.29
3.35
0.89
0.81
0.12
second approximation γ2 is selected and the cycle of calculations as described above is repeated. The iterations are continued until two sufficiently close values γi and γi+1 are found such that the following is satisfied:
γi + eµi +σi Φ
−1 [ 1 ] n+1
< (>)y1 < (>)γi+1 + eµi+1 +σi+1 Φ
−1 [ 1 ] n+1
.
Thus, the final MMLE γ˜ of γ is obtained using which the final MMLEs µ ˜ , σ˜ of µ and σ are then obtained. 6. Illustrative examples To illustrate the computational methods presented in this article, we use three examples given in Cohen and Whitten (1988). We modified these original examples to consider progressively censored data. We then assume that each datum comes from a three-parameter lognormal distribution, and use them to carry out all the estimation procedures discussed in the preceding sections in order to see whether they produce similar results or not. Also, in Section 7, we compare the coverage probabilities of confidence intervals based on LMLE1 and LMLE2 using a simulation study. We used same γ0 value for obtaining LMLE1 and LMLE2. For the determination of MMLE, as discussed in Section 5, Φ −1 [ n+1 1 ] = −1.67 for Examples 1 and 3 since n = 20 in both these cases. For Example 2, Φ −1 [ n+1 1 ] = −1.34 since n = 10 here. Proceeding as explained in Section 5, we get two versions of MMLE, viz., MMLE1 and MMLE2 (for the Newton–Raphson algorithm and EM algorithm, respectively). Example 1. The maximum flood levels (in millions of cubic feet per second) for 20 four-year periods from 1890 to 1969 in the Susquehanna river at Harrisburg, Pennsylvania are given in Cohen and Whitten (1988) and were also used by Dumonceaux and Antle (1973). We modified these data to make it progressively censored with m = 12 stages, and these progressively censored data are presented in Table 1. The LMLE2 and their standard deviations are presented in Table 2. For the data in Table 1, the LMLE1 are given by γˆ = 0.1822, µ ˆ = −1.291 and σˆ = 0.362, the MMLE1 are given by γ˜ = 0.1581, µ ˜ = −1.103 and σ˜ = 0.279, and the MMLE2 are given by γ˜ = 0.1563, µ ˜ = −1.1725 and σ˜ = 0.273. Example 2. This example was used by McCool (1974) and is also given in Cohen and Whitten (1988). The data are fatigue lives (in hours) of 10 bearings of a certain type. We modified these data to make it progressively censored with m = 7 stages, and these progressively censored data are presented in Table 3. LMLE1 did not converge for this sample. The LMLE2 and their standard deviations are presented in Table 4. The MMLE1 are given by γ˜ = 112.7815, µ ˜ = 3.9231 and σ˜ = 0.8932, while MMLE2 are given by γ˜ = 109.7232, µ ˜ = 3.2763 and σ˜ = 0.8036.
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
3589
Table 5 Progressively censored data for Example 3. i
1
2
3
4
5
6
7
8
9
10
11
12
Yi Ri
127.2 0
128.7 0
131.4 0
133.0 2
133.1 0
135.9 0
144.3 2
145.8 0
148.3 0
153.1 0
157.3 0
164.3 4
Table 6 Estimates and standard deviations of estimates for γ , µ and σ based on data in Table 5.
γ Estimator LMLE2
µ
σ
γˆ
σγˆ
µ ˆ
σµˆ
σˆ
σσˆ
108.25
11.27
3.82
0.86
0.39
0.01
Table 7 Progressively censored data for Example 3 with PCS-1. i
1
2
3
4
5
6
7
8
9
10
11
12
Yi Ri
127.2 0
128.7 0
131.4 0
133.0 0
133.1 0
135.9 0
137.3 0
144.3 0
145.8 0
148.3 2
153.1 2
157.2 4
Table 8 Progressively censored data for Example 3 with PCS-2. i
1
2
3
4
5
6
7
8
9
10
11
12
Yi Ri
127.2 4
128.7 2
133.0 2
135.9 0
144.3 0
148.3 0
153.1 0
157.2 0
166.5 0
174.8 0
184.1 0
201.4 0
Table 9 Estimates and standard deviations of estimates for γ , µ and σ for data in Table 7.
γ Estimator LMLE2
µ
σ
γˆ
σγˆ
µ ˆ
σµˆ
σˆ
σσˆ
102.11
11.79
4.01
1.38
0.36
0.03
Table 10 Estimates and standard deviations of estimates for γ , µ and σ for data in Table 8.
γ Estimator LMLE2
µ
σ
γˆ
σγˆ
µ ˆ
σµˆ
σˆ
σσˆ
108.37
11.99
4.12
0.6
0.39
0.02
Example 3. In this example, the complete data are given by Cohen and Whitten (1988) and were also used by Cohen (1951). The complete data consists of 20 observations from a three-parameter lognormal distribution with γ = 100, µ = 3.912023 and σ = 0.4. We modified these data to make it progressively censored with m = 12 stages, and these progressively censored data are presented in Table 5. Although Yi values are reported up to one decimal place in Tables 5, 7 and 8, we used all three decimal places for these Yi values in Cohen (1951) for our computations. The LMLE2 and their standard deviations are presented in Table 6. For the data in Table 5, the LMLE1 are given by γˆ = 110.3112, µ ˆ = 3.573 and σˆ = 0.422, the MMLE1 are given by γ˜ = 109.1563, µ ˜ = 4.419 and σ˜ = 0.295, and the MMLE2 are given by γ˜ = 107.8036, µ ˜ = 3.927 and σ˜ = 0.378. In order to examine the effect of delayed censoring, we considered two schemes of progressive censoring. In one scheme, the censoring was delayed than the other. Keeping the number of stages m = 12 the same as in Table 4, for the delayed censoring scheme, we considered PCS-1: R1 = 0, R2 = 0, R3 = 0, R4 = 0, R5 = 0, R6 = 0, R7 = 0, R8 = 0, R9 = 0, R10 = 2, R11 = 2, R12 = 4. For the other scheme, we took PCS-2: R1 = 4, R2 = 2, R3 = 2, R4 = 0, R5 = 0, R6 = 0, R7 = 0, R8 = 0, R9 = 0, R10 = 0, R11 = 0, R12 = 0. The progressively censored data obtained under these two schemes are presented in Tables 7 and 8, respectively. The LMLE2 estimates and their standard deviations for the data in Tables 7 and 8 are presented in Tables 9 and 10, respectively. For the data in Table 7, the LMLE1 are given by γˆ = 104.4622, µ ˆ = 3.992 and σˆ = 0.373, the MMLE1 are given by γ˜ = 103.5689, µ ˜ = 4.151 and σ˜ = 0.298, and the MMLE2 are given by γ˜ = 103.8039, µ ˜ = 4.107 and σ˜ = 0.315. For the data in Table 8, the LMLE1 are given by γˆ = 109.3784, µ ˆ = 4.041 and σˆ = 0.413, the MMLE1 are given by γ˜ = 107.2501, µ ˜ = 4.123 and σ˜ = 0.292, and the MMLE2 are given by γ˜ = 107.9231, µ ˜ = 4.137 and σ˜ = 0.368. It is
3590
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
Table 11 95% confidence intervals for Example 1. Estimators
γ
µ
σ
LMLE1 LMLE2 MMLE1 MMLE2
(−.255, .647) (−.071, .415) (−.250, .574) (−.165, .423)
(−2.165, −.477) (−1.644, −.816) (−1.855, −.327) (−1.651, −.593)
(.085, .673) (.137, .423) (.047, .595) (.085, .437)
Table 12 95% confidence intervals for Example 2. Estimators
γ
µ
σ
LMLE2 MMLE1 MMLE2
(66.682, 152.942) (69.882, 157.815) (65.721, 139.061)
(1.689, 5.535) (1.284, 6.968) (1.232, 4.956)
(.323, 1.303) (.423, 1.482) (.403, 1.030)
Table 13 Coverage probabilities of 95% confidence intervals for γ based on Monte Carlo simulations. n
m
R
LMLE1
LMLE2
20 20 20 20 20 20 20 20 40 40 40 40 40 40 40 40
12 12 14 14 16 16 18 18 24 24 28 28 32 32 36 36
(7,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,7) (5,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,5) (3,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,3) (1,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,1) (14,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,14) (10,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,10) (6,2,0, . . . ,0,0,0) 0,0,0, . . . ,0,2,6) (2,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,2)
90.71 90.92 91.02 91.12 93.34 93.78 93.92 94.35 92.10 92.94 93.21 94.28 94.37 93.48 93.92 94.73
92.25 92.81 93.01 93.82 94.71 94.29 94.35 95.72 93.01 94.15 93.91 94.23 94.41 94.87 94.91 94.28
clear from Tables 6, 9 and 10 that the estimates of γ are affected not only by censoring but also by the pattern of censoring, while the estimates of µ and σ remain more or less the same in all the situations. 7. Simulation study In this section, we report results from two simulation studies that we carried out. In the first simulation study, we constructed 95% confidence intervals based on 10,000 randomly generated progressively censored samples corresponding to different estimates of γ , µ and σ for each of the estimation methods discussed in Section 6. We used the same progressive censoring schemes as given in those examples. We provide the confidence intervals from this simulation study in Tables 11 and 12 for Examples 1 and 2, respectively, for different estimators. In the second simulation study, we compared the performance of LMLE1 and LMLE2 in terms of coverage probabilities of 95% confidence intervals for the parameters γ , µ and σ for different sample sizes and different degrees of censoring. 10,000 samples were simulated from the lognormal distribution with γ = 100, µ = 3.912023 and σ = 0.4 with sample size n = 20, 40. For sample size n = 20, we considered m = 12, 14, 16, 18 stages of censoring, and for sample size n = 40 we considered m = 24, 28, 32, 36 stages. In each case, we considered two censoring schemes with one reflecting comparatively delayed censoring than the other. The coverage probabilities of 95% confidence intervals from this simulation study for the parameters γ , µ and σ are presented in Tables 13–15, respectively. It is observed from these tables that the coverage probabilities are better when the proportion of uncensored data is larger. The coverage probabilities seem to be almost the same for the two different censoring schemes. The confidence intervals for the parameter γ seems to be most sensitive to the censoring pattern, while the confidence intervals for the parameters µ and σ seem to be stable and quite satisfactory and close to the nominal level of 95%. In a similar tone, it is known that different censoring schemes do not change the estimates for the two-parameter exponential distribution (see for example, Balakrishnan and Sandhu (1996)). The sensitivity of estimation of γ may be due to the fact that γ is the threshold parameter of the lognormal distribution and the pattern of censoring applied to the sample might have affected the nature of the sample observations as far as the estimation of the threshold parameter is involved. Also, LMLE1 and LMLE2 seem to have nearly the same coverage probabilities.
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
3591
Table 14 Coverage probabilities of 95% confidence intervals for µ based on Monte Carlo simulations. n
m
R
LMLE1
LMLE2
20 20 20 20 20 20 20 20 40 40 40 40 40 40 40 40
12 12 14 14 16 16 18 18 24 24 28 28 32 32 36 36
(7,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,7) (5,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,5) (3,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,3) (1,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,1) (14,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,14) (10,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,10) (6,2,0, . . . ,0,0,0) 0,0,0, . . . ,0,2,6) (2,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,2)
93.12 93.28 93.27 94.02 94.96 95.02 94.81 95.38 93.24 93.38 94.11 94.43 95.02 95.15 95.39 95.58
93.94 93.82 93.73 93.61 95.11 95.25 95.42 95.46 94.12 94.22 94.09 94.42 95.36 95.37 95.38 96.01
Table 15 Coverage probabilities of 95% confidence intervals for σ based on Monte Carlo simulations. n
m
R
PC for LMLE1
PC for LMLE2
20 20 20 20 20 20 20 20 40 40 40 40 40 40 40 40
12 12 14 14 16 16 18 18 24 24 28 28 32 32 36 36
(7,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,7) (5,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,5) (3,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,3) (1,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,1) (14,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,14) (10,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,10) (6,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,6) (2,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,2)
92.94 93.17 93.25 93.74 94.82 94.29 95.03 95.41 93.15 93.21 93.36 93.81 95.01 95.21 95.11 95.14
93.66 93.79 93.73 93.85 94.21 94.62 95.43 95.44 93.62 93.79 93.36 93.76 95.21 95.01 95.35 95.91
8. Concluding remarks In this article, we have discussed the EM algorithm for the maximum likelihood estimation based on progressively Type-II censored samples from a three-parameter lognormal distribution. We have also considered the traditional Newton–Raphson method for this purpose. Additionally, we have discussed two versions of modified maximum likelihood estimators. The EM algorithm and the Newton–Raphson method produced similar results in all three examples considered, as well as in terms of coverage probabilities determined from a Monte Carlo simulation study. The coverage probabilities for both methods are better and closer to the nominal level of 95% when the proportion of uncensored data is larger. Acknowledgments The authors are thankful to the referees for their valuable comments which led to a considerable improvement in the presentation of this article. References Balakrishnan, N., 2007. Progressive censoring methodology: An appraisal (with discussions). Test 20, 211–296. Balakrishnan, N., Aggarwala, R., 2000. Progressive Censoring: Theory, Methods, and Applications. Birkhäuser, Boston. Balakrishnan, N., Cohen, A.C., 1991. Order Statistics and Inference: Estimation Methods. Academic Press, San Diego. Balakrishnan, N., Kannan, N., Lin, C.T., Ng, H.K.T., 2003. Point and interval estimation for Gaussian distribution based on progressively Type-II censored samples. IEEE Transactions on Reliability 52, 90–95. Balakrishnan, N., Sandhu, R.A., 1996. Best linear unbiased and maximum likelihood estimation for exponential distributions under general progressive type-II censored samples. Sankhya, Series B 58, 1–9. Calitz, F., 1973. Maximum likelihood estimation of the parameters of the three-parameter lognormal distribution — A reconsideration. Australian Journal of Statistics 15, 185–199. Cohen, A.C., 1951. Estimating parameters of logarithmic-normal distributions by maximum likelihood. Journal of the American Statistical Association 46, 206–212.
3592
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
Cohen, A.C., 1988. Estimation in the lognormal distributions. In: Crow, E.L., Shimizu, K. (Eds.), Lognormal Distributions: Theory and Applications. Marcel Dekker, New York. Cohen, A.C., Whitten, B.J., 1988. Parameter Estimation in Reliability and Life Span Models. Marcel Dekker, New York. Crow, E.L., Shimizu, K. (Eds.), 1988. Lognormal Distributions: Theory and Applications. Marcel Dekker, New York. Dempster, A.P., Laird, N.M., Rubin, D.B., 1977. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B 39, 1–38. Dumonceaux, R., Antle, C.E., 1973. Discrimination between the log-normal and the Weibull distributions. Technometrics 15, 923–926. Harter, H.L., Moore, A.H., 1966. Local-maximum-likelihood estimation of the parameters of the three-parameter lognormal populations from complete and censored samples. Journal of the American Statistical Association 61, 842–851. Corrections: 61 (1966) 1247; 62 (1967) 1519–1520; 63 (1968) 1549. Hill, B.M., 1963. The three-parameter lognormal distribution and Bayesian analysis of a point-source epidemic. Journal of the American Statistical Association 58, 72–84. Johnson, N.L., Kotz, S., Balakrishnan, N., 1994. Continuous Univariate Distributions – vol. 1, second ed. John Wiley & Sons, New York. Louis, T.A., 1982. Finding the observed information matrix when using the EM algorithm. Journal of the Royal Statistical Society, Series B 44, 226–233. McCool, J.I., 1974. Inferential techniques for Weibull populations. Aerospace Research Laboratories Report, ARL TR 74-0180, Wright-Patterson AFB, Ohio. McLachlan, G.J., Krishnan, T., 1997. The EM Algorithm and Extensions. Marcel Dekker, New York. Mi, J., Balakrishnan, N., 2003. Existence and uniqueness of the MLEs for normal distribution based on general progressively Type-II censored samples. Statistics & Probability Letters 64, 407–414. Munro, A.H., Wixley, R.A.J., 1970. Estimators based on order statistics of small samples from a three-parameter lognormal distribution. Journal of the American Statistical Association 65, 212–225. Nelson, W., 1982. Applied Life Data Analysis. John Wiley & Sons, New York. 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. Rukhin, A.L., 1984. Improved estimation in lognormal models. Technical Report No. 84-38, Department of Statistics, Purdue University, West Lafayette, Indiana. Tanner, M.A., 1993. Tools for Statistical Inference: Observed Data and Data Augmentation Methods, second ed. Springer-Verlag, New York.
Contents lists available at ScienceDirect
Computational Statistics and Data Analysis journal homepage: www.elsevier.com/locate/csda
Estimation for the three-parameter lognormal distribution based on progressively censored data Prasanta Basak a , Indrani Basak a,∗ , N. Balakrishnan b a b
Penn State Altoona, Altoona, PA, United States McMaster University, Hamilton, Canada
article
info
Article history: Received 27 February 2009 Accepted 27 March 2009 Available online 1 April 2009
a b s t r a c t Some work has been done in the past on the estimation of parameters of the threeparameter lognormal distribution based on complete and censored samples. In this article, we develop inferential methods based on progressively Type-II censored samples from a three-parameter lognormal distribution. In particular, we use the EM algorithm as well as some other numerical methods to determine maximum likelihood estimates (MLEs) of parameters. The asymptotic variances and covariances of the MLEs from the EM algorithm are computed by using the missing information principle. An alternative estimator, which is a modification of the MLE, is also proposed. The methodology developed here is then illustrated with some numerical examples. Finally, we also discuss the interval estimation based on large-sample theory and examine the actual coverage probabilities of these confidence intervals in case of small samples by means of a Monte Carlo simulation study. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Lognormal distribution is one of the distributions commonly used for modeling lifetimes or reaction-times, and is particularly useful for modeling data which are long-tailed and positively skewed. It has been discussed extensively by many authors including Cohen (1951, 1988), Hill (1963), Harter and Moore (1966), Crow and Shimizu (1988), Johnson et al. (1994), Munro and Wixley (1970) and Rukhin (1984). It is well-known that there is a close relationship between normal and lognormal distributions. If X = log(Y − γ ) is normally distributed with mean µ and standard deviation σ , then the distribution of Y becomes a three-parameter lognormal distribution with parameter θ = (γ , µ, σ ). The probability density function of such a three-parameter lognormal distribution is
[log(y − γ ) − µ]2 √ 1 exp − , f (y; θ) = σ 2π (y − γ ) 2σ 2 0, otherwise.
γ < y < ∞, σ > 0, −∞ < µ < ∞
(1.1)
In (1.1), σ 2 and µ are the variance and mean of the underlying normal variable X , but become the shape and scale parameters of the lognormal variable Y . It is more convenient to use β = exp(µ) and w = exp(σ 2 ) as the scale and shape parameters of the lognormal variable Y , respectively. Also, γ is the threshold (location) parameter of the lognormal random variable.
∗
Corresponding author. Tel.: +1 814 9495263; fax: +1 814 949 5456. E-mail address: [email protected] (I. Basak).
0167-9473/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.csda.2009.03.015
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
3581
When the threshold parameter γ is known, the parameter estimation can be done using the well-known results for normal distribution by making the transformation from Y to X . Estimation methods become more complex when γ is unknown. Adding to that complexity is the fact that quite often the data on lifetimes or reaction-times come with ‘‘censoring.’’ Censoring occurs when exact survival times are known only for a portion of the individuals or items under study. The complete survival times may not have been observed by the experimenter either intentionally or unintentionally, and there are numerous examples for each kind; see, for example, Nelson (1982) and Balakrishnan and Cohen (1991). In this article, we consider a general scheme of progressive Type-II right censoring. Under this scheme, n units are placed on a life-testing experiment and only m (< n) are completely observed until failure. The censoring occurs progressively in m stages. These m stages offer failure times of the m completely observed units. At the time of the first failure (the first stage), R1 of the n − 1 surviving units are randomly withdrawn (censored intentionally) from the experiment, R2 of the n − 2 − R1 surviving units are withdrawn at the time of the second failure (the second stage), and so on. Finally, at the time of the m-th failure (the m-th stage), all the remaining Rm = n − m − R1 − · · · − Rm−1 surviving units are withdrawn. This scheme (R1 , R2 , . . . , Rm ) is referred to as progressive Type-II right censoring scheme. It is clear that this scheme includes the conventional Type-II right censoring scheme (when R1 = R2 = · · · = Rm−1 = 0 and Rm = n − m) and complete sampling scheme (when n = m and R1 = R2 = · · · = Rm = 0). The ordered lifetime data which arise from such a progressive Type-II right censoring scheme are called progressively Type-II right censored order statistics. For theory, methods and applications of progressive censoring, readers are referred to the book by Balakrishnan and Aggarwala (2000) and the recent discussion paper by Balakrishnan (2007). Suppose n independent lognormally distributed units are placed on a life-testing experiment. Let Y1:m:n ≤ · · · ≤ Ym:m:n denote the above mentioned m progressively Type-II right censored order statistics. For ease in notation, let us use Yj (j = 1, . . . , m) to denote these Yj:m:n ’s. Note that we observe only Y = (Y1 , . . . , Ym ). The purpose of this article is to discuss different estimation procedures for the parameter θ based on the progressively Type-II right censored order statistics Y. Recently, Ng et al. (2002) considered two-parameter lognormal distribution (which does not include the threshold parameter) and discussed Newton–Raphson algorithm as well as EM algorithm for finding the MLEs. Inclusion of the third parameter γ introduces an unusual feature in the likelihood function. Hill (1963) has shown that there exist paths along which the likelihood function of any ordered sample y1 , . . . , yn from the three-parameter lognormal distribution tends to ∞ as (γ , µ, σ ) approaches (y1 , −∞, ∞). Thus, the global maximization of the likelihood leads to the unreasonable estimate (y1 , −∞, ∞) although, in fact, the likelihood at the point is zero. But Hill (1963) also retained the idea that reasonable estimates could be obtained by solving the likelihood equations. To get these reasonable estimates, Cohen (1951), Cohen and Whitten (1988) and Harter and Moore (1966) equated partial derivatives of the likelihood function to zero and solved the resulting equations. These estimates are called local maximum likelihood estimates (LMLEs). Harter and Moore (1966) and Calitz (1973) noted that although these LMLEs are not true MLEs according to the usual definition, they are reasonable estimates and appear to possess most of the desirable properties associated with MLEs. However, it is noted in Harter and Moore (1966) and Cohen and Whitten (1988) that sometimes the likelihood function may have no clearly defined local maximum for small samples and so LMLEs fail to produce estimates in that case. In the past, some work has been done on estimation methods for the three-parameter lognormal distribution based on complete and censored samples; see, for example, the books by Cohen and Whitten (1988) and Balakrishnan and Cohen (1991). In Section 2, we describe the Newton–Raphson algorithm for determining the MLEs of the parameter θ based on a progressively censored sample. The second derivatives of the log-likelihood are required in order to use the algorithm. These computations are complicated when data are progressively censored. Another viable alternative to Newton–Raphson algorithm is the well-known EM algorithm and in Section 3 we discuss how that can be used to determine the MLEs in this case. Asymptotic variances and covariances of the maximum likelihood estimates generated through the EM algorithm are given in Section 4. Section 5 describes some simplified estimation methods which yields simple alternative estimators. All these methods discussed in this article are then illustrated with some numerical examples in Section 6. With these illustrative examples, we also discuss the interval estimation based on large-sample theory and then examine the actual coverage probabilities in case of small samples through Monte Carlo simulations. 2. Newton–Raphson algorithm One of the standard methods of determining the maximum likelihood estimates is the Newton–Raphson algorithm. In this section, we describe the Newton–Raphson algorithm for finding the MLEs numerically when life-times are distributed as a three-parameter lognormal distribution with parameter θ . These MLEs are local MLEs, as mentioned earlier, corresponding to the Newton–Raphson algorithm, and would be denoted here by LMLE1. The log-likelihood function log L(θ) = log L based on the progressively Type-II right censored order statistics Y, is log L(θ) = const . − m log σ −
m X j =1
log(y −γ )−µ
log(yj − γ ) −
m 1X
2 j =1
Ψ0j2 +
m X
¯ Ψ0j Rj log Φ
,
(2.1)
j =1
j ¯ = 1−Φ where Ψ0j = ; see Balakrishnan and Aggarwala (2000). In (2.1) and throughout this article, φ and Φ σ denote the probability density and survival function of the standard normal distribution, respectively. Three likelihood
3582
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
equations which need to be solved simultaneously for the required estimate θˆ = (µ, ˆ σˆ , γˆ ) are as follows: m m X 1 ∂ log L 1 X 1 = Ψ0j + Rj L0j + =0 ∂γ σ j=1 yj − γ y − γ j j =1 " # m m X 1 X ∂ log L = Ψ0j + Rj L0j = 0 ∂µ σ j =1 j =1 " # m m X X ∂ log L 1 2 = −m + Ψ0j + Rj Ψ0j L0j = 0 ∂σ σ j =1 j =1
(2.2)
¯ (Ψ0j ) is the hazard function of the standard normal distribution at Ψ0j . in which L0j = φ(Ψ0j )/Φ In the Newton–Raphson algorithm, the simultaneous solution is obtained through an iterative procedure. In each iterative step, the corrections a, b, c to the previous estimates γ0 , µ0 , σ0 produce new estimates γˆ , µ, ˆ σˆ as γˆ = γ0 + a,
µ ˆ = µ0 + b and σˆ = σ0 + c .
The iteration method is based on Taylor series expansions of the estimating equations in (2.2) in the neighborhood of the previous simultaneous estimates. Neglecting powers of a, b and c above the first order and using Taylor’s theorem, we get the following equations which need to be solved for a, b, c:
∂ 2 log L ∂ 2 log L + b +c ∂γ 2 0 ∂γ ∂µ 0 ∂ 2 log L ∂ 2 log L +b +c a ∂µ∂γ 0 ∂µ2 0 ∂ 2 log L ∂ 2 log L a + b +c ∂σ ∂γ ∂σ ∂µ
a
0
0
∂ 2 log L ∂ log L = − ∂γ ∂σ 0 ∂γ 0 ∂ 2 log L ∂ log L =− ∂µ∂σ 0 ∂µ 0 ∂ 2 log L ∂ log L =− , ∂σ 2 0 ∂σ 0
(2.3)
where the notation A|0 , for any partial derivative A, means the partial derivative evaluated at (γ0 , µ0 , σ0 ). The second derivatives needed in (2.3) are as follows: m m m ∂ 2 log L X 1 1 X 1 1 X 1 = + ( Ψ + R L ) − (1 − Rj Ψ0j L0j + Rj L20j ) 0j j 0j 2 2 2 ∂γ 2 ( y − γ ) σ ( y − γ ) σ ( y − γ )2 j j j j=1 j =1 j =1 " # m X ∂ 2 log L 1 = − m + Rj L0j (L0j − Ψ0j ) ∂µ2 σ2 j =1 " # m X 2 ∂ 2 log L 1 3 2 2 = 2 m− 3Ψ0j + Rj (2Ψ0j L0j − Ψ0j L0j + Ψ0j L0j ) ∂σ 2 σ j =1 m 1 1 X ∂ 2 log L =− 2 1 − Rj Ψ0j L0j + Rj L20j ∂γ ∂µ σ j =1 yj − γ m 1 ∂ 2 log L 1 X =− 2 2Ψ0j + Rj (L0j − Ψ0j2 L0j + Ψ0j L20j ) ∂γ ∂σ σ j =1 yj − γ m ∂ 2 log L 1 X =− 2 2Ψ0j + Rj (L0j − Ψ0j2 L0j + Ψ0j L20j ) . ∂µ∂σ σ j =1
(2.4)
Using the Newton–Raphson method, Ng et al. (2002) discussed methods of finding MLEs while Balakrishnan et al. (2003) discussed the construction of confidence intervals of µ and σ for the two-parameter lognormal distribution. Mi and Balakrishnan (2003) have made use of the fact that the lognormal density is log-concave in order to establish that the MLEs of µ and σ do exist and are unique. Hence, in that case, the EM algorithm and the Newton–Raphson algorithm will converge to the same values. In order to obtain LMLE1 for the three-parameter lognormal distribution, we suggest a variation of the above method, which was used by Cohen (1951). Calitz (1973) found this method of Cohen (1951) produced better convergence in the estimation process. Instead of solving all three equations in (2.2) simultaneously, the suggested method is to start with a γ(0) which is less than y1 and use the second and third equations of (2.2) to get µ(0) = µ(γ(0) ) and σ(0) = σ (γ(0) ). The first equation is used as the test equation. If the left hand side of this equation equals zero when (γ(0) , µ(0) , σ(0) ) substituted there, then no further iteration is required. Otherwise, another γ(1) is chosen which is less than y1 yielding µ(1) = µ(γ(1) )
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
3583
and σ(1) = σ (γ(1) ) and look for the sign change in the left hand side of the first equation of (2.2). Finally, interpolating on γ values, first find the final estimate γˆ and then the final estimates µ ˆ = µ(γˆ ) and σˆ = σ (γˆ ) are obtained using the second and third equations of (2.2). The final estimates (γˆ , µ, ˆ σˆ ) when substituted in the first equation of (2.2) should produce zero or close to zero. As noted by Cohen and Whitten (1988), usually a single root γˆ will be found. In the event that multiple roots of γˆ occur, Cohen and Whitten (1988) suggested to use the root which results in the closest agreement between y¯ and 2
ˆ σˆ2 E[ [Y ] = γˆ + eµ+ .
3. EM Algorithm In the case of progressively censored samples, an alternative to the Newton–Raphson algorithm is the use of EM algorithm for numerically finding the MLEs. One advantage of the EM algorithm is that asymptotic variances and covariances of the EM algorithm estimates can be computed, which is discussed in Section 4. EM algorithm, introduced by Dempster et al. (1977), is a very popular tool to handle any missing or incomplete data situation; readers are referred to the book by McLachlan and Krishnan (1997) for a detailed discussion on EM algorithm and its applications. This algorithm is an iterative method which has two steps. In the E-step, it replaces any missing data by its expected value and in the M-step the log-likelihood function is maximized with the observed data and expected value of the incomplete data, producing an update of the parameter estimates. The MLEs of the parameters are obtained by repeating the E- and M-steps until convergence occurs. Since progressive censoring model can be viewed as a missing data problem, EM algorithm can be applied to obtain the MLEs of the parameters in this case. Let us denote the censored data vector as Z = (Z1 , Z2 , . . . , Zm ), where the j-th stage censored data vector Zj is a 1 ×Rj vector, Zj = (Zj1 , Zj2 , . . . , ZjRj ), for j = 1, 2, . . . , m. The complete data set is then obtained by combining the observed data Y and the censored data Z. E-step of the algorithm requires the computation of the conditional expectation of functions of censored data vector Z, conditional on the observed data vector Y and the current value of the parameters. In particular, one computes the conditional expectation of the log-likelihood E [log L(Y, Z, θ )|Y = y] as E [log L(Y, Z, θ )|Y = y] = const . − n log σ −
m X
log(yj − γ ) −
m 1X
log(yj − γ ) − µ
2
2 j =1 σ " # 2 Rj Rj m X m X X 1X log(Zjk − γ ) − µ E log(Zjk − γ ) | Zjk > yj − E − | Zjk > yj . 2 j=1 k=1 σ j=1 k=1 j =1
(3.1)
The above conditional expectations are obtained using the result that given Yj = yj , Zj ’s have a left-truncated distribution F , truncated at yj . More specifically, the conditional probability density of Z, given Y, is given by (see Balakrishnan and Aggarwala (2000))
fZ|Y (z|y; θ) =
Rj m Y Y
fZjk |Yj (zjk |yj ; θ),
(3.2)
j=1 k=1
where fZjk |Yj (zjk |yj ; θ) =
f (zjk ; θ)
(3.3)
1 − F (yj ; θ)
and f (zjk ; θ) is given by (1.1) and F denotes the corresponding cumulative distribution function. In the M-step of the (h + 1)-th iteration, we will denote the updated estimates of the parameter θ as θ (h+1) . This θ (h+1) maximizes the loglikelihood function involving the observed data Y, conditional expectation of the log-likelihood function of censored data vector Z given the observed data vector Y, and the h-th iteration value of the parameter θ (h) . As a starting value θ (0) , one can use a γ(0) < y1 and µ(0) and σ(0) are computed on the basis of the so-called ‘‘pseudo-complete’’ sample which involves observed data Y and the censored observations at the j-th step Zj all taken to be yj . Thus, µ(0) = µ(γ(0) ) and σ(0) = σ (γ(0) ) are then given by
µ(0) = σ(0)
m 1X
n
(Rj + 1) log yj − γ(0)
1 v j=" # . u m u1 X 2 2 t = (Rj + 1) log yj − γ(0) − µ(0) n j=1
(3.4)
3584
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
Starting with the initial estimates in (3.4), the (h + 1)-th iteration value θ (h+1) is obtained using the h-th iteration value θ h as follows:
µ(h+1) = σ(h+1)
1
( m X
n
log yj − γ(h) +
m X
Rj E log(Z − γ(h) )|Z > yj ; θ (h)
)
j =1
j=1
v ( ) u m m u1 X X 2 2 t = Rj E log (Z − γ(h) )|Z > yj ; µ(h+1) , σ(h) , γ(h) − µ2(h+1) log yj − γ(h) + n
j =1
j =1
.
(3.5)
The conditional expectations in the above expression (3.5) can be obtained as follows:
)
E log(Z − γ(h) )|Z > yj ; θ (h) = σ(h) L0j(h) + µ(h)
E log2 (Z − γ(h) )|Z > yj ; µ(h+1) , σ(h) , γ(h) = σ(2h) 1 + Ψ0j∗ L∗0j + 2σ(h) µ(h+1) L∗0j + µ2(h+1)
where log yj − γ(h) − µ(h)
Ψ0j(h) = Ψ0j (θ (h) ) =
σ(h)
,
φ(Ψ0j(h) ) , ¯ (Ψ0j(h) ) Φ log yj − γ(h) − µ(h+1) Ψ0j∗ = , σ(h) φ(Ψ0j∗ ) L∗0j = . ¯ (Ψ0j∗ ) Φ L0j(h) = L0j (θ (h) ) =
γ(h+1) is obtained by solving the following equation for γ : m X µ(h+1) − σ(2h+1)
−
m X
Rj E
j =1 ∗ θ (h+1) =
m X + µ(h+1) − σ(2h+1) Rj E
− yj − γ j =1 log(Z − γ ) Z > yj ; γ , θ ∗ (h+1) = 0, Z −γ j =1
m X log(yj − γ )
1
yj − γ
j =1
Z > yj ; γ , θ ∗ (h+1) Z −γ 1
(3.6)
(µ(h+1) , σ(h+1) ). The conditional expectations in the above expression (3.6) can be obtained as follows: σ2 (h+1) 1 ∗ −µ(h+1) 2 P0j(h+1) (γ ) E Z > yj ; γ , θ (h+1) = e Z −γ 2 σ (h+1) log(Z − γ ) −µ(h+1) Z > yj ; γ , θ ∗ 2 E σ(h+1) P0j(h+1) (γ ) + (µ(h+1) − σ(2h+1) )P0j(h+1) (γ ) , (h+1) = e Z −γ
where
where log yj − γ − µ(h+1)
Ψ0j(h+1) (γ ) =
Ψ0j (γ , θ ∗(h+1) )
=
P0j(h+1) (γ ) = P0j (γ , θ ∗(h+1) ) =
σ(h+1)
,
¯ (Ψ0j(h+1) (γ ) + σˆ (h+1) ) Φ . ¯ (Ψ0j(h+1) (γ )) Φ
The final MLEs of the parameters are local MLEs corresponding to the EM algorithm, and would be denoted here by LMLE2. In order to obtain LMLE2, we suggest a variation of the above method. As mentioned before, one can start with a γ(0) which is less than y1 and get µ(0) = µ(γ(0) ) and σ(0) = σ (γ(0) ) by using (3.4) and then one gets µ(1) and σ(1) by using (3.5). The suggestion is to treat (3.6) as the test equation. If the left hand side of that equation equals zero when (γ(0) , µ(1) , σ(1) ) substituted there, then no further iteration is needed. Otherwise, the suggested method involves choosing another γ(1) which is less than y1 and repeating the procedure looking for the sign change in (3.6). Finally, interpolating on γ values, one gets the final estimate (γˆ , µ, ˆ σˆ ) which when substituted in (3.6) should produce zero or close to zero. 4. Asymptotic variances and covariances of the EM algorithm estimates Asymptotic variances and covariances of the MLEs when the EM algorithm is used can be obtained by using the missing information principle of Louis (1982) and Tanner (1993). This principle is basically Observed information = Complete information − Missing information. Based on this principle, Louis (1982) developed the procedure of finding the observed information matrix when EM algorithm is used to find MLEs in an incomplete data situation. We adopt this principle in the situation of progressive
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
3585
Type II right censoring. We will denote the complete, observed and missing (expected) information by I (θ), IY (θ) and IZ |Y (θ), respectively. The complete information I (θ) is given by I (θ) = −E
∂ 2 E [log L(Y, Z, θ)|Y = y]
(4.1)
∂θ 2
in which E [log L(Y, Z, θ)|Y = y] is given by (3.1). In (4.1), the expectation is taken with respect to both Y and Z. The Fisher information matrix for a single observation which is censored at the time of the j-th failure is given by (j) IZ |Y (θ)
" = −E
∂ 2 log fZjk |Yj (zjk |yj ; θ)
# (4.2)
∂θ 2 (j)
in which fZjk |Yj (zjk |yj ; θ) is given by (3.3). In (4.2), the expectation is taken with respect to zjk so that IZ |Y (θ) is a function of yj and θ . Then, the missing (expected) information is simply IZ|Y (θ) =
m X
(j)
Rj IZ |Y (θ),
(4.3)
j =1
(j)
where IZ |Y (θ) is given by (4.2). The observed information IY (θ) is then obtained as follows: IY (θ) = I (θ) − IZ|Y (θ),
(4.4)
where I (θ) and IZ|Y (θ) are given by (4.1) and (4.3), respectively, and are derived below. Finally, upon inverting the observed information matrix IY (θ) in (4.4), one gets the asymptotic variances and covariances of the MLEs when EM algorithm is used. 4.1. Complete information matrix I (θ) The log-likelihood function log L∗ (θ) based on n uncensored observations yi , i = 1, . . . , n, is given by log L∗ (θ) =
n X
log f (yi ; θ),
(4.5)
i=1
where f (yi ; θ) is as given in (1.1). On differentiating the log-likelihood in (4.5) and equating to zero, one obtains the estimating equations given by (2.2) with m = n and Rj = 0; j = 1, 2, . . . , m. Negative of the second derivatives of the log-likelihood function log L∗ (θ) are obtained by appropriately differentiating the first derivatives in (2.2) with m = n and Rj = 0, j = 1, 2, . . . , m, and are given by (2.4) with m = n and Rj = 0, j = 1, 2, . . . , m. For the three-parameter lognormal distribution, it can be shown that
1
σ2
= e 2 −µ Yi − γ 1 2 E = e2(σ −µ) (Yi − γ )2 log(Yi − γ ) σ2 E = (µ − σ 2 )e 2 −µ Yi − γ log(Yi − γ ) 2 E = (µ − 2σ 2 )e2(σ −µ) . (Yi − γ )2
E
(4.6)
One gets the complete information matrix I (θ) in (4.1), by using (2.4) with m = n, Rj = 0, j = 1, 2, . . . , m, and (4.6), as follows:
1 + σ2
2(σ 2 −µ)
σ2 e 1 σ 2 −µ I (θ) = n e 2 σ2 2 2 σ − e 2 −µ σ
1
σ2
e
σ 2 −µ 2
1
2
− e σ
σ 2 −µ
2
σ2
0
0
σ2
2
.
(4.7)
Expressions for the asymptotic variances and covariances of the MLE of θ in the uncensored case are obtained by inverting
−1 2 the matrix I (θ) in (4.7). Denoting β = eµ , w = eσ and H = w(1 + σ 2 ) − (1 + 2σ 2 ) , the asymptotic variances and
3586
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
covariances are:
σ 2 β2 H n w σ2 V (µ) ˆ = (1 + H )
n 2 σ 2 V (σˆ ) = (1 + 2σ H ) 2n 2 σ β Cov(γˆ , µ) ˆ = − √ H n w σ3 β Cov(γˆ , σˆ ) = √ H n w 3 σ Cov(µ, ˆ σˆ ) = − H . V (γˆ ) =
n It is worthwhile to mention here that, for the uncensored case, Cohen (1951) has obtained the asymptotic variances and covariances of the MLEs of the parameters (γ , β, σ ), while those corresponding to the parameters (γ , µ, σ 2 ) were obtained by Hill (1963). 4.2. Missing information matrix IZ|Y (θ) The logarithm of the density function of an observation zjk = z censored at yj , the time of the j-th failure, is given by [see Eq. (3.2)]
log fz |yj (z |yj ; θ) = const . − log σ − log(z − γ ) − log 1 − Φ
log(yj − γ ) − µ
σ
−
1 2σ 2
[log(z − γ ) − µ]2 . (4.8)
Differentiating (4.8) with respect to γ , µ and σ , one gets
∂ log fz |yj ∂γ ∂ log fz |yj ∂µ ∂ log fz |yj ∂σ
= = =
1 log(z − γ ) − µ
σ 2 1
z−γ log(z − γ ) − µ
σ
σ
1
σ
"
+
1 z −γ
−
1
1
σ yj − γ
L0j ,
− L0j ,
log(z − γ ) − µ
2
σ
− 1 + ψ0j L0j
#
.
(4.9)
By using the properties of the left-truncated log-normal distribution, it can then be shown that E log(Z − γ ) − µ|Z > yj ; θ = σ L0j ,
E {log(Z − γ ) − µ}2 |Z > yj ; θ = σ 2 1 + Ψ0j L0j , E {log(Z − γ ) − µ}3 |Z > yj ; θ = σ 3 2 + Ψ0j2 , E {log(Z − γ ) − µ}4 |Z > yj ; θ = σ 4 3 1 + Ψ0j L0j + Ψ0j3 L0j , 1 σ2 E |Z > yj ; θ = e 2 −µ LL1j , Z −γ " # 2 1 2 E Z > yj ; θ = e2(σ −µ) LL2j , Z −γ 2 log(Z − γ ) Z > yj ; θ = e σ2 −µ σ L1j + (µ − σ 2 )LL1j , E Z −γ " # 2 log(Z − γ ) 2 E Z > yj ; θ = e2(σ −µ) σ 2 Ψ2j + 2σ (µ − 2σ 2 ) L2j + σ 2 + (µ − 2σ 2 )2 LL2j , Z −γ log(Z − γ ) Z > yj ; θ = e2(σ 2 −µ) σ L2j + (µ − 2σ 2 )LL2j , E 2 (Z − γ ) σ2 (log(Z − γ ))2 E Z > y ; θ = e 2 −µ σ 2 Ψ1j + 2σ (µ − σ 2 ) L2j + σ 2 + (µ − σ 2 )2 LL2j j Z −γ σ2 (log(Z − γ ))3 E Z > y ; θ = e 2 −µ σ 2 Ψ1j2 + 2σ 3 − 3σ (µ − σ 2 )2 L1j j Z −γ + 3(µ − σ 2 )σ 2 Ψ1j + 6σ (µ − σ 2 )2 L2j − 2(µ − σ 2 )3 LL1j + 3(µ − σ 2 ) σ 2 + (µ − σ 2 )2 LL2j ,
(4.10)
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
3587
in which Ψ0j , Ψ1j , Ψ2j , L0j , L1j , L2j , LL1j and LL2j are given by log yj − γ − µ
Ψ0j =
σ
log(yj − γ ) − (µ − σ 2 )
, σ 2 log(yj − γ ) − (µ − 2σ ) Ψ2j = , σ φ(Ψ1j ) φ(Ψ0j ) , L1j = , L0j = ¯ (Ψ0j ) ¯ (Ψ0j ) Φ Φ ¯ (Ψ1j ) ¯ (Ψ2j ) Φ Φ LL1j = , LL2j = . ¯ (Ψ0j ) ¯ (Ψ0j ) Φ Φ Ψ1j =
,
φ(Ψ2j ) L2j = , ¯ (Ψ0j ) Φ
(4.11)
(j)
The Fisher information matrix based on one observation zjk censored at yj , IZ |Y (θ), in (4.2), can be obtained using (4.9)–(4.11) as follows:
∂ log fZ |yj
2 Ψ e2(σ −µ) 2 + σ2j2
2
1
2(σ 2 −µ)
1
2e
σ 2 −µ 2
− L2j + e 1 + 2 LL2j − 2 L0j L1j σ 2 (yj − γ ) σ σ σ (yj − γ ) ∂ log fZ |yj 2 1 E = 2 1 + Ψ0j L0j − L20j ∂µ σ ∂ log fZ |yj 2 1 E = 2 2 + Ψ0j L0j 1 − Ψ0j L0j + Ψ0j2 ∂σ σ 2 ∂ log fZ |yj ∂ log fZ |yj 1 σ = 3 e 2 −µ σ (σ + 1)L1j − σ 2 (σ + 1)LL1j E ∂γ ∂µ σ σ L0j + 1 − σ 2 −σ L0j LL1j − L0j yj − γ " ∂ log fZ |yj ∂ log fZ |yj 1 2µ (σ − 1)L0j + µ = L0j E ∂γ ∂σ σ σ 2 (yj − γ ) ) ( 2µ2 + Ψ1j2 1 − Ψ0j L0j + 4µ 2µ2 σ2 + e 2 −µ − 3 + + − 3σ L1j 2 σ σ σ n µ µ o µ µ σ2 − − 1 − 4 − 1 σ + µ − + e 2 −µ 2Ψ1j L2j σ 2 σ 2 σ2 σ 2 2 2 µ − σ 2 − µσ µ µ µ σ σ −µ 2 2 2 2 2 + 2e 2 −µ − 1 − − σ + µ LL + 2e σ + (µ − σ ) LL 1j 2j σ2 σ σ2 σ4 ∂ log fZ |yj ∂ log fZ |yj 1 = 2 L0j + Ψ0j L0j (Ψ0j − L0j ) . E ∂µ ∂σ σ E
∂γ
=
L2 2 0j
(4.12)
Using IZ|Y (θ) from (4.12) and I (θ) from (4.7), the observed information can be obtained from (4.4). Finally, one can get the asymptotic variances and covariances of the MLEs, when the EM algorithm is used, by inverting this observed information matrix IY (θ). 5. Alternative estimators As an alternative estimator, one can use modified MLE (MMLE). MMLEs use first order statistic for estimating γ and are therefore easier to compute than the MLEs. For small samples when LMLEs fail to converge, these MMLEs (which always exist) produce reasonable estimates. For progressively censored data, one can have two versions of MMLEs. One version corresponds to the Newton–Raphson algorithm and the other version corresponds to the EM algorithm, and they will be ∂ log L denoted by MMLE1 and MMLE2, respectively. In both versions, the likelihood equation ∂γ = 0 is replaced by
γ + eµ+σ Φ
−1 [ 1 ] n+1
Pm
= y1
(5.1)
in which n = m + j=1 Rj is the total sample size. Eq. (5.1) is used as a test equation and it replaces the first equation in (2.2) for the Newton–Raphson algorithm version and Eq. (3.6) for the EM algorithm version. One starts with a first approximation γ1 and make the transformation log(yi − γ1 ). One then calculates conditional estimates µ1 = µ(γ1 ) and σ1 = σ (γ1 ) by using the second and third equations of (2.2) for the Newton–Raphson algorithm and (3.5) for the EM algorithm. The values (γ1 , µ1 , σ1 ) are substituted in Eq. (5.1). If the test equation is satisfied, then no further iteration is required. Otherwise, a
3588
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
Table 1 Progressively censored data for Example 1. i
1
2
3
4
5
6
7
8
9
10
11
12
Yi Ri
0.265 0
0.269 0
0.297 2
0.315 0
0.338 0
0.379 2
0.392 0
0.402 0
0.412 0
0.416 0
0.418 0
0.484 4
Table 2 Estimates and standard deviations of estimates for γ , µ and σ based on data in Table 1.
γ Estimator LMLE2
µ
γˆ
σγˆ
0.16
0.09
σ
µ ˆ −1.21
σµˆ
σˆ
σσˆ
0.19
0.29
0.05
Table 3 Progressively censored data for Example 2. i
1
2
3
4
5
6
7
Yi Ri
152.7 0
172.0 0
172.5 0
173.3 1
193.0 0
204.7 0
234.9 2
Table 4 Estimates and standard deviations of estimates for γ , µ and σ based on data in Table 3.
γ Estimator LMLE2
µ
σ
γˆ
σγˆ
µ ˆ
σµˆ
σˆ
σσˆ
110.30
10.29
3.35
0.89
0.81
0.12
second approximation γ2 is selected and the cycle of calculations as described above is repeated. The iterations are continued until two sufficiently close values γi and γi+1 are found such that the following is satisfied:
γi + eµi +σi Φ
−1 [ 1 ] n+1
< (>)y1 < (>)γi+1 + eµi+1 +σi+1 Φ
−1 [ 1 ] n+1
.
Thus, the final MMLE γ˜ of γ is obtained using which the final MMLEs µ ˜ , σ˜ of µ and σ are then obtained. 6. Illustrative examples To illustrate the computational methods presented in this article, we use three examples given in Cohen and Whitten (1988). We modified these original examples to consider progressively censored data. We then assume that each datum comes from a three-parameter lognormal distribution, and use them to carry out all the estimation procedures discussed in the preceding sections in order to see whether they produce similar results or not. Also, in Section 7, we compare the coverage probabilities of confidence intervals based on LMLE1 and LMLE2 using a simulation study. We used same γ0 value for obtaining LMLE1 and LMLE2. For the determination of MMLE, as discussed in Section 5, Φ −1 [ n+1 1 ] = −1.67 for Examples 1 and 3 since n = 20 in both these cases. For Example 2, Φ −1 [ n+1 1 ] = −1.34 since n = 10 here. Proceeding as explained in Section 5, we get two versions of MMLE, viz., MMLE1 and MMLE2 (for the Newton–Raphson algorithm and EM algorithm, respectively). Example 1. The maximum flood levels (in millions of cubic feet per second) for 20 four-year periods from 1890 to 1969 in the Susquehanna river at Harrisburg, Pennsylvania are given in Cohen and Whitten (1988) and were also used by Dumonceaux and Antle (1973). We modified these data to make it progressively censored with m = 12 stages, and these progressively censored data are presented in Table 1. The LMLE2 and their standard deviations are presented in Table 2. For the data in Table 1, the LMLE1 are given by γˆ = 0.1822, µ ˆ = −1.291 and σˆ = 0.362, the MMLE1 are given by γ˜ = 0.1581, µ ˜ = −1.103 and σ˜ = 0.279, and the MMLE2 are given by γ˜ = 0.1563, µ ˜ = −1.1725 and σ˜ = 0.273. Example 2. This example was used by McCool (1974) and is also given in Cohen and Whitten (1988). The data are fatigue lives (in hours) of 10 bearings of a certain type. We modified these data to make it progressively censored with m = 7 stages, and these progressively censored data are presented in Table 3. LMLE1 did not converge for this sample. The LMLE2 and their standard deviations are presented in Table 4. The MMLE1 are given by γ˜ = 112.7815, µ ˜ = 3.9231 and σ˜ = 0.8932, while MMLE2 are given by γ˜ = 109.7232, µ ˜ = 3.2763 and σ˜ = 0.8036.
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
3589
Table 5 Progressively censored data for Example 3. i
1
2
3
4
5
6
7
8
9
10
11
12
Yi Ri
127.2 0
128.7 0
131.4 0
133.0 2
133.1 0
135.9 0
144.3 2
145.8 0
148.3 0
153.1 0
157.3 0
164.3 4
Table 6 Estimates and standard deviations of estimates for γ , µ and σ based on data in Table 5.
γ Estimator LMLE2
µ
σ
γˆ
σγˆ
µ ˆ
σµˆ
σˆ
σσˆ
108.25
11.27
3.82
0.86
0.39
0.01
Table 7 Progressively censored data for Example 3 with PCS-1. i
1
2
3
4
5
6
7
8
9
10
11
12
Yi Ri
127.2 0
128.7 0
131.4 0
133.0 0
133.1 0
135.9 0
137.3 0
144.3 0
145.8 0
148.3 2
153.1 2
157.2 4
Table 8 Progressively censored data for Example 3 with PCS-2. i
1
2
3
4
5
6
7
8
9
10
11
12
Yi Ri
127.2 4
128.7 2
133.0 2
135.9 0
144.3 0
148.3 0
153.1 0
157.2 0
166.5 0
174.8 0
184.1 0
201.4 0
Table 9 Estimates and standard deviations of estimates for γ , µ and σ for data in Table 7.
γ Estimator LMLE2
µ
σ
γˆ
σγˆ
µ ˆ
σµˆ
σˆ
σσˆ
102.11
11.79
4.01
1.38
0.36
0.03
Table 10 Estimates and standard deviations of estimates for γ , µ and σ for data in Table 8.
γ Estimator LMLE2
µ
σ
γˆ
σγˆ
µ ˆ
σµˆ
σˆ
σσˆ
108.37
11.99
4.12
0.6
0.39
0.02
Example 3. In this example, the complete data are given by Cohen and Whitten (1988) and were also used by Cohen (1951). The complete data consists of 20 observations from a three-parameter lognormal distribution with γ = 100, µ = 3.912023 and σ = 0.4. We modified these data to make it progressively censored with m = 12 stages, and these progressively censored data are presented in Table 5. Although Yi values are reported up to one decimal place in Tables 5, 7 and 8, we used all three decimal places for these Yi values in Cohen (1951) for our computations. The LMLE2 and their standard deviations are presented in Table 6. For the data in Table 5, the LMLE1 are given by γˆ = 110.3112, µ ˆ = 3.573 and σˆ = 0.422, the MMLE1 are given by γ˜ = 109.1563, µ ˜ = 4.419 and σ˜ = 0.295, and the MMLE2 are given by γ˜ = 107.8036, µ ˜ = 3.927 and σ˜ = 0.378. In order to examine the effect of delayed censoring, we considered two schemes of progressive censoring. In one scheme, the censoring was delayed than the other. Keeping the number of stages m = 12 the same as in Table 4, for the delayed censoring scheme, we considered PCS-1: R1 = 0, R2 = 0, R3 = 0, R4 = 0, R5 = 0, R6 = 0, R7 = 0, R8 = 0, R9 = 0, R10 = 2, R11 = 2, R12 = 4. For the other scheme, we took PCS-2: R1 = 4, R2 = 2, R3 = 2, R4 = 0, R5 = 0, R6 = 0, R7 = 0, R8 = 0, R9 = 0, R10 = 0, R11 = 0, R12 = 0. The progressively censored data obtained under these two schemes are presented in Tables 7 and 8, respectively. The LMLE2 estimates and their standard deviations for the data in Tables 7 and 8 are presented in Tables 9 and 10, respectively. For the data in Table 7, the LMLE1 are given by γˆ = 104.4622, µ ˆ = 3.992 and σˆ = 0.373, the MMLE1 are given by γ˜ = 103.5689, µ ˜ = 4.151 and σ˜ = 0.298, and the MMLE2 are given by γ˜ = 103.8039, µ ˜ = 4.107 and σ˜ = 0.315. For the data in Table 8, the LMLE1 are given by γˆ = 109.3784, µ ˆ = 4.041 and σˆ = 0.413, the MMLE1 are given by γ˜ = 107.2501, µ ˜ = 4.123 and σ˜ = 0.292, and the MMLE2 are given by γ˜ = 107.9231, µ ˜ = 4.137 and σ˜ = 0.368. It is
3590
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
Table 11 95% confidence intervals for Example 1. Estimators
γ
µ
σ
LMLE1 LMLE2 MMLE1 MMLE2
(−.255, .647) (−.071, .415) (−.250, .574) (−.165, .423)
(−2.165, −.477) (−1.644, −.816) (−1.855, −.327) (−1.651, −.593)
(.085, .673) (.137, .423) (.047, .595) (.085, .437)
Table 12 95% confidence intervals for Example 2. Estimators
γ
µ
σ
LMLE2 MMLE1 MMLE2
(66.682, 152.942) (69.882, 157.815) (65.721, 139.061)
(1.689, 5.535) (1.284, 6.968) (1.232, 4.956)
(.323, 1.303) (.423, 1.482) (.403, 1.030)
Table 13 Coverage probabilities of 95% confidence intervals for γ based on Monte Carlo simulations. n
m
R
LMLE1
LMLE2
20 20 20 20 20 20 20 20 40 40 40 40 40 40 40 40
12 12 14 14 16 16 18 18 24 24 28 28 32 32 36 36
(7,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,7) (5,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,5) (3,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,3) (1,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,1) (14,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,14) (10,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,10) (6,2,0, . . . ,0,0,0) 0,0,0, . . . ,0,2,6) (2,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,2)
90.71 90.92 91.02 91.12 93.34 93.78 93.92 94.35 92.10 92.94 93.21 94.28 94.37 93.48 93.92 94.73
92.25 92.81 93.01 93.82 94.71 94.29 94.35 95.72 93.01 94.15 93.91 94.23 94.41 94.87 94.91 94.28
clear from Tables 6, 9 and 10 that the estimates of γ are affected not only by censoring but also by the pattern of censoring, while the estimates of µ and σ remain more or less the same in all the situations. 7. Simulation study In this section, we report results from two simulation studies that we carried out. In the first simulation study, we constructed 95% confidence intervals based on 10,000 randomly generated progressively censored samples corresponding to different estimates of γ , µ and σ for each of the estimation methods discussed in Section 6. We used the same progressive censoring schemes as given in those examples. We provide the confidence intervals from this simulation study in Tables 11 and 12 for Examples 1 and 2, respectively, for different estimators. In the second simulation study, we compared the performance of LMLE1 and LMLE2 in terms of coverage probabilities of 95% confidence intervals for the parameters γ , µ and σ for different sample sizes and different degrees of censoring. 10,000 samples were simulated from the lognormal distribution with γ = 100, µ = 3.912023 and σ = 0.4 with sample size n = 20, 40. For sample size n = 20, we considered m = 12, 14, 16, 18 stages of censoring, and for sample size n = 40 we considered m = 24, 28, 32, 36 stages. In each case, we considered two censoring schemes with one reflecting comparatively delayed censoring than the other. The coverage probabilities of 95% confidence intervals from this simulation study for the parameters γ , µ and σ are presented in Tables 13–15, respectively. It is observed from these tables that the coverage probabilities are better when the proportion of uncensored data is larger. The coverage probabilities seem to be almost the same for the two different censoring schemes. The confidence intervals for the parameter γ seems to be most sensitive to the censoring pattern, while the confidence intervals for the parameters µ and σ seem to be stable and quite satisfactory and close to the nominal level of 95%. In a similar tone, it is known that different censoring schemes do not change the estimates for the two-parameter exponential distribution (see for example, Balakrishnan and Sandhu (1996)). The sensitivity of estimation of γ may be due to the fact that γ is the threshold parameter of the lognormal distribution and the pattern of censoring applied to the sample might have affected the nature of the sample observations as far as the estimation of the threshold parameter is involved. Also, LMLE1 and LMLE2 seem to have nearly the same coverage probabilities.
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
3591
Table 14 Coverage probabilities of 95% confidence intervals for µ based on Monte Carlo simulations. n
m
R
LMLE1
LMLE2
20 20 20 20 20 20 20 20 40 40 40 40 40 40 40 40
12 12 14 14 16 16 18 18 24 24 28 28 32 32 36 36
(7,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,7) (5,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,5) (3,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,3) (1,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,1) (14,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,14) (10,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,10) (6,2,0, . . . ,0,0,0) 0,0,0, . . . ,0,2,6) (2,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,2)
93.12 93.28 93.27 94.02 94.96 95.02 94.81 95.38 93.24 93.38 94.11 94.43 95.02 95.15 95.39 95.58
93.94 93.82 93.73 93.61 95.11 95.25 95.42 95.46 94.12 94.22 94.09 94.42 95.36 95.37 95.38 96.01
Table 15 Coverage probabilities of 95% confidence intervals for σ based on Monte Carlo simulations. n
m
R
PC for LMLE1
PC for LMLE2
20 20 20 20 20 20 20 20 40 40 40 40 40 40 40 40
12 12 14 14 16 16 18 18 24 24 28 28 32 32 36 36
(7,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,7) (5,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,5) (3,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,3) (1,1,0, . . . ,0,0,0) (0,0,0, . . . ,0,1,1) (14,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,14) (10,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,10) (6,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,6) (2,2,0, . . . ,0,0,0) (0,0,0, . . . ,0,2,2)
92.94 93.17 93.25 93.74 94.82 94.29 95.03 95.41 93.15 93.21 93.36 93.81 95.01 95.21 95.11 95.14
93.66 93.79 93.73 93.85 94.21 94.62 95.43 95.44 93.62 93.79 93.36 93.76 95.21 95.01 95.35 95.91
8. Concluding remarks In this article, we have discussed the EM algorithm for the maximum likelihood estimation based on progressively Type-II censored samples from a three-parameter lognormal distribution. We have also considered the traditional Newton–Raphson method for this purpose. Additionally, we have discussed two versions of modified maximum likelihood estimators. The EM algorithm and the Newton–Raphson method produced similar results in all three examples considered, as well as in terms of coverage probabilities determined from a Monte Carlo simulation study. The coverage probabilities for both methods are better and closer to the nominal level of 95% when the proportion of uncensored data is larger. Acknowledgments The authors are thankful to the referees for their valuable comments which led to a considerable improvement in the presentation of this article. References Balakrishnan, N., 2007. Progressive censoring methodology: An appraisal (with discussions). Test 20, 211–296. Balakrishnan, N., Aggarwala, R., 2000. Progressive Censoring: Theory, Methods, and Applications. Birkhäuser, Boston. Balakrishnan, N., Cohen, A.C., 1991. Order Statistics and Inference: Estimation Methods. Academic Press, San Diego. Balakrishnan, N., Kannan, N., Lin, C.T., Ng, H.K.T., 2003. Point and interval estimation for Gaussian distribution based on progressively Type-II censored samples. IEEE Transactions on Reliability 52, 90–95. Balakrishnan, N., Sandhu, R.A., 1996. Best linear unbiased and maximum likelihood estimation for exponential distributions under general progressive type-II censored samples. Sankhya, Series B 58, 1–9. Calitz, F., 1973. Maximum likelihood estimation of the parameters of the three-parameter lognormal distribution — A reconsideration. Australian Journal of Statistics 15, 185–199. Cohen, A.C., 1951. Estimating parameters of logarithmic-normal distributions by maximum likelihood. Journal of the American Statistical Association 46, 206–212.
3592
P. Basak et al. / Computational Statistics and Data Analysis 53 (2009) 3580–3592
Cohen, A.C., 1988. Estimation in the lognormal distributions. In: Crow, E.L., Shimizu, K. (Eds.), Lognormal Distributions: Theory and Applications. Marcel Dekker, New York. Cohen, A.C., Whitten, B.J., 1988. Parameter Estimation in Reliability and Life Span Models. Marcel Dekker, New York. Crow, E.L., Shimizu, K. (Eds.), 1988. Lognormal Distributions: Theory and Applications. Marcel Dekker, New York. Dempster, A.P., Laird, N.M., Rubin, D.B., 1977. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B 39, 1–38. Dumonceaux, R., Antle, C.E., 1973. Discrimination between the log-normal and the Weibull distributions. Technometrics 15, 923–926. Harter, H.L., Moore, A.H., 1966. Local-maximum-likelihood estimation of the parameters of the three-parameter lognormal populations from complete and censored samples. Journal of the American Statistical Association 61, 842–851. Corrections: 61 (1966) 1247; 62 (1967) 1519–1520; 63 (1968) 1549. Hill, B.M., 1963. The three-parameter lognormal distribution and Bayesian analysis of a point-source epidemic. Journal of the American Statistical Association 58, 72–84. Johnson, N.L., Kotz, S., Balakrishnan, N., 1994. Continuous Univariate Distributions – vol. 1, second ed. John Wiley & Sons, New York. Louis, T.A., 1982. Finding the observed information matrix when using the EM algorithm. Journal of the Royal Statistical Society, Series B 44, 226–233. McCool, J.I., 1974. Inferential techniques for Weibull populations. Aerospace Research Laboratories Report, ARL TR 74-0180, Wright-Patterson AFB, Ohio. McLachlan, G.J., Krishnan, T., 1997. The EM Algorithm and Extensions. Marcel Dekker, New York. Mi, J., Balakrishnan, N., 2003. Existence and uniqueness of the MLEs for normal distribution based on general progressively Type-II censored samples. Statistics & Probability Letters 64, 407–414. Munro, A.H., Wixley, R.A.J., 1970. Estimators based on order statistics of small samples from a three-parameter lognormal distribution. Journal of the American Statistical Association 65, 212–225. Nelson, W., 1982. Applied Life Data Analysis. John Wiley & Sons, New York. 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. Rukhin, A.L., 1984. Improved estimation in lognormal models. Technical Report No. 84-38, Department of Statistics, Purdue University, West Lafayette, Indiana. Tanner, M.A., 1993. Tools for Statistical Inference: Observed Data and Data Augmentation Methods, second ed. Springer-Verlag, New York.