- Email: [email protected]
n-fold Weibull multiplicative model
Reliability Engineering and System Safety 74 (2001) 211±219
www.elsevier.com/locate/ress
n-fold Weibull multiplicative model R. Jiang a,*, D.N.P. Murthy b, P. Ji a a
Department of Manufacturing Engineering, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong, China b Department of Mechanical Engineering, The University of Queensland, Q 4072, Australia Received 10 June 2001; accepted 16 August 2001
Abstract In this paper we study the n-fold multiplicative model involving Weibull distributions and examine some properties of the model. These include the shapes for the density and failure rate functions and the WPP plot. These allow one to decide if a given data set can be adequately modelled by the model. We also discuss the estimation of model parameters based on the WPP plot. q 2001 Elsevier Science Ltd. All rights reserved. Keywords: Weibull distribution; Multiplicative model; WPP plot; Density function; Failure rate function; Parameter estimation
1. Introduction The CDF of the n-fold multiplicative model is given by G t Gn t
n Y i1
Fi t:
1
where Fi t; i 1; 2; ¼; n, are probability distribution functions. If Fi (t) are 2-parameter Weibull distributions given by h i 2 Fi t 1 2 exp 2 t=ai bi then Eq. (1) is an n-fold Weibull multiplicative model. We can assume, without loss of generality, that b1 # b2 # ¼ # bn , and ai $ aj for i , j if bi bj . Basu and Klein [1] call the Weibull multiplicative model as the Weibull complementary risk model. When the Fi(t) are identical, the multiplicative model becomes a special case (u n) of the exponentiated Weibull model [2] given by " b !#u t G t 1 2 exp 2 : 3 a An application of the n-fold multiplicative model is as follows. Consider an unreliable system comprising of n components connected in parallel. The failure time (Ti ) for component i is distributed according to Fi t and the failures are statistically independent. The failure time for the system is given by the maximum of T1 ; T2 ; ¼; Tn * Corresponding author. Fax: 1852-2362-5267. E-mail address: [email protected] (R. Jiang).
with distribution function given by Eq. (1). The above relationship can be used to generate random data. A review of the relevant literature is as follows. Jiang and Murthy [3,4] study the 2-fold Weibull multiplicative model. The study includes the parametric characterisation of the shapes of the density and failure rate functions, WPP plot to determine if a given data set can be modelled by a 2-fold Weibull multiplicative model and estimation of model parameters based on the WPP plot. Jiang and Murthy [2] deal with a similar study for the exponentiated Weibull model. Jiang et al. [5] deal with a similar study of the 2-fold inverse Weibull multiplicative model. In this paper, we study the n-fold Weibull multiplicative model. Our study deals with the shapes for the density and failure rate functions and the WPP plot. These allow one to determine the suitability of the model to model a given data set. We examine the problem of estimating the model parameters based on the WPP plot and suggest a method for estimation. The outline of the paper is as follows. Section 2 deals with model analysis. Following this, we study the WPP plot of the model in Section 3. We propose a graphical parameter estimation method in Section 4. The shapes for the density and failure rate functions are examined in Sections 5 and 6, respectively. Finally, we conclude the paper with a brief summary in Section 7. 2. Model analysis The model is mathematically intractable. We use the following three-part approximation for G(t).
0951-8320/01/$ - see front matter q 2001 Elsevier Science Ltd. All rights reserved. PII: S 0951-832 0(01)00108-9
212
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
Approximation 1.
When G(t) is small, then
and
G t < 1 2 e2G t
4
and for small t we have Fi t < zi ; zi t=ai bi :
5
Using these in Eq. (1) yields P bi n Y ti zi Q bi t=a0 b0 z0 < 1 2 exp 2z0 G t < i ai i1 F0 t; 6 where F0 t is the CDF of the Weibull distribution with parameters b0 and a0 are given by
b0
n X i1
bi ; a0
n Y i1
abi i =b0 :
7
Approximation 2.
When G(t) is large, then " # n n X X 1 2 G t < 1 2 1 2 Ri t Ri t; G t i1
G t
i1
1 2 Ri t 1 2
n X i1
Ri t 1
nX 21
n X
i1 ji 1 1
Ri tRj t
and the ®rst summation dominates all the subsequent summations as the Ri ts are all very small. For very large t, we have the following approximation: < kR1 t G t
8
where k is the number of sub-populations with Ri t ; R1 t; i 2; ¼; k: This follows since
i1
(
Ri t R1 t
1
if zi ; z1
0
if zi ± z1
n X i1
exp z1 2 zi ; lim exp z1 2 zi
Example 1. (the well-mixed case): (a 1, a 2, a 3) (1.0, 2.6, 3.6). This yields (a , b ) (3.0076, 2.8168). Example 2. (the normal case): (a 1, a 2, a 3) (2.5, 2.6, 5.5). This yields (a , b ) (4.4427, 2.4010).
i1
n X
Case (a) Well-mixed case: G(t) is large (.0.9) when t tw. Case (b) Well-separated case: G(t) is small (,0.1) when t tw. Case (c) Normal case: G(t) [ (0.1, 0.9) when t tw.
i1
n Y Ri t 1 ¼ 1 21n
< G t
In the above three-part approximation, the ®rst two approximations are for large and small G(t) and the last one is for intermediate range. These approximations will be used later in the analysis of the model. When the sub-populations are not identical, we de®ne the weakest sub-population as the sub-population which ®rst reaches a large value (for example, 0.99) and let this occur at t tw. Let Fw t denote the distribution function for this sub-population. Based on this we have the following three different cases.
Three examples with n 3 and (b 1, b 2, b 3) (0.8, 2.5, 4.3) are as follows.
where Ri t 1 2 Fi t: This follows since n Y
bi zi t0:5 = ezi t0:5 2 1 : b i1 1=n 2 n 2 2 1 ln 1 2 221=n n X
t!1
) :
Approximation 3. For intermediate value of G(t), the model can be approximated by the exponentiated Weibull model given by Eq. (3) with u n. The parameters a and b are determined by the relationships: h i n G t0:5 1 2 exp 2 t0:5 =a b 0:5;
Example 3. (the well-separated case): (a 1, a 2, a 3) (5.0, 2.6, 7.5). This yields (a , b ) (6.2309, 2.2297). In these examples, the weakest sub-population is always sub-population 2 with tw 4.7893. G tw 0.9285 for Example 1, 0.3417 for Example 2 and 0.0830 for Example 3. Note that the weakest sub-population is not necessarily the one with the smallest scale parameter. Fig. 1 shows the plots of G(t) and the three-part approximation. The following observations can be made. 1. The three-part approximation is quite accurate. We measure the accuracy by the absolute average (maximum) errors between the approximation and the true values over the interval (t0.001, t0.999). They are 0.00083 (0.0079) for Example 1, 0.00047 (0.0065) for Example 2 and 0.00048 (0.0167) for Example 3. Moreover, the accuracy is insensitive to n. For example, if we add a new sub-population with the scale and shape parameters 2.0 and 1.5 to the populations of Examples 1±3, we obtain three 4-fold models with G(tw) 0.9057, 0.3333, and 0.0809, respectively. The absolute average (maximum) errors now are 0.00076 (0.0082) for Example 1,
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
213
1 0.9 0.8
Example 1
Example 3
0.7 0.6 0.5
Approximation 1
0.4
Approximation 2 Approximation 3
0.3
G(t)
Example 2
0.2 0.1 0 0
2
4
6
8
10
12
14
16
18
20
t
Fig. 1. G(t) and its approximations.
0.00050 (0.0076) for Example 2, and 0.00052 (0.0156) for Example 3. 2. In the case of the well-mixed case, the distribution can be effectively approximated by Approximation 3 over the whole interval. 3. Approximation 3 can also be used when G(t) is small, but this is not as good as Approximation 1. For our analysis of the model, we approximate the model by these approximations over different intervals of the independent variable t. Since G(t) and its approximations are smooth functions, we assume that their ®rst and second derivatives are close to each other and this has been veri®ed for Examples 1±3. 3. WPP plot Under the Weibull transformations: x ln t; y ln 2ln G t or exp 2ey t ex ; G t
the multiplicative model given by Eq. (1) gets transformed into y y x ln 2ln 1 2 G ex 9 and is called the WPP (Weibull Probability Paper) plot. We look at some of the properties of this WPP plot. Left asymptote. From Eq. (6), as x ! 21 (or t ! 0 the left asymptote is given by: yL
n X i1
yi
with
yi bi x 2 ln ai :
10
This can be rewritten as yL b0 x 2 ln a0
11
with a0 and b0 given by Eq. (7). Right asymptote. From Eq. (8) and noting that uln ku p ln R1 when t is very large, as x ! 1 (or t ! 1) the right asymptote is given by: yR ln 2ln kR1 < ln 2ln R1 y1 :
Shape of WPP plot. For the case n 2 the WPP plot is concave. It is not possible to prove this analytically for general n. Using the three-part approximation, we show that the WPP plot is concave. We need to consider the three cases (Cases (a)±(c)) separately. Case (a) [Well-mixed case]: If the sub-populations are identical or the population is well-mixed, then the model is or can be approximated by Approximation 3 over the whole interval. This implies that the model behaves similar to the exponentiated Weibull model given by Eq. (3) (with u n . 1), for which Jiang and Murthy [2] show that the WPP plot is concave. Case (b) [Well-separated case]: We use induction to prove that the WPP plot is concave. The WPP plot is a straight line for n 1 and is concave for n 2, see Jiang and Murthy [4]. Assume that the WPP plot is concave for n (k 2 1). We now show that it is also true for n k. We need to consider t , tw and t . tw separately.When t , tw, the model can be effectively approximated by Approximation 1 since G(tw) is small. As a result, the WPP plot for x , ln(tw) is similar to the WPP plot for the approximation, which is concave because using Eq. (5) in Eq. (9) and differentiating yields y 00
n X i1
b2i zi zi e 1 2 zi 2 1 , 0 z i e 2 1
214
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
When t . tw, the model can be effectively viewed as an (n 2 1)-fold model by excluding the weakest sub-population since Fw(t) < 1. The WPP plot for x . ln(tw) is concave from the assumption that the WPP plot is concave for n (k 2 1). As a result, the conjecture is true for n k. Case (c) [Normal case]: The proof in this case is similar to that for the well-separated case. The difference is that for t , tw the model can be piecewise approximated by Approximation 1 and Approximation 3, respectively. The WPP plots for these two approximations both are concave. As a result, the WPP plot is concave. Fig. 2 shows the exact WPP plots and the WPP plots obtained using the three-part approximation for Examples 1±3. As can be seen, the ®t between the two is very good and indicates that the WPP plots are concave. 4. Modelling data set Given a data set, one can plot the data on the WPP plot. The plotting procedure depends on the type of data (complete, censored, grouped etc.) and can be found in Ref. [6]. If the WPP plot of the data is roughly concave with asymptote at either end, then an n-fold Weibull multiplicative model can be viewed as potential model for the data set. The WPP plot also provides crude estimates of the model parameters and these are used as starting values for more re®ned statistical estimation methods. We now discuss the estimation of model parameters based on the WPP plot assuming that the model is appropriate. For the 2-fold Weibull multiplicative model, Jiang and Murthy [4] propose a graphical method to estimate the parameters. The basic steps are as follows. From the left asymptote of the WPP plot one obtains estimates of b0 and a 0, and from the right asymptote one obtains estimates of b1 (or b 2) and a 1 (or a 2). Using these in Eq. (7) yields estimates of b 2 (or b1 ) and a 2 (or a 1).
For n . 2, we propose the following three part approach. 4.1. Part 1 [estimating n] The basis for this is that an exponentiated Weibull model (with parameters a; b and u n) has the left asymptote given by YE nbx 2 ln a:
12
On the other hand, the WPP plot of the model has the left asymptote given by Eq. (11). If the selected value of n is appropriate, then the two straight lines should be close to each other. Note that a 0 and b 0 in Eq. (11) can be estimated from the slope and intercept of the ®t to the left asymptote of the WPP plot. Alternatively, for a given n the parameters b and a in Eq. (12) can be obtained from the linear regression equation for Eq. (3) which is given by n n oo ln 2ln 1 2 F t1=n b x 2 ln a: The following quantity measures how close the lines yL and yE are: Z x2 D n yL x 2 yE x 2 dx: 13 x1
Since yL and yE are the left asymptote (for the model) and its approximation (based on the exponentiated Weibull model), respectively, we evaluate D(n) for small values of x1 and x2. Let x1 ln(tmin) 2 0.5 and x2 ln(tmin) 1 0.5, where tmin is the minimum of a given data set. Using Eqs. (11) and (12) in Eq. (13) yields: D n A2 =12 1 Aln tmin 1 B 2 ; 14 where
A b0 2 nb; B ln anb =ab0 0 : D n is small if the selected value of n equals the true value. As such, we can use the following procedure to determine n.
3 2
Example 1
1 0 y
0
0.5
1
1.5
2
2.5
3
3.5
-1 -2
Example 2
y(x) Approximation 1
-3 -4
Approximation 2
Example 3
-5 x
Fig. 2. WPP plots for Examples 1±3.
Approximation 3
4
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
215
2 1 0
y
-1 0
0.5
1
1.5
2
2.5
3
3.5
-2 -3 -4 -5 -6 x Fig. 3. WPP plots for Data-sets 1±3.
Step 1: Fit a tangent line to the left end of the WPP plot of the data and obtain a^ 0 and b^ 0 using Eq. (11). Step 2: For a speci®ed value of n ( 2, 3, 4, ¼) and using data points (ti, G(ti)) such that G(ti) # 0.8 (as Assumption 3 is good only for small and intermediate values of G(t)) obtain the estimates (a^ , b^ ), which are functions of n, from the following regression: Y i b X i 2 ln a where X i ln ti ; Y i ln 2 ln 1 2 G^ ti 1=n : Step 3: Compute the distance function D n given by Eq. (14) for n 2, 3, 4, ¼, using the estimates from Steps 1 and 2. The value of n which yields the minimum for D(n) is the estimate of n. To illustrate this procedure, we look at data generated for the following three models. Models 1 and 3 have n 3 and the model parameters are the same as in Examples 1 and 3 of Section 2. Model 2 has n 4 with the parameters for the three sub-populations being the same as in Example 2 of Section 2 and the scale and shape parameters for the fourth sub-population being 2.0 and 1.5, respectively. For each model, M ( 50) data points were generated using a random generator and these are given below. Data-set I (from Model 1). 1.36, 1.97, 1.98, 2.00, 2.30, 2.63, 2.65, 2.67, 2.68, 2.71, 2.75, 2.77, 2.78, 2.82, 2.86, 2.86, 2.91, 2.98, 2.98, 3.07, 3.10, 3.23, 3.30, 3.31, 3.31, 3.43, 3.66, 3.66, 3.77, 3.78, 3.80, 3.89, 3.99, 4.03, 4.15, 4.21, 4.23, 4.31, 4.33, 4.36,4.50, 4.62, 5.05, 5.33, 7.37
2.02, 2.80, 3.24, 3.89, 4.34,
Data-set 2 (from Model 2). 3.49, 3.67, 3.71, 4.22, 4.30, 4.49, 4.62, 4.68, 4.74, 4.75, 5.27, 5.27, 5.28, 5.34, 5.52,
3.07, 4.46, 5.22, 5.76,
2.81, 2.89, 2.90, 2.92, 4.30, 4.32, 4.33, 4.44, 4.81, 4.96, 4.96, 5.06, 5.52, 5.56, 5.69, 5.70,
5.76, 5.81, 5.97, 6.39, 6.42, 6.42, 6.57, 6.66, 6.67, 6.75, 6.80, 8.02, 8.13, 9.88, 14.87 Data-set 3 (from Model 3). 4.45, 4.79, 5.00, 5.31, 5.45, 5.63, 5.70, 5.72, 5.83, 5.88, 6.06, 6.22, 6.35, 6.45, 6.66, 6.68, 6.94, 6.97, 7.02, 7.15, 7.18, 7.23, 7.26, 7.32, 7.39, 7.57, 7.64, 7.87, 8.00, 8.08, 8.21, 8.40, 8.49, 8.66, 8.70, 8.87, 8.77, 8.89, 9.06, 9.10, 9.11, 9.66, 9.99, 10.76, 13.39, 15.12, 15.44, 16.46, 19.22, 26.18 Fig. 3 shows the WPP plots for the data sets. As can be seen, they are roughly concave. Table 1 shows the computed D(n) (shown in the rows marked (a)) as a function of n for the three data sets. As can be seen, the values of n corresponding to minimum D(n) are 3 for all three models. For Models 1 and 3 this is the same as the true value and for Model 2 it differs by one. 4.2. Part 2 [estimating parameters of sub-populations] ^ a^ 0 and b^ 0 (the outputs from The inputs to Part 2 are n, Part 1). This part deals with the problem of estimating Table 1 Estimating n
a^ 0
b^ 0
n
a^ (n) b^ (n) (a) 3.00 6.24 D(n) (b) 2.524 7.782 D(n) Model 2 a^ (n) b^ (n) (a) 3.97 8.77 D(n) (b) 3.354 10.40 D(n) Model 3 a^ (n) b^ (n) (a) 6.50 10.52 D(n) (b) 6.096 9.954 D(n) Model 1
2
3
4
5
6
3.119 2.748 0.187 0.498 4.831 2.939 0.721 3.511 7.210 3.524 1.347 0.776
2.719 2.175 0.180 0.217 4.249 2.329 0.285 2.070 6.480 2.800 1.062 0.202
2.420 1.875 0.510 0.247 3.812 2.009 0.381 0.838 5.922 2.419 1.550 0.142
2.189 1.686 1.254 0.675 3.471 1.807 1.265 0.159 5.480 2.178 2.965 0.822
2.003 1.553 2.556 1.647 3.195 1.665 3.180 0.324 5.118 2.009 5.487 2.461
216
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
parameters of sub-populations using these inputs and the WPP plots. The basic idea is as follows. We estimate the parameters of one sub-population from the right tangent line of the WPP plot for the data. If the sample size is suf®ciently large, this sub-population would be the sub-population 1. However, it can be any other when the sample size is not very large. Denote this sub-population with subscript `(1)' so as not to confuse it with the subpopulation 1 of the model. Denote the estimated parameters (for the sub-population `(1)') as (a^ 1 ; b^ 1 ). The data associated with this sub-population are removed from the original data set and the remaining data are viewed as coming from an (n 2 1)-fold Weibull multiplicative model. The process is repeated in an iterative manner until the estimates (a^ n21 ; b^ n21 ) have been obtained. Finally, Eq. (7) is used to compute the parameters of the last sub-population and these are given by " #1=b^ n nX 21 nY 21 b^ i b^ 0 ^ ^ ^ b n b 0 2 b i ; a^ n a^ = a^ i1
0
i1
i
, where a^ i and b^ i are the estimates of (i)-th sub-population's parameters obtained at stage i. We now look at how to remove the data associated with the (i)-th sub-population from the (n 2 i 1 1)-fold population. We ®rst introduce the following notations. (i): Gn2i :
The sub-population whose parameters are estimated at stage i, (1 # i # n) Distribution function for the (n 2 i)-fold model excluding the sub-populations `(1)'±`(i)'. This implies that: Gn2i11 t Gn2i tF i t
mn2i :
Number of data points available for estimating the parameters of the (n 2 i)-fold model after the data from the sub-populations `(1)'±`(i)' are deleted
The procedure involves the following two steps: 1. Step 1: Compute the transform data set for the (n 2 i 1 1)-fold population by the following relationship: Gn2i t Gn2i11 t=F^ i t:
15
2. Step 2: Identify the data to be retained (for the (n 2 i)fold population) by ensuring that the data satisfy the following conditions: 1 2 0:3 m 2 0:3 ; n2i , for i [ the set mn2i 1 0:4 mn2i 1 0:4 retained; and ± Gn2i t versus t is roughly monotonic.
± Gn2i ti [
The results of this procedure for the three data sets are given in the column under the heading `Graphical' in Table 2. The
Table 2 Estimates of sub-populations' parameters
Model 1 Model 2
Model 3
(a^ 1 , b^ 1 ) (a^ 2 , b^ 2 ) (a^ 3 , b^ 3 ) (a^ 1 , b^ 1 ) (a^ 2 , b^ 2 ) (a^ 3 , b^ 3 ) (a^ 4 , b^ 4 ) (a^ 1 , b^ 1 ) (a^ 2 , b^ 2 ) (a^ 3 , b^ 3 )
Graphical
MLE
True
(1.23, 0.84) (3.52, 1.99) (3.41, 3.41) (2.11, 0.81) (5.15, 4.01) (3.57, 2.32) (3.34, 1.62) (4.38, 0.86) (6.78, 3.33) (6.70, 6.33)
(1.11, 0.84) (2.06, 3.03) (3.52, 3.91) (2.39, 0.88) (5.21, 4.48) (2.39, 3.59) (2.45, 1.45) (4.85, 0.90) (4.84, 3.88) (7.54, 5.18)
(1.0, 0.8) (2.6, 2.5) (3.6, 4.3) (2.5, 0.8) (5.5, 4.3) (2.6, 2.5) (2.0, 1.5) (5.0, 0.8) (2.6, 2.5) (7.5, 4.3)
values of n associated with the results are the true values. The above estimates can be used as initial values for a more re®ned statistical approach such as the maximum likelihood (ML) method. The ML estimates and the true parameter values for the three data sets are also shown in Table 2. Some observations are as follows. 1. The graphical estimates are close to the estimates obtained from the ML method and hence can be used as good initial values for a numerical iterative process. 2. The proposed graphical estimation method is suitable for all three casesÐwell-mixed, well-separated and normal. 3. The ML estimate values are close to the true values except the estimates for the weakest sub-populations.
4.3. Part 3 [re-estimating n] The value of n has a signi®cant in¯uence for the goodness-of-®t between the model and the data because it determines the number of the model parameters that need to be estimated. From Part 2 we have the estimates of the sub-population parameters. These can be used in Eq. (7) to compute new estimates of a 0 and b 0 (which are generally more accurate than the earlier estimates obtained from the graphical method). As a result, the parameter n can be re-estimated based on the new a^ 0 and b^ 0 (which are functions of the n used in Part 2). We illustrate this by looking at the three data sets again. The re-estimated (a 0, b 0) and the corresponding D(n) are shown in the rows marked `(b)' in Table 1. The new esti^ are 3 for Data Set 1, 5 for Data Set mates of n (denote it by n) 2, and 4 for Data Set 3. As can be seen, the estimate of n for Model 1 is the same as the true value. However, for Models 2 and 3, the estimates differ from the true values by one. One explanation for this phenomenon is that the asymptotic line (yE) of Approximation 3 is not accurately equal to the left asymptote (yL) of the model. When the model is well-mixed, for example Model 1, the proposed method can provide good estimate of n since Approximation 3 is very close to the model. When the model is not well-mixed, for example
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
Model 2 or 3, the proposed method overestimates n. In fact, it can be identi®ed to which case the ®tted model belongs after Part 2. As such, the ®nal estimate for n is n^ if the ®tted model is well-mixed, or (n^ 2 1) if the ®tted model is not well-mixed. 5. Probability density function When n 1, the probability density function (pdf) can have two different shapes: decreasing and unimodal. When n 2, Jiang and Murthy [3] have proved that the pdf can have three different shapes: decreasing, unimodal, and bimodal. For the n-fold Weibull multiplicative model, we show that the possible shapes for the pdf are k-modal (k 1, 2, 3, ¼, n) or decreasing. The proof is in two parts. Part (1): The pdf is either decreasing or k-modal with the upper limit of k not speci®ed. Part (2): The upper limit for k is n. First, note that the pdf for the n-fold model can be written as g t Fn tgn21 t 1 Gn21 tfn t;
16
where Gn21 t and gn21 t are the cdf and pdf of the (n 2 1)fold model excluding the n-th sub-population and that the WPP plot y(x) is concave, i.e. y 0 is decreasing and y 00 , 0: For Part (1), differentiating Eq. (9) with respect to x (noting that y is a function of x, which in turn is a function of t) yields: f t R tey y 0 =t: Differentiating Eq. (17) with respect to t yields: h i f 0 t f t y 02 1 2 ey 2 1=y 0 1 y 00 = ty 0 :
17
2. the shape parameter bn is suf®ciently large. In this case, Fn tgn21 t < 0 and g(t) < Gn21 tfn t is unimodal over the interval (0, tw). On the other hand, g(t)< gn21 t is (n 2 1)-modal over the interval tw ; 1 : As a result, the pdf is n-modal over the interval (0, 1). This proves the result. 6. Failure rate function From Eq. (1), the failure rate function of the n-fold Weibull multiplicative model is given by: r t g t=G t
We need to consider b 0 # 1 and b 0 . 1 separately. 1. When b 0 # 1, we have y 0 , 1. As a result, from Eq. (18) we have f 0 t , 0, implying that the pdf is decreasing. 2. When b 0 . 1, from Eq. (6) we have f 0 t 01 f 00 t 01 . 0. Since f t is decreasing as t ! 1, this implies that the pdf is k-modal (with k $ 1). We now prove Part (2) by induction. Note that it is true for the cases n 1 and 2. Assume that it is true for the (n 2 1)-fold model. We now prove that it is also true for the n-fold model. Consider an (n 2 1)-fold model whose pdf shape is (n 2 1)-modal. We now construct the n-fold model by adding the n-th sub-population to this (n 2 1)-fold model so that: 1. the new n-fold model is well-separated and the n-th subpopulation is the weakest one; and
n G t X f t=Fi t: 1 2 G t i1 i
From Eqs. (6) and (8) we have: 8 b0 21 b0 > > ; as t ! 0 > < r0 t a0 t=a0 r t < : > b1 21 b1 > > r t t= a ; as t ! 1 : 1 1 a1
19
20
This implies that as t ! 0, the failure rate is similar to that for a Weibull distribution with the parameters b0 and a0 , and as t ! 1, it is similar to the failure rate of the subpopulation with the smallest shape parameter. From Eq. (17) we have r t ey y 0 =t
21
and differentiating this with respect to t yields r 0 t
18
217
r t 0 0 00 0 y y 2 1 2 2y : ty
22
We discuss the different possible shapes for r(t) using Eqs. (20) and (22). We need to consider the following three cases. Case (1): b0 # 1. Case (2): b0 . 1 and b1 # 1: Case (3): b0 . 1 and b1 . 1:
Case (1). In this case, r 0 t , 0 since y 0 , 1 and y 00 , 0, implying that the failure rate is decreasing. Case (2). In this case, from Eq. (20) we see that r(t) is increasing for small t and decreasing for large t. This implies that the shape of r(t) is multi-modal. We show that the maximum possible number of peaks for r(t) is (n 2 1). To prove this, we need to study the interrelationship between the functions y 0 y 0 2 1 and 2y 00 : The failure rate r(t) has the maximum possible number of peaks only when G(t) is well-separated and any sub-model formed by different sub-populations is also well- separated. For such a
218
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
case, we have:
Y 00n21 has the following important properties:
0 , tw1 p tw2 p ¼ p tw n21 ;
1. Y 00n2i , 0. This con®rms that the WPP plot is concave for the well-separated case under consideration; and 2. 2 Y 00n2i is a unimodal curve with the mode occurring at location of the modal at tm;n2i given by zn2i tm;n2i 1:8605, which corresponds to Fn2i tm;n2i 0:8444:
where twi (i 1, 2, ¼, n 2 1) is the time instant when the ith weakest sub-population's cdf reaches a certain large value. This implies G t < Gn2i t; t [ twi ; tw i11 ; i 0; 1; 2; ¼; n 2 1;
23
where Gn2i(t) is the cdf of the (n 2 i)-fold model excluding the ®rst i weakest sub-populations, tw0 0 and twn 1. Without loss of generality, we can assume that the subpopulations are numbered such that they form a sequence where the j-th weakest sub-population is Fn2j11 (t). Since Gn2i(t) is well-separated, we can use Approximation 1 and we have Yn2i ln 2ln 1 2 Gn2i t < ln Gn2i t ; 24
For the well-separated case under consideration, we have twi , tm;n2i , tw(i11). This implies that 2y 00 x is (n 2 1)modal over (0, 1). On the other hand, y 0 y 0 2 1 is decreasing since y 0 (x) is decreasing. Fig. 4 illustrates these for a well-separated 3-fold model with model parameters given by Example 3. From Eq. (22) the zeros of r 0 (t) correspond to the intersections of the curves y 0 (y 0 2 1) and 2y 00 : An extreme case is that y 0 (y 0 2 1) intersects with all the peaks of 2y 00 : Note that:
where Yn2i is the WPP plot of Gn2i(t), and it approximates y(x) over the interval x [ (ln(twi), ln(tw(i11))). For subpopulations 1 through (n 2 i 2 1), Fj t < zj . As a result, we have: 25 Yn2i < b 00 x 2 ln a 00 1 ln 1 2 exp 2zn2i ;
1. there are two intersections if y 0 (y 0 2 1) intersects with any one of the ®rst (n 2 2) peaks of 2y 00 ; and 2. only one intersection if it intersects with the (n 2 1)-th peak, since y 0 (y 0 2 1) tends to b 1(b 1 2 1) as t ! 1, which is zero or a negative constant when b1 # 1.
where
b 00
n 2X i21 j1
bj ; a 00
n2 i21 Y j1
b =b 00
aj j
; zn2i t=an2i
bn2i
As such, there are totally 2 £ (n 2 2) 1 1 intersections. As a result, r(t) is (n 2 1)-modal. Jiang and Murthy [3] show that the failure rate for 2-fold model can be unimodal and this agrees with the above conclusion.
: 26
Differentiating Eq. (25) with respect to x yields y 0 x < Y 0n2i < b 00 1 and 00
y x <
Y 00n2i
bn2i zn2i exp zn2i 2 1
27
Case (3). In this case, y 0 y 0 2 1) is decreasing and tends to a positive constant as t ! 1. On the other hand, 2y 00 is unimodal, or bimodal, ¼, or (n 2 1)-modal and tends to zero as t ! 1. As a result, y 0 (y 0 2 1) does not intersect with 2y 00 or intersects with 2 £ 1, or 2 £ 2, ¼, or 2 £ (n 2 1) intersections, implying that r(t) is either mono-
" # zn2i exp zn2i b2n2i zn2i 12 : 28 < exp zn2i 2 1 exp zn2i 2 1 14
y’(y’-1)
12 10 8
y’
6 4 -y"
2 0 -1
-0.5
0
0.5
1
1.5
x=ln(t) Fig. 4. Curves y 0 ; y 00 ; and r 0 t.
2
2.5
3
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
tonically increasing or k-modal-followed-by-increasing with k 1, 2, ¼, n 2 1. For n 2, Jiang and Murthy [3] show that the failure rate can be either monotonically increasing or unimodalfollowed-by-increasing and this agrees with the above conclusion. In summary, the shape types of the failure rate can be grouped into four groups: ² ² ² ²
monotonically decreasing; monotonically increasing; k-modal-followed-by-increasing, k 1, 2, ¼, n 2 1; and k-modal, k 1, 2, ¼, n 2 1.
Finally, it is worth noting that failure rate can never be decreasing for small t and increasing for large t. This rules out the possibility for the failure rate to have a bathtub shape. 7. Conclusions In this paper, we have carried out a detailed characterisation of the n-fold Weibull multiplicative model in terms of the shapes for the WPP plot, density and failure rate functions. This understanding is useful for determining whether or not a given set of data can be modelled by an n-fold
219
Weibull multiplicative model. We have also looked at a graphical method to estimate the model parameters.
Acknowledgements The ®rst and third authors wish to acknowledge the Hong Kong Polytechnic University for the ®nancial support of the project (No. G-YW41).
References [1] A.P. Basu, J.P. Klein. Some recent results in competing risks theory, in: Survival analysis, Vol. 2, Proceedings of the special topics meeting sponsored by the Institute of Mathematical Statistics. Columbus, Ohio, October 26±28, 1981, pp. 216±229. [2] Jiang R, Murthy DNP. The exponentiated Weibull family: a graphical approach. IEEE Transactions on Reliability 1999;48(1):68±72. [3] Jiang R, Murthy DNP. Parametric study of multiplicative model involving two Weibull distributions. Reliability Engineering and System Safety 1997;55:217±26. [4] Jiang R, Murthy DNP. Reliability modelling involving two Weibull distributions. Reliability Engineering and System Safety 1995;47:187± 98. [5] Jiang R, Murthy DNP, Ji P. Models involving two inverse Weibull distributions. Reliability Engineering and System Safety 2001;73(1):73±81. [6] Nelson W. Applied Life Data Analysis. New York: John Wiley, 1982.
www.elsevier.com/locate/ress
n-fold Weibull multiplicative model R. Jiang a,*, D.N.P. Murthy b, P. Ji a a
Department of Manufacturing Engineering, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong, China b Department of Mechanical Engineering, The University of Queensland, Q 4072, Australia Received 10 June 2001; accepted 16 August 2001
Abstract In this paper we study the n-fold multiplicative model involving Weibull distributions and examine some properties of the model. These include the shapes for the density and failure rate functions and the WPP plot. These allow one to decide if a given data set can be adequately modelled by the model. We also discuss the estimation of model parameters based on the WPP plot. q 2001 Elsevier Science Ltd. All rights reserved. Keywords: Weibull distribution; Multiplicative model; WPP plot; Density function; Failure rate function; Parameter estimation
1. Introduction The CDF of the n-fold multiplicative model is given by G t Gn t
n Y i1
Fi t:
1
where Fi t; i 1; 2; ¼; n, are probability distribution functions. If Fi (t) are 2-parameter Weibull distributions given by h i 2 Fi t 1 2 exp 2 t=ai bi then Eq. (1) is an n-fold Weibull multiplicative model. We can assume, without loss of generality, that b1 # b2 # ¼ # bn , and ai $ aj for i , j if bi bj . Basu and Klein [1] call the Weibull multiplicative model as the Weibull complementary risk model. When the Fi(t) are identical, the multiplicative model becomes a special case (u n) of the exponentiated Weibull model [2] given by " b !#u t G t 1 2 exp 2 : 3 a An application of the n-fold multiplicative model is as follows. Consider an unreliable system comprising of n components connected in parallel. The failure time (Ti ) for component i is distributed according to Fi t and the failures are statistically independent. The failure time for the system is given by the maximum of T1 ; T2 ; ¼; Tn * Corresponding author. Fax: 1852-2362-5267. E-mail address: [email protected] (R. Jiang).
with distribution function given by Eq. (1). The above relationship can be used to generate random data. A review of the relevant literature is as follows. Jiang and Murthy [3,4] study the 2-fold Weibull multiplicative model. The study includes the parametric characterisation of the shapes of the density and failure rate functions, WPP plot to determine if a given data set can be modelled by a 2-fold Weibull multiplicative model and estimation of model parameters based on the WPP plot. Jiang and Murthy [2] deal with a similar study for the exponentiated Weibull model. Jiang et al. [5] deal with a similar study of the 2-fold inverse Weibull multiplicative model. In this paper, we study the n-fold Weibull multiplicative model. Our study deals with the shapes for the density and failure rate functions and the WPP plot. These allow one to determine the suitability of the model to model a given data set. We examine the problem of estimating the model parameters based on the WPP plot and suggest a method for estimation. The outline of the paper is as follows. Section 2 deals with model analysis. Following this, we study the WPP plot of the model in Section 3. We propose a graphical parameter estimation method in Section 4. The shapes for the density and failure rate functions are examined in Sections 5 and 6, respectively. Finally, we conclude the paper with a brief summary in Section 7. 2. Model analysis The model is mathematically intractable. We use the following three-part approximation for G(t).
0951-8320/01/$ - see front matter q 2001 Elsevier Science Ltd. All rights reserved. PII: S 0951-832 0(01)00108-9
212
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
Approximation 1.
When G(t) is small, then
and
G t < 1 2 e2G t
4
and for small t we have Fi t < zi ; zi t=ai bi :
5
Using these in Eq. (1) yields P bi n Y ti zi Q bi t=a0 b0 z0 < 1 2 exp 2z0 G t < i ai i1 F0 t; 6 where F0 t is the CDF of the Weibull distribution with parameters b0 and a0 are given by
b0
n X i1
bi ; a0
n Y i1
abi i =b0 :
7
Approximation 2.
When G(t) is large, then " # n n X X 1 2 G t < 1 2 1 2 Ri t Ri t; G t i1
G t
i1
1 2 Ri t 1 2
n X i1
Ri t 1
nX 21
n X
i1 ji 1 1
Ri tRj t
and the ®rst summation dominates all the subsequent summations as the Ri ts are all very small. For very large t, we have the following approximation: < kR1 t G t
8
where k is the number of sub-populations with Ri t ; R1 t; i 2; ¼; k: This follows since
i1
(
Ri t R1 t
1
if zi ; z1
0
if zi ± z1
n X i1
exp z1 2 zi ; lim exp z1 2 zi
Example 1. (the well-mixed case): (a 1, a 2, a 3) (1.0, 2.6, 3.6). This yields (a , b ) (3.0076, 2.8168). Example 2. (the normal case): (a 1, a 2, a 3) (2.5, 2.6, 5.5). This yields (a , b ) (4.4427, 2.4010).
i1
n X
Case (a) Well-mixed case: G(t) is large (.0.9) when t tw. Case (b) Well-separated case: G(t) is small (,0.1) when t tw. Case (c) Normal case: G(t) [ (0.1, 0.9) when t tw.
i1
n Y Ri t 1 ¼ 1 21n
< G t
In the above three-part approximation, the ®rst two approximations are for large and small G(t) and the last one is for intermediate range. These approximations will be used later in the analysis of the model. When the sub-populations are not identical, we de®ne the weakest sub-population as the sub-population which ®rst reaches a large value (for example, 0.99) and let this occur at t tw. Let Fw t denote the distribution function for this sub-population. Based on this we have the following three different cases.
Three examples with n 3 and (b 1, b 2, b 3) (0.8, 2.5, 4.3) are as follows.
where Ri t 1 2 Fi t: This follows since n Y
bi zi t0:5 = ezi t0:5 2 1 : b i1 1=n 2 n 2 2 1 ln 1 2 221=n n X
t!1
) :
Approximation 3. For intermediate value of G(t), the model can be approximated by the exponentiated Weibull model given by Eq. (3) with u n. The parameters a and b are determined by the relationships: h i n G t0:5 1 2 exp 2 t0:5 =a b 0:5;
Example 3. (the well-separated case): (a 1, a 2, a 3) (5.0, 2.6, 7.5). This yields (a , b ) (6.2309, 2.2297). In these examples, the weakest sub-population is always sub-population 2 with tw 4.7893. G tw 0.9285 for Example 1, 0.3417 for Example 2 and 0.0830 for Example 3. Note that the weakest sub-population is not necessarily the one with the smallest scale parameter. Fig. 1 shows the plots of G(t) and the three-part approximation. The following observations can be made. 1. The three-part approximation is quite accurate. We measure the accuracy by the absolute average (maximum) errors between the approximation and the true values over the interval (t0.001, t0.999). They are 0.00083 (0.0079) for Example 1, 0.00047 (0.0065) for Example 2 and 0.00048 (0.0167) for Example 3. Moreover, the accuracy is insensitive to n. For example, if we add a new sub-population with the scale and shape parameters 2.0 and 1.5 to the populations of Examples 1±3, we obtain three 4-fold models with G(tw) 0.9057, 0.3333, and 0.0809, respectively. The absolute average (maximum) errors now are 0.00076 (0.0082) for Example 1,
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
213
1 0.9 0.8
Example 1
Example 3
0.7 0.6 0.5
Approximation 1
0.4
Approximation 2 Approximation 3
0.3
G(t)
Example 2
0.2 0.1 0 0
2
4
6
8
10
12
14
16
18
20
t
Fig. 1. G(t) and its approximations.
0.00050 (0.0076) for Example 2, and 0.00052 (0.0156) for Example 3. 2. In the case of the well-mixed case, the distribution can be effectively approximated by Approximation 3 over the whole interval. 3. Approximation 3 can also be used when G(t) is small, but this is not as good as Approximation 1. For our analysis of the model, we approximate the model by these approximations over different intervals of the independent variable t. Since G(t) and its approximations are smooth functions, we assume that their ®rst and second derivatives are close to each other and this has been veri®ed for Examples 1±3. 3. WPP plot Under the Weibull transformations: x ln t; y ln 2ln G t or exp 2ey t ex ; G t
the multiplicative model given by Eq. (1) gets transformed into y y x ln 2ln 1 2 G ex 9 and is called the WPP (Weibull Probability Paper) plot. We look at some of the properties of this WPP plot. Left asymptote. From Eq. (6), as x ! 21 (or t ! 0 the left asymptote is given by: yL
n X i1
yi
with
yi bi x 2 ln ai :
10
This can be rewritten as yL b0 x 2 ln a0
11
with a0 and b0 given by Eq. (7). Right asymptote. From Eq. (8) and noting that uln ku p ln R1 when t is very large, as x ! 1 (or t ! 1) the right asymptote is given by: yR ln 2ln kR1 < ln 2ln R1 y1 :
Shape of WPP plot. For the case n 2 the WPP plot is concave. It is not possible to prove this analytically for general n. Using the three-part approximation, we show that the WPP plot is concave. We need to consider the three cases (Cases (a)±(c)) separately. Case (a) [Well-mixed case]: If the sub-populations are identical or the population is well-mixed, then the model is or can be approximated by Approximation 3 over the whole interval. This implies that the model behaves similar to the exponentiated Weibull model given by Eq. (3) (with u n . 1), for which Jiang and Murthy [2] show that the WPP plot is concave. Case (b) [Well-separated case]: We use induction to prove that the WPP plot is concave. The WPP plot is a straight line for n 1 and is concave for n 2, see Jiang and Murthy [4]. Assume that the WPP plot is concave for n (k 2 1). We now show that it is also true for n k. We need to consider t , tw and t . tw separately.When t , tw, the model can be effectively approximated by Approximation 1 since G(tw) is small. As a result, the WPP plot for x , ln(tw) is similar to the WPP plot for the approximation, which is concave because using Eq. (5) in Eq. (9) and differentiating yields y 00
n X i1
b2i zi zi e 1 2 zi 2 1 , 0 z i e 2 1
214
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
When t . tw, the model can be effectively viewed as an (n 2 1)-fold model by excluding the weakest sub-population since Fw(t) < 1. The WPP plot for x . ln(tw) is concave from the assumption that the WPP plot is concave for n (k 2 1). As a result, the conjecture is true for n k. Case (c) [Normal case]: The proof in this case is similar to that for the well-separated case. The difference is that for t , tw the model can be piecewise approximated by Approximation 1 and Approximation 3, respectively. The WPP plots for these two approximations both are concave. As a result, the WPP plot is concave. Fig. 2 shows the exact WPP plots and the WPP plots obtained using the three-part approximation for Examples 1±3. As can be seen, the ®t between the two is very good and indicates that the WPP plots are concave. 4. Modelling data set Given a data set, one can plot the data on the WPP plot. The plotting procedure depends on the type of data (complete, censored, grouped etc.) and can be found in Ref. [6]. If the WPP plot of the data is roughly concave with asymptote at either end, then an n-fold Weibull multiplicative model can be viewed as potential model for the data set. The WPP plot also provides crude estimates of the model parameters and these are used as starting values for more re®ned statistical estimation methods. We now discuss the estimation of model parameters based on the WPP plot assuming that the model is appropriate. For the 2-fold Weibull multiplicative model, Jiang and Murthy [4] propose a graphical method to estimate the parameters. The basic steps are as follows. From the left asymptote of the WPP plot one obtains estimates of b0 and a 0, and from the right asymptote one obtains estimates of b1 (or b 2) and a 1 (or a 2). Using these in Eq. (7) yields estimates of b 2 (or b1 ) and a 2 (or a 1).
For n . 2, we propose the following three part approach. 4.1. Part 1 [estimating n] The basis for this is that an exponentiated Weibull model (with parameters a; b and u n) has the left asymptote given by YE nbx 2 ln a:
12
On the other hand, the WPP plot of the model has the left asymptote given by Eq. (11). If the selected value of n is appropriate, then the two straight lines should be close to each other. Note that a 0 and b 0 in Eq. (11) can be estimated from the slope and intercept of the ®t to the left asymptote of the WPP plot. Alternatively, for a given n the parameters b and a in Eq. (12) can be obtained from the linear regression equation for Eq. (3) which is given by n n oo ln 2ln 1 2 F t1=n b x 2 ln a: The following quantity measures how close the lines yL and yE are: Z x2 D n yL x 2 yE x 2 dx: 13 x1
Since yL and yE are the left asymptote (for the model) and its approximation (based on the exponentiated Weibull model), respectively, we evaluate D(n) for small values of x1 and x2. Let x1 ln(tmin) 2 0.5 and x2 ln(tmin) 1 0.5, where tmin is the minimum of a given data set. Using Eqs. (11) and (12) in Eq. (13) yields: D n A2 =12 1 Aln tmin 1 B 2 ; 14 where
A b0 2 nb; B ln anb =ab0 0 : D n is small if the selected value of n equals the true value. As such, we can use the following procedure to determine n.
3 2
Example 1
1 0 y
0
0.5
1
1.5
2
2.5
3
3.5
-1 -2
Example 2
y(x) Approximation 1
-3 -4
Approximation 2
Example 3
-5 x
Fig. 2. WPP plots for Examples 1±3.
Approximation 3
4
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
215
2 1 0
y
-1 0
0.5
1
1.5
2
2.5
3
3.5
-2 -3 -4 -5 -6 x Fig. 3. WPP plots for Data-sets 1±3.
Step 1: Fit a tangent line to the left end of the WPP plot of the data and obtain a^ 0 and b^ 0 using Eq. (11). Step 2: For a speci®ed value of n ( 2, 3, 4, ¼) and using data points (ti, G(ti)) such that G(ti) # 0.8 (as Assumption 3 is good only for small and intermediate values of G(t)) obtain the estimates (a^ , b^ ), which are functions of n, from the following regression: Y i b X i 2 ln a where X i ln ti ; Y i ln 2 ln 1 2 G^ ti 1=n : Step 3: Compute the distance function D n given by Eq. (14) for n 2, 3, 4, ¼, using the estimates from Steps 1 and 2. The value of n which yields the minimum for D(n) is the estimate of n. To illustrate this procedure, we look at data generated for the following three models. Models 1 and 3 have n 3 and the model parameters are the same as in Examples 1 and 3 of Section 2. Model 2 has n 4 with the parameters for the three sub-populations being the same as in Example 2 of Section 2 and the scale and shape parameters for the fourth sub-population being 2.0 and 1.5, respectively. For each model, M ( 50) data points were generated using a random generator and these are given below. Data-set I (from Model 1). 1.36, 1.97, 1.98, 2.00, 2.30, 2.63, 2.65, 2.67, 2.68, 2.71, 2.75, 2.77, 2.78, 2.82, 2.86, 2.86, 2.91, 2.98, 2.98, 3.07, 3.10, 3.23, 3.30, 3.31, 3.31, 3.43, 3.66, 3.66, 3.77, 3.78, 3.80, 3.89, 3.99, 4.03, 4.15, 4.21, 4.23, 4.31, 4.33, 4.36,4.50, 4.62, 5.05, 5.33, 7.37
2.02, 2.80, 3.24, 3.89, 4.34,
Data-set 2 (from Model 2). 3.49, 3.67, 3.71, 4.22, 4.30, 4.49, 4.62, 4.68, 4.74, 4.75, 5.27, 5.27, 5.28, 5.34, 5.52,
3.07, 4.46, 5.22, 5.76,
2.81, 2.89, 2.90, 2.92, 4.30, 4.32, 4.33, 4.44, 4.81, 4.96, 4.96, 5.06, 5.52, 5.56, 5.69, 5.70,
5.76, 5.81, 5.97, 6.39, 6.42, 6.42, 6.57, 6.66, 6.67, 6.75, 6.80, 8.02, 8.13, 9.88, 14.87 Data-set 3 (from Model 3). 4.45, 4.79, 5.00, 5.31, 5.45, 5.63, 5.70, 5.72, 5.83, 5.88, 6.06, 6.22, 6.35, 6.45, 6.66, 6.68, 6.94, 6.97, 7.02, 7.15, 7.18, 7.23, 7.26, 7.32, 7.39, 7.57, 7.64, 7.87, 8.00, 8.08, 8.21, 8.40, 8.49, 8.66, 8.70, 8.87, 8.77, 8.89, 9.06, 9.10, 9.11, 9.66, 9.99, 10.76, 13.39, 15.12, 15.44, 16.46, 19.22, 26.18 Fig. 3 shows the WPP plots for the data sets. As can be seen, they are roughly concave. Table 1 shows the computed D(n) (shown in the rows marked (a)) as a function of n for the three data sets. As can be seen, the values of n corresponding to minimum D(n) are 3 for all three models. For Models 1 and 3 this is the same as the true value and for Model 2 it differs by one. 4.2. Part 2 [estimating parameters of sub-populations] ^ a^ 0 and b^ 0 (the outputs from The inputs to Part 2 are n, Part 1). This part deals with the problem of estimating Table 1 Estimating n
a^ 0
b^ 0
n
a^ (n) b^ (n) (a) 3.00 6.24 D(n) (b) 2.524 7.782 D(n) Model 2 a^ (n) b^ (n) (a) 3.97 8.77 D(n) (b) 3.354 10.40 D(n) Model 3 a^ (n) b^ (n) (a) 6.50 10.52 D(n) (b) 6.096 9.954 D(n) Model 1
2
3
4
5
6
3.119 2.748 0.187 0.498 4.831 2.939 0.721 3.511 7.210 3.524 1.347 0.776
2.719 2.175 0.180 0.217 4.249 2.329 0.285 2.070 6.480 2.800 1.062 0.202
2.420 1.875 0.510 0.247 3.812 2.009 0.381 0.838 5.922 2.419 1.550 0.142
2.189 1.686 1.254 0.675 3.471 1.807 1.265 0.159 5.480 2.178 2.965 0.822
2.003 1.553 2.556 1.647 3.195 1.665 3.180 0.324 5.118 2.009 5.487 2.461
216
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
parameters of sub-populations using these inputs and the WPP plots. The basic idea is as follows. We estimate the parameters of one sub-population from the right tangent line of the WPP plot for the data. If the sample size is suf®ciently large, this sub-population would be the sub-population 1. However, it can be any other when the sample size is not very large. Denote this sub-population with subscript `(1)' so as not to confuse it with the subpopulation 1 of the model. Denote the estimated parameters (for the sub-population `(1)') as (a^ 1 ; b^ 1 ). The data associated with this sub-population are removed from the original data set and the remaining data are viewed as coming from an (n 2 1)-fold Weibull multiplicative model. The process is repeated in an iterative manner until the estimates (a^ n21 ; b^ n21 ) have been obtained. Finally, Eq. (7) is used to compute the parameters of the last sub-population and these are given by " #1=b^ n nX 21 nY 21 b^ i b^ 0 ^ ^ ^ b n b 0 2 b i ; a^ n a^ = a^ i1
0
i1
i
, where a^ i and b^ i are the estimates of (i)-th sub-population's parameters obtained at stage i. We now look at how to remove the data associated with the (i)-th sub-population from the (n 2 i 1 1)-fold population. We ®rst introduce the following notations. (i): Gn2i :
The sub-population whose parameters are estimated at stage i, (1 # i # n) Distribution function for the (n 2 i)-fold model excluding the sub-populations `(1)'±`(i)'. This implies that: Gn2i11 t Gn2i tF i t
mn2i :
Number of data points available for estimating the parameters of the (n 2 i)-fold model after the data from the sub-populations `(1)'±`(i)' are deleted
The procedure involves the following two steps: 1. Step 1: Compute the transform data set for the (n 2 i 1 1)-fold population by the following relationship: Gn2i t Gn2i11 t=F^ i t:
15
2. Step 2: Identify the data to be retained (for the (n 2 i)fold population) by ensuring that the data satisfy the following conditions: 1 2 0:3 m 2 0:3 ; n2i , for i [ the set mn2i 1 0:4 mn2i 1 0:4 retained; and ± Gn2i t versus t is roughly monotonic.
± Gn2i ti [
The results of this procedure for the three data sets are given in the column under the heading `Graphical' in Table 2. The
Table 2 Estimates of sub-populations' parameters
Model 1 Model 2
Model 3
(a^ 1 , b^ 1 ) (a^ 2 , b^ 2 ) (a^ 3 , b^ 3 ) (a^ 1 , b^ 1 ) (a^ 2 , b^ 2 ) (a^ 3 , b^ 3 ) (a^ 4 , b^ 4 ) (a^ 1 , b^ 1 ) (a^ 2 , b^ 2 ) (a^ 3 , b^ 3 )
Graphical
MLE
True
(1.23, 0.84) (3.52, 1.99) (3.41, 3.41) (2.11, 0.81) (5.15, 4.01) (3.57, 2.32) (3.34, 1.62) (4.38, 0.86) (6.78, 3.33) (6.70, 6.33)
(1.11, 0.84) (2.06, 3.03) (3.52, 3.91) (2.39, 0.88) (5.21, 4.48) (2.39, 3.59) (2.45, 1.45) (4.85, 0.90) (4.84, 3.88) (7.54, 5.18)
(1.0, 0.8) (2.6, 2.5) (3.6, 4.3) (2.5, 0.8) (5.5, 4.3) (2.6, 2.5) (2.0, 1.5) (5.0, 0.8) (2.6, 2.5) (7.5, 4.3)
values of n associated with the results are the true values. The above estimates can be used as initial values for a more re®ned statistical approach such as the maximum likelihood (ML) method. The ML estimates and the true parameter values for the three data sets are also shown in Table 2. Some observations are as follows. 1. The graphical estimates are close to the estimates obtained from the ML method and hence can be used as good initial values for a numerical iterative process. 2. The proposed graphical estimation method is suitable for all three casesÐwell-mixed, well-separated and normal. 3. The ML estimate values are close to the true values except the estimates for the weakest sub-populations.
4.3. Part 3 [re-estimating n] The value of n has a signi®cant in¯uence for the goodness-of-®t between the model and the data because it determines the number of the model parameters that need to be estimated. From Part 2 we have the estimates of the sub-population parameters. These can be used in Eq. (7) to compute new estimates of a 0 and b 0 (which are generally more accurate than the earlier estimates obtained from the graphical method). As a result, the parameter n can be re-estimated based on the new a^ 0 and b^ 0 (which are functions of the n used in Part 2). We illustrate this by looking at the three data sets again. The re-estimated (a 0, b 0) and the corresponding D(n) are shown in the rows marked `(b)' in Table 1. The new esti^ are 3 for Data Set 1, 5 for Data Set mates of n (denote it by n) 2, and 4 for Data Set 3. As can be seen, the estimate of n for Model 1 is the same as the true value. However, for Models 2 and 3, the estimates differ from the true values by one. One explanation for this phenomenon is that the asymptotic line (yE) of Approximation 3 is not accurately equal to the left asymptote (yL) of the model. When the model is well-mixed, for example Model 1, the proposed method can provide good estimate of n since Approximation 3 is very close to the model. When the model is not well-mixed, for example
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
Model 2 or 3, the proposed method overestimates n. In fact, it can be identi®ed to which case the ®tted model belongs after Part 2. As such, the ®nal estimate for n is n^ if the ®tted model is well-mixed, or (n^ 2 1) if the ®tted model is not well-mixed. 5. Probability density function When n 1, the probability density function (pdf) can have two different shapes: decreasing and unimodal. When n 2, Jiang and Murthy [3] have proved that the pdf can have three different shapes: decreasing, unimodal, and bimodal. For the n-fold Weibull multiplicative model, we show that the possible shapes for the pdf are k-modal (k 1, 2, 3, ¼, n) or decreasing. The proof is in two parts. Part (1): The pdf is either decreasing or k-modal with the upper limit of k not speci®ed. Part (2): The upper limit for k is n. First, note that the pdf for the n-fold model can be written as g t Fn tgn21 t 1 Gn21 tfn t;
16
where Gn21 t and gn21 t are the cdf and pdf of the (n 2 1)fold model excluding the n-th sub-population and that the WPP plot y(x) is concave, i.e. y 0 is decreasing and y 00 , 0: For Part (1), differentiating Eq. (9) with respect to x (noting that y is a function of x, which in turn is a function of t) yields: f t R tey y 0 =t: Differentiating Eq. (17) with respect to t yields: h i f 0 t f t y 02 1 2 ey 2 1=y 0 1 y 00 = ty 0 :
17
2. the shape parameter bn is suf®ciently large. In this case, Fn tgn21 t < 0 and g(t) < Gn21 tfn t is unimodal over the interval (0, tw). On the other hand, g(t)< gn21 t is (n 2 1)-modal over the interval tw ; 1 : As a result, the pdf is n-modal over the interval (0, 1). This proves the result. 6. Failure rate function From Eq. (1), the failure rate function of the n-fold Weibull multiplicative model is given by: r t g t=G t
We need to consider b 0 # 1 and b 0 . 1 separately. 1. When b 0 # 1, we have y 0 , 1. As a result, from Eq. (18) we have f 0 t , 0, implying that the pdf is decreasing. 2. When b 0 . 1, from Eq. (6) we have f 0 t 01 f 00 t 01 . 0. Since f t is decreasing as t ! 1, this implies that the pdf is k-modal (with k $ 1). We now prove Part (2) by induction. Note that it is true for the cases n 1 and 2. Assume that it is true for the (n 2 1)-fold model. We now prove that it is also true for the n-fold model. Consider an (n 2 1)-fold model whose pdf shape is (n 2 1)-modal. We now construct the n-fold model by adding the n-th sub-population to this (n 2 1)-fold model so that: 1. the new n-fold model is well-separated and the n-th subpopulation is the weakest one; and
n G t X f t=Fi t: 1 2 G t i1 i
From Eqs. (6) and (8) we have: 8 b0 21 b0 > > ; as t ! 0 > < r0 t a0 t=a0 r t < : > b1 21 b1 > > r t t= a ; as t ! 1 : 1 1 a1
19
20
This implies that as t ! 0, the failure rate is similar to that for a Weibull distribution with the parameters b0 and a0 , and as t ! 1, it is similar to the failure rate of the subpopulation with the smallest shape parameter. From Eq. (17) we have r t ey y 0 =t
21
and differentiating this with respect to t yields r 0 t
18
217
r t 0 0 00 0 y y 2 1 2 2y : ty
22
We discuss the different possible shapes for r(t) using Eqs. (20) and (22). We need to consider the following three cases. Case (1): b0 # 1. Case (2): b0 . 1 and b1 # 1: Case (3): b0 . 1 and b1 . 1:
Case (1). In this case, r 0 t , 0 since y 0 , 1 and y 00 , 0, implying that the failure rate is decreasing. Case (2). In this case, from Eq. (20) we see that r(t) is increasing for small t and decreasing for large t. This implies that the shape of r(t) is multi-modal. We show that the maximum possible number of peaks for r(t) is (n 2 1). To prove this, we need to study the interrelationship between the functions y 0 y 0 2 1 and 2y 00 : The failure rate r(t) has the maximum possible number of peaks only when G(t) is well-separated and any sub-model formed by different sub-populations is also well- separated. For such a
218
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
case, we have:
Y 00n21 has the following important properties:
0 , tw1 p tw2 p ¼ p tw n21 ;
1. Y 00n2i , 0. This con®rms that the WPP plot is concave for the well-separated case under consideration; and 2. 2 Y 00n2i is a unimodal curve with the mode occurring at location of the modal at tm;n2i given by zn2i tm;n2i 1:8605, which corresponds to Fn2i tm;n2i 0:8444:
where twi (i 1, 2, ¼, n 2 1) is the time instant when the ith weakest sub-population's cdf reaches a certain large value. This implies G t < Gn2i t; t [ twi ; tw i11 ; i 0; 1; 2; ¼; n 2 1;
23
where Gn2i(t) is the cdf of the (n 2 i)-fold model excluding the ®rst i weakest sub-populations, tw0 0 and twn 1. Without loss of generality, we can assume that the subpopulations are numbered such that they form a sequence where the j-th weakest sub-population is Fn2j11 (t). Since Gn2i(t) is well-separated, we can use Approximation 1 and we have Yn2i ln 2ln 1 2 Gn2i t < ln Gn2i t ; 24
For the well-separated case under consideration, we have twi , tm;n2i , tw(i11). This implies that 2y 00 x is (n 2 1)modal over (0, 1). On the other hand, y 0 y 0 2 1 is decreasing since y 0 (x) is decreasing. Fig. 4 illustrates these for a well-separated 3-fold model with model parameters given by Example 3. From Eq. (22) the zeros of r 0 (t) correspond to the intersections of the curves y 0 (y 0 2 1) and 2y 00 : An extreme case is that y 0 (y 0 2 1) intersects with all the peaks of 2y 00 : Note that:
where Yn2i is the WPP plot of Gn2i(t), and it approximates y(x) over the interval x [ (ln(twi), ln(tw(i11))). For subpopulations 1 through (n 2 i 2 1), Fj t < zj . As a result, we have: 25 Yn2i < b 00 x 2 ln a 00 1 ln 1 2 exp 2zn2i ;
1. there are two intersections if y 0 (y 0 2 1) intersects with any one of the ®rst (n 2 2) peaks of 2y 00 ; and 2. only one intersection if it intersects with the (n 2 1)-th peak, since y 0 (y 0 2 1) tends to b 1(b 1 2 1) as t ! 1, which is zero or a negative constant when b1 # 1.
where
b 00
n 2X i21 j1
bj ; a 00
n2 i21 Y j1
b =b 00
aj j
; zn2i t=an2i
bn2i
As such, there are totally 2 £ (n 2 2) 1 1 intersections. As a result, r(t) is (n 2 1)-modal. Jiang and Murthy [3] show that the failure rate for 2-fold model can be unimodal and this agrees with the above conclusion.
: 26
Differentiating Eq. (25) with respect to x yields y 0 x < Y 0n2i < b 00 1 and 00
y x <
Y 00n2i
bn2i zn2i exp zn2i 2 1
27
Case (3). In this case, y 0 y 0 2 1) is decreasing and tends to a positive constant as t ! 1. On the other hand, 2y 00 is unimodal, or bimodal, ¼, or (n 2 1)-modal and tends to zero as t ! 1. As a result, y 0 (y 0 2 1) does not intersect with 2y 00 or intersects with 2 £ 1, or 2 £ 2, ¼, or 2 £ (n 2 1) intersections, implying that r(t) is either mono-
" # zn2i exp zn2i b2n2i zn2i 12 : 28 < exp zn2i 2 1 exp zn2i 2 1 14
y’(y’-1)
12 10 8
y’
6 4 -y"
2 0 -1
-0.5
0
0.5
1
1.5
x=ln(t) Fig. 4. Curves y 0 ; y 00 ; and r 0 t.
2
2.5
3
R. Jiang et al. / Reliability Engineering and System Safety 74 (2001) 211±219
tonically increasing or k-modal-followed-by-increasing with k 1, 2, ¼, n 2 1. For n 2, Jiang and Murthy [3] show that the failure rate can be either monotonically increasing or unimodalfollowed-by-increasing and this agrees with the above conclusion. In summary, the shape types of the failure rate can be grouped into four groups: ² ² ² ²
monotonically decreasing; monotonically increasing; k-modal-followed-by-increasing, k 1, 2, ¼, n 2 1; and k-modal, k 1, 2, ¼, n 2 1.
Finally, it is worth noting that failure rate can never be decreasing for small t and increasing for large t. This rules out the possibility for the failure rate to have a bathtub shape. 7. Conclusions In this paper, we have carried out a detailed characterisation of the n-fold Weibull multiplicative model in terms of the shapes for the WPP plot, density and failure rate functions. This understanding is useful for determining whether or not a given set of data can be modelled by an n-fold
219
Weibull multiplicative model. We have also looked at a graphical method to estimate the model parameters.
Acknowledgements The ®rst and third authors wish to acknowledge the Hong Kong Polytechnic University for the ®nancial support of the project (No. G-YW41).
References [1] A.P. Basu, J.P. Klein. Some recent results in competing risks theory, in: Survival analysis, Vol. 2, Proceedings of the special topics meeting sponsored by the Institute of Mathematical Statistics. Columbus, Ohio, October 26±28, 1981, pp. 216±229. [2] Jiang R, Murthy DNP. The exponentiated Weibull family: a graphical approach. IEEE Transactions on Reliability 1999;48(1):68±72. [3] Jiang R, Murthy DNP. Parametric study of multiplicative model involving two Weibull distributions. Reliability Engineering and System Safety 1997;55:217±26. [4] Jiang R, Murthy DNP. Reliability modelling involving two Weibull distributions. Reliability Engineering and System Safety 1995;47:187± 98. [5] Jiang R, Murthy DNP, Ji P. Models involving two inverse Weibull distributions. Reliability Engineering and System Safety 2001;73(1):73±81. [6] Nelson W. Applied Life Data Analysis. New York: John Wiley, 1982.