Discrete competing risk model with application to modeling bus-motor failure data

ARTICLE IN PRESS Reliability Engineering and System Safety 95 (2010) 981–988

Contents lists available at ScienceDirect

Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress

Discrete competing risk model with application to modeling bus-motor failure data R. Jiang n Faculty of Automotive and Mechanical Engineering, Changsha University of Science and Technology, Changsha, Hunan 410114, China

a r t i c l e in fo

abstract

Article history: Received 27 September 2009 Received in revised form 15 April 2010 Accepted 18 April 2010 Available online 24 April 2010

Failure data are often modeled using continuous distributions. However, a discrete distribution can be appropriate for modeling interval or grouped data. When failure data come from a complex system, a simple discrete model can be inappropriate for modeling such data. This paper presents two types of discrete distributions. One is formed by exponentiating an underlying distribution, and the other is a two-fold competing risk model. The paper focuses on two special distributions: (a) exponentiated Poisson distribution and (b) competing risk model involving a geometric distribution and an exponentiated Poisson distribution. The competing risk model has a decreasing-followed-by-unimodal mass function and a bathtub-shaped failure rate. Five classical data sets on bus-motor failures can be simultaneously and appropriately fitted by a general 5-parameter competing risk model with the parameters being functions of the number of successive failures. The lifetime and aging characteristics of the fitted distribution are analyzed. & 2010 Elsevier Ltd. All rights reserved.

Keywords: Discrete distribution Exponentiated geometric distribution Exponentiated Poisson distribution Competing risk model Bathtub failure rate

1. Introduction Davis [1] presents five sets of typical repair–reuse data on busmotor major failures. The data are widely cited (e.g., see [2,3]), and have been fitted to various continuous distribution models, e.g., the normal and exponential distributions [1], the exponential-Rayleigh competing risk model [4], the exponential-power distribution [5], the exponentiated Weibull model [6], the multiplicative intensities model [7], the Weibull and mixed Weibull distributions [8], and the Gompertz-sinh model for the 2nd data set [9]. In addition, various statistical procedures are developed to test whether the 3rd–5th data sets follow an exponential distribution (e.g., see [10,11]). The previous works focus on seeking the best model to each of data sets and the considered models are not favored due to the following two main reasons: (a) the models lack of a plausible physical interpretation [7] and (b) the models cannot be used to predict the reliability after the jth repair for j Z5, which is needed in maintenance decision analysis. This necessitates developing more plausible and flexible models to fit lifetime data with complex features. Since the bus-motor failure data are given in the number of failures in an equispaced usage period, it is possible to model the n

Tel.: + 86 731 82309122; fax: +86 731 85258646. E-mail address: [email protected]

0951-8320/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ress.2010.04.009

data using a discrete distribution. It is desired that all the five data sets are modeled by a distribution family and the distributional parameters can be functions of the number of successive failures. This will considerably facilitate the maintenance decision analysis. Many discrete distributions have been developed in the literature (e.g., see [12–14]). However, the common discrete distributions (e.g., Poisson, geometric and negative binomial) can be unsuitable for modeling the bus-motor failure data due to the facts:

 A preliminary analysis to the data (to be presented in Section



2.2) shows that the empirical mass function (see Fig. 1) has a decreasing-followed-by-unimodal (DFU for short) shape and the empirical failure rate (see Fig. 2) is either increasing or bathtub-shaped. The number of intervals for the 4th and 5th data sets is 5 so that the number of the parameters of an appropriate distribution must be smaller than 5–1 ¼4 to ensure that the degree of freedom of the w2 test is not smaller than 1.

As a result, a possible model for modeling the bus-motor failure data must be sufficiently flexible and structurally simple. In this paper, we first present an exponentiated Poisson distribution (EPD for short). The EPD has only two distributional parameters and its failure rate can be increasing or bathtub shaped. We then show that a competing risk model involving a geometric distribution and an EPD has three parameters and can

ARTICLE IN PRESS 982

R. Jiang / Reliability Engineering and System Safety 95 (2010) 981–988

0.4

Table 1 Bus-motor failure data to the jth major failure (in 103 miles).

5th

fx

0.3

4th

0.2 1st 0.1

3rd 2nd

0 0

2

4

6

8

x Fig. 1. Empirical mass functions of the bus-motor failure data.

Interval

1st

Interval

2nd

Interval

3rd

Interval

4th

5th

0–20 20–40 40–60 60–80 80–100 100–120 120–140 140–160 160–180 180–up Size

6 11 16 25 34 46 33 16 2 2 191

0–20 20–40 40–60 60–80 80–100 100–120 120–140 140–up

19 13 13 15 15 18 7 4

0–20 20–40 40–60 60–80 80–100 100–up

27 16 18 13 11 16

0–20 20–40 40–60 60–80 80–up

34 20 15 15 12

29 27 14 8 7

96

85

104

101

0.9 0.8 0.7

The bus-motor failure data have the following problems (e.g., see [2,7]):

2nd

rx

0.6 5th

0.5 0.4

 The sample size varies, implying that there are some censored

3rd

4th

observations are not included.

0.3

 The individual sequences of failure times for each bus-motor

0.2

1st

0.1 0 0

2

4

6

8

10

x Fig. 2. Empirical failure rates of the bus-motor failure data.

have a bathtub-shaped failure rate and a DFU shaped mass function. Finally, we develop a 6-parameter competing risk model (with the parameters are functions of the number of successive failures) to fit all the five bus-motor failure data sets simultaneously. The competing risk model has a plausible physical interpretation and provides appropriate fit to the five data sets. After a comprehensive analysis, one sub-model (with 5 parameters) of the 6-parameter competing risk model is identified as the final model for all the five data sets. The paper is organized as follows. We present the bus-motor failure data and carry out a preliminary data analysis in Section 2. The proposed models are presented and their properties are discussed in Section 3. Section 4 deals with modeling of the data and analysis of the various fitted models. The paper is concluded in Section 5.

2. Bus-motor failure data and preliminary data analysis 2.1. Data Field data for bus-motor major failures (abrupt breakdown or an event in which the motor maximum power is smaller than a certain value) come from Davis [1] and are reproduced in Table 1. The description of bus-motor failure implies two failure modes (also see [7,8]): (a) Serious accidents (usually involved worn cylinders, pistons, piston rings, valves, camshafts or connecting rod or crankshaft bearings). (b) Performance deterioration (i.e., the maximum power produced, measured by dynamometer, fell below a specified proportion of the normal value). The former can be viewed as a random failure mode, and the latter can be viewed as a wear-out failure mode.



were not supplied. This makes it impossible to consider the dependences between the successive failures of each busmotor. Repair actions are not clearly stated whilst repair actions can have a significant effect on types and times of succeeding failures.

To make the modeling possible, the previous works explicitly or implicitly assume: (a) The failures of each bus-motor are independent but not identically distributed. (b) For a fixed j (1rjr5), the buses, no matter whether or not the times to the jth failure were observed by the time that the data were collected, are statistically similar in the sense that they all have the same failure distribution Fj(t) and hence the samples given in Table 1 represent the population’s behavior.

2.2. Empirical probability mass and failure rate functions Let Dt denote the interval length, tx ¼xDt, x¼0, 1, 2,y,K 1, denote the observed times, tK ¼N, and nx denote the failure number in the interval (tx, tx + 1). A typical grouped data set can be represented by fnx ,0 rx r K1g:

ð1Þ

The sample size or the total number of observations is given by N¼

K 1 X

nx :

ð2Þ

x¼0

The empirical probability mass function (pmf) is given by fx ¼ nx =N,

0 r x rK2:

The empirical failure rate function (frf) is given by , K 1 X ni , 0 r x r K2: rx ¼ nx

ð3Þ

ð4Þ

i¼x

Figs. 1 and 2 show the empirical pmf and frf of the bus-motor failure data, respectively. As can be seen from the figures, both the pmf and frf have complex shapes (e.g., DFU pmf and bathtub-shaped frf). This implies that a discrete distribution that

ARTICLE IN PRESS R. Jiang / Reliability Engineering and System Safety 95 (2010) 981–988

is suitable for modeling the data must meet the following requirements:

 the pmf can be decreasing, unimodal or DFU;  the frf can be bathtub-shaped or increasing and  the number of the parameters is not larger than 3.

983

The frf of the PD is increasing (e.g., see [16]). This property implies that the PD is appropriate for representing the wear-out failure. 3.2. Distribution derived by exponentiating an underlying distribution

3. Proposed models and their properties

Mudholkar and Srivasta [17] define the exponentiated Weibull model (EWM for short) as below:

To have sufficient flexibility, the potential model for modeling the bus-motor data should not be a simple model. A relatively complex model can be built in the following two ways: (a) modifying a simple model by introducing an additional parameter and (b) combining two simple models in the ways used in [15] (e.g., competing risk and mixture).

g

FðtÞ ¼ FW ðt; b, ZÞ,

g 40,

ð10Þ

where FW (t; b, Z) is the Weibull cdf with the shape and scale parameters b and Z, respectively. Generally, let F0(x) [f0(x)] denote the cdf [pmf] of an underlying discrete distribution, then one can define three exponentiated models by introducing an additional parameter g ( 40) as follows: 1. Cdf-based exponentiated model:

We derive some new models in such two ways as follows.

g

FðxÞ ¼ F0 ðxÞ:

ð11Þ g

2. Rf-based exponentiated model: R(x)¼[1  F0(x)] and 1 P g g f0 ðxÞ. 3. Pmf-based exponentiated model: f ðxÞ ¼ f0 ðxÞ=

3.1. Geometric and Poisson distributions Consider a discrete random variable X that assumes integer values from the set {0, 1, 2,y}. Let f(x), F(x) and R(x) denote the probabilities of X¼x, Xrx and X4x, respectively, and we call them probability mass, cumulative distribution and reliability functions (pmf, cdf and Rf), respectively. The relations between them are given by

Since the last two exponentiated models sometimes suffer from numerical problems (e.g., overflow), we focus on the model given by (11) in the remainder of the paper. 3.2.1. Properties of the exponentiated model given by (11)

f ð0Þ ¼ Fð0Þ ¼ 1Rð0Þ RðxÞ ¼ 1FðxÞ,

i¼0

f ðxÞ ¼ FðxÞFðx1Þ ¼ Rðx1ÞRðxÞ,

x Z 1:

ð5Þ Property 1. When x is large, we have f0(x)-0, F0(x 1)-1 and

The discrete failure rate function, r(x), is given by f ðxÞ , rðxÞ ¼ Rðx1Þ

rð0Þ ¼ f ð0Þ,

x Z 1:

g

ð6Þ

3.1.1. Geometric distribution The geometric distribution (GD for short) can be the simplest discrete distribution. The pmf, cdf and frf are given respectively by f ðxÞ ¼ pqx ,

FðxÞ ¼ 1qx þ 1 ,

rðxÞ ¼ p,

0 op o1,q ¼ 1p:

g

f ðxÞ ¼ ½F0 ðx1Þ þ f0 ðxÞg F0 ðx1Þ  gf0 ðxÞF0 1ðx1Þ  gf0 ðxÞ:

ð7Þ

Clearly, the pmf is decreasing and the frf is constant. These properties are similar to those for the exponential distribution. Actually, if T ( Z0) follows the exponential distribution and the observed times are at tx ¼xDt, then X ( ¼0, 1, 2,y) follows the GD. The proof is as follows. Let r denote the failure rate of the exponential distribution. We have f ðxÞ ¼ PrðX ¼ xÞ ¼ expðrtx Þexpðrtx þ 1 Þ ¼ ð1erDt ÞexpðrxDtÞ:

ð12Þ

This implies that the exponentiated model and the underlying distribution have the same pmf shape for large x.

When x-N, Property 2. Let z ¼F0(x  1). limgð1zÞ=ð1zg Þ ¼ 1 . As a result, we have:

z-1

and

z-1

rðxÞ 

gf0 ðxÞ g½1F0 ðx1Þ ¼ r0 ðxÞ g g 1F 0 ðx1Þ 1F 0 ðx1Þ

ð13Þ

This implies that the exponentiated model and the underlying distribution have the same frf shape for large x (also see Figs. 3 and 4).

When the underlying distribution F0(x) is the GD and PD, we obtain two new two-parameter discrete distributions as follows.

ð8Þ Letting p¼ 1 e  rDt and q¼e  rDt, (8) becomes (7). Namely, X follows the GD. This property of the GD implies that it is appropriate for representing the random failure.

0.4 0.35

q = 0.8, γ = 0.6

0.3

x l

f ðxÞ ¼

le

x!

,

l 4 0:

ð9Þ

The Microsoft Excel has a standard function Poisson (x, l, false or true) to directly evaluate its pmf or cdf. In this sense, the PD is simple.

r (x)

0.25

3.1.2. Poisson distribution The Poisson distribution (PD for short) is one of widely used distributions in reliability analysis. The pmf is given by

0.2 0.15 q = 0.8, γ = 1.5

0.1 0.05 0 0

2

4

6

8

10

x Fig. 3. Failure rates of the EGD.

12

14

ARTICLE IN PRESS 984

R. Jiang / Reliability Engineering and System Safety 95 (2010) 981–988

0.45

0.9 0.8

0.4

0.7

0.35

λ = 2.5, γ = 0.25

0.3

0.5

f (x)

r (x)

0.6 0.4 0.3

0.25

q = 0.8

0.2 0.15

λ = 2.5, γ = 0.75

0.2

q = 0.6

q = 0.95

0.1 0.05

0.1

0

0 0

5

10

15

0

1

2

3

x

3.2.2. Exponentiated geometric distribution The exponentiated geometric distribution (EGD) is given by ð14Þ

Different from the GD, the pmf of the EGD can be either decreasing or unimodal. Noting that f(0)¼pg and f(1)¼pg[(1+q)g  1], the pmf is unimodal when f(0)of(1) (i.e., g 4ln(2)/ln(1+q)); otherwise, decreasing. The frf of the EGD can be either decreasing or increasing (see Fig. 3). The frf is increasing when r(0) or(1) or ð1þ qÞg þ ð1qÞg 4 2

ð15Þ

otherwise, decreasing. 3.2.3. Exponentiated Poisson distribution The exponenstiated Poisson distribution (i.e., EPD) is given by !g x X lx lg FðxÞ ¼ e : ð16Þ x! i¼0 The pmf of the EPD can be decreasing or unimodal; the frf can be either increasing or bathtub-shaped (see Fig. 4). Noting the property given by (13), the frf is bathtub shaped when r(0) 4r(1), which holds when ð1þ lÞg þ elg o 2:

ð17Þ

Eq. (17) can be met when l is relatively large and/or g is relatively small. 3.3. Two-fold competing risk model The item failure can occur due to one of two competing causes. The time to failure (T1) due to cause 1 is distributed according to F1(t), and the time to failure (T2) due to cause 2 is distributed according to F2(t). The item failure is given by the minimum of T1 and T2, and F(t) is given by (see [15]) FðtÞ ¼ 1½1F1 ðtÞ½1F2 ðtÞ:

ð18Þ

Eq. (18) is a continuous distribution of a two-fold competing risk model. A discrete counterpart can be obtained by replacing t with x. In this case, the pmf is given by f ðxÞ ¼ f1 ðxÞR2 ðxÞ þ f2 ðxÞR1 ðxÞ þ f1 ðxÞf2 ðxÞ

ð20Þ

ð21Þ

and rðxÞ ¼ p þqr2 ðxÞ:

7

8

Eq. (22) implies that r(x) and r2(x) have the same failure rate shape. For example, r(x) can be bathtub-shaped if F2(x) is an EPD. According to the above discussion, in the competing risk model used to analyze the bus-motor data, we take F1(x) as the GD and F2(x) as the EPD. The model has three parameters; the pmf has all the shapes shown in Fig. 1; and the frf has all the shapes shown in Fig. 2. Furthermore, the GD can represent the random failure and the EPD can represent the wear-out failure, implying that the physical interpretation of the model is consistent with the background of the bus-motor failure data. As a result, the competing risk model can be appropriate for modeling the bus-motor failure data. Finally, the competing risk model (written as ‘‘CR model’’ hereafter) reduces into the EPD when q¼1, and the PD when q¼ g ¼1. 3.4. Life characteristics The probability that an item fails in (tx, tx + 1) is f(x). Assume that the time to failure for the item that fails in (tx, tx + 1) is uniformly distributed with a probability density function f(x)/Dt. The MTBF (m) is given by 1 Z tx þ 1 1 X X f ðxÞ dt ¼ m¼ t ðx þ0:5Þf ðxÞDt ¼ ½EðXÞ þ0:5Dt: ð23Þ D t x ¼ 0 tx x¼0 The lifetime variance (s2) is given by   1 Z tx þ 1 X f ðxÞ 1 dt ¼ VðXÞ þ s2 ¼ ðtmÞ2 ðDtÞ2 : Dt 12 x ¼ 0 tx

ð24Þ

The relative dispersion of the lifetime can be measured by the coefficient of variation s/m. The larger it is, the more dispersive the lifetime is. When s/m is close to 1, the distribution of time between two successive failures is equivalent to an exponential distribution [18].

4. Modeling the bus-motor failure data In this section, we use the CR model developed above to fit the bus-motor failure data. 4.1. Parameter estimation method and goodness-of-fit

Specially, when F1(x) is a GD, from (19) and (20) we have f ðxÞ ¼ pqx ½R2 ðxÞ þð1 þq=pÞf2 ðxÞ

6

ð19Þ

and the frf is given by rðxÞ ¼ r1 ðxÞ þ r2 ðxÞr1 ðxÞr2 ðxÞ:

5

Fig. 5. Mass functions of the competing risk model (l ¼ 2, g ¼3.5).

Fig. 4. Failure rates of the EPD.

FðxÞ ¼ ð1qx þ 1 Þg :

4 x

We use the maximum likelihood (ML) method to estimate the model parameters. For the interval data given by (1), the loglikelihood function is given by lnðLÞ ¼

ð22Þ

Eq. (21) implies that f(x) can be decreasing (e.g., when f2(x) is decreasing), unimodal or DFU (when f2(x) is unimodal), see Fig. 5.

K1 X

nx ln½Fðtx þ 1 ÞFðtx Þ:

ð25Þ

x¼0

The parameters can be obtained by directly maximizing ln(L) using the Solver of Microsoft Excel.

ARTICLE IN PRESS R. Jiang / Reliability Engineering and System Safety 95 (2010) 981–988

The goodness-of-fit is evaluated using the w2 statistic given by

w2 ¼

K 1 X ðnx Ex Þ2 , Ex x¼0

  Ex ¼ N Fðtx þ 1 ÞFðtx Þ :

ð26Þ

shown in Table 2. From the table, we have the following observations:

 As anticipated, the CR model provides appropriate fits to all the data sets, and is the unique model among the five models that P can appropriately fit the first data set. In terms of 5j ¼ 1 w2ðjÞ , the P5 2 CR model (with j ¼ 1 wðjÞ ¼ 11:4382) is superior to the EWM P5 (with j ¼ 1 w2ðjÞ ¼ 15:6293).

2

The smaller the w , the better the fitted model. Optionally, the goodness-of-fit can be measured by the p-value given by pv ¼ Pr{Y 4 w2}, where Y follows the chi-square distribution with a degree of freedom K 1  m, and m is the number of the estimated parameters (e.g., m¼3 for the CR model). The larger the p-value, the better the goodness-of-fit. A p-value that is not smaller than 0.2 is desired for a good model [2]. Since the PD model is nested in the EPD one, and the EPD is in turn nested in the CR model, we use the likelihood ratio procedure to test the null hypothesis that a data set follows the PD model or the EPD one against the alternative hypothesis of the CR model. The test statistic is _ _ D ¼ 2lnð L 0 = L 1 Þ

ð27Þ

_ _ where L 0 and L 1 are the estimated likelihoods (i.e., the likelihood functions evaluated at the ML estimates of the model parameters) associated with the null hypothesis (i.e., a simpler model, denoted as H0) and the alternative hypothesis (i.e., a more complex model, denoted as H1), respectively. The D statistic is approximately a chi-square variate with n degrees of freedom (n ¼n1  n0, where n1 and n0 are the numbers of parameters of the complex and simple models, respectively) when the null hypothesis H0 holds. Let d1  a denote the upper 100 (1  a) percentile point of the chi-square distribution with n degrees of freedom. We reject H0 if D4 d1  a. In this paper, we perform the likelihood ratio test at an a ¼0.05 significance level.

 The EPD provides appropriate fits for the 3rd–5th data sets, and is the best model in terms of p-value among the five models for the 3rd data set, whilst the PD cannot appropriately fit any one of these data sets. The GD provides fairly good fits to the 3rd and 5th data sets.



These illustrate the usefulness of the proposed exponentiated and CR models. A likelihood ratio test has been performed to assess whether or not the null hypothesis of a sub-model of the CR model has to be rejected, and the test results are shown in Table 3. As can be seen from the table, the null hypothesis for all the sub-models is Table 3 Likelihood ratio test d0.95(n ¼ 2)¼5.9915). Data set

1 2 3

4.2. Fitted models

4

Using the approach outlined above, we fitted the data to the CR model, EWM, EPD, PD and GD. The estimated model parameters, corresponding log-likelihood values and p-values are

985

5

H1:CR Reject H0? H1:CR Reject H0? H1:CR Reject H0? H1:CR Reject H0? H1:CR Reject H0? H1:EPD

of

the

fitted

models

(d0.95(n ¼ 1)¼3.8415,

H0:EPD

H0:PD

H0:GD

Preferred model

19.1208 Yes 9.1476 Yes 0 No 1.2306 No 0.2578 No

19.1510 Yes 34.8002 Yes 32.7382 Yes 20.0102 Yes 5.0790 No 4.8212 (Reject H0)

201.7944 Yes 34.7626 Yes 3.8360 No 4.3966 No 3.5518 No

CR CR EPD EPD EPD

Table 2 Parameters and goodness-of-fit. Data set

Model

1st

CR EWM EPD PD GD CR EWM EPD PD GD CR EWM EPD PD GD CR EWM EPD PD GD CR EWM EPD PD GD

2nd

3rd

4th

5th

q (or b) 0.9477 5.3565

0.8142 0.8265 13.2293

0.7525 1 3.9366

0.7167 0.6582 16.5284

0.6300 0.7578 0.9584

0.5784

l (or Z)

g

3.1630 133.89 4.3007 4.3418

4.2620 0.3916 1.0246

2.6820 140.38 4.5705 2.9423

7.9348 0.0688 0.4052

4.3432 118.03 4.3432 2.2053

0.3074 0.1909 0.3074

1.3040 98.09 2.9364 1.5370

24.2245 0.0395 0.3582

2.8518 24.83 1.8919 1.2840

0.6701 1.8311 0.5563

 ln(L)

w2

pv

387.4005 387.5464 396.9609 396.9760 488.2977 208.6969 209.9150 213.2707 226.0970 226.0782 176.9042 176.9874 176.9042 193.2733 178.8222 147.3900 147.3744 148.0053 157.3951 149.5883 124.1740 123.8257 124.3029 126.7135 125.9499

8.5235 10.9021 28.1422 27.7580 234.9502 1.4902 3.8617 11.4295 49.2667 41.9020 0.4696 0.6438 0.4696 41.0946 4.4811 0.1652 0.1349 1.4926 21.1634 5.2355 0.7897 0.0868 1.0414 6.5414 4.3673

0.2022 0.0914 0.0002 0.0005 0.0000 0.8284 0.4250 0.0435 6.6  10  9 1.9  10  7 0.7907 0.7248 0.9255 2.610  8 0.3448 0.6844 0.7134 0.4741 9.710  5 0.1553 0.3742 0.7682 0.5941 0.0880 0.2244

ARTICLE IN PRESS 986

R. Jiang / Reliability Engineering and System Safety 95 (2010) 981–988

rejected for the first two data sets, and the null hypothesis for EPD and GD models cannot be rejected for the last three data sets. After comprehensively considering the likelihood ratio test results, model simplicity (in terms of number of model parameters) and goodness-of-fit (given in Table 2), the preferred model for each data set is shown in the last column of Table 3, i.e., the CR model for the first two data sets and the EPD for the last three data sets. According to the above discussion, we have two candidate models for all the five data sets: (a) the CR model with totally 15 parameters (simply denoted as 15-CR model hereafter) and (b) the model with the CR model for the first two data sets and the EPD for the other three data sets (simply denote it as CR-EPD hereafter). Discussion: Lindsey [7] gives the failure data grouped into intervals of 10 (in 103 miles). Using the 15-CR model to fit the grouped data with Dt¼10 and 20, respectively, we obtained the estimates of the life mean and standard deviation shown in Table 4. As can be seen from the table, the estimates are fairly consistent. Since the grouped data with Dt ¼10 contain a number of failure observations that are smaller than 5 or equal to zero, we use the grouped data given in Table 1 in this paper. 4.3. A general model Lindsey [7] uses the multiplicative intensities model (MIM for short) 2

ln½rðtÞ ¼ b0 þ b1 ln ðtÞ

ð28Þ

to fit the five data sets. In the final MIM, b1 is the same for the 2nd–5th data sets, and hence the model totally has 7 parameters for the five data sets. Clearly, the fitted model cannot be used to predict the reliability after the jth repair (j Z5), while this is needed in maintenance decision analysis. Purohit [19] proposes a model for the rate of occurrence of failures, in which the number of prior failures is considered as a covariate. In this section we develop a general CR model with the parameters being functions of the number of successive failures. Let j (¼ 1, 2,y) denote the number of successive failures. Assume that the time between the (j 1)st and the jth failures follows the CR model with the parameters: qðjÞ ¼ q1 =ja , lðjÞ ¼ l1 =jb , g ¼ g0 ecj

Table 4 ML estimates of the life mean and standard deviation for two kinds of grouped data under the 15-CR model.

m 96.28 70.38 53.68 41.45 36.10 Average

4.4. Model analysis of the general model In this subsection, we confine our attention on the fitted 5-CR model. 4.4.1. Comparison of goodness-of-fit To illustrate the appropriateness of the 5-CR model, we compare its goodness-of-fit with those obtained from the EWM, MIM and the 15-CR model. The results are shown in Table 7. As can be seen from the table, the 15-CR model is superior to the EWM in terms of the p-value; and the 5-CR model is superior to the MIM in terms of the model simplicity and p-value. Particularly, the 5-CR model has the prediction ability for j Z5. As a result, the 5-CR model is the best among these models. It is noted that for the MIM given by (28) r(0)¼N and r0 (1)¼0. This implies that the failure rate is always bathtub-shaped and the change point is always at t¼1 no matter what the time scale is. Table 5 Parameters and goodness-of-fit of the 6-CR and 5-CR models.

q1 a

l1 b

g0 c  ln(L)

w2 pv

a, b, c a0

a¼ 0

b ¼0

c¼0

0.9444 0.2090 2.9544 -0.2406 9.8115 0.5516 1047.7090 18.2849 0.7419

0.9579 0 4.4988 0.2597 1.8087 0.3778 1068.1644 61.9337 3.39  10  5

0.9450 0.2154 3.1523 0 6.4562 0.2951 1049.4956 21.8525 0.5881

0.9509 0.2502 3.1937 0.0154 4.1512 0 1050.5860 23.8011 0.4730

ð29Þ

where q1, l1, g0, a, b and c are the model parameters to be specified. We simply denote it as 6-CR model hereafter. The model includes three 5-parameter models (5-CR model for short) as its special cases: a¼ 0 (i.e., q ¼q1), b¼0 (i.e., l ¼ l1) and c¼0 (i.e., g ¼ g0).

Dt ¼10 (  103)

We fitted these models to all the five data sets simultaneously, and obtained the results shown in Table 5. In terms of p-value, the 6-CR and 5-CR models with b¼0 and c ¼0 fit the data sets well. In terms of likelihood and p-value, the 5-CR model with b¼0 is superior to the 5-CR model with c ¼0. Hereafter, the 5-CR model means the one with b¼0. A likelihood ratio test has been performed to assess whether or not the null hypothesis of a sub-model of a CR model has to be rejected, and the results are shown in Table 6. As can be seen from the table, the null hypothesis of the CR-EPD or 6-CR [5-CR] model against the alternative hypothesis of the 15-CR [6-CR] model is not rejected. After comprehensively considering the model simplicity, applicability and goodness-of-fit, the 5-CR model is chosen as the final model.

Dt ¼20 (  103)

Relative error (%)

s

m

s

38.47 44.92 39.04 30.58 30.11

96.42 68.72 55.48 40.55 36.06

39.95 42.63 41.78 30.20 27.80

in m 0.15 2.36 3.35 2.17 0.12 1.63

Table 6 Likelihood ratio test of the nested models. H0

H1

D

n

d0.95(v)

Reject H0?

CR-EPD 6-CR 5-CR

15-CR 15-CR 6-CR

1.4884 6.2868 3.5732

3 9 1

7.8147 16.9190 3.8415

No No No

Table 7 Comparison of goodness-of-fit of different models.

in s 3.84 5.10 7.03 1.24 7.69 4.98

5-CR 15-CR EWM, from Table 2 MIM (Dt ¼ 20, [7])

Number of parameters

w2

Degree of freedom

p-Value

5 15 15 7

21.8525 11.4382 15.6293 23.4

24 14 14 18

0.5881 0.6513 0.3365 0.1757

ARTICLE IN PRESS R. Jiang / Reliability Engineering and System Safety 95 (2010) 981–988

Table 8 Life mean and standard deviation for different j.

0.4

EWM

0.3 f (x)

987

1st

5th

5-CR

j

m

s

E(X)

pffiffiffiffiffiffiffiffiffiffi V ðXÞ

m

s

s/m

1 2 3 4 5 6 7 8 9 10

95.85 66.39 54.05 39.91 36.43

36.44 39.43 41.20 28.72 29.05

4.3813 2.7020 2.0724 1.7117 1.4570 1.2507 1.0675 0.8979 0.7406 0.5985

2.0148 2.0355 1.8122 1.6265 1.4738 1.3450 1.2336 1.1336 1.0394 0.9472

97.6258 64.0409 51.4482 44.2334 39.1404 35.0136 31.3506 27.9575 24.8126 21.9692

40.7085 41.1168 36.7016 33.0378 30.0359 27.5132 25.3391 23.3956 21.5749 19.8048

0.4170 0.6420 0.7134 0.7469 0.7674 0.7858 0.8082 0.8368 0.8695 0.9015

0.2

0.1

0 0

1

2

3

4

5

6

7

8

9

x 0.4 0.35 0.3 4th

f (x)

0.25

0.8

0.2

j = 1, 2, …, 5, 10

2nd

0.15 0.6

0.1 r (x)

0.05 0 0

1

2

3

4

5

6

7

8

0.4

9 0.2

x 0.3

0 0

2

4

0.25 0.2 f (x)

6

8

10

x Fig. 7. Failure rates of the 5-CR model as functions of x and j.

0.15 0.1 Table 9 Aging characteristics of the 5-CR model.

0.05 0 0

1

2

3

4

5

6

7

8

9

x Fig. 6. Empirical (dotted lines) and estimated pmf: (a) pmf for the 1st and 5th data sets; (b) pmf for the 2nd and 4th data sets and (c) pmf for the 3rd data set.

4.4.2. Graphical comparison between the empirical and estimated pmf To further illustrate the appropriateness of the 5-CR model, we show the empirical and estimated pmf for j ¼1, 2,y,5 in Fig. 6(a)–(c). As can be seen from the figures, the empirical and estimated pmf values are fairly close to each other. 4.4.3. Lifetime characteristics The lifetime characteristics of the fitted 5-CR model are summarized in Table 8. The 2nd and 3rd columns show that the estimates of m and s obtained from the EWM. As can be seen, the estimates of m and s from the two models are fairly close to each other, and hence the appropriateness of the (23) and (24) are illustrated. The further remarks can be made from the table:

j

1

2

3

4

5

p r(0) r(0)  p

0.0550 0.0550 0.0000

0.1860 0.1860 0.0000

0.2541 0.2543 0.0002

0.2989 0.3003 0.0014

0.3318 0.3382 0.0064

failure rate is increasing for small j (¼ 1, 2,y,9) and bathtubshaped for large j ( Z10). This is clearly shown in Fig. 7. Another observation from the figure is that the failure rate for large j is close to a constant, implying that X (observed at tx ¼ xDt) tends to the GD or T tends to the exponential distribution. This phenomenon can be explained by the interaction between the random and wear-out failures based on (22). Table 9 shows the GD parameter p ( ¼1 q) and initial failure rate (or the failure probability in the first usage interval) r(0) as functions of j. As can be seen from the table,

 p increases as j increases, implying that the random failure mode becomes a main failure mode after several repairs and

 The MTBF quickly decreases for small j and slowly decreases

 r(0) Ep for small j, implying that the early failure is mainly due

for large j, implying that the repairs after major failures are imperfect and s/m quickly increases and tends to 1 as j increases, implying that the time between two successive failures is approximately equivalent to an exponential distribution when j is large.

to the random failure mode but may be due to the wear-out failure mode when j is large (see the last row in Table 7).



4.4.4. Aging characteristics In the fitted 5-CR model, both q(j) and g(j) decreases as j increases. According to the analysis in Sections 3.2 and 3.3, the

Fig. 8 shows the failure rate of the EPD nested in the 5-CR model, r2(x), as functions of x and j. As can be seen, the failure rate curve becomes flatter (though it can be non-monotonic) as j increases, implying that the wear-out failure occurs earlier for larger j.

ARTICLE IN PRESS 988

R. Jiang / Reliability Engineering and System Safety 95 (2010) 981–988

fields for modeling and analyzing various grouped or count data (e.g., failure, defects and accident data).

0.8

r2 (x)

0.6 j = 1, 2, …, 5, 10

Acknowledgements

0.4

The author would like to thank the editor and anonymous referees for their useful comments and suggestions for improvement on earlier versions of this paper. This research was supported by the National Natural Science Foundation (No. 70771015).

0.2

0 0

2

4

6

8

10

x

References

Fig. 8. Failure rate of the EPD nested in the 5-CR model as functions of x and j.

The two trends (increasing p and advancing wear-out failure) interact and result in the outcome that the failure rate tends to a constant as j increases.

5. Conclusions In this paper, we have examined the features of the bus-motor failure data and identified the requirements for a discrete distribution that can appropriately model the bus-motor failure data. We have developed two discrete distribution families with the bathtub-shaped failure rate. It has been shown that the competing risk model involving a GD and an EPD can provide appropriate fits to all the five sets of bus-motor failure data. The model has a physical interpretation: the GD represents the random failure and the EPD represents the wear-out failure. We have shown the equivalence between the discrete GD and the continuous exponential distribution. A 6-parameter general model has been developed by considering the three parameters of the competing risk model as functions of the number of successive failures. After comprehensively considering the likelihood ratio test results, goodness-of-fit and model simplicity, one of its sub-models has been chosen as the final model. The 5-parameter model not only fits the observed data sets well but also can be used to predict the reliability after j (jZ5) repairs. This considerably facilitates the maintenance decision analysis. Finally, the dynamic aging characteristics of the bus-motor have been analyzed in detail using the fitted 5-parameter model. The proposed models (EGD, EPD and the competing risk model involving GD and EPG) are attractive due to their simplicity and flexibility, and hence can be used in reliability, quality and safety

[1] Davis DJ. An analysis of some failure data. Journal of the American Statistical Association 1952;47:113–50. [2] Blischke WR, Murthy DNP. Reliability. New York: Wiley; 2000. [3] Lai CD, Xie M. Stochastic Ageing and Dependence for Reliability. New York: Springer; 2006. [4] Bain LJ. Analysis for the linear failure-rate life-testing distribution. Technometrics 1974;16:551–9. [5] Smith RM, Bain LJ. An exponential power life-testing distribution. Communications in Statistics, Part A—Theory and Methods 1975;4:469–81. [6] Mudholkar GS, Srivasta DK. Exponentiated Weibull family: a reanalysis of the bus-motor-failure data. Technometrics 1995;37(4):436–45. [7] Lindsey JK. Parametric multiplicative intensities models fitted to bus motor failure data. Applied Statistics 1997;46(245–252):542. [8] Wang KS, Po HJ, Hsu FS, Liu CS. Analysis of equivalent dynamic reliability with repairs under partial information. Reliability Engineering and System Safety 2002;76:29–42. [9] Cooray K, Ananda MAM. Analyzing survival data with highly negatively skewed distribution: The Gompertz-sinh family. Journal of Applied Statistics 2010;37:1–11. [10] Gulati S, Neus J. Goodness of fit statistics for the exponential distribution when the data are grouped. Communications in Statistics—Theory and Methods 2003;32:681–700. [11] Best DJ, Rayner JCW. Chi-squared components for tests of fit and improved models for the grouped exponential distribution. Computational Statistics & Data Analysis 2007;51(8):3946–54. [12] Khalique A. On discrete failure-time distributions. Reliability Engineering and System Safety 1989;25(2):99–107. [13] Johnson NL, Kotz S, Kemp AW. Univariate Discrete Distributions, 2nd ed. New York: Wiley; 1992. [14] Bracquemond C, Gaudoin O. A survey on discrete lifetime distributions. International Journal of Reliability, Quality and Safety Engineering 2003;10:69–98. [15] Jiang R, Murthy DNP. Reliability modeling involving two Weibull distributions. Reliability Engineering and System Safety 1995;47:187–98. [16] Gupta PL, Gupta RC, Tripathi RC. On the monotonic properties of discrete failure rates. Journal of Statistical Planning and Inference 1997;65(2):255–68. [17] Mudholkar GS, Srivastava DK. Exponentiated Weibull family for analyzing bathtub failure-rate data. IEEE Transactions on Reliability 1993;42(2): 299–302. [18] Jiang R, Ji P, Xiao X. Aging property of unimodal failure rate models. Reliability Engineering and System Safety 2003;79(1):113–6. [19] Purohit SG. Recurring event data: a general model and related inference procedures. Journal of Statistical Planning and Inference 1997;59:31–44.