Confidence Interval Procedures for System ... - Semantic Scholar
Confidence Interval Procedures for System Reliability and Applications to Competing-Risk Models Yili Hong
William Q. Meeker
Department of Statistics
Department of Statistics
Virginia Tech
Iowa State University
Blacksburg, VA 24061
Ames, IA 50011 Abstract
System reliability depends on the reliability of the system’s components and the structure of the system. For example, in a competing-risk model, the system fails when the weakest component fails. The reliability function and the quantile function of a complicated system are two important metrics for characterizing the system’s reliability. When there are data available at the component level, the system reliability can be estimated by using the component level information. Confidence intervals (CI) are needed to quantify the statistical uncertainty in the estimation. Obtaining system reliability CI procedures with good properties is not straightforward, especially when the system structure is complicated. In this paper, we develop a general procedure for constructing a CI for the system failure-time quantile function by using the implicit delta method. We also develop general procedures for constructing a CI for the cdf of the system. We show that the recommended procedures are asymptotically valid and have good statistical properties. We conduct simulations to study the finitesample coverage properties of the proposed procedures and compare them with existing procedures. We apply the proposed procedures to three applications; two applications in competing-risk models and an application with a k-out-of-s system. The paper concludes with some discussion and an outline of areas for future research. Key Words: Asymptotic approximation, failure-time data, k-out-of-s system, maximum likelihood, normal approximation, series system.
1
1
Introduction
1.1
Background
For complicated systems, failures occurring at the component level can lead to system failures. For example, in a competing-risk setting (i.e., series system), a system fails when the first component fails. Meeker and Escobar (1998, page 383) described an application of Device G, where the system fails due to one of two competing failure modes. Hong and Meeker (2010) described another application of Product D, where the system fails due to one of four competing failure modes. Thus the reliability of a complicated system depends on the reliability of its components. The relationship between system reliability and component reliability is often of interest. In the competing risk setting, system reliability is the product of reliability of individual components. The relationship can be complicated in some systems such as the k-out-of-s system, where a system of s components operates when k or more of the system’s components operate. Nevertheless, with the relationship between system reliability and component reliability and information on component reliability, one can obtain information on system reliability. The reliability function and the quantile function of a system are two important metrics for characterizing system reliability. The reliability function is defined by R(t) = 1 − F (t) where F (t) is the cumulative distribution function (cdf) of the system failure-time random variable T . The cdf F (t) is interpreted as the probability that a unit will fail by time t, or the proportion of units in the population that will fail by time t. The quantile function tp is the inverse of the cdf, and corresponds to the time at which a specified proportion p of the population of systems will fail. When there are data available at the component level, these metrics can be estimated by using the component level data. Confidence intervals (CIs) are needed to quantify the statistical uncertainty in the estimation of F (t) and tp . In practice, it is common to plot on one graph an estimate and a set of pointwise CIs for F (t) over a range of t values. These sets of pointwise CIs are referred to as a pointwise confidence band (CB). Obtaining a CI procedure for the system reliability with good statistical properties is not straightforward, especially when the system structure is complicated. The following properties of the CI procedure should be considered. • The coverage probability of the CI procedure should be approximately equal to the nominal coverage probability for small or moderate sample size and converge to the nominal coverage probability when the sample size goes to infinity. • The endpoints of the CI should fall inside the parameter space, because the cdf is bounded by 0 and 1. Some normal approximation CI procedures may produce CIs 2
with endpoints that fall outside [0, 1]. • A transformation function can help restrict the CI endpoints to be within [0, 1] but may adversely affect the monotonicity of the boundary of the pointwise CBs. Hong, Meeker, and Escobar (2008a) considered the abnormal behavior of pointwise CBs caused by the use of an inappropriate transformation function, resulting in a boundary of the CB is not monotone increasing. That abnormal behavior was called the “bend-back” behavior. The bend-back behavior is not desirable in applications and needs to be avoided where possible. • In many applications, especially in software implementations, computational efficiency of the CI procedure is also important. Although CI procedures based on likelihood or simulation may have important advantages, normal approximation CI procedures are widely used due to their computational efficiency. The normal approximation CI procedures for F (t) of systems with simple structures are widely implemented in popular software packages. These software package includes JMP (2011), Minitab (2011), WEIBULL++ 7.5 by ReliaSoft (2011), PROC RELIABILITY in SAS (2011), and SPLIDA/RSPLIDA described in Meeker and Escobar (2008). In this paper, we develop a general procedure for constructing CIs for the quantile function by using the implicit delta method. We also develop a general procedure constructing CIs for the system cdf. We show that the proposed procedures are asymptotically valid and have good statistical properties. The proposed procedures can also be easily implemented in most software packages.
1.2
Related Literature
Rausand and Høyland (2004) provides a comprehensive discussion of system reliability. In the area of system reliability CIs, some earlier work includes Myhre and Saunders (1968) and Easterling (1972), which are based on binary data and a simple series systems. Coit (1997) used a log transformation to construct CI for the reliability function using binary or lifetime data collected at component level. Meeker and Escobar (1998) used the logit transformation to construct CI for the cdf of the system. Ramirez-Marquez and Jiang (2006) considered a CI procedure by using a recursive approach. Spall (2009) constructed a CI for system reliability by using the information from both full system and component level tests. At the component level, a failure-time distribution is needed for the component reliability. Meeker and Escobar (1998) is a useful resource for statistical methods for reliability data. Hamada et al. (2008) provide statistical methods for analyzing reliability data from a 3
Bayesian perspective. For systems failing due to competing risks, David and Moeschberger (1978), and Crowder (2001) are useful books on the subject of competing-risk models.
1.3
Overview
The rest of the paper is organized as follows. Section 2 introduces the model for the system and maximum likelihood (ML) estimation for unknown parameters. Section 3 presents the proposed CI procedures for the cdf and the quantile function, based on normal approximations. Section 4 presents likelihood based CI procedures for the cdf and the quantile function. Section 5 describes the results of a simulation study to evaluate the performance of several CI procedures, using various criteria. Section 6 applies the proposed CI procedures to three different applications. Section 7 contains some conclusions and areas for future research.
2
Model and ML Estimation
2.1
System Models
We consider a coherent system with s components. A coherent systems is a system that meets two requirements: 1) the improvement of any component will improve the system reliability, and 2) every component in the system makes some nonzero contribution to the system’s reliability. Let T be the failure time of the system. The system cdf is defined by F (t; θ) = Pr(T ≤ t) where θ is the vector for the parameters in the system cdf. The cdf for the failure time of component i is denoted by Fi (t; θ i ), i = 1, . . . , s. Here θ i is a vector for the parameter in the cdf of component i. The relationship between a system’s cdf and components’ cdfs is characterized by a system structure function g. The cdf for failure-time distribution of the system is expressed as F (t; θ) = g[F1 (t; θ 1 ), . . . , Fs (t; θ s )].
(1)
Here θ = (θ ′1 , . . . , θ ′s )′ and g = g(x1 , · · · , xs ) is a function mapping from [ 0, 1 ]s (an s dimensional unit hypercube) to [ 0, 1 ]. We assume the function g is a strictly increasing and twice differentiable function with respect to each of its arguments, which can be met for most common system structures. The form of the function g is determined by system’s structure. Some common examples include the following (see Rausand and Høyland 2004 for a comprehensive discussion).
4
• For a series system with s components, g(x1 , . . . , xs ) = 1 −
s Y
(1 − xi ).
(2)
i=1
• For a parallel system with s components, g(x1 , . . . , xs ) =
Qs
i=1
xi .
• For a 2 × 2 series-parallel system with component-level redundancy, g(x1 , x2 , x3 , x4 ) = 1 − (1 − x1 x2 )(1 − x3 x4 ). • For a k-out-of-s system with s components, the system fails when there are s − k + 1 or more components fail. The g function is ( ) s s 1 X exp[−iωl(s − k + 1)] − exp(−iωl) Y g(x1 , . . . , xs ) = [1 − xi + xi exp(iωl)] . s + 1 l=0 1 − exp(−iωl) i=1 (3)
The details for the derivation of equation (3) can be found in Hong (2011). Depending on the system structure, the form of g may be more complicated than these examples. The CI procedures proposed in this paper works for a general form of g. Note that we use θ = (θ ′1 , . . . , θ ′s )′ to denote the unknown parameters in the system. For some system structures (e.g., when a system has several copies of the same component in the system structure), it is possible that the components share the same parameters. Thus there may be some redundancy in parameter θ. The extension of our proposed procedure to handle this situation is straightforward. For notational simplicity, we assume that each component has its own set of parameters. In addition to the cdf F (t; θ), the quantile function of the system failure-time distribution is another quantity that is often of interest. The p quantile is defined by the inverse of the cdf tp = F −1 (p). Numerical methods may be needed to solve tp from F (tp ; θ) = p. The log of the failure time T is denoted by Y = log(T ). The quantile function of Y is denoted by yp = log(tp ). Note that both tp and yp depend on the unknown parameter vector θ. That is tp = tp (θ) and yp = yp (θ).
2.2
Component Failure-time Model
Here we use the log-location-scale family of distributions to model the failure-time distribution of the components. The Weibull and lognormal distributions, which are the most commonly-used distributions for product/system failure times, are members of this family.
5
The cdf and probability density function (pdf) of a log-location-scale distribution can be expressed as
log(t) − µ F (t; µ, σ) = Φ σ
1 log(t) − µ f (t; µ, σ) = φ , σt σ
and
where Φ is the standard cdf for the location-scale family of distributions (location 0 and scale 1), µ is the location parameter, and σ is the scale parameter. For the lognormal distribution, replace Φ and φ above with Φnor and φnor , the standard normal cdf and pdf, respectively. For the Weibull distribution, replace Φ and φ above with Φsev (z) = 1 − exp[− exp(z)] and φsev (z) = exp[z − exp(z)], respectively. The cdf of the component failure time is denoted by Fi (t; θ i ) = Φi (zi ), i = 1, 2, · · · , s (as the form of the distribution need not be the same for all components) where θ i = (µi , σi )′ and zi = [log(t) − µi ]/σi . Here Φi (·) is a standard location-scale cdf. Although the loglocation-scale family of distribution is used to model the failure time of the components, the CI procedures developed in this paper are general. We provide specific results for the log-location-scale family of distributions because this family of distributions is widely used.
2.3
Data and Maximum Likelihood Estimation
Suppose that data are collected at the component level from either a life test and/or the field. The data for component i is denoted by Di . These data are possibly censored. A loglocation-scale distribution is used to describe the failure time of component i. The maximum likelihood (ML) estimators of parameters of component i based on data Di , denoted by bi , i = 1, . . . , s. From large-sample theory, θ bi ∼ θ ˙ N(θ i , Σi ). Here Σi is the variance-covariance b To estimate Σi , we use the observed information matrix which can be written matrix of θ. as
bi2 bi = σ Σ ni
b11i λ b12i λ
b12i λ b22i λ
!
, i = 1, 2, · · · , s,
b11i , b where ni is the number of observations in the data Di , and λ λ12i and b λ22i are elements in the scaled variance-covariance matrix.
b = (θ b′ , . . . , θ b′ )′ . We assume The ML estimator for the system parameters is obtained by θ 1 s bi based on Di are independent from each other. Thus θ b is asymptotically the estimators θ multivariate normal (MVN) distributed. That is,
b∼ θ ˙ MVN [θ, Diag (Σ1 , Σ2 , · · · , Σs )] .
By the invariance property of ML estimators, the ML estimator of the distribution probability at te is b = g[F1 (te ; θ b1 ), . . . , Fs (te ; θ bs )]. Fb(te ) = F (te ; θ) 6
b Similarly, the ML estimator of tp is b tp = exp(b yp ) where ybp = yp (θ).
3
Normal Approximation Confidence Interval Procedures
3.1
Normal Approximation Confidence Interval Procedure for tp
An estimate of the variance of ybp based on the delta method and observed information is ∂b yp yp d (b b b b ∂b Var yp ) = . ′ Diag(Σ1 , Σ2 , · · · , Σs ) b b ∂θ ∂θ
(4)
In many applications, there is no closed-form expression for yp = yp (θ). Thus we use the b Note that implicit derivative technique to find ∂b yp /∂ θ. p = g {F1 [exp(b yp ); θ 1 ], . . . , Fs [exp(b yp ); θ s ]} .
(5)
Taking derivatives on both sides of (5) with respect to µ bi , gives X ∂b yp ∂g ∂b yp ∂g 0= fi (b yp ) −1 + fj (b yp ) , ∂Fi ∂b µi ∂Fj ∂b µi j6=i yp ) = φi (b zip )/b σi . Here φi is the standard location-scale pdf of component i and where fi (b zbip = (b yp − µ bi )/b σi , i = 1, 2, · · · , s. Hence, ∂b yp ∂b µi
∂g fi (b yp ) ∂Fi = , Ps ∂g fi (b yp ) i=1 ∂Fi
i = 1, · · · , s,
P provided that si=1 (∂g/∂Fi )fi (b yp ) 6= 0. Because g is a strictly increasing function, Ps yp ) 6= 0. Similarly, i=1 (∂g/∂Fi )fi (b ∂b yp ∂b σi From (4),
d (b Var yp ) =
∂g fi (b yp )b zip ∂Fi = , Ps ∂g fi (b yp ) i=1 ∂Fi
Ps
i=1
2 2 ∂g σ bi fi (b yp ) ∂Fi ni Ps i=1
7
i = 1, 2, · · · , s.
2 b b b λ11i + 2λ12i zbip + λ22i zbip . 2 ∂g fi (b yp ) ∂Fi
Using the large-sample result that yb − yp qp ∼ ˙ N(0, 1), d Var (b yp )
we obtain the following CI procedure for yp .
Result 1 A large-sample approximate 100(1 − α)% CI for yp is s 2 Ps ∂g 2 b b12i zbip + λ b22i zb2 f (b y ) σ b γ λ + 2 λ i p i 11i i ip i=1 h i ∂Fi y p , yep = ybp ∓ , Ps ∂g e f (b y ) i p i=1 ∂Fi
(6)
2 where γi = z1−α/2 /ni , i = 1, 2, · · · , s. A large-sample approximate CI for tp is i h tp , e tp = exp(y p ), exp(e yp ) . e e Such approximate CIs based on Wald statistics may have a positive probability of bend-
back behavior (i.e., the boundary of the pointwise CB is not monotonically increasing, see Hong, Meeker, and Escobar 2008a for more details). We will study the probability of bending back for this CI procedure in Section 5.
3.2
Normal Approximation Confidence Interval Procedure for F (te )
A CI for the cdf at a particular time point t = te is often needed. Let p = F (te ; θ) and b Using the large-sample result that (b pb = F (te ; θ). p − p)/se b pb ∼ ˙ N(0, 1), a direct approach to construct a CI for p is h
i p, pe = pb − z1−α/2 se b pb, pb + z1−α/2 se b pb , (7) e where se b pb is the standard error which is obtained by using the delta method. In particular, v u s 2 2 uX ∂g σ bi b 2 t b b fi (ye ) λ11i + 2λ12i zbie + λ22i zbie . se b pb = ∂Fi ni i=1
The CI procedure in (7) is referred to as the Fb procedure. The Fb procedure may produce a CI with endpoints falling outside of [0, 1], which is not desirable. A logit or some other similar
transformations can be applied to restrict the CI endpoints to fall inside [0, 1] (e.g., page
190 of Meeker and Escobar 1998). Such a procedure, however, can have a high probability of the bend-back, as shown in Hong, Meeker, and Escobar (2008a) for a system with a single component. 8
The following alternative procedure, that we call ye procedure, will have CI endpoints that always fall inside the interval [0, 1] and has a much smaller probability of bending back. The construction of the CI, which is illustrated in Figure 1, is based, indirectly, on the endpoints obtained by using the Fb procedure. As shown in Figure 1, the endpoints of the CI obtained by the Fb procedure, reflected by the tangent line of Fb[exp(ye )] at ye = log(te ), are denoted by y e and yee , respectively. In particular, e s 2 h i Ps ∂g 2 2 b b b f (y ) σ b γ λ + 2 λ z b + λ z b i e 11i 12i ie 22i ie i i i=1 ∂Fi z1−α/2 se b pb y e = ye − = ye − , Ps ∂g ∂ Fb[exp(ye )] e fi (ye ) i=1 ∂Fi ∂ye s 2 h i Ps ∂g 2 b11i + 2λ b12i zbie + λ b22i zb2 f (y ) σ b γ λ i e i i ie i=1 ∂Fi z1−α/2 se b pb yee = ye + , = ye + Ps ∂g ∂ Fb[exp(ye )] fi (ye ) i=1 ∂Fi ∂ye
where fi (ye ) = φi (b zie )/b σi and zbie = (ye − µ bi )/b σi . The lower and upper endpoints of the b at exp(y e ) and exp(e proposed CI are obtained by evaluating F (t; θ) ye ), respectively. e Result 2 A large-sample approximate 100(1 − α)% CI for p = F (te ; θ) is h i h i p, pe = Fb{exp(y e )}, Fb{exp(e ye )} . e e
(8)
The proof of this result is given in Appendix A. The CI procedure in Result 2 may have a positive probability of bend-back behavior. We will study the probability of bending back for this CI procedure in Section 5. Hong, Meeker, and Escobar (2008a) studied a CI procedure for the cdf of the failuretime distribution of a system with a single component, which was called the zb procedure.
In particular, the zb procedure is based on the result that (b z1e − z1e )/se b zb1e ∼ ˙ N(0, 1), where
z1e = (ye − µ b)/b σ. The zb procedure has a small probability of bending back.
Here we show that the ye procedure in (8) is a generalization of the zb procedure. When
g(x) = x and s = 1 (i.e., there is only one component), the ye procedure (8) is identical to the zb procedure. Consider that ye e
r h i 1 2 2 2 b b b = ye − [f1 (ye )] σ b1 γ1 (λ111 + 2λ121 zb1e + λ221 zb1e ) f1 (ye ) r b111 + 2λ b121 zb1e + λ b221 zb2 = ye − σ b12 γ1 λ 1e 9
F (y)
Fb procedure
p=1
pb
p e
pe
z1−α/2 se b pb
t z1−α/2 se b pb gen n ta
e lin
pb
(ye , pb )
z1−α/2 se b pb
estimated cdf
p e
z1−α/2 se b pb
∂F ∂y
ye procedure
pe
∂F ∂y
p=0
ye e
ye
yee
y
Figure 1: Illustration of the proposed ye procedure to obtain a CI for F (te ). and q b111 + 2λ b121 zb1e + λ b221 zb2 ) − µ ye − σ b12 γ1 (λ b 1 1e p = Fb[exp(y e )] = Φ1 σ b 1 e e q b111 + 2λ b121 zb1e + λ b221 zb2 ) = Φ1 zb1e − γ1 (λ 1e
= Φ1 (z 1e ). e
Note that Φ1 (z 1e ) is the lower endpoint obtained by the zb procedure. Hence the lower e endpoints are identical. In a similar manner, one can show that the upper endpoints are identical. Thus the zb procedure is a special case of the ye procedure. 10
4
Likelihood Based Confidence Interval Procedures
In this section, we present likelihood based CI procedures for tp and F (te ). In general, there are no closed-form expressions for likelihood based CI procedures. Thus numerical methods are needed to compute the CIs and these can be computationally intensive. For this reason, in some applications, such as producing dynamic graphics (like those used in the profiler on the JMP reliability platform), the likelihood based CI procedures may not be suitable with currently available computer hardware.
4.1
Likelihood Based Confidence Interval Procedure for tp
Let L(θ) be the log-likelihood function of the data introduced in Section 2.3. Standard large-sample theory for an ML estimator, under standard regularity conditions (e.g., page 241 of Pawitan 2001, and Chapter 6 of Lehmann and Casella 1998), provides the result that ) ( b − e 2 L(θ) max L(θ) ∼ ˙ χ21 (9) e F (tp ; θ)=p} e {θ:
for fixed 0 < p < 1 and tp = tp (θ). The result in (9) provides the basis for a likelihood ratio test for tp . A 100(1 − α)% likelihood based CI for the tp , obtained by inverting the likelihood ratio test, is
where
tp , e
tp = min t : t satisfying e e tp = max t : t satisfying
e tp
(10)
e ≥ L(θ) b − 1 χ2 L(θ) e F (t; θ)=p} e 2 1;1−α {θ: e ≥ L(θ) b − 1 χ2 max L(θ) 1;1−α . e F (t; θ)=p} e 2 {θ: max
Here χ21;1−α is the 1 − α quantile of the χ21 distribution.
4.2
Likelihood Based Confidence Interval Procedure for F (te )
Standard large-sample theory for ML estimator also provides the result that ( ) b − e 2 L(θ) max L(θ) ∼ ˙ χ21 e (te ; θ)=p} e {θ:F
(11)
for fixed te > 0 and p = F (te ; θ). Using the result in (11), a 100(1 − α)% likelihood based CI for p is [ p, e
pe ]
11
(12)
where p = min u : u satisfying e pe = max u : u satisfying
e ≥ L(θ) b − 1 χ2 L(θ) 1;1−α e (te ; θ)=u} e 2 {θ:F e ≥ L(θ) b − 1 χ2 max L(θ) . e (te ; θ)=u} e 2 1;1−α {θ:F max
Interestingly, it can be shown that the CB for cdf obtained by using (12) and the CB for the quantile function obtained by using (10) are identical. Hong, Meeker, and Escobar (2008b) presented similar equivalence result but only for the cdf and quantile function of a single distribution.
5
Simulation Studies
In this section, we conduct a small simulation study to investigate several properties of the proposed procedure and other existing procedures. Our simulation is illustrative rather than comprehensive.
5.1
Simulation Setup and Evaluation Criteria
We consider a scenario that is similar to the Device G application. The details of this application are given in Section 6.1. In particular, we consider a system that fails due to competing risks. The system has two failure modes which are denoted by FM1 and FM2. Each failure mode has a failure-time distribution that follows a Weibull distribution. The parameters for FM1 and FM2 are θ 1 = (6.1, 1.5)′ and θ 2 = (5.8, 0.23)′, respectively. Failure-time data are available at the component level and the ML method is used to estimate θ 1 and θ 2 . In the simulation set up, the sample size for each failure mode is n. We chose n = 10, 20, 50, 100, 200, 500, 1,000 to see the effect of sample size. Type I censoring is used in the simulation. We examine the following combinations of censoring proportions for failure modes FM1 and FM2, respectively: (75%, 25%), (50%, 50%), and (25%, 75%). We consider the following procedures to construct CIs for F (te ) in the simulation; the ye procedure, the Fb procedure, logit transform procedure, and the likelihood based procedure. We use the following criteria in the comparisons; coverage probability, percentage of bending back, and percentage of CI with endpoints outside [0, 1].
5.2
Simulation Results
Table 1 shows the estimated probability of bending back and being outside [0,1] for various CI procedures, based on 10,000 trials. The bend back behavior occurs for the ye procedure, the 12
Fb procedure, and logit transform procedure when the sample size is small. The ye procedure,
however, has the smallest probability of bending back among these normal approximation CI procedures. The logit transformation procedure suffers seriously with bend back behavior.
When the sample size is more than 200, the bend-back behavior does not occur with the ye procedure. Theoretically we know the probability of bending back for the likelihood procedure is zero. Only the Fb procedure results in CIs with endpoints outside [0, 1] and this phenomena
happens with probability one. This abnormal behavior usually happens when the estimated failure probabilities are close to zero or one. Of course, one can always correct those values less than zero to zero and those values larger than one to one. Figure 2 shows the estimated coverage probability for those four procedures with different configurations of sample size n and censoring proportions, based on 10,000 trials. For each sub-figure, the upper, middle, and lower panels show the results when the censoring proportions are (75%, 25%), (50%, 50%), and (25%, 75%), respectively. The Fb procedure
has poor performance in term of coverage probability. When the sample size is small (e.g., n ≤ 100), the performance of the ye procedure and the logit procedure are reasonably good, while the likelihood based procedure has the best performance. When the sample size is large (e.g., n ≥ 500), the coverage probabilities of the ye procedure, the logit procedure, and the likelihood based procedure are all close to the nominal coverage probability. Although the likelihood based procedure performs best, the computational efficiency is of concern in practice. The computing time of the ye procedure, the Fb procedure, and logit transform procedure is nearly negligible (i.e., less than 0.1 second). The likelihood based
procedure is computationally intensive because numerical optimizations are involved. The
computing time of the likelihood based procedure also depends the number of time points where the CIs for F (te ) are needed, and the size of the dataset. For example, it takes about half hour to compute the likelihood CIs in Figure 4. In software implementations of these kinds of computations, users expect the results come out quickly. Thus, normal approximation CI procedures with good properties are still widely used in practice.
6
Applications
In this section, we apply to the proposed procedures to three applications; two applications in competing-risk (series system) models and an application for a k-out-of-s system.
13
Table 1: Estimated probability of bending back and being outside [0,1] for various CI procedures, based on 10,000 trials. n
20
50
100
200
500
1000
censoring proportion
bend back
FM1
FM2
ye
0.75
0.25
0.719
0.50
0.50
0.25
Fb
outside [ 0, 1 ]
logit
likelihood
ye
1
0.996
0
0
0.368
1
0.999
0
0.75
0.521
1
0.995
0.75
0.25
0.264
1
0.50
0.50
0.023
0.25
0.75
0.75
Fb
logit likelihood
1
0
0
0
1
0
0
0
0
1
0
0
0.314
0
0
1
0
0
1
0.996
0
0
1
0
0
0.036
1
0.999
0
0
1
0
0
0.25
0.012
1
0
0
0
1
0
0
0.50
0.50
0
1
0.533
0
0
1
0
0
0.25
0.75
0
1
0.999
0
0
1
0
0
0.75
0.25
0
1
0
0
0
1
0
0
0.50
0.50
0
1
0
0
0
1
0
0
0.25
0.75
0
1
0.967
0
0
1
0
0
0.75
0.25
0
1
0
0
0
1
0
0
0.50
0.50
0
1
0
0
0
1
0
0
0.25
0.75
0
1
0.001
0
0
1
0
0
0.75
0.25
0
1
0
0
0
1
0
0
0.50
0.50
0
1
0
0
0
1
0
0
0.25
0.75
0
1
0
0
0
1
0
0
14
1.00
1.00
0.95
0.95
0.90
0.90
0.85
0.85
Estimated Coverage Probability
Estimated Coverage Probability
replacemen
1.00
0.95
0.90
0.85
1.00
0.95
0.90
0.85
1.00
1.00
0.95
0.95
0.90
0.90
0.85
0.85 1
2
5
10
20
50
100
200
500
1
2
5
10
Time
50
100
200
500
50
100
200
500
200
500
(b) n = 50
1.00
1.00
0.95
0.95
0.90
0.90
0.85
0.85
Estimated Coverage Probability
Estimated Coverage Probability
(a) n = 20
1.00
0.95
0.90
0.85
1.00
0.95
0.90
0.85
1.00
1.00
0.95
0.95
0.90
0.90
0.85
0.85 1
2
5
10
20
50
100
200
500
1
2
5
10
Time
20
Time
(c) n = 100
(d) n = 200
1.00
1.00
0.95
0.95
0.90
0.90
0.85
0.85
Estimated Coverage Probability
Estimated Coverage Probability
20
Time
1.00
0.95
0.90
0.85
1.00
0.95
0.90
0.85
1.00
1.00
0.95
0.95
0.90
0.90
0.85
ye procedure ^ F procedure
logit transformation likelihood procedure
0.85 1
2
5
10
20
Time
(e) n = 500
50
100
200
500
1
2
5
10
20
50
100
Time
(f) n = 1000
Figure 2: Estimated coverage probabilities for the four procedures, based on 10,000 trials. For each sub-figure, the upper, middle, and lower panels show the results when the censoring proportions are (75%, 25%), (50%, 50%), and (25%, 75%), respectively. 15
Table 2: ML Estimates for the parameters in the failure-time distribution for the components in Device G. Failure Mode Surge Wearout
6.1
n
r
Distribution
µ b
σ b
b11 λ
b12 λ
b22 λ
30 15
Weibull
6.108 1.490 2.462 0.876 1.659
30
Weibull
5.830 0.231 6.365 2.642 3.358
7
Device G Data
Meeker and Escobar (1998, page 383) described the Device G data, providing failure-time information for a sample of 30 units from a field tracking study. Device G has two failure modes: surge and wearout. The surge failures correspond to the failure of a particular unprotected electronic component, which results from an accumulation of randomly occurring damage from power-line voltage spikes during electric storms. The wearout failures were caused by normal product wear over the time. Because the surge failure mode is caused by random damage events, it is reasonable to assume that the two failure modes are independent. The device fails due to the competing risks of surges and wearout. One can only observe the failure time for one failure mode while the failure time of the other failure modes are treated as censored observations, which is typical in competing-risk data setting. Thus there are data for each failure mode (i.e., data are available at component level). Table 2 shows the sample size, the number of failures from each failure mode, and the ML estimates for the parameters in the failure-time distribution for the components. This is an example where the sample size is relatively small. Because Device G is a series system, the g function is given in (2) with s = 2. Figure 3 shows the CBs for Device-G data computed by using various procedures. The bend-back behavior occurs for lower boundary of the pointwise CB for the logit transformation procedure. The endpoints of the CIs computed by the Fb procedure are outside [0, 1]. Both the ye procedure and the likelihood based procedure behave normally.
6.2
Product D Data
The background for Product D application is given in Hong and Meeker (2010). Product D is an electro-mechanical system used in offices or residences. The data were obtained from the manufacturer’s warranty database (simulated data were, however, used in Hong and Meeker 2010 to protect sensitive information). Product D can be thought of as a copying machine where the failure can be due to electronic components or mechanical components. Product D fails from causes that can be categorized into one of four major failure mode groups, which 16
1.0 0.6 0.0
0.2
0.4
Probability
0.8
estimated cdf ye procedure ^ F procedure logit transformation likelihood procedure
1
2
5
10
20
50
100
200
500
time
Figure 3: CBs for Device G data computed by using various procedures.
17
Table 3: ML Estimates for the parameters in the failure-time distribution for the components in Product D.
58.426
33.821
1.109 161.645
b12 λ
46.454
15.985
Lognormal
20.871 7.677 156.381
72.217
35.065
Lognormal
10.628 1.872 421.395 169.679 72.456
µ b
σ b
b11 λ
Failure Mode
n
r
Distribution
FM1
986
9
Lognormal
8.222
0.743 129.137
FM2
986 37
Weibull
9.401
FM3
986 22 986
FMOther
6
b22 λ
are coded as failure mode 1 (FM1), failure mode 2 (FM2), failure mode 3 (FM3), and all other failure modes (FMOther). Thus Product D also fails due to competing risks of the four failure modes. One can only observe the failure time for the first failure mode while the failure times of the rest failure modes are treated as censored observations. Similar to Device G, there are data are available at the component level. Table 3 shows the sample size, the number of failures from each failure mode, the ML estimates for the parameters in the failure-time distribution for the components. This is an example where the sample size is relatively large. Because Product D is also a series system, the g function is given in (2) with s = 4. Figure 4 shows the CBs for Product D data computed by using various procedures. The endpoints of the CIs computed by the Fb procedure are outside [0, 1]. All other procedures
behave normally because the sample size and the total number of failures are relative large in this example. Note that the likelihood and the ye procedure have excellent agreement.
6.3
k-out-of-s System
For a k-out-of-s system with s components, the system operates when there are k or more than k components operate. For example, an airplane that has four engines requires at least two engines to be functioning for the aircraft to function properly. This is an example 2-out-of-4 system. For another example, consider a system of eight pumps of which requires at least four pumps function properly for the system to success. Thus this is an example of a 4-out-of-8 system. The relationship between the system failure-time cdf and the component failure-time cdf is given in (3). This is an example where the relationship function g is complicated. More details on k-out-of-s systems are available in Kuo and Zuo (2003, Chpater 7). For illustration, we consider a 4-out-of-8 system. We simulated failure-time data at the component level. The distribution and sample size used in the simulation can be found 18
1.0 0.6 0.0
0.2
0.4
Probability
0.8
estimated cdf ye procedure ^ F procedure logit transformation likelihood procedure
10
50
100
500
5000
50000
time
Figure 4: CBs for Product D data computed by using various procedures.
19
Table 4: ML Estimates for the parameters in the failure-time distribution for the components 4-out-of-8 system. µ b
σ b
b λ11
b12 λ
b22 λ
Failure Mode
n
r
Distribution
FM1
20
11
Weibull
2.664 1.517 2.101 0.660 1.541
FM2
20
17
Lognormal
2.474 0.499 1.038 0.070 0.635
FM3
30
13
Weibull
4.002 0.581 3.301 1.407 1.993
FM4
30
6
Lognormal
3.057 0.879 5.681 3.604 3.425
FM5
50
30
Lognormal
0.230 8.530 1.274 0.359 0.992
FM6
50
24
Weibull
2.344 0.918 2.727 1.080 1.810
FM7
100 41
Weibull
4.846 6.485 3.755 1.691 2.173
FM8
100 18
Lognormal
2.979 0.298 6.935 4.473 4.036
in Table 4. To illustrate the general applicability of the proposed procedure, we allow components to have different failure-time distributions with different values of parameters (e.g., pumps made by different manufacturers may have different distributions in the pump system). Table 4 shows the ML estimates for the parameters in the failure-time distribution for the components in the 4-out-of-8 system. Figure 5 shows the CBs for the 4-out-of-8 system computed by using various procedures. Again, the endpoints of the CIs computed by the Fb procedure are outside [0, 1] while all other procedures behave normally.
7
Conclusions and Areas for Future Research
In this paper, we develop a general procedure for constructing the CI for the quantile function by using the implicit delta method. We also develop a general procedure constructing the CI for the system cdf. We show that the proposed procedures are asymptotically valid. Simulation studies also show that the proposed procedure has desirable properties. We also illustrate the proposed procedures with applications in competing-risk models and k-out-of-s system. There are several possible areas for future research. It is straightforward to derive specific results for non-location-scale distribution using the ye procedure. In this paper, we assume the data from component levels are independent. In the future, one can consider a situation in which the data from the component levels are dependent. It would also be useful to have CI procedure for nonparametric estimation of the system cdf by using approach similar to the ye procedure. 20
1.0 0.6 0.4 0.0
0.2
Probability
0.8
estimated cdf ye procedure ^ F procedure logit transformation likelihood procedure
2
5
10
20
50
100
time
Figure 5: CBs for 4-out-of-8 cdf computed by using various procedures.
21
A
Proof
In this appendix, we provide a proof for Result 2. When we consider the asymptotic behavior, we assume n1 , . . . , ns are of the same order and each of the ri /ni is held as constant where ri is the number of failures or the expected number of failures for the ith component. To show that the CI procedure in (8) is asymptotically valid, we need to prove that the true coverage probability (TCP)
Let
n o TCP = Pr Fb[exp(y e )] ≤ F (te ) ≤ Fb[exp(e ye )] e n o n o = Pr F (te ) ≤ Fb[exp(e ye )] − Pr F (te ) ≤ Fb[exp(y e )] e → 1 − α as n → ∞. s ∂ Fb[exp(ye )] X ∂g = fi (ye ) ∂ye ∂F i i=1 v u s 2 2 uX ∂g σ bi = se b pb = t fi (ye ) ∂Fi ni i=1
∆1 =
∆2
and b11i + 2λ b12i zbie + λ b22i zb2 . λ ie
Here ∆2 is op (1) (i.e., ∆2 can be small enough when the sample size n is large enough). Then, yee = ye + z1−α/2 ∆2 /∆1 . Consider the Taylor series expansion of Fb[exp(e ye )] at ye
2 ∂ Fb[exp(ye )] ∆2 ∂ 2 Fb[exp(ye )] ∆2 b b F [exp(e ye )] = F [exp(ye )] + z1−α/2 + (ξ1e ) z1−α/2 ∂ye ∆1 ∂ye2 ∆1 = Fb(te ) + z1−α/2 ∆2 + D1 ∆2 . (13)
Here ξ1e is between ye and ye + z1−α/2 ∆2 /∆1 , and
∂ 2 Fb[exp(ye )] ∆2 (ξ1e ) z1−α/2 2 . D1 = ∂ye2 ∆1
Similarly, consider the Taylor series expansion of Fb[exp(y e )] at ye e Fb[exp(y e )] = Fb(te ) − z1−α/2 ∆2 + D2 ∆2 , e where
∂ 2 Fb[exp(ye )] ∆2 D2 = (ξ2e ) z1−α/2 2 , ∂ye2 ∆1
22
(14)
and ξ2e is between ye − z1−α/2 ∆2 /∆1 and ye . From (13), n o h i b b Pr F (te ) ≤ F [exp(e ye )] = Pr F (te ) ≤ F (te ) + z1−α/2 ∆2 + D1 ∆2 " # F (te ) − Fb(te ) = Pr ≤ z1−α/2 + D1 ∆2 " # F (te ) − Fb(te ) = Pr − D1 ≤ z1−α/2 . se b pb From (14),
n o h i Pr F (te ) ≤ Fb[exp(y e )] = Pr F (te ) ≤ Fb(te ) − z1−α/2 ∆2 + D2 ∆2 e " # F (te ) − Fb(te ) = Pr ≤ −z1−α/2 + D2 ∆2 " # b F (te ) − F (te ) − D2 ≤ −z1−α/2 . = Pr se b pb
Assuming that ∂ 2 g/∂Fi2 and ∂fi /∂ye are bounded in the neighborhood of ye (conditions that can be easily satisfied for most structure functions g and failure-time distributions for components), D1 and D2 are op (1). Because F (te ) − Fb(te ) d −→ N(0, 1), as n → ∞, se b pb
and D1 and D2 are op (1), by Slutsky’s theorem (e.g., page 290 of Athreya and Lahiri 2006), " # " # F (te ) − Fb(te ) F (te ) − Fb(te ) TCP = Pr − D1 ≤ z1−α/2 − Pr − D2 ≤ −z1−α/2 se b pb se b pb α α → 1 − − = 1 − α, as n → ∞. 2 2 This proves that, in large samples, the CI procedure in (8) has an approximate coverage probability 1 − α.
References Athreya, K. B. and S. N. Lahiri (2006). Measure Theory and Probability Theory. New York: Springer. Coit, D. W. (1997). System-reliability confidence-intervals for complex-systems with estimated component-reliability. IEEE Transactions on Reliability 46, 487–493.
23
Crowder, M. J. (2001). Classical Competing Risks. Florida, Boca Raton: Chapman & Hall/CRC. David, H. A. and M. L. Moeschberger (1978). The Theory of Competing Risks. London: Griffin. Easterling, R. G. (1972). Approximate confidence limits for system reliability. Journal of the American Statistical Association 67, 220–222. Hamada, M. S., A. Wilson, C. S. Reese, and H. Martz (2008). Bayesian Reliability. New York: Springer. Hong, Y. (2011). On computing the distribution function for the sum of independent and non-identically distributed random indicators. Technical Report No. 112: Department of Statitics, Virginia Tech, Blacksburg, Virginia 24061. Available at http://www.stat.vt.edu/tech reports/index.html. Hong, Y. and W. Q. Meeker (2010). Field-failure and warranty prediction based on auxiliary use-rate information. Technometrics 52, 148–159. Hong, Y., W. Q. Meeker, and L. A. Escobar (2008a). Avoiding problems with normal approximation confidence intervals for probabilities. Technometrics 50, 64–68. Hong, Y., W. Q. Meeker, and L. A. Escobar (2008b). The relationship between confidence intervals for failure probabilities and life time quantile. IEEE Transactions on Reliability 57, 260–266. JMP (2011). JMP 9. http://www.jmp.com/. Kuo, W. and M. Zuo (2003). Optimal Reliability Modeling: Principles and Applications. Hoboken, NJ: John Wiley & Sons, Inc. Lehmann, E. L. and G. Casella (1998). Theory of Point Estimation (Second ed.). New York: Springer-Verlag. Meeker, W. Q. and L. A. Escobar (1998). Statistical Methods for Reliability Data. New York: John Wiley & Sons, Inc. Meeker, W. Q. and L. A. Escobar (2008). SPLIDA (S-PLUS Life Data Analysis). www.public.iastate.edu/∼ splida. Minitab (2011). MINITAB Statistical Software, Release 16. State College, Pennsylvania: Minitab Inc. Myhre, J. M. and S. C. Saunders (1968). On confidence limits for the reliability of systems. Annals of Mathematical Statistics 39, 1463–1472.
24
Pawitan, Y. (2001). In All Likelihood: Statistical Modelling and Inference Using Likelihood. USA: Oxford University Press. Ramirez-Marquez, J. E. and W. Jiang (2006). On improved confidence bounds for system reliability. IEEE Transactions on Reliability 55, 26–36. Rausand, M. and A. Høyland (2004). System Reliability Theory: Models, Statistical Methods, and Applications (Second ed.). Hoboken, New Jersey: John Wiley & Sons, Inc. ReliaSoft (2011). Weibull++ 7.5. http://weibull.reliasoft.com/. SAS (2011). SAS/QC User’s Guide, Version 9.3. NC: SAS Institute Inc. Spall, J. C. (2009). System reliability estimation and confidence regions from subsystem and full system tests. In American Control Conference, 2009, St. Louis, MO, pp. 5067– 5072.
25
William Q. Meeker
Department of Statistics
Department of Statistics
Virginia Tech
Iowa State University
Blacksburg, VA 24061
Ames, IA 50011 Abstract
System reliability depends on the reliability of the system’s components and the structure of the system. For example, in a competing-risk model, the system fails when the weakest component fails. The reliability function and the quantile function of a complicated system are two important metrics for characterizing the system’s reliability. When there are data available at the component level, the system reliability can be estimated by using the component level information. Confidence intervals (CI) are needed to quantify the statistical uncertainty in the estimation. Obtaining system reliability CI procedures with good properties is not straightforward, especially when the system structure is complicated. In this paper, we develop a general procedure for constructing a CI for the system failure-time quantile function by using the implicit delta method. We also develop general procedures for constructing a CI for the cdf of the system. We show that the recommended procedures are asymptotically valid and have good statistical properties. We conduct simulations to study the finitesample coverage properties of the proposed procedures and compare them with existing procedures. We apply the proposed procedures to three applications; two applications in competing-risk models and an application with a k-out-of-s system. The paper concludes with some discussion and an outline of areas for future research. Key Words: Asymptotic approximation, failure-time data, k-out-of-s system, maximum likelihood, normal approximation, series system.
1
1
Introduction
1.1
Background
For complicated systems, failures occurring at the component level can lead to system failures. For example, in a competing-risk setting (i.e., series system), a system fails when the first component fails. Meeker and Escobar (1998, page 383) described an application of Device G, where the system fails due to one of two competing failure modes. Hong and Meeker (2010) described another application of Product D, where the system fails due to one of four competing failure modes. Thus the reliability of a complicated system depends on the reliability of its components. The relationship between system reliability and component reliability is often of interest. In the competing risk setting, system reliability is the product of reliability of individual components. The relationship can be complicated in some systems such as the k-out-of-s system, where a system of s components operates when k or more of the system’s components operate. Nevertheless, with the relationship between system reliability and component reliability and information on component reliability, one can obtain information on system reliability. The reliability function and the quantile function of a system are two important metrics for characterizing system reliability. The reliability function is defined by R(t) = 1 − F (t) where F (t) is the cumulative distribution function (cdf) of the system failure-time random variable T . The cdf F (t) is interpreted as the probability that a unit will fail by time t, or the proportion of units in the population that will fail by time t. The quantile function tp is the inverse of the cdf, and corresponds to the time at which a specified proportion p of the population of systems will fail. When there are data available at the component level, these metrics can be estimated by using the component level data. Confidence intervals (CIs) are needed to quantify the statistical uncertainty in the estimation of F (t) and tp . In practice, it is common to plot on one graph an estimate and a set of pointwise CIs for F (t) over a range of t values. These sets of pointwise CIs are referred to as a pointwise confidence band (CB). Obtaining a CI procedure for the system reliability with good statistical properties is not straightforward, especially when the system structure is complicated. The following properties of the CI procedure should be considered. • The coverage probability of the CI procedure should be approximately equal to the nominal coverage probability for small or moderate sample size and converge to the nominal coverage probability when the sample size goes to infinity. • The endpoints of the CI should fall inside the parameter space, because the cdf is bounded by 0 and 1. Some normal approximation CI procedures may produce CIs 2
with endpoints that fall outside [0, 1]. • A transformation function can help restrict the CI endpoints to be within [0, 1] but may adversely affect the monotonicity of the boundary of the pointwise CBs. Hong, Meeker, and Escobar (2008a) considered the abnormal behavior of pointwise CBs caused by the use of an inappropriate transformation function, resulting in a boundary of the CB is not monotone increasing. That abnormal behavior was called the “bend-back” behavior. The bend-back behavior is not desirable in applications and needs to be avoided where possible. • In many applications, especially in software implementations, computational efficiency of the CI procedure is also important. Although CI procedures based on likelihood or simulation may have important advantages, normal approximation CI procedures are widely used due to their computational efficiency. The normal approximation CI procedures for F (t) of systems with simple structures are widely implemented in popular software packages. These software package includes JMP (2011), Minitab (2011), WEIBULL++ 7.5 by ReliaSoft (2011), PROC RELIABILITY in SAS (2011), and SPLIDA/RSPLIDA described in Meeker and Escobar (2008). In this paper, we develop a general procedure for constructing CIs for the quantile function by using the implicit delta method. We also develop a general procedure constructing CIs for the system cdf. We show that the proposed procedures are asymptotically valid and have good statistical properties. The proposed procedures can also be easily implemented in most software packages.
1.2
Related Literature
Rausand and Høyland (2004) provides a comprehensive discussion of system reliability. In the area of system reliability CIs, some earlier work includes Myhre and Saunders (1968) and Easterling (1972), which are based on binary data and a simple series systems. Coit (1997) used a log transformation to construct CI for the reliability function using binary or lifetime data collected at component level. Meeker and Escobar (1998) used the logit transformation to construct CI for the cdf of the system. Ramirez-Marquez and Jiang (2006) considered a CI procedure by using a recursive approach. Spall (2009) constructed a CI for system reliability by using the information from both full system and component level tests. At the component level, a failure-time distribution is needed for the component reliability. Meeker and Escobar (1998) is a useful resource for statistical methods for reliability data. Hamada et al. (2008) provide statistical methods for analyzing reliability data from a 3
Bayesian perspective. For systems failing due to competing risks, David and Moeschberger (1978), and Crowder (2001) are useful books on the subject of competing-risk models.
1.3
Overview
The rest of the paper is organized as follows. Section 2 introduces the model for the system and maximum likelihood (ML) estimation for unknown parameters. Section 3 presents the proposed CI procedures for the cdf and the quantile function, based on normal approximations. Section 4 presents likelihood based CI procedures for the cdf and the quantile function. Section 5 describes the results of a simulation study to evaluate the performance of several CI procedures, using various criteria. Section 6 applies the proposed CI procedures to three different applications. Section 7 contains some conclusions and areas for future research.
2
Model and ML Estimation
2.1
System Models
We consider a coherent system with s components. A coherent systems is a system that meets two requirements: 1) the improvement of any component will improve the system reliability, and 2) every component in the system makes some nonzero contribution to the system’s reliability. Let T be the failure time of the system. The system cdf is defined by F (t; θ) = Pr(T ≤ t) where θ is the vector for the parameters in the system cdf. The cdf for the failure time of component i is denoted by Fi (t; θ i ), i = 1, . . . , s. Here θ i is a vector for the parameter in the cdf of component i. The relationship between a system’s cdf and components’ cdfs is characterized by a system structure function g. The cdf for failure-time distribution of the system is expressed as F (t; θ) = g[F1 (t; θ 1 ), . . . , Fs (t; θ s )].
(1)
Here θ = (θ ′1 , . . . , θ ′s )′ and g = g(x1 , · · · , xs ) is a function mapping from [ 0, 1 ]s (an s dimensional unit hypercube) to [ 0, 1 ]. We assume the function g is a strictly increasing and twice differentiable function with respect to each of its arguments, which can be met for most common system structures. The form of the function g is determined by system’s structure. Some common examples include the following (see Rausand and Høyland 2004 for a comprehensive discussion).
4
• For a series system with s components, g(x1 , . . . , xs ) = 1 −
s Y
(1 − xi ).
(2)
i=1
• For a parallel system with s components, g(x1 , . . . , xs ) =
Qs
i=1
xi .
• For a 2 × 2 series-parallel system with component-level redundancy, g(x1 , x2 , x3 , x4 ) = 1 − (1 − x1 x2 )(1 − x3 x4 ). • For a k-out-of-s system with s components, the system fails when there are s − k + 1 or more components fail. The g function is ( ) s s 1 X exp[−iωl(s − k + 1)] − exp(−iωl) Y g(x1 , . . . , xs ) = [1 − xi + xi exp(iωl)] . s + 1 l=0 1 − exp(−iωl) i=1 (3)
The details for the derivation of equation (3) can be found in Hong (2011). Depending on the system structure, the form of g may be more complicated than these examples. The CI procedures proposed in this paper works for a general form of g. Note that we use θ = (θ ′1 , . . . , θ ′s )′ to denote the unknown parameters in the system. For some system structures (e.g., when a system has several copies of the same component in the system structure), it is possible that the components share the same parameters. Thus there may be some redundancy in parameter θ. The extension of our proposed procedure to handle this situation is straightforward. For notational simplicity, we assume that each component has its own set of parameters. In addition to the cdf F (t; θ), the quantile function of the system failure-time distribution is another quantity that is often of interest. The p quantile is defined by the inverse of the cdf tp = F −1 (p). Numerical methods may be needed to solve tp from F (tp ; θ) = p. The log of the failure time T is denoted by Y = log(T ). The quantile function of Y is denoted by yp = log(tp ). Note that both tp and yp depend on the unknown parameter vector θ. That is tp = tp (θ) and yp = yp (θ).
2.2
Component Failure-time Model
Here we use the log-location-scale family of distributions to model the failure-time distribution of the components. The Weibull and lognormal distributions, which are the most commonly-used distributions for product/system failure times, are members of this family.
5
The cdf and probability density function (pdf) of a log-location-scale distribution can be expressed as
log(t) − µ F (t; µ, σ) = Φ σ
1 log(t) − µ f (t; µ, σ) = φ , σt σ
and
where Φ is the standard cdf for the location-scale family of distributions (location 0 and scale 1), µ is the location parameter, and σ is the scale parameter. For the lognormal distribution, replace Φ and φ above with Φnor and φnor , the standard normal cdf and pdf, respectively. For the Weibull distribution, replace Φ and φ above with Φsev (z) = 1 − exp[− exp(z)] and φsev (z) = exp[z − exp(z)], respectively. The cdf of the component failure time is denoted by Fi (t; θ i ) = Φi (zi ), i = 1, 2, · · · , s (as the form of the distribution need not be the same for all components) where θ i = (µi , σi )′ and zi = [log(t) − µi ]/σi . Here Φi (·) is a standard location-scale cdf. Although the loglocation-scale family of distribution is used to model the failure time of the components, the CI procedures developed in this paper are general. We provide specific results for the log-location-scale family of distributions because this family of distributions is widely used.
2.3
Data and Maximum Likelihood Estimation
Suppose that data are collected at the component level from either a life test and/or the field. The data for component i is denoted by Di . These data are possibly censored. A loglocation-scale distribution is used to describe the failure time of component i. The maximum likelihood (ML) estimators of parameters of component i based on data Di , denoted by bi , i = 1, . . . , s. From large-sample theory, θ bi ∼ θ ˙ N(θ i , Σi ). Here Σi is the variance-covariance b To estimate Σi , we use the observed information matrix which can be written matrix of θ. as
bi2 bi = σ Σ ni
b11i λ b12i λ
b12i λ b22i λ
!
, i = 1, 2, · · · , s,
b11i , b where ni is the number of observations in the data Di , and λ λ12i and b λ22i are elements in the scaled variance-covariance matrix.
b = (θ b′ , . . . , θ b′ )′ . We assume The ML estimator for the system parameters is obtained by θ 1 s bi based on Di are independent from each other. Thus θ b is asymptotically the estimators θ multivariate normal (MVN) distributed. That is,
b∼ θ ˙ MVN [θ, Diag (Σ1 , Σ2 , · · · , Σs )] .
By the invariance property of ML estimators, the ML estimator of the distribution probability at te is b = g[F1 (te ; θ b1 ), . . . , Fs (te ; θ bs )]. Fb(te ) = F (te ; θ) 6
b Similarly, the ML estimator of tp is b tp = exp(b yp ) where ybp = yp (θ).
3
Normal Approximation Confidence Interval Procedures
3.1
Normal Approximation Confidence Interval Procedure for tp
An estimate of the variance of ybp based on the delta method and observed information is ∂b yp yp d (b b b b ∂b Var yp ) = . ′ Diag(Σ1 , Σ2 , · · · , Σs ) b b ∂θ ∂θ
(4)
In many applications, there is no closed-form expression for yp = yp (θ). Thus we use the b Note that implicit derivative technique to find ∂b yp /∂ θ. p = g {F1 [exp(b yp ); θ 1 ], . . . , Fs [exp(b yp ); θ s ]} .
(5)
Taking derivatives on both sides of (5) with respect to µ bi , gives X ∂b yp ∂g ∂b yp ∂g 0= fi (b yp ) −1 + fj (b yp ) , ∂Fi ∂b µi ∂Fj ∂b µi j6=i yp ) = φi (b zip )/b σi . Here φi is the standard location-scale pdf of component i and where fi (b zbip = (b yp − µ bi )/b σi , i = 1, 2, · · · , s. Hence, ∂b yp ∂b µi
∂g fi (b yp ) ∂Fi = , Ps ∂g fi (b yp ) i=1 ∂Fi
i = 1, · · · , s,
P provided that si=1 (∂g/∂Fi )fi (b yp ) 6= 0. Because g is a strictly increasing function, Ps yp ) 6= 0. Similarly, i=1 (∂g/∂Fi )fi (b ∂b yp ∂b σi From (4),
d (b Var yp ) =
∂g fi (b yp )b zip ∂Fi = , Ps ∂g fi (b yp ) i=1 ∂Fi
Ps
i=1
2 2 ∂g σ bi fi (b yp ) ∂Fi ni Ps i=1
7
i = 1, 2, · · · , s.
2 b b b λ11i + 2λ12i zbip + λ22i zbip . 2 ∂g fi (b yp ) ∂Fi
Using the large-sample result that yb − yp qp ∼ ˙ N(0, 1), d Var (b yp )
we obtain the following CI procedure for yp .
Result 1 A large-sample approximate 100(1 − α)% CI for yp is s 2 Ps ∂g 2 b b12i zbip + λ b22i zb2 f (b y ) σ b γ λ + 2 λ i p i 11i i ip i=1 h i ∂Fi y p , yep = ybp ∓ , Ps ∂g e f (b y ) i p i=1 ∂Fi
(6)
2 where γi = z1−α/2 /ni , i = 1, 2, · · · , s. A large-sample approximate CI for tp is i h tp , e tp = exp(y p ), exp(e yp ) . e e Such approximate CIs based on Wald statistics may have a positive probability of bend-
back behavior (i.e., the boundary of the pointwise CB is not monotonically increasing, see Hong, Meeker, and Escobar 2008a for more details). We will study the probability of bending back for this CI procedure in Section 5.
3.2
Normal Approximation Confidence Interval Procedure for F (te )
A CI for the cdf at a particular time point t = te is often needed. Let p = F (te ; θ) and b Using the large-sample result that (b pb = F (te ; θ). p − p)/se b pb ∼ ˙ N(0, 1), a direct approach to construct a CI for p is h
i p, pe = pb − z1−α/2 se b pb, pb + z1−α/2 se b pb , (7) e where se b pb is the standard error which is obtained by using the delta method. In particular, v u s 2 2 uX ∂g σ bi b 2 t b b fi (ye ) λ11i + 2λ12i zbie + λ22i zbie . se b pb = ∂Fi ni i=1
The CI procedure in (7) is referred to as the Fb procedure. The Fb procedure may produce a CI with endpoints falling outside of [0, 1], which is not desirable. A logit or some other similar
transformations can be applied to restrict the CI endpoints to fall inside [0, 1] (e.g., page
190 of Meeker and Escobar 1998). Such a procedure, however, can have a high probability of the bend-back, as shown in Hong, Meeker, and Escobar (2008a) for a system with a single component. 8
The following alternative procedure, that we call ye procedure, will have CI endpoints that always fall inside the interval [0, 1] and has a much smaller probability of bending back. The construction of the CI, which is illustrated in Figure 1, is based, indirectly, on the endpoints obtained by using the Fb procedure. As shown in Figure 1, the endpoints of the CI obtained by the Fb procedure, reflected by the tangent line of Fb[exp(ye )] at ye = log(te ), are denoted by y e and yee , respectively. In particular, e s 2 h i Ps ∂g 2 2 b b b f (y ) σ b γ λ + 2 λ z b + λ z b i e 11i 12i ie 22i ie i i i=1 ∂Fi z1−α/2 se b pb y e = ye − = ye − , Ps ∂g ∂ Fb[exp(ye )] e fi (ye ) i=1 ∂Fi ∂ye s 2 h i Ps ∂g 2 b11i + 2λ b12i zbie + λ b22i zb2 f (y ) σ b γ λ i e i i ie i=1 ∂Fi z1−α/2 se b pb yee = ye + , = ye + Ps ∂g ∂ Fb[exp(ye )] fi (ye ) i=1 ∂Fi ∂ye
where fi (ye ) = φi (b zie )/b σi and zbie = (ye − µ bi )/b σi . The lower and upper endpoints of the b at exp(y e ) and exp(e proposed CI are obtained by evaluating F (t; θ) ye ), respectively. e Result 2 A large-sample approximate 100(1 − α)% CI for p = F (te ; θ) is h i h i p, pe = Fb{exp(y e )}, Fb{exp(e ye )} . e e
(8)
The proof of this result is given in Appendix A. The CI procedure in Result 2 may have a positive probability of bend-back behavior. We will study the probability of bending back for this CI procedure in Section 5. Hong, Meeker, and Escobar (2008a) studied a CI procedure for the cdf of the failuretime distribution of a system with a single component, which was called the zb procedure.
In particular, the zb procedure is based on the result that (b z1e − z1e )/se b zb1e ∼ ˙ N(0, 1), where
z1e = (ye − µ b)/b σ. The zb procedure has a small probability of bending back.
Here we show that the ye procedure in (8) is a generalization of the zb procedure. When
g(x) = x and s = 1 (i.e., there is only one component), the ye procedure (8) is identical to the zb procedure. Consider that ye e
r h i 1 2 2 2 b b b = ye − [f1 (ye )] σ b1 γ1 (λ111 + 2λ121 zb1e + λ221 zb1e ) f1 (ye ) r b111 + 2λ b121 zb1e + λ b221 zb2 = ye − σ b12 γ1 λ 1e 9
F (y)
Fb procedure
p=1
pb
p e
pe
z1−α/2 se b pb
t z1−α/2 se b pb gen n ta
e lin
pb
(ye , pb )
z1−α/2 se b pb
estimated cdf
p e
z1−α/2 se b pb
∂F ∂y
ye procedure
pe
∂F ∂y
p=0
ye e
ye
yee
y
Figure 1: Illustration of the proposed ye procedure to obtain a CI for F (te ). and q b111 + 2λ b121 zb1e + λ b221 zb2 ) − µ ye − σ b12 γ1 (λ b 1 1e p = Fb[exp(y e )] = Φ1 σ b 1 e e q b111 + 2λ b121 zb1e + λ b221 zb2 ) = Φ1 zb1e − γ1 (λ 1e
= Φ1 (z 1e ). e
Note that Φ1 (z 1e ) is the lower endpoint obtained by the zb procedure. Hence the lower e endpoints are identical. In a similar manner, one can show that the upper endpoints are identical. Thus the zb procedure is a special case of the ye procedure. 10
4
Likelihood Based Confidence Interval Procedures
In this section, we present likelihood based CI procedures for tp and F (te ). In general, there are no closed-form expressions for likelihood based CI procedures. Thus numerical methods are needed to compute the CIs and these can be computationally intensive. For this reason, in some applications, such as producing dynamic graphics (like those used in the profiler on the JMP reliability platform), the likelihood based CI procedures may not be suitable with currently available computer hardware.
4.1
Likelihood Based Confidence Interval Procedure for tp
Let L(θ) be the log-likelihood function of the data introduced in Section 2.3. Standard large-sample theory for an ML estimator, under standard regularity conditions (e.g., page 241 of Pawitan 2001, and Chapter 6 of Lehmann and Casella 1998), provides the result that ) ( b − e 2 L(θ) max L(θ) ∼ ˙ χ21 (9) e F (tp ; θ)=p} e {θ:
for fixed 0 < p < 1 and tp = tp (θ). The result in (9) provides the basis for a likelihood ratio test for tp . A 100(1 − α)% likelihood based CI for the tp , obtained by inverting the likelihood ratio test, is
where
tp , e
tp = min t : t satisfying e e tp = max t : t satisfying
e tp
(10)
e ≥ L(θ) b − 1 χ2 L(θ) e F (t; θ)=p} e 2 1;1−α {θ: e ≥ L(θ) b − 1 χ2 max L(θ) 1;1−α . e F (t; θ)=p} e 2 {θ: max
Here χ21;1−α is the 1 − α quantile of the χ21 distribution.
4.2
Likelihood Based Confidence Interval Procedure for F (te )
Standard large-sample theory for ML estimator also provides the result that ( ) b − e 2 L(θ) max L(θ) ∼ ˙ χ21 e (te ; θ)=p} e {θ:F
(11)
for fixed te > 0 and p = F (te ; θ). Using the result in (11), a 100(1 − α)% likelihood based CI for p is [ p, e
pe ]
11
(12)
where p = min u : u satisfying e pe = max u : u satisfying
e ≥ L(θ) b − 1 χ2 L(θ) 1;1−α e (te ; θ)=u} e 2 {θ:F e ≥ L(θ) b − 1 χ2 max L(θ) . e (te ; θ)=u} e 2 1;1−α {θ:F max
Interestingly, it can be shown that the CB for cdf obtained by using (12) and the CB for the quantile function obtained by using (10) are identical. Hong, Meeker, and Escobar (2008b) presented similar equivalence result but only for the cdf and quantile function of a single distribution.
5
Simulation Studies
In this section, we conduct a small simulation study to investigate several properties of the proposed procedure and other existing procedures. Our simulation is illustrative rather than comprehensive.
5.1
Simulation Setup and Evaluation Criteria
We consider a scenario that is similar to the Device G application. The details of this application are given in Section 6.1. In particular, we consider a system that fails due to competing risks. The system has two failure modes which are denoted by FM1 and FM2. Each failure mode has a failure-time distribution that follows a Weibull distribution. The parameters for FM1 and FM2 are θ 1 = (6.1, 1.5)′ and θ 2 = (5.8, 0.23)′, respectively. Failure-time data are available at the component level and the ML method is used to estimate θ 1 and θ 2 . In the simulation set up, the sample size for each failure mode is n. We chose n = 10, 20, 50, 100, 200, 500, 1,000 to see the effect of sample size. Type I censoring is used in the simulation. We examine the following combinations of censoring proportions for failure modes FM1 and FM2, respectively: (75%, 25%), (50%, 50%), and (25%, 75%). We consider the following procedures to construct CIs for F (te ) in the simulation; the ye procedure, the Fb procedure, logit transform procedure, and the likelihood based procedure. We use the following criteria in the comparisons; coverage probability, percentage of bending back, and percentage of CI with endpoints outside [0, 1].
5.2
Simulation Results
Table 1 shows the estimated probability of bending back and being outside [0,1] for various CI procedures, based on 10,000 trials. The bend back behavior occurs for the ye procedure, the 12
Fb procedure, and logit transform procedure when the sample size is small. The ye procedure,
however, has the smallest probability of bending back among these normal approximation CI procedures. The logit transformation procedure suffers seriously with bend back behavior.
When the sample size is more than 200, the bend-back behavior does not occur with the ye procedure. Theoretically we know the probability of bending back for the likelihood procedure is zero. Only the Fb procedure results in CIs with endpoints outside [0, 1] and this phenomena
happens with probability one. This abnormal behavior usually happens when the estimated failure probabilities are close to zero or one. Of course, one can always correct those values less than zero to zero and those values larger than one to one. Figure 2 shows the estimated coverage probability for those four procedures with different configurations of sample size n and censoring proportions, based on 10,000 trials. For each sub-figure, the upper, middle, and lower panels show the results when the censoring proportions are (75%, 25%), (50%, 50%), and (25%, 75%), respectively. The Fb procedure
has poor performance in term of coverage probability. When the sample size is small (e.g., n ≤ 100), the performance of the ye procedure and the logit procedure are reasonably good, while the likelihood based procedure has the best performance. When the sample size is large (e.g., n ≥ 500), the coverage probabilities of the ye procedure, the logit procedure, and the likelihood based procedure are all close to the nominal coverage probability. Although the likelihood based procedure performs best, the computational efficiency is of concern in practice. The computing time of the ye procedure, the Fb procedure, and logit transform procedure is nearly negligible (i.e., less than 0.1 second). The likelihood based
procedure is computationally intensive because numerical optimizations are involved. The
computing time of the likelihood based procedure also depends the number of time points where the CIs for F (te ) are needed, and the size of the dataset. For example, it takes about half hour to compute the likelihood CIs in Figure 4. In software implementations of these kinds of computations, users expect the results come out quickly. Thus, normal approximation CI procedures with good properties are still widely used in practice.
6
Applications
In this section, we apply to the proposed procedures to three applications; two applications in competing-risk (series system) models and an application for a k-out-of-s system.
13
Table 1: Estimated probability of bending back and being outside [0,1] for various CI procedures, based on 10,000 trials. n
20
50
100
200
500
1000
censoring proportion
bend back
FM1
FM2
ye
0.75
0.25
0.719
0.50
0.50
0.25
Fb
outside [ 0, 1 ]
logit
likelihood
ye
1
0.996
0
0
0.368
1
0.999
0
0.75
0.521
1
0.995
0.75
0.25
0.264
1
0.50
0.50
0.023
0.25
0.75
0.75
Fb
logit likelihood
1
0
0
0
1
0
0
0
0
1
0
0
0.314
0
0
1
0
0
1
0.996
0
0
1
0
0
0.036
1
0.999
0
0
1
0
0
0.25
0.012
1
0
0
0
1
0
0
0.50
0.50
0
1
0.533
0
0
1
0
0
0.25
0.75
0
1
0.999
0
0
1
0
0
0.75
0.25
0
1
0
0
0
1
0
0
0.50
0.50
0
1
0
0
0
1
0
0
0.25
0.75
0
1
0.967
0
0
1
0
0
0.75
0.25
0
1
0
0
0
1
0
0
0.50
0.50
0
1
0
0
0
1
0
0
0.25
0.75
0
1
0.001
0
0
1
0
0
0.75
0.25
0
1
0
0
0
1
0
0
0.50
0.50
0
1
0
0
0
1
0
0
0.25
0.75
0
1
0
0
0
1
0
0
14
1.00
1.00
0.95
0.95
0.90
0.90
0.85
0.85
Estimated Coverage Probability
Estimated Coverage Probability
replacemen
1.00
0.95
0.90
0.85
1.00
0.95
0.90
0.85
1.00
1.00
0.95
0.95
0.90
0.90
0.85
0.85 1
2
5
10
20
50
100
200
500
1
2
5
10
Time
50
100
200
500
50
100
200
500
200
500
(b) n = 50
1.00
1.00
0.95
0.95
0.90
0.90
0.85
0.85
Estimated Coverage Probability
Estimated Coverage Probability
(a) n = 20
1.00
0.95
0.90
0.85
1.00
0.95
0.90
0.85
1.00
1.00
0.95
0.95
0.90
0.90
0.85
0.85 1
2
5
10
20
50
100
200
500
1
2
5
10
Time
20
Time
(c) n = 100
(d) n = 200
1.00
1.00
0.95
0.95
0.90
0.90
0.85
0.85
Estimated Coverage Probability
Estimated Coverage Probability
20
Time
1.00
0.95
0.90
0.85
1.00
0.95
0.90
0.85
1.00
1.00
0.95
0.95
0.90
0.90
0.85
ye procedure ^ F procedure
logit transformation likelihood procedure
0.85 1
2
5
10
20
Time
(e) n = 500
50
100
200
500
1
2
5
10
20
50
100
Time
(f) n = 1000
Figure 2: Estimated coverage probabilities for the four procedures, based on 10,000 trials. For each sub-figure, the upper, middle, and lower panels show the results when the censoring proportions are (75%, 25%), (50%, 50%), and (25%, 75%), respectively. 15
Table 2: ML Estimates for the parameters in the failure-time distribution for the components in Device G. Failure Mode Surge Wearout
6.1
n
r
Distribution
µ b
σ b
b11 λ
b12 λ
b22 λ
30 15
Weibull
6.108 1.490 2.462 0.876 1.659
30
Weibull
5.830 0.231 6.365 2.642 3.358
7
Device G Data
Meeker and Escobar (1998, page 383) described the Device G data, providing failure-time information for a sample of 30 units from a field tracking study. Device G has two failure modes: surge and wearout. The surge failures correspond to the failure of a particular unprotected electronic component, which results from an accumulation of randomly occurring damage from power-line voltage spikes during electric storms. The wearout failures were caused by normal product wear over the time. Because the surge failure mode is caused by random damage events, it is reasonable to assume that the two failure modes are independent. The device fails due to the competing risks of surges and wearout. One can only observe the failure time for one failure mode while the failure time of the other failure modes are treated as censored observations, which is typical in competing-risk data setting. Thus there are data for each failure mode (i.e., data are available at component level). Table 2 shows the sample size, the number of failures from each failure mode, and the ML estimates for the parameters in the failure-time distribution for the components. This is an example where the sample size is relatively small. Because Device G is a series system, the g function is given in (2) with s = 2. Figure 3 shows the CBs for Device-G data computed by using various procedures. The bend-back behavior occurs for lower boundary of the pointwise CB for the logit transformation procedure. The endpoints of the CIs computed by the Fb procedure are outside [0, 1]. Both the ye procedure and the likelihood based procedure behave normally.
6.2
Product D Data
The background for Product D application is given in Hong and Meeker (2010). Product D is an electro-mechanical system used in offices or residences. The data were obtained from the manufacturer’s warranty database (simulated data were, however, used in Hong and Meeker 2010 to protect sensitive information). Product D can be thought of as a copying machine where the failure can be due to electronic components or mechanical components. Product D fails from causes that can be categorized into one of four major failure mode groups, which 16
1.0 0.6 0.0
0.2
0.4
Probability
0.8
estimated cdf ye procedure ^ F procedure logit transformation likelihood procedure
1
2
5
10
20
50
100
200
500
time
Figure 3: CBs for Device G data computed by using various procedures.
17
Table 3: ML Estimates for the parameters in the failure-time distribution for the components in Product D.
58.426
33.821
1.109 161.645
b12 λ
46.454
15.985
Lognormal
20.871 7.677 156.381
72.217
35.065
Lognormal
10.628 1.872 421.395 169.679 72.456
µ b
σ b
b11 λ
Failure Mode
n
r
Distribution
FM1
986
9
Lognormal
8.222
0.743 129.137
FM2
986 37
Weibull
9.401
FM3
986 22 986
FMOther
6
b22 λ
are coded as failure mode 1 (FM1), failure mode 2 (FM2), failure mode 3 (FM3), and all other failure modes (FMOther). Thus Product D also fails due to competing risks of the four failure modes. One can only observe the failure time for the first failure mode while the failure times of the rest failure modes are treated as censored observations. Similar to Device G, there are data are available at the component level. Table 3 shows the sample size, the number of failures from each failure mode, the ML estimates for the parameters in the failure-time distribution for the components. This is an example where the sample size is relatively large. Because Product D is also a series system, the g function is given in (2) with s = 4. Figure 4 shows the CBs for Product D data computed by using various procedures. The endpoints of the CIs computed by the Fb procedure are outside [0, 1]. All other procedures
behave normally because the sample size and the total number of failures are relative large in this example. Note that the likelihood and the ye procedure have excellent agreement.
6.3
k-out-of-s System
For a k-out-of-s system with s components, the system operates when there are k or more than k components operate. For example, an airplane that has four engines requires at least two engines to be functioning for the aircraft to function properly. This is an example 2-out-of-4 system. For another example, consider a system of eight pumps of which requires at least four pumps function properly for the system to success. Thus this is an example of a 4-out-of-8 system. The relationship between the system failure-time cdf and the component failure-time cdf is given in (3). This is an example where the relationship function g is complicated. More details on k-out-of-s systems are available in Kuo and Zuo (2003, Chpater 7). For illustration, we consider a 4-out-of-8 system. We simulated failure-time data at the component level. The distribution and sample size used in the simulation can be found 18
1.0 0.6 0.0
0.2
0.4
Probability
0.8
estimated cdf ye procedure ^ F procedure logit transformation likelihood procedure
10
50
100
500
5000
50000
time
Figure 4: CBs for Product D data computed by using various procedures.
19
Table 4: ML Estimates for the parameters in the failure-time distribution for the components 4-out-of-8 system. µ b
σ b
b λ11
b12 λ
b22 λ
Failure Mode
n
r
Distribution
FM1
20
11
Weibull
2.664 1.517 2.101 0.660 1.541
FM2
20
17
Lognormal
2.474 0.499 1.038 0.070 0.635
FM3
30
13
Weibull
4.002 0.581 3.301 1.407 1.993
FM4
30
6
Lognormal
3.057 0.879 5.681 3.604 3.425
FM5
50
30
Lognormal
0.230 8.530 1.274 0.359 0.992
FM6
50
24
Weibull
2.344 0.918 2.727 1.080 1.810
FM7
100 41
Weibull
4.846 6.485 3.755 1.691 2.173
FM8
100 18
Lognormal
2.979 0.298 6.935 4.473 4.036
in Table 4. To illustrate the general applicability of the proposed procedure, we allow components to have different failure-time distributions with different values of parameters (e.g., pumps made by different manufacturers may have different distributions in the pump system). Table 4 shows the ML estimates for the parameters in the failure-time distribution for the components in the 4-out-of-8 system. Figure 5 shows the CBs for the 4-out-of-8 system computed by using various procedures. Again, the endpoints of the CIs computed by the Fb procedure are outside [0, 1] while all other procedures behave normally.
7
Conclusions and Areas for Future Research
In this paper, we develop a general procedure for constructing the CI for the quantile function by using the implicit delta method. We also develop a general procedure constructing the CI for the system cdf. We show that the proposed procedures are asymptotically valid. Simulation studies also show that the proposed procedure has desirable properties. We also illustrate the proposed procedures with applications in competing-risk models and k-out-of-s system. There are several possible areas for future research. It is straightforward to derive specific results for non-location-scale distribution using the ye procedure. In this paper, we assume the data from component levels are independent. In the future, one can consider a situation in which the data from the component levels are dependent. It would also be useful to have CI procedure for nonparametric estimation of the system cdf by using approach similar to the ye procedure. 20
1.0 0.6 0.4 0.0
0.2
Probability
0.8
estimated cdf ye procedure ^ F procedure logit transformation likelihood procedure
2
5
10
20
50
100
time
Figure 5: CBs for 4-out-of-8 cdf computed by using various procedures.
21
A
Proof
In this appendix, we provide a proof for Result 2. When we consider the asymptotic behavior, we assume n1 , . . . , ns are of the same order and each of the ri /ni is held as constant where ri is the number of failures or the expected number of failures for the ith component. To show that the CI procedure in (8) is asymptotically valid, we need to prove that the true coverage probability (TCP)
Let
n o TCP = Pr Fb[exp(y e )] ≤ F (te ) ≤ Fb[exp(e ye )] e n o n o = Pr F (te ) ≤ Fb[exp(e ye )] − Pr F (te ) ≤ Fb[exp(y e )] e → 1 − α as n → ∞. s ∂ Fb[exp(ye )] X ∂g = fi (ye ) ∂ye ∂F i i=1 v u s 2 2 uX ∂g σ bi = se b pb = t fi (ye ) ∂Fi ni i=1
∆1 =
∆2
and b11i + 2λ b12i zbie + λ b22i zb2 . λ ie
Here ∆2 is op (1) (i.e., ∆2 can be small enough when the sample size n is large enough). Then, yee = ye + z1−α/2 ∆2 /∆1 . Consider the Taylor series expansion of Fb[exp(e ye )] at ye
2 ∂ Fb[exp(ye )] ∆2 ∂ 2 Fb[exp(ye )] ∆2 b b F [exp(e ye )] = F [exp(ye )] + z1−α/2 + (ξ1e ) z1−α/2 ∂ye ∆1 ∂ye2 ∆1 = Fb(te ) + z1−α/2 ∆2 + D1 ∆2 . (13)
Here ξ1e is between ye and ye + z1−α/2 ∆2 /∆1 , and
∂ 2 Fb[exp(ye )] ∆2 (ξ1e ) z1−α/2 2 . D1 = ∂ye2 ∆1
Similarly, consider the Taylor series expansion of Fb[exp(y e )] at ye e Fb[exp(y e )] = Fb(te ) − z1−α/2 ∆2 + D2 ∆2 , e where
∂ 2 Fb[exp(ye )] ∆2 D2 = (ξ2e ) z1−α/2 2 , ∂ye2 ∆1
22
(14)
and ξ2e is between ye − z1−α/2 ∆2 /∆1 and ye . From (13), n o h i b b Pr F (te ) ≤ F [exp(e ye )] = Pr F (te ) ≤ F (te ) + z1−α/2 ∆2 + D1 ∆2 " # F (te ) − Fb(te ) = Pr ≤ z1−α/2 + D1 ∆2 " # F (te ) − Fb(te ) = Pr − D1 ≤ z1−α/2 . se b pb From (14),
n o h i Pr F (te ) ≤ Fb[exp(y e )] = Pr F (te ) ≤ Fb(te ) − z1−α/2 ∆2 + D2 ∆2 e " # F (te ) − Fb(te ) = Pr ≤ −z1−α/2 + D2 ∆2 " # b F (te ) − F (te ) − D2 ≤ −z1−α/2 . = Pr se b pb
Assuming that ∂ 2 g/∂Fi2 and ∂fi /∂ye are bounded in the neighborhood of ye (conditions that can be easily satisfied for most structure functions g and failure-time distributions for components), D1 and D2 are op (1). Because F (te ) − Fb(te ) d −→ N(0, 1), as n → ∞, se b pb
and D1 and D2 are op (1), by Slutsky’s theorem (e.g., page 290 of Athreya and Lahiri 2006), " # " # F (te ) − Fb(te ) F (te ) − Fb(te ) TCP = Pr − D1 ≤ z1−α/2 − Pr − D2 ≤ −z1−α/2 se b pb se b pb α α → 1 − − = 1 − α, as n → ∞. 2 2 This proves that, in large samples, the CI procedure in (8) has an approximate coverage probability 1 − α.
References Athreya, K. B. and S. N. Lahiri (2006). Measure Theory and Probability Theory. New York: Springer. Coit, D. W. (1997). System-reliability confidence-intervals for complex-systems with estimated component-reliability. IEEE Transactions on Reliability 46, 487–493.
23
Crowder, M. J. (2001). Classical Competing Risks. Florida, Boca Raton: Chapman & Hall/CRC. David, H. A. and M. L. Moeschberger (1978). The Theory of Competing Risks. London: Griffin. Easterling, R. G. (1972). Approximate confidence limits for system reliability. Journal of the American Statistical Association 67, 220–222. Hamada, M. S., A. Wilson, C. S. Reese, and H. Martz (2008). Bayesian Reliability. New York: Springer. Hong, Y. (2011). On computing the distribution function for the sum of independent and non-identically distributed random indicators. Technical Report No. 112: Department of Statitics, Virginia Tech, Blacksburg, Virginia 24061. Available at http://www.stat.vt.edu/tech reports/index.html. Hong, Y. and W. Q. Meeker (2010). Field-failure and warranty prediction based on auxiliary use-rate information. Technometrics 52, 148–159. Hong, Y., W. Q. Meeker, and L. A. Escobar (2008a). Avoiding problems with normal approximation confidence intervals for probabilities. Technometrics 50, 64–68. Hong, Y., W. Q. Meeker, and L. A. Escobar (2008b). The relationship between confidence intervals for failure probabilities and life time quantile. IEEE Transactions on Reliability 57, 260–266. JMP (2011). JMP 9. http://www.jmp.com/. Kuo, W. and M. Zuo (2003). Optimal Reliability Modeling: Principles and Applications. Hoboken, NJ: John Wiley & Sons, Inc. Lehmann, E. L. and G. Casella (1998). Theory of Point Estimation (Second ed.). New York: Springer-Verlag. Meeker, W. Q. and L. A. Escobar (1998). Statistical Methods for Reliability Data. New York: John Wiley & Sons, Inc. Meeker, W. Q. and L. A. Escobar (2008). SPLIDA (S-PLUS Life Data Analysis). www.public.iastate.edu/∼ splida. Minitab (2011). MINITAB Statistical Software, Release 16. State College, Pennsylvania: Minitab Inc. Myhre, J. M. and S. C. Saunders (1968). On confidence limits for the reliability of systems. Annals of Mathematical Statistics 39, 1463–1472.
24
Pawitan, Y. (2001). In All Likelihood: Statistical Modelling and Inference Using Likelihood. USA: Oxford University Press. Ramirez-Marquez, J. E. and W. Jiang (2006). On improved confidence bounds for system reliability. IEEE Transactions on Reliability 55, 26–36. Rausand, M. and A. Høyland (2004). System Reliability Theory: Models, Statistical Methods, and Applications (Second ed.). Hoboken, New Jersey: John Wiley & Sons, Inc. ReliaSoft (2011). Weibull++ 7.5. http://weibull.reliasoft.com/. SAS (2011). SAS/QC User’s Guide, Version 9.3. NC: SAS Institute Inc. Spall, J. C. (2009). System reliability estimation and confidence regions from subsystem and full system tests. In American Control Conference, 2009, St. Louis, MO, pp. 5067– 5072.
25
