Study on segmented distribution for reliability evaluation

CJA 741 27 December 2016 Chinese Journal of Aeronautics, (2016), xxx(xx): xxx–xxx

No. of Pages 21

1

Chinese Society of Aeronautics and Astronautics & Beihang University

Chinese Journal of Aeronautics [email protected] www.sciencedirect.com

4

Study on segmented distribution for reliability evaluation

5

Li Huaiyuan a,b, Zuo Hongfu b,*, Su Yan b, Xu Juan b, Yin Yibing b

3

6 7

8

a b

School of Software Engineering, Jinling Institute of Technology, Nanjing 211169, China College of Civil Aviation, Nanjing University of Aeronautics and Astronautics, Nanjing 211106, China

Received 17 December 2015; revised 31 January 2016; accepted 18 October 2016

9

11 12

KEYWORDS

13

Bayesian estimates; Gibbs sampling method; Maximum likelihood estimation; Proximity sensors in leadingedge flap in aircraft; Reliability; Segmented distribution; Weibull distribution

14 15 16 17 18 19 20 21

22

Abstract In practice, the failure rate of most equipment exhibits different tendencies at different stages and even its failure rate curve behaves a multimodal trace during its life cycle. As a result, traditionally evaluating the reliability of equipment with a single model may lead to severer errors. However, if lifetime is divided into several different intervals according to the characteristics of its failure rate, piecewise fitting can more accurately approximate the failure rate of equipment. Therefore, in this paper, failure rate is regarded as a piecewise function, and two kinds of segmented distribution are put forward to evaluate reliability. In order to estimate parameters in the segmented reliability function, Bayesian estimation and maximum likelihood estimation (MLE) of the segmented distribution are discussed in this paper. Since traditional information criterion is not suitable for the segmented distribution, an improved information criterion is proposed to test and evaluate the segmented reliability model in this paper. After a great deal of testing and verification, the segmented reliability model and its estimation methods presented in this paper are proven more efficient and accurate than the traditional non-segmented single model, especially when the change of the failure rate is time-phased or multimodal. The significant performance of the segmented reliability model in evaluating reliability of proximity sensors of leading-edge flap in civil aircraft indicates that the segmented distribution and its estimation method in this paper could be useful and accurate. Ó 2016 Production and hosting by Elsevier Ltd. on behalf of Chinese Society of Aeronautics and Astronautics. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/ licenses/by-nc-nd/4.0/).

* Corresponding author. E-mail addresses: [email protected] (H. Li), [email protected], [email protected] (H. Zuo), [email protected] (J. Xu). Peer review under responsibility of Editorial Committee of CJA.

Production and hosting by Elsevier

1. Introduction

23

As one of the essential devices affecting aircraft safety, the proximity sensors in leading-edge flap is an important part of aircraft operating system. However, proximity sensors in leading-edge flap are also one of the devices with higher failure rate in aircraft.1 In order to ensure aircraft safety and reduce its maintenance costs, it’s necessary to accurately evaluate the reliability of proximity sensors in leading-edge flap.2

24

http://dx.doi.org/10.1016/j.cja.2016.12.008 1000-9361 Ó 2016 Production and hosting by Elsevier Ltd. on behalf of Chinese Society of Aeronautics and Astronautics. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

25 26 27 28 29 30

CJA 741 27 December 2016

No. of Pages 21

2 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92

However, like most other equipment, the failure rate of the proximity sensors in leading-edge flap is complicated, exhibiting different changes in different operation phases.3 In order to accurately evaluate the reliability of devices such as proximity sensors in leading-edge flap, whose failure rate is multimodal or multi-stage, only segmented distribution is applicable and inevitable.4 In this paper, two segmented reliability models are proposed: (1) failure rate is a piecewise linear function about lifetime, namely, segmented linear failure rate mode and (2) failure rate in any stage is of the same formula as the failure rate function of Weibull distribution, called segmented Weibull distribution. Currently, the application of segmented distribution to reliability evaluation is mainly classified into three categories: describing the degradation process of devices; creating reliability model; and testing and detecting change-points.5 (1) Now, how to fit the equipment degradation process by segmented distribution is discussed in a lot of literature. For example, Bae et al. 6 proposed a hierarchical Bayesian change-point regression model to fit two-phase degradation patterns and derived the failure-time distribution of a unit randomly selected from its population. Kvam 7 introduced a log-linear model with random coefficients and a change-point to describe the nonlinear degradation path. (2) Since the failure rates of many products may perform different trends at different stages, an application of segmented distribution to evaluating product reliability has been widely concerned. Li et al. 8 estimated the change-point for a piecewise hazard regression model in the presence of right censoring and long-term survivors. In a study by He et al. 9, a sequential testing approach to detect multiple change-points in the hazard function by likelihood ratio statistics and resampling was proposed, which was applicable to both right-censored and interval-censored data. Uhm et al. 10 studied a weighted least squares estimator for Aalen’s additive risk model with right-censored survival data which may allow for very flexible handling of covariates. (3) In addition, another important application of segmented distribution is to test whether there are change-points in product failure rate or reliability. For a segmented regression system with an unknown change-point over two domains of a predictor, a new empirical likelihood ratio statistic was proposed to test the null hypothesis of no change by Liu and Qian 11. Goodman et al. 12 expanded the set of alternatives to allow for multiple change-points, and proposed a model selection algorithm using sequential testing for the piecewise constant hazard model. Suresh 13 derived a test statistic and its asymptotic distribution, and compared the power of the test with other existing tests such as likelihood ratio, Weibull, and log Gamma tests. Nosek and Szkutnik 14 discussed a regression model with a possible structural change and a small number of measurements. Recently, segmented distribution and its estimation method have caught more attentions. To minimize the expected sum of manufacturing cost, bum-in cost, and warranty cost of failed items found during their warranty period, a cost model was formulated to find the optimal bum-in time based on Weibull hyperexponential distribution by Chou and Tang 15. Considering different prior densities for parameters and censored survival data, Achcar and Loibel 16 discussed Bayesian analysis on segmented constant hazard function models and put forward their inference methods by Metropolis algorithms. Weibull-exponential distribution, a common segmented

H. Li et al. distribution, was discussed and its accurate calculation formula by Bayes estimation was proposed by Boukai 17. Patra and Dey 18 proposed a general class of change-point hazard models for survival data, which included and extended many different types of segmented distribution. However, the existent segmented distribution is generally considered as a two-segment distribution, or their failure rate is supposed as a constant in recent literature 6–18. Compared with many other current segmented distributions, the contributions of this paper are as follows. (1) In this paper, not only two kinds of universal n-segment distribution have been discussed, but also two general methods to estimate parameters in the segmented distributions, namely Bayesian estimation and maximum likelihood estimation (MLE), are given. (2) An information criterion adequate for the n-segment distribution is proposed. It can not only test whether a change-point exists or not, but also measure the appropriateness and correctness of the segmented reliability model. (3) A new bathtub curve model is put forward to verify the effects of the segmented distributions.

93

2. Segmented distribution

113

In fact, the failure rate of equipment not only is changeable but also often manifests different trends at different phases. Reliability-centered maintenance (RCM) recommended six failure models as shown in Fig. 1.19 As Fig. 1 shows, apart from failure models C and E, the failure rates of the other failure models show different patterns at different phases, and change-points obviously exist in failure rate curves. Therefore, if one and same distribution is utilized to evaluate a failure whose failure rate is characterized by multiple change-points and multi-peak, evaluation may be rather difficult and inaccurate.20 Due to unevenness and inconsistency of failure rate, it is not appropriate to describe the change of reliability by one and same analytical function in the entire lifetime cycle. However, if we divide lifetime into some time intervals according to failure rate trend, the failure rate in every different interval can be described by a different function. Just as piecewise interpolation with a simple function can improve approximation precision, we can more accurately approximate the true change of failure rate. For example, as for failure model A, i.e., bathtub curve, if three different failure rate functions are in correspondence to each phase, then the change of failure rate in the

114

Fig. 1

Six kinds of failure model.

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112

115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135

CJA 741 27 December 2016

No. of Pages 21

Study on segmented distribution 136 137 138 139 140 141 142 143

145 146 147 148 149 150 151 152 153 154 155 156 157 158 160

entire lifetime cycle of equipment can be accurately fitted by segmented distribution. Supposing that lifetime is divided into n intervals and the failure rate function in the ith interval ½si ; siþ1 Þ is expressed by hi ðtÞ, then the segmented failure rate function hðtÞ can be expressed by Eq. (1) below, where hi ðtÞ > 0, t 2 ½si ; siþ1 Þ is requisite. 8 0 t < s1 > > > > > > > > h1 ðtÞ s1 6 t < s2 > > > > > > > > . > > < .. hðtÞ ¼ ð1Þ > > > hi ðtÞ si 6 t < siþ1 > > > > > > > .. > > > . > > > > > : hn ðtÞ sn 6 t < 1 From Eq. (1), devices have no failure probability before s1 , but may fail after s1 . Usually, s1 is a time point when devices are likely to begin degeneration, such as the first use of devices, namely, s1 ¼ 0.17,18 However, s1 –0 is possible in rare particular cases. If hi1 ðsi Þ–hi ðsi Þ, then Eq. (1) is discontinuous at si , or Eq. (1) is continuous. In most cases, si may be a changepoint of failure rate. For convenience, notation Ji and Ki(t) are introduced, expressed by Eqs. (2) and (3). According to Rt Rt fðtÞ ¼ hðtÞ expð s1 hðuÞduÞ and RðtÞ ¼ expð s1 hðuÞduÞ, we can get the probability density function fi ðtÞ and the reliability function Ri ðtÞ in the ith interval as Eqs. (4) and (5), respectively. Z siþ1 Ji ¼ hi ðuÞdu ð2Þ si

161

163

Z Ki ðtÞ ¼

si

t

hi ðuÞdu;

t 2 ½si ; siþ1 Þ

164 i1 X fi ðtÞ ¼ hi ðtÞ exp  Jj  Ki ðtÞ

166

i1 X Ri ðtÞ ¼ exp  Jj  Ki ðtÞ

170 171 172 173 174

176

ð3Þ

si 6 t < siþ1

ð7Þ

.. . sn 6 t < 1 179

2.1. Segmented linear failure rate model

180

2.1.1. n-segment linear failure rate model

181

As calculation of piecewise linear interpolation is easy and can meet most engineering requirements, we divide lifetime into n intervals and set the failure rate function hi ðtÞ in the ith interval ½si ; siþ1 Þ as Eq. (8), where li 6 si is required. Usually, li ¼ 0 or li ¼ si . According to Eqs. (2) and (3), we get Ji and Ki ðtÞ shown as Eqs. (9) and (10), respectively.

182

hi ðtÞ ¼ ai þ bi ðt  li Þ

li 6 si ;

t 2 ½si ; siþ1 Þ

ð8Þ

183 184 185 186 187 188 190 191

1 1 Ji ¼ ai siþ1  bi li siþ1 þ bi s2iþ1  ai si þ bi li si  bi s2i 2 2

ð9Þ

193 194

1 1 Ki ðtÞ ¼ ai t  bi li t þ bi li si þ bi t2  ai si  bi s2i 2 2

t 2 ½si ; siþ1 Þ ð10Þ

Furthermore, according to Eqs. (6)–(10), we can obtain the probability density function fi ðtÞ and the reliability function Ri ðtÞ in the ith interval, respectively, shown as Eqs. (11) and (12). fi ðtÞ ¼ hi ðtÞe



i1 P

196 197 198 199 200 201

ðaj sjþ1 bj lj sjþ1 þ12bj s2jþ1 aj sj þbj lj sj 12bj s2j Þðai tbi li tþ12bi t2 ai si þbi li si 12bi s2i Þ

j¼1

t 2 ½si ;siþ1 Þ

ð11Þ

203 204

t 2 ½si ; siþ1 Þ

ð4Þ

! t 2 ½si ; siþ1 Þ

ð5Þ

j¼1

Furthermore, fi ðtÞ and Ri ðtÞ can piecewise compose the segmented probability density function fðtÞ and the reliability function RðtÞ in the entire lifetime cycle, shown as Eqs. (6) and (7), respectively. 8 h1 ðtÞ expðK1 ðtÞÞ s1 6 t < s2 > > > > > h2 ðtÞ expðJ1  K2 ðtÞÞ s2 6 t < s3 > > > > .. .. > > > . > > ! . > i1 > X < si 6 t < siþ1 fðtÞ ¼ hi ðtÞ exp  Jj  Ki ðtÞ ð6Þ > j¼1 > > > > > .. .. > > . > > ! . > > n1 X > > > > sn 6 t < 1 : hn ðtÞ exp  Jj  Kn ðtÞ j¼1

177

s1 6 t < s2 s2 6 t < s3 .. .

j¼1

!

j¼1

167

169

3 8 expðK1 ðtÞÞ > > > > > expðJ1  K2 ðtÞÞ > > > > .. > > > . > > ! > > i1 X < RðtÞ ¼ exp  Jj  Ki ðtÞ > j¼1 > > > > > . > .. > > > ! > > n1 X > > > > : exp  Jj  Kn ðtÞ



i1 X

ðaj sjþ1 bj lj sjþ1 þ12bj s2jþ1 aj sj þbj lj sj 12bj s2j Þðai tbi li tþ12bi t2 ai si þbi li si 12bi s2i Þ

Ri ðtÞ ¼ e j¼1 t 2 ½si ;siþ1 Þ

ð12Þ

If bi > 0, then the failure rate in the ith interval will linearly increase; if bi ¼ 0, the failure rate in the ith interval will be constantly kept as ai ; if bi < 0, then the failure rate in the ith interval will linearly decrease. If t 2 ½tn ; 1Þ, then bn > 0 is required in order to ensure hn ðtÞ > 0; t 2 ½tn ; 1Þ, but if hðtÞ is defined in finite support on some special occasions, bn 6 0 may be acceptable.

206 207 208 209 210 211 212 213

2.1.2. Linear failure rate model

214

Particularly, if lifetime isn’t segmented, namely, nonsegmented, then the failure rate will linearly change in the entire lifetime cycle. hðtÞ, fðtÞ and RðtÞ of the linear failure rate mode are expressed as Eqs. (13)–(15), respectively. Fig. 2 shows some typical function curves of this distribution, from which we can see that l1 and s1 can adjust the curve shape and its starting point. Specially, when b1 ¼ 0 and l1 ¼ s1 , the distribution will become an exponential distribution.

215

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

216 217 218 219 220 221 222 223 224

CJA 741 27 December 2016

No. of Pages 21

4

H. Li et al.

Fig. 2

226

Some curves of linear failure rate model.

hðtÞ ¼ a1 þ b1 ðt  l1 Þ s1 6 t < 1

ð13Þ

227

fðtÞ ¼ ða1 þ b1 ðt  l1 ÞÞeða1 tb1 l1 tþ2b1 t 1

2 a s þb l s 1b s2 Þ 1 1 1 1 1 2 1 1

232

233 234 235 236 237 238 239 240 241 242 243 244 245 246

248

eða1 tb1 l1 tþ2b1 t 1

2 a

1 2 1 s1 þb1 l1 s1 2b1 s1 Þ

s1 6 t < s2

eða1 s2 b1 l1 s2 þ2b1 s2 a1 s1 þb1 l1 s1 2b1 s1 Þða2 tb2 l2 tþ2b2 t 1

2

1

2

1

2 a

1 2 2 s2 þb2 l2 s2 2b2 s2 Þ

s1 6 t < 1

s2 6 t < 1 ð17Þ

252

ð14Þ

229 230

( RðtÞ ¼

RðtÞ ¼ eða1 tb1 l1 tþ2b1 t 1

2 a

1 2 1 s1 þb1 l1 s1 2b1 s1

Þ

s1 6 t < 1

ð15Þ

2.1.3. Two-segment linear failure rate model In general, n ¼ 2 or n ¼ 3, namely two- or three-segment distribution can meet most of reliability evaluation. For twosegment distribution, hðtÞ and RðtÞ will be Eqs. (16) and (17), respectively. Supposing that l1 ¼ s1 ¼ 0, l2 ¼ s2 ¼ 1:5, some curves of the two-segment linear failure rate model are shown in Fig. 3. From Fig. 3, it can be seen that the segmented distribution may degenerate into a non-segmented mode when the failure rate changes consistently and continuously in the entire lifetime cycle. When the failure rate changes in an obviously phased manner, this segmented distribution can also precisely reflect the different changes of the failure rate in its every interval.  a1 þ b1 ðt  l1 Þ s1 6 t < s2 hðtÞ ¼ ð16Þ a2 þ b2 ðt  l2 Þ s2 6 t < 1

2.1.4. Three-segment linear failure rate model

253

When n=3, hðtÞ and R(t) of the three-segment linear failure rate model are respectively shown as Eqs. (18) and (19). When l1 ¼ s1 ¼ 0, l2 ¼ s2 and l3 ¼ s3 ¼ 2:5, some function curves of the three-segment linear failure rate model are shown in Fig. 4. From Fig. 4, it’s clear that the adaptability of segmented distribution is extensive, as it can express both the case where the failure rate changes sharply and the case where the failure rate has no obvious change when ai  aiþ1 ; bi  biþ1 . From Figs. 3 and 4, it’s apparent that fðtÞ seems more fluctuating when the trend of the failure rate is obviously in a staged or piecewise manner.

254

8 > < a1 þ b1 ðt  l1 Þ s1 6 t < s2 hðtÞ ¼ a2 þ b2 ðt  l2 Þ s2 6 t < s3 > : a3 þ b3 ðt  l3 Þ s3 6 t < 1

265

255 256 257 258 259 260 261 262 263 264

ð18Þ 267

249 250

268 269

8 ða1 tb1 l1 tþ12b1 t2 a1 s1 þb1 l1 s1 12b1 s21 Þ > : ða1 s2 b1 l1 s2 þ1b1 s2 a1 s1 þb1 l1 s1 1b1 s2 Þða2 s3 b2 l2 s3 þ1b2 s2 a2 s2 þb2 l2 s2 1b2 s2 Þða3 tb3 l3 tþ1b3 t2 a3 s3 þb3 l3 s3 1b3 s2 Þ 2 2 2 2 2 2 2 1 3 2 3 e

s1 6 t < s2 s2 6 t < s3 s3 6 t < 1

ð19Þ

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

CJA 741 27 December 2016

No. of Pages 21

Study on segmented distribution

270

271 272 273 274 275 276 277 278 279 280 281

5

Fig. 3

Some curves of two-segment linear failure rate model.

Fig. 4

Some curves of three-segment linear failure rate model.

2.2. Segmented Weibull distribution Although the n-segment linear failure rate function hðtÞ can express the change of failure rate in a staged manner, it can only express a linear relation. Meanwhile, linear fitting is not precise enough and its adaptability is not good enough. Nevertheless, by adjusting shape parameters, Weibull distribution can not only express a distribution with an increasing failure rate, but also model a distribution whose failure rate is decreasing.21,22 Therefore, the failure rate in each interval is fitted by Weibull failure rate functions with different parameters, and better precision and adaptability can be achieved.23,24

2.2.1. n-segment Weibull distribution

282

Let the failure rate function hi ðtÞ in the ith interval ½si ; siþ1 Þ be Eq. (20). From Eqs. (2) and (3), we obtain Ji and Ki ðtÞ, respectively shown as Eqs. (21) and (22). hi ðtÞ ¼ ai bi ðai ðt  li ÞÞbi 1

t 2 ½si ; siþ1 Þ

ð20Þ

283 284 285 286 288 289

Ji ¼ ðai ðsiþ1  li ÞÞbi  ðai ðsi  li ÞÞbi

ð21Þ

291 292

Ki ðtÞ ¼ ðai ðt  li ÞÞbi  ðai ðsi  li ÞÞbi

t 2 ½si ; siþ1 Þ

ð22Þ

Here li 6 si ; ai > 0; bi > 0. Furthermore, according to Eqs. (6), (7) and (20)–(22), the probability density function fi ðtÞ

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

294 295 296

CJA 741 27 December 2016

No. of Pages 21

6 297 298 299

H. Li et al.

and the reliability function Ri ðtÞ in the ith interval are respectively obtained as Eqs. (23) and (24). 

301

fi ðtÞ ¼ ai bi ðai ðt  li ÞÞ t 2 ½si ;siþ1 Þ

bi 1

e

i1 X

ððaj ðsjþ1 lj ÞÞbj ðaj ðsj lj ÞÞbj Þððai ðtli ÞÞbi ðai ðsi li ÞÞbi Þ

j¼1

ð23Þ

Let l1 ¼ s1 ¼ 0, some curves of the two-segment Weibull distribution are drawn in Fig. 5. ( a1 b1 ða1 ðt  l1 ÞÞb1 1 s1 6 t < s2 ð25Þ hðtÞ ¼ a2 b2 ða2 ðt  l2 ÞÞb2 1 s2 6 t < 1

302

318 319 320

322 323 324

( fðtÞ ¼

b

b

a1 b1 ða1 ðt  l1 ÞÞb1 1 eðða1 ðtl1 ÞÞ 1 ða1 ðs1 l1 ÞÞ 1 Þ

s1 6 t < s2

b2 1 ðða1 ðs2 l1 ÞÞb1 ða1 ðs1 l1 ÞÞb1 þða2 ðtl2 ÞÞb2 ða2 ðs2 l2 ÞÞb2 Þ

a2 b2 ða2 ðt  l2 ÞÞ

e

i1 X

ððaj ðsjþ1 lj ÞÞbj ðaj ðsj lj ÞÞbj Þððai ðtli ÞÞbi ðai ðsi li ÞÞbi Þ



Ri ðtÞ ¼ e t 2 ½si ; siþ1 Þ

RðtÞ ¼

j¼1

304

ð24Þ

313

If bi > 1, then hi ðtÞ increases in the interval ½si ; siþ1 Þ; if bi ¼ 2, then hi ðtÞ becomes a linearly increasing function; if bi ¼ 1, then hi ðtÞ is fixed at ai in the interval ½si ; siþ1 Þ; if bi < 1, then hi ðtÞ decreases in the interval ½si ; siþ1 Þ. In this way, the segmented Weibull distribution can flexibly approximate the changeful failure rate. If lifetime is non-segmented, it will become a classical Weibull distribution. Usually, the segment number n 6 3 can satisfy the requirement of reliability evaluation in most cases.

314

2.2.2. Two-segment Weibull distribution

315

In particular, if lifetime is divided into two intervals, according to Eqs. (23) and (24), hðtÞ; fðtÞ and R(t) of the two-segment Weibull distribution are respectively shown as Eqs. (25)–(27).

305 306 307 308 309 310 311 312

316 317

Fig. 5

ð26Þ

s2 6 t < 1

325

8 ððk ðtl ÞÞb1 ðk ðs l ÞÞb1 Þ 1 1 1 1 1 > :

b

b

s 1 6 t < s2 ð27Þ b

b

eððk1 ðs2 l1 ÞÞ 1 ðk1 ðs1 l1 ÞÞ 1 þðk2 ðtl2 ÞÞ 2 ðk2 ðs2 l2 ÞÞ 2 Þ s2 6 t < 1

327

2.2.3. Three-segment Weibull distribution

328

Similarly, when n = 3, then hðtÞ and RðtÞ of the three-segment Weibull distribution can be respectively written as Eqs. (28) and (29). When l1 ¼ s1 ¼ 0, some three-segment curves of the Weibull distribution are shown in Fig. 6.

329

8 > a1 b1 ða1 ðt  l1 ÞÞb1 1 > > > > > < hðtÞ ¼ a2 b2 ða2 ðt  l2 ÞÞb2 1 > > > > > > : a3 b3 ða3 ðt  l3 ÞÞb3 1

330 331 332 333 334

s1 6 t < s2 s2 6 t < s3

ð28Þ

s3 6 t < 1

Some curves of two-segment Weibull distribution.

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

336

CJA 741 27 December 2016

No. of Pages 21

Study on segmented distribution

7

8 ðða ðtl ÞÞb1 ða1 ðs1 l1 ÞÞb1 Þ > : ðða1 ðs2 l1 ÞÞb1 ða1 ðs1 l1 ÞÞb1 þða2 ðs3 l2 ÞÞb2 ða2 ðs2 l2 ÞÞb2 þða3 ðtl3 ÞÞb3 ða3 ðs3 l3 ÞÞb3 Þ e

337 338 339 340 341 342 343 344 345 346 347 348 349 350

351

From Figs. 5 and 6, the segmented Weibull distribution can precisely approximate the sharp fluctuation of failure rate. As for multimodal distribution or varying failure rate, such as bathtub curve, the segmented Weibull distribution is very appropriate. However, if the parameters in the segmented distribution are not appropriate, fðtÞ of the segmented distribution at t ¼ si may incur a discontinuity point. For example, when li ¼ si , hðtÞ and fðtÞ of the segmented Weibull distribution at t ¼ si are discontinued, but this discontinuity may not affect the measurement of reliability, because RðtÞ is continuous. As a whole, the segmented distribution in this paper can be fit to all kinds of failure mode recommended by RCM. 3. Estimation of segmented distribution

s1 6 t < s2 s2 6 t < s3 s3 6 t < 1

3.1. Bayesian regression estimation of segmented distribution

361

Usually, a piecewise-linear regression model is defined as Eq. (30). Namely, the response of any input xj 2 ½si ; siþ1 Þ contains two terms: the linear function ai þ bi xj and the stochastic disturbance ei that obeys normal distribution Nð0; r2i Þ. The aim of piecewise-linear regression is to infer parameters ai and bi under the circumstance that m input samples xj and their correspondingly output samples yj are known. Denoted by a = (ai)n1, b = (bi)n1, r2 = (r2i )n1, s = (si)n1, X = (xi)m1, Y = (yi)m1, under the condition that parameters a; b; r2 and s are given, then the maximum likelihood function of nsegment distribution is written as Eq. (31). In addition, the residual sum of squares Q in Eq. (32) can reflect precision of fitting.

362

y ¼ ai þ bi x þ ei 352 353 354 355 356 357 358 359 360

If the sample set S ¼ ft1 ; t2 ; . . . ; tk ; . . . ; tm g of segmented distribution is known, in order to estimate parameters in segmented distribution, we should, first of all, estimate change-points s1 ; s2 ; . . . ; sn and location parameters l1 ; l2 ; . . . ; ln .25 Usually, let li ¼ 0ði ¼ 1; 2; . . . ; nÞ or li ¼ si ði ¼ 1; 2; . . . ; nÞ. Of course, there are many other more accurate methods to estimate li and si . In this paper, Bayesian estimation and MLE are utilized to preliminarily estimate parameters in segmented distribution.26

Fig. 6

ð29Þ

ei 

Nð0; r2i Þ;

ð30Þ

fðX; Yja; b; r ; sÞ ¼ ð2pÞ

m2

n Y

364 365 366 367 368 369 370 371 372 373 374 375

x 2 ½si ; siþ1 Þ;

i ¼ 1; 2; . . . ; n 2

363

mi

 ðr2i Þ 2



e

1 2r2 i

P xj 2½si ;siþ1 Þ

377 378

ðyj ai bi xj Þ2

i¼1

ð31Þ Q¼

n X X

380 381

ðyj  ðai þ bi xj ÞÞ

2

ð32Þ

i¼1 xj 2½si ;siþ1 Þ

Some curves of three-segment Weibull distribution.

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

383

CJA 741 27 December 2016

No. of Pages 21

8

H. Li et al. n Y

384

3.1.1. Prior distribution of Bayesian regression

385

Usually, properly determining prior distribution is a key of Bayesian estimation. Although ai , bi , r2i , si in most Bayesian estimation are supposed to be independent27,28, in fact, there are close relations among ai , bi , r2i , si .

386 387 388

389 390 391 392 393 394

395 396 397 398

400 401 402

404 405 406 407 408

Theorem 1. In the ith interval ½si ; siþ1 Þ of the piecewise-linear regression model in Eq. (30), under the condition that r2i and s are given, parameters bi jr2i ; s obey normal distribution pffiffiffi 2 NðBi =A ði ¼ 1; 2; . . . ; nÞ, P where i =Ai Þ Þ Pi ; ð 2mi rP Ai ¼ xj 2½si ;siþ1 Þ xj 2½si ;siþ1 Þ ðxj1  xj2 Þ, Bi ¼ xj 2½si ;siþ1 Þ 1 2 1 P ðy  y Þ. j2 xj 2½si ;siþ1 Þ j1 2

Proof. From Eq. (30), if xj1 2 ½si ; siþ1 Þ; xj2 2 ½si ; siþ1 Þ, then their stochastic disturbances ej1 and ej2 will also obey normal distribution Nð0; r2i Þ, namely, (

ej1 ¼ yj1  ðai þ bi xj1 Þ  Nð0; r2i Þ ej2 ¼ yj2  ðai þ bi xj2 Þ  Nð0; r2i Þ

:

Thus, we get yj1  yj2  bi ðxj1  xj2 Þ pffiffiffi  Nð0; 1Þ: 2ri Furthermore, if any two samples are combined, then summation of all combinations will also obey normal distribution, i.e., X

yj1  yj2  bi ðxj1  xj2 Þ pffiffiffi  Nð0; m2i Þ: 2ri 2½si ;siþ1 Þ

xj1 2½si ;siþ1 Þ xj2

411

After deduction and simplification, it’s concluded that bi obeys normal distribution when r2i ; s are given, namely,

413 415

416 417 418 419 420

421 422 423 424 425 426 427

428 429 430 431 432 433

pffiffiffi 2 bi jr2i ; s  NðBi =Ai ; ð 2mi ri =Ai Þ Þ i ¼ 1; 2; . . . ; n

1 Ai pffiffiffiffiffiffiffiffi 2 e 2ri pða; b; r2 ; sÞ ¼ r 2p 2m i i i¼1

mi ðai Yi þbi Xi Þ þ12 2

434

A b B 2  i i mi

i

pðr2i jsÞpðsi Þ ð33Þ

436

Proof. From the topic supposes, as the parameters in different segments are independent of each other, we get

437

n Y pða; b; r2 ; sÞ ¼ pðsÞpða; b; r2 jsÞ ¼ pðsÞ pðai ; bi ; r2i jsÞ

438 439

i¼1

¼ pðsÞ

n Y

pðai jbi ; r2i ; sÞpðbi ; r2i jsÞ

i¼1

¼ pðsÞ

n Y

pðai jbi ; r2i ; sÞpðbi jr2i ; sÞpðr2i jsÞ

i¼1

¼

n Y

pðai jbi ; r2i ; sÞpðbi jr2i ; sÞpðr2i jsÞpðsi Þ:

i¼1

From Theorems 1 and 2, the prior distribution of the model in Eq. (30) is obtained as Eq. (33). Unlike other linear regression models, the above prior distribution in Theorem 3 can express relations among parameters ai ; bi ; r2i ; si . Usually, inverse Gamma distribution IGðui ; vi Þ or uniform distribution Uðui ; vi Þ or non-information distribution can be taken as a prior distribution of r2i js. If samples xj and yj have been correctly divided, namely s is known, the following Algorithm 1 can estimate hyper-parameters ui and vi of the ith segment. h

441 442 443 444 445 446 447 448 449 450 451 452 454

X

410

412



:

Theorem 2. In the ith interval ½si ; siþ1 Þ of the piecewise-linear regression model in Eq. (30), under the condition that bi , r2i , s are given, parameters ai jbi ; r2i ; sði ¼ 1; 2; . . . ; nÞ will obey normal P distribution NðYi  bi Xi ; r2i =mi Þ, where X ¼ m1i xj 2½si ;siþ1 Þ xj , P Y ¼ m1i xj 2½si ;siþ1 Þ yj . Proof. From Eq. (30), the stochastic disturbance ej of the jth sample xj and yj obeys normal distribution Nð0; r2i Þ, so it’s easy to know that ej ¼ yj  ðai þ bi xj Þ  Nð0; r2i Þ. Furthermore, we get Yi  ai  bi Xi  Nð0; r2i =mi Þ. bi ; r2i

Therefore, when and s are given, parameters ai will obviously obey normal distribution, denoted as ai jbi ; r2i ; s  NðYi  bi Xi ; r2i =mi Þi ¼ 1; 2; . . . ; n: h Theorem 3. Assuming that pðsi Þ is the prior distribution of si ; the distribution of r2i js is denoted as pðr2i jsÞ under the condition that s are given; and parameters ai , bi , r2i , si in the ith segment are independent from parameters aj , bj , r2j , sj in the jth segment, then the prior distribution of the piecewise-linear regression model in Eq. (30) can be written as

Algorithm 1

455

Step 1 Ascendingly sort all samples in the ith segment according to xj þ yj ðj ¼ 1; 2; . . . ; mi Þ, and obtain ðxð1Þ ; yð1Þ Þ; ðxð2Þ ; yð2Þ Þ; . . . ; ðxðjÞ ; yðjÞ Þ; . . . ; ðxðmi Þ ; yðmi Þ Þ ðj ¼ 1; 2; . . . ; mi Þ. Step 2 Determine empirical samples of r2i . For l ¼ 1 : l

456

 For m ¼ 1 : m y yðlÞ y yðlÞ Set alm ¼ yðlÞ  xðNmþ1Þ xðlÞ , blm ¼ xðNmþ1Þ . ðNmþ1Þ xðlÞ ðNmþ1Þ xðlÞ

Set elmj ¼ yj  ðalm þ blm xj Þ ðj ¼ 1; 2; . . . ; mi Þ. Calculate the variance of elmj denoted as rlm i , namely, Pmi lmj 1 PN lmj 2 lm 1 ri ¼ mi 1 j¼1 ðe  mi j¼1 e Þ . End End 2 Step 3 If rlm i can be taken as the samples of ri , then estimate hyper-parameters ui and vi by moment estimation.

457 458 459 460 461 462 463 464 465 466 467 468 469 470 471

 lines, where land m  are pre-set conAlgorithm 1 constructs lm stants. These lines are supposed as an approximation of optimal fitting, and then residuals caused by the samples relative to these lines can be regarded as an approximation of the stochastic disturbance ei . Therefore, rlm can be taken as an i approximation of r2i js, which can be utilized to estimate its hyper-parameters ui and vi of r2i js. The prior distribution of change-points si , denoted as pðsi Þ, is often considered as a discrete distribution defined in samples xi ði ¼ 1; 2; . . . ; mÞ.16 The greater the first- and second-order differences of ðxj ; yj Þ are, the higher the probability that change-points occur at xj will be, and thus it should be at xj

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

472 473 474 475 476 477 478 479 480 481 482 483

CJA 741 27 December 2016

No. of Pages 21

Study on segmented distribution 484 485

9

that samples are divided. In this paper, the flow of determination of pðsi Þ is shown as the following Algorithm 2.

486 488 489

Algorithm 2

494

Step 1 Ascendingly sort all samples ðxj ; yj Þ according to xj ðj ¼ 1; 2; . . . ; mÞ and get ðxð1Þ ; yð1Þ Þ; ðxð2Þ ; yð2Þ Þ; . . . ; ðxðjÞ ; yðjÞ Þ; . . . ; ðxðmÞ ; yðmÞ Þ ðj ¼ 1; 2; . . . ; mÞ. y yðj1Þ Step 2 Set @ 0j ¼ xðjþ1Þ ðj ¼ 2; 3; . . . ; m  1Þ, @ 01 ¼ 0, @ 0m ¼ 0, ðjþ1Þ xðj1Þ

495

and then normalize @ 0j , namely set p0j ¼ Pm j

490 491 492 493

j@ 0 j

496 497 498

499

500 501 502

503 504

Step 3 Set

@ 00j

¼

@ 0jþ1 @ 0j1 xðjþ1Þ xðj1Þ

and then normalize

@ 00j ,

j@ 0 j i¼1 i

ðj ¼ 2; 3; . . . ; m  1Þ,

namely set

p00j

. @ 001

j@ 00j j

¼ Pm

j@ 00 j i¼1 i

¼ 0,

@ 00m

Eq. (35) merely includes hyper-parameters ui , vi ði ¼ 1; 2; . . . ; nÞ; however, the posteriori distribution of parameters is difficult to be explicitly expressed. Therefore, only the Gibbs sampling method can find out the numerical solution.

525 526 527

Theorem 4. In the joint distribution in Eq. (35), if parameters al , b, r2 , s and samples X; Y are given, then the full conditional distribution of parameters ai ; i–l will obey normal distribution, namely,

528



r2 ai jal ; b; r ; s; X; Y  N Yi  bi Xi ; i ; 2mi

532

2

¼ 0,

523 524

l–i:

529 530 531

534

.

Proof. if al ; b; r2 ; s; X; Y; l–i in Eq. (35) are fixed, then we obtain

Step 4 Calculate wj . When m is an even number, set ( 4ðj1Þ j ¼ 1; 2; . . . ; m2 wj ¼ mðm2Þ ; when m is an odd 4ðmjÞ m m mðm2Þ j ¼ 2 þ 1; 2 þ 2; . . . ; m 8 < 4ðj1Þ2 j ¼ 1; 2; . . . ; mþ1 2 . number, set wj ¼ ðm1Þ : 4ðmjÞ2 j ¼ mþ3 ; mþ5 ; . . . ; m 2 2 ðm1Þ  m   m Step 5 Set lpi ¼ max 1; ði  2Þ n , rpi ¼ min m; i n , and get the prior distribution 8 < Prwj ðp0j þp00j Þ j ¼ lpi ; lpi þ 1; . . . ; rpi pi w ðp0 þp00i Þ pðsi ¼ xj Þ ¼ i¼lpi i i : 0 Others j ¼ 1; 2; . . . ; m; i ¼ 2; 3; . . . ; n.

0 

fðai jal ; b; r ; s; X; YÞ / e 2

1 2r2 i

1

X

@

ðyj ai bi xj Þ2 þmi ðai Yi þbi Xi Þ

2

537

:

Assuming that Ci ¼ m1i xj 2½si ;siþ1 Þ xj yj , Di ¼ m1i P Ei ¼ m1i xj 2½si ;siþ1 Þ y2j , then we get

536

A

xj 2½si ;siþ1 Þ

P

535

539

P

2 xj 2½si ;siþ1 Þ xj ,

540 541 542

mi X

ðyj  ðai þ bi xj ÞÞ2 ¼ mi a2i  2ai mi Yi þ 2ai bi mi Xi  2bi mi Ci

j¼1

þ b2i mi Di þ mi Ei :

505 506 507 508 509 510 511

513

In Algorithm 2, Step 4 ensures that the middle samples are more likely to be regarded as si , which is inclined to make all segmentation as equalized as possible. In addition, pðs1 Þ and pðsnþ1 Þ are deemed as constants, for example, its distribution may be often written as the following formula.   1 s1 ¼ 0 1 s1 ¼ minðxi Þ pðs1 Þ ¼ or pðs1 Þ ¼ : 0 Others 0 Others

Pmi 2 2 If j¼1 ðyj  ðai þ bi xj ÞÞ þ mi ðai  Yi þ bi Xi Þ is transformed into a perfect square expression about ai , then fðai jal ;b;r2 ;s;X;YÞ / e

 12 2r i



2mi ðai ðYi bi Xi ÞÞ

2

515

If

r2i js

obeys inverse Gamma distribution IGðui ; vi Þ, namely, v

516 517 518 519 520

ðu þ1Þ 

ui

vi 2

pðr2i jsÞ ¼ Cðui i Þ ðr2i Þ i e ri , from Eqs. (31) and (33), it can be seen that the prior distribution and joint distribution of the piecewise-linear regression model are written as respectively Eqs. (34) and (35). n 1 Y Ai vui i pðsi Þ ðui þ2Þ 2ri pffiffiffiffiffiffiffiffi pða; b; r ; sÞ ¼ e ðr2i Þ 2p 2m Cðu Þ i i i¼1 2



mi ðai Yi þbi Xi Þ þ12 2

A b B 2 i i mi

i



þ2vi

ð34Þ

522

Þ2mi ðYi bi XÞ

2



is obtained. Furthermore, we get ð

ai  Yi bi Xi ri

pffiffiffiffi

fðai jal ; b; r2 ; s; X; YÞ / e

3.1.2. Bayesian piecewise-linear regression estimation

ð

þmi 2bi Xi Yi 2bi Ci þb2i X2i þb2i Di þEi þY2i

545 546 547 549

12

514

544

Þ

2mi

550 551

!2 :

553

Obviously, it is a kernel of normal distribution. Thus we get

r2i 2 : ai jal ; b; r ; s; X; Y  N Yi  bi Xi ; 2mi

554 555



557

Theorem 5. In the joint distribution in Eq. (35), if parameters a, bl , r2 , s and samples X, Y are given, then the full conditional distribution of parameters bi , i–l will obey normal distribution, namely,

558 559 560 561 562

0 n Y



m Ai vui i pðsi Þ ð i þui þ2Þ ðr2i Þ 2 e fðX; Y; a; b; r ; sÞ ¼ m pffiffiffiffiffiffiffiffi 2 2mi Cðui Þ i¼1 ð2pÞ

2

1 2r2 i

@

X xj 2½si ;siþ1 Þ

ðyj ai bi xj Þ2 þmi ðai Yi þbi Xi Þ þ12 2

A b B 2 i i mi

i

1 þ2vi A

ð35Þ

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

CJA 741 27 December 2016

No. of Pages 21

10

564

H. Li et al.

bi ja; bl ; r2 ; s; X; Y  N 3

2mi Ci þ 2m3i Xi Yi þ Ai Bi  4m3i ai Xi 2m2i r2i ; ; 2m3i Di þ 2m3i X2i þ A2i 2m3i Di þ 2m3i X2i þ A2i

bi ja; bl ; r2 ; s; X; Y  N

565

where

566

Ei ¼

1 mi

P

Di ¼ m1i

xj 2½si ;siþ1 Þ xj yj ,

586 587 588

3

2mi Ci þ 2m3i Xi Yi þ Ai Bi  4m3i ai Xi 2m2i r2i : ; 2m3i Di þ 2m3i X2i þ A2i 2m3i Di þ 2m3i X2i þ A2i

P

Ci ¼ m1i

is obtained. Obviously, it is also a kernel of normal distribution. Therefore, we get

P

2 xj 2½si ;siþ1 Þ xj ,

2 xj 2½si ;siþ1 Þ yj .



Theorem 6. In the joint distribution in Eq. (35), if parameters a, b, r2l , s and samples X, Y are given, then the full conditional distribution of parameters r2i , i–l will obey inverse Gamma distribution, namely,

589 590 591 592 593

 0 1  2  2 3 2 þ ðAi bi  Bi Þ2 þ 4m2i vi 2m a  2a Y þ 2a b X  2b C þ b D þ E þ a  Y þ b X i i i i i i i i i i i i i i i i m i 2 2 A: ri ja; b; rl ; s; X; Y  IG@ þ ui þ 1; 2 4m2i

599 567 568

Proof. Supposing that X, Y, a, bl , r2 , s, l–i in Eq. (35) are constants, then we get

569

0 

1 2r2 i

@

X

ðyj ai bi xj Þ2 þmi ðai Yi þbi Xi Þ þ12 2

i

576

After derivation and simplification, we get

A2 mi Di þmi X2i þ i2 2m i

1 2r2 i

0

A :

½fðbi ja;bl ;r2 ;s;X;YÞ / e



i i mi

xj 2½si ;siþ1 Þ

575 574 573 572 571

577 578

A b B 2

1

b2i þ

fðbi ja;bl ;r2 ;s;X;YÞ / e

581

In addition, because it’s easy to know

fðr2 ja;b;r2l ;s;X;YÞ / ðr2i Þ

ð 2i þui þ2Þ



1 2r2 i

m

X

@

e

ðyj ai bi xj Þ þmi ðai Yi þbi Xi Þ þ12 2

2

A b B 2 i i mi

xj 2½si ;siþ1 Þ

i

601 602 603

1 þ2vi A

:

605 606



580

600 Proof. Supposing that a, b, r2l , s, X, Y in Eq. (35) are constants, then we get

607



4ai mi Xi 2mi Ci 2mi Xi Yi 

Ai Bi bi m2 i

bi

608

After further derivation and simplification, we get

:

609 610





A2 Ai B i mi Di þ mi X2i þ i 2 b2i þ 4ai mi Xi  2mi Ci  2mi Xi Yi  2 bi 2mi mi 00 12 0 12 1 2 mi Ci þ mi Xi Yi þ A2mi B2i  2ai mi Xi mi Ci þ mi Xi Yi þ A2mi B2i  2ai mi Xi A i i A @ A A; ¼ mi Di þ mi X2i þ i 2 @@bi  A2 A2 2mi mi Di þ mi X2i þ 2mi 2 mi Di þ mi X2 þ 2mi 2 i

582 583

0 B

12@

585

fðbi ja; bl ; r2 ; s; X; YÞ / e

bi 

2m3 Ci þ2m3 Xi Yi þAi Bi 4m3 ai Xi i i i 2m3 Di þ2m3 X2 þA2 i p i i i 2mi ri

ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 3 2 2

12 C A

i

fðr

2

ja; b; r2l ; s; X; YÞ

/

ð þui þ2Þ ðr2i Þ 2 e mi



2m3 i



ð

a2 2ai Yi þ2ai bi Xi 2bi Ci þb2 Di þEi þ ai Yi þbi Xi i i 4m2 r2 i i

Þ

2



þðAi bi Bi Þ2 þ4m2 vi i

Therefore, it’s tenable that

2m Di þ2m X þA i i i i

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

:

612 613 614

CJA 741 27 December 2016

No. of Pages 21

Study on segmented distribution Table 1

11

Results of Bayesian estimation in simulation.

Prior distribution ai jbi ; r2i ; s

bi jr2i ; s

NðYi  bi Xi ; mii Þ

 pffiffi 2 N ABii ; ð 2Ami i ri Þ  pffiffi 2  N ABii ; 2Ami i ri

Uðuai ; vai Þ

Uðubi ; vbi Þ

r2

NðYi  bi Xi ; mii Þ r2

EðkÞ

EðtÞ

EðQÞ

IGðui ; vi Þ

455.14

0.098106

0.559669

Uðui ; vi Þ

476.13

0.098392

0.559673

Uðui ; vi Þ

672.85

0.098240

0.560546

r2i js

Fig. 7

Illustration of a sampling result.

 0 1  2  þ ðAi bi  Bi Þ2 þ 4m2i vi 2m3i a2i  2ai Yi þ 2ai bi Xi  2bi Ci þ b2i Di þ Ei þ ai  Yi þ bi Xi m i A: r2 ja; b; r2l ; s; X; Y  IG@ þ ui þ 1; 2 4m2i

618 619 620 621

Theorem 7. In the joint distribution in Eq. (35), if parameters a, b, r2 , sl and samples X, Y are given, then the full conditional distribution of si ja; b; r2 ; sl ; X; Y about parameters si , i–l is same as the prior distribution pðsi Þ, namely,

fðsi ja; b; r2 ; sl ; X; YÞ / pðsi Þ:

fðsi ja; b; r2 ; sl ; X; YÞ ¼ pðsi Þ:

633

fðsi ja; b; r2 ; sl ; X; YÞ ¼ pðsi Þ:

According to the above full conditional distribution, if the model in Eq. (30) is inferred by Bayesian estimation, then the Gibbs sampling method is shown as the following Algorithm 3. h

634

622 624

625 626



Proof. Similarly, if parameters a, b, r , sl , X, Y in Eq. (35) are deemed as constants, then we obtain 2

Because pðsi Þ is a discrete distribution, we get

629 630 631

627

635 636 637 638

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

CJA 741 27 December 2016

No. of Pages 21

12

H. Li et al.

640 641

Algorithm 3

642

Step 1 Infer pðsi Þ by Algorithm 2, and set k ¼ 0. Step 2 According to Theorem 7, generate n1 discrete random numbers obeying the full conditional distribution

643 644

ðkÞ

ðkÞ

ðkÞ ; sl ; X; YÞ,

645

fðsi ja ; b ; r

646

ði ¼ 1; 2; . . . ; nÞ, which

647 648 649 650 651 652 653 654 655

2ðkÞ

ðkþ1Þ

ðkþ1Þ ðkþ1Þ and satisfy si 6 siþ1 ðkþ1Þ can be taken as si ði ¼ 2; 3; . . . ; nÞ. ðkþ1Þ 0, snþ1 ¼ maxðxi Þ, and then sðkþ1Þ is

Moreover, set s1 ¼ obtained. Step 3 The prior distribution of r2i ði ¼ 1; 2; . . . ; nÞ is estimated according to Algorithm 1. Step 4 Compute Ai , Bi , Ci , Di , Ei , Xi , Yi , mi , li ði ¼ 1; 2; . . . ; nÞ. ð0Þ

ð0Þ

2ð0Þ

Step 5 Determine initial values ai , bi , ri ði ¼ 1; 2; . . . ; nÞ, and set s ¼ 0. Step 6 According to Theorem 4, generate n normal random numbers obeying the full conditional distribution ðsÞ

659

fðai jal ; bðsÞ ; r2ðsÞ ; sðkþ1Þ ; X; YÞ ði ¼ 1; 2; . . . ; nÞ, which are taken as aðsþ1Þ . Step 7 According to Theorem 5, generate n normal random numbers obeying the full conditional distribution

660

fðbi jaðsþ1Þ ; bl ; r2ðsÞ ; sðkþ1Þ ; X; YÞ ði ¼ 1; 2; . . . ; nÞ, which are taken

656 657 658

661 662 663 664 665 666 667 668 669 670 671 672 673 674

ðsÞ

as bðsþ1Þ . Step 8 According to Theorem 6, generate n inverse Gamma random numbers obeying the full conditional distribution

fðr2i jaðsþ1Þ ; bðsþ1Þ ; rl ; sðkþ1Þ ; X; YÞ ði ¼ 1; 2; . . . ; nÞ, which are taken as r2ðsþ1Þ . Step 9 If aðsþ1Þ , bðsþ1Þ , r2ðsþ1Þ tend to converge to a standstill, then set aðkþ1Þ ¼ aðsþ1Þ , bðkþ1Þ ¼ bðsþ1Þ , r2ðkþ1Þ ¼ r2ðsþ1Þ and go to Step 10; otherwise, set s ¼ s þ 1 and go to Step 6. Step 10 If aðkþ1Þ , bðkþ1Þ , r2ðkþ1Þ , sðkþ1Þ satisfy the termination condition, then go to Step 11; otherwise, set k ¼ k þ 1 and go to Step 2. Step 11 Calculate expectations of random sequences aðkÞ , bðkÞ , r2ðkÞ , sðkÞ , denoted as a, b, r 2 , s, which are taken as estimation 2 ^ values ^ a, b, r ^ , ^s. 2ðsÞ

675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699

Assuming that parameters and their estimated values are ^2 Þ, then respectively denoted as h ¼ ða; b; r2 Þ and ^ h ¼ ð^ a; ^ b; r t ¼ k^ h  hk=khk can reflect estimation accuracy. Q in Eq. (32) can measure fitting precision. In addition, iterations k can imply algorithm efficiency. Therefore, the averages of k, t, Q from all simulating calculations, denoted as EðkÞ, EðtÞ, EðQÞ, can indicate the effectiveness and efficiency of Algorithm 3. From Table 1, the effectiveness and efficiency of the algorithm in this paper are satisfying, because the relations among parameters are fully considered.

700

3.1.3. Bayesian regression estimation of segmented distribution

710

In this paper, estimation of segmented distribution by Bayesian regression is based on the following hypotheses: (1) li are given or known and (2) parameters ai , bi , r2i , si in the ith segment are independent of parameters aj , bj , r2j , sj in the jth segment. Usually, these hypotheses can be met in most practical cases. Assuming that failure rate is linear, shown as Eq. (8), according to Eq. (36), Eq. (8) can be transformed into a linear relationship y ¼ ai þ bi x, which can meet the requirements of regression. If the estimation values of ai , bi with the Bayesian ^i ¼ b^i are regression model are denoted as a^i , b^i , then ^ ai ¼ a^i , b achieved.

711

yj ¼ hðtj Þ;

xj ¼ tj  li ;

tj 2 ½si ; siþ1 Þ;

ai ¼ ai ;

¼ bi

ð36Þ

Under the circumstance that samples follow segmented Weibull distribution, if we transform Eq. (20) according to Eq. (37), then the failure rate of samples in the ith interval can be expressed as Eq. (30). Therefore, after all samples are classified, if li , si are known and estimation values a^i , b^i have already been obtained by Bayesian estimation, moreover we  ^ ^ 1 ai bi þ1 ^ i ¼ b^i þ 1, ^ get the estimate values of ai , bi , i.e., b ai ¼ b^eþ1 . yj ¼ ln hi ðtj Þ; ¼

xj ¼ lnðtj  li Þ;

lnðabi i bi Þ;

tj 2 ½si ; siþ1 Þ;

702 703 704 705 706 707 708 709

712 713 714 715 716 717 718 719 720 721 722 723

bi

i

In Algorithm 3, there are many methods to estimate parameters li . In brief, we can set li ¼ 0 or li ¼ si . If n ¼ 1, Step 2 can be skipped, what’s more, Algorithm 3 only needs one loop. In Step 5, the expectations of ai , bi , r2i , which are worked out according to their full conditional distribution, can be ð0Þ ð0Þ 2ð0Þ regarded as initial values ai , bi , ri . In addition, if the dis2 tribution of ri js is a uniform or non-information distribution, then we can get a conclusion similar to Theorems 4–6. Under the same parameters circumstance, compute 1024 groups of simulation sample set. The average results are compared with those of other methods so as to verify the effects of the algorithms in this paper. In simulation verification, set n ¼ 1, a ¼ 2, b ¼ 2, r2 ¼ 0:1, s ¼ 0, l ¼ 0 in Eq. (30). Then set m ¼ 16, namely, every 16 pairs random number ðxj ; yj Þ compose one group sample set. Simulation results are listed in Table 1. The results of non-information and uniform prior distributions are also tabulated in Table 1. Fig. 7 is a scatter diagram illustrating all sampling samples in a certain simulating calculation, which shows that the sampling samples are concentratedly distributed in the neighborhood of the true values of parameters. Fig. 8 illustrates sampling samples in course of iteration, which implies the algorithm begins to converge only at a cost of small quantity of iteration.

701

725 726 727 728 729 730 731 732 733

ai

bi ¼ ðbi  1Þ

ð37Þ

Thus, the failure rate hðtj Þ of sample tj is inferred firstly, and then ðtj ; hðtj ÞÞ is transformed into sample ðxj ; yj Þ according to Eqs. (36) and (37). Next, regression coefficients can be estimated by Algorithm 3. In the end, the regression coefficients are transformed into parameters of segmented distribution.

735 736 737 738 739 740

3.2. MLE of segmented distribution

741

The unknown parameters in the ith interval are written as hi ¼ ½ai bi li T . Supposing that n, si ði ¼ 1; 2; . . . ; nÞ are known, then the likelihood function LðhÞ of n-segment distribution is written as Eq. (38) according to Eq. (6). Taking the logarithm of Eq. (38), we get its logarithmic likelihood function ln L shown as Eq. (39). We take the partial derivative of Eq. (39) with respect to hi , which is the parameter in the ith interval, and get Eq. (40).

742

!

n1 n n X X X X LðhÞ ¼ exp  Jj Ni  Ki ðtk Þ j¼1

i¼jþ1

i¼1 si 6tk
! n Y Y

743 744 745 746 747 748 749 750

hi ðtk Þ

i¼1 si 6tk
ð38Þ

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

752

CJA 741 27 December 2016

No. of Pages 21

Study on segmented distribution

13

Fig. 8

!

753

ln L ¼ 

n1 X j¼1

þ 755 756

758

Jj

n X

Ni

i¼jþ1

n X X



n X X

Parameters changing with iteration in course of sampling.

Ki ðtk Þ

i¼1 si 6tk
ln hi ðtk Þ

ð39Þ

i¼1 si 6tk
ð40Þ

759

3.2.1. Test change-points by improved information criterion

760

Usually, only when the sample number of each interval is approximately equal or the difference among the lengths of different intervals is slight, partitioning is most reasonable.29,30 In addition, the lower the number of interval is, the more concise the model will be.31 In order to measure the appropriateness of segmented distribution, learning from the idea that the information criterion SIC penalizes the complexity of a model, the information criterion SIC is improved. Thus a novel information criterion ICSD appropriated for segmented distribution is proposed in this paper, which is shown as Eq. (43).32,33 In Eq. (43), ln LðhÞ means the logarithm likelihood function shown as Eq. (39); h are the unknown parameters to be estimated; and dh is the dimensions of unknown parameters h; m is the total number of samples; n is the number of intervals, i.e., the segmentation number; cpenalty is the penalty

761 762 763 764 765 766 767 768 769 770 771 772 773 774

item expressed by the standard deviation, which is calculated by Eq. (41) or (42). In Eq. (41), mi indicates the sample number of the ith interval; snþ1 in Eq. (42) is usually supposed as the maximum sample, and of course, there are other methods to estimate snþ1 . In the special case of n ¼ 1, Eq. (43) is exactly the same as SIC information criterion.34 Therefore, when n and si are given, the MLE of the n-segment distribution can be taken as an optimization problem, which is shown as Eq. (44). sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

2 1 Xn 1 Xn n cpenalty ¼ m  m ; snþ1 ¼ maxðti Þ ð41Þ i j i¼1 j¼1 i¼1 n n

cpenalty

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2ffi u n u1 Xn X 1 ¼t ðsiþ1  si Þ  ðsjþ1  sj Þ ; i¼1 n n j¼1

775 776 777 778 779 780 781 782 783 784

786 787

snþ1

n

¼ maxðti Þ i¼1

ð42Þ

789 790

ICSDðhÞ ¼ 2 ln LðhÞ þ ðdh þ ðn  1Þ2 þ cpenalty Þ ln m

ð43Þ

792 793

min ICSDsi ;n ðhÞ ¼ min 2 ln LðhÞ þ ðdh þ ðn  1Þ2 þ cpenalty Þ ln m

ð44Þ

Eq. (44) can also be utilized to test whether change-points exist or not. Namely, min ICSDsi ;n ðhÞ < SIC indicates that there obviously exist some change-points and estimating by the seg-

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

795 796 797 798

CJA 741 27 December 2016

No. of Pages 21

14

H. Li et al.

805

mented distribution is appropriate. Let Eq. (40) be 0, we can get likelihood equations. Although the solution of likelihood equations can be regarded as an estimation of unknown parameters, usually we may get no solution at all or multiple solutions due to the too high complexity and dimensions of the likelihood equations, especially in the case of the estimation of li .

806

3.2.2. MLE of segmented linear failure rate model

807

When failure rate is expressed by Eq. (8), the derivatives of Eqs. (8)–(10) are substituted into Eq. (40), and we can get the likelihood equation with respect to ai ; bi in the ith interval of the linear failure rate model, shown as Eq. (45). As the likeln L lihood equation @@l ¼ 0 is an inconsistent equation, only i other methods are utilized to estimate li . In this paper, we set li ¼ minðtk Þ.

799 800 801 802 803 804

808 809 810 811 812 813 814

816

8 n X X X > 1 > > ðsiþ1  si Þ mj þ ðtk  si Þ  ¼0 > ai þbi ðtk li Þ > > si 6tk > > > > > > > < n X 1 X   1 mj þ ðtk  si Þ 12 tk þ 12 si  li > ðsiþ1  si Þ 2 siþ1 þ 2 si  li > si 6tk j¼iþ1 > > > > > > > X > > tk li > ¼0 > ai þbi ðtk li Þ : si 6tk
Fig. 9

In particular, when n = 2, the likelihood equation about a1 ; b1 ; a2 ; b2 of the two-segment linear failure rate model is expressed as Eq. (46). X X 8 1 > ðs2  s1 Þm2 þ ðtk  s1 Þ  ¼0 > a1 þb1 ðtk l1 Þ > > s1 6tk > > > > > >   > > > l1 s2 þ 12 s22 þ l1 s1  12 s21 m2 > > > > > > > X X   > > tk l1 <þ l1 tk þ 12 t2k þ l1 s1  12 s21  ¼0 a1 þb1 ðtk l1 Þ s 6t
819 820

s 6t
1 2 1 2 k k > > > > > X X > > 1 > > ðtk  s2 Þ  ¼0 > a2 þb2 ðtk l2 Þ > > s2 6tk s2 6tk > > > > > > >  X tk l2 > > X > l2 tk þ 12 t2k þ l2 s2  12 s22  ¼0 > a2 þb2 ðtk l2 Þ :

s2 6tk

s2 6tk

ð46Þ ð45Þ

817 818

822

3.2.3. MLE of segmented Weibull distribution

823

For n-segment Weibull distribution, we can substitute the derivative of Eqs. (20)–(22) into Eq. (40), and thus we can get the likelihood equation with respect to ai ; bi ; li in the ith interval of segmented Weibull distribution, shown as Eq. (47).

824

Some distribution curves of bathtub mode.

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

825 826 827 828

CJA 741 27 December 2016

No. of Pages 21

Study on segmented distribution Table 2

15

Verification results of segmented distribution in bathtub model.

Distribution

n

Bayesian estimation E(ICSD)

E(ln L)

E(Q)

E(ICSD)

E(ln L)

Linear failure rate model

1 2 3

0.763215 0.013991 0.005517

712.025561 1272.742802 1237.180250

351.368390 626.278686 598.496703

0.771646 0.017703 0.009220

586.574952 1136.994683 1236.104312

288.776617 558.576941 597.986081

Segmented Weibull distribution

1 2 3

50.175355 15.869131 2.610704

1210.214449 1227.338421 1251.915052

600.103278 601.092377 600.921539

78.316891 43.976807 48.027381

1082.648733 1113.451129 1137.161508

536.320420 544.134123 543.476008

8 n X > b b > mj > ððsiþ1  li Þ i  ðsi  li Þ i Þ > > > j¼iþ1 > > X > > i > þ ðtk  li Þbi  ðab þ ðsi  li Þbi Þmi ¼ 0 > i > > > si 6tk > > > n X > > > > ððsiþ1  li Þbi lnðai ðsiþ1  li ÞÞ  ðsi  li Þbi lnðai ðsi  li ÞÞÞ mj > > > j¼iþ1 < X i þ ððtk  li Þbi  ab i Þ lnðai ðtk  li ÞÞ > > > si 6tk > > > > > mi ðsi  li Þbi lnðai ðsi  li ÞÞ  bmi i ¼ 0 > ai bi > > > > n > X X > bi 1 bi 1 > > ððs  l Þ  ðs  l Þ Þ mj þ ððtk  li Þbi 1 > iþ1 i i i > > > s 6t > > > :  b bi 1 Þ  mi ðsi  li Þbi 1 ¼ 0 i ai bi ðtk li Þ

ð47Þ

830 831 832 833

In particular, if n = 2, we get the likelihood equation of two-segment Weibull distribution as Eq. (48). 8 X b1 b1 > ððs ðtk  l1 Þb1 2  l1 Þ  ðs1  l1 Þ Þm2 þ > > > s1 6tk > > > > b1 b1 > ða1 þ ðs1  l1 Þ Þm1 ¼ 0 > > > > > > ððs2  l1 Þb1 lnða1 ðs2  l1 ÞÞ  ðs1  l1 Þb1 lnða1 ðs1  l1 ÞÞÞm2 > > > X > > 1 > þ ððtk  l1 Þb1  ab > 1 Þ lnða1 ðtk  l1 ÞÞ > > > s1 6tk > > > > > m1 ðs1  l1 Þb1 lnða1 ðs1  l1 ÞÞ  bm1 1 ¼ 0 > a1 b1 > > > > > b1 1 b1 1 >  ðs1  l1 Þ Þm2 < ððs2  l1 Þ

X b1 1 b1 1 > >  m1 ðs1  l1 Þb1 1 ¼ 0 þ ðt  k  l1 Þ b1 > > ðtk l1 Þa1 b1 > s 6t 1 2 k > > > >X 2 > ðtk  l2 Þb2  m2 ðab þ ðs2  l2 Þb2 Þ ¼ 0 > 2 > > > s 6t 2 k > > > >X > 2 > ððtk  l2 Þb2  ab > 2 Þ lnða2 ðtk  l2 ÞÞ > > s2 6tk > > > > > > m2 ðs2  l2 Þb2 lnða2 ðs2  l2 ÞÞ  bm2 2 ¼ 0 > > a2 b2 > > > X > b2 1 b 1 > 2 > ððtk  l2 Þ  Þ  m2 ðs2  l2 Þb2 1 ¼ 0 > b : ðtk l2 Þa2 2 b2 s2 6tk

835

836

837 838

MLE

E(Q)

and MLE, the detailed steps of the method to estimate parameters are shown as the following Algorithm 4.

839 840 841

Algorithm 4 Step 1 Ascendingly sort lifetime samples to form sample set S ¼ ft1 ; t2 ; . . . ; tk ; . . . ; tm g and then estimate the failure rate hðtk Þ of the kth sample tk . Step 2 Determine the feasible region of si ði ¼ 1; 2; . . . ; nÞ and initial values of parameters. There are usually two methods to estimate the feasible region of si . (1) The parameters of segmented distribution are estimated by Algorithm 3. If the results of Bayesian estimation are ^si ¼ t^i , then set Wi ¼ ft^in0 ; . . . ; t^i1 ; t^i ; t^iþ1 ; . . . ; t^iþn0 g, namely, there are n0 samples around t^i which are regarded as the feasible region of si . (2) ðtk ; hðtk ÞÞ ðk ¼ 1; 2; . . . ; mÞ are categorized by a certain clustering method. If the best number of categories is n, then the distribution is also divided into n intervals. If the center of the ith class is denoted as ci , then we let Wi ¼ ft; t 2 S ^ ci1 < t < ci g ði ¼ 2; 3; . . . ; nÞ. In addition, we usually set W1 ¼ f0g or W1 ¼ fmint2S ðtÞg. We take Wi as the feasible solutions of si , namely si 2 Wi ði ¼ 1; 2; . . . ; nÞ. Step 3 Estimate the distribution parameters by MLE. Step 3.1 Initialize ICSDbest , abest , bbest , sbest , lbest . Step 3.2 Choose a group si ði ¼ 1; 2; . . . ; nÞ from Wi ði ¼ 1; 2; . . . ; nÞ. Step 3.3 The sample set S is divided into n intervals according to si ði ¼ 1; 2; . . . ; nÞ. Step 3.4 Solve the optimization problem in Eq. (44) or the likelihood equation in Eq. (45) or (47), and the result is denoted as ICSDsi ;n . If ICSDsi ;n < ICSDbest , set ICSDbest ¼ ICSDs;n , abest ¼ a, bbest ¼ b, sbest ¼ s, lbest ¼ l. Step 3.5 Go to Step 3.2 and compute repeatedly until all si ði ¼ 1; 2; . . . ; nÞ have been searched. Step 4 Output the solutions abest , bbest , sbest , lbest . 843

Steps 3.2–3.5 construct a loop, which can search optimal si ði ¼ 1; 2; . . . ; nÞ from their feasible region Wi . In order to estimate the parameters of segmented distribution, Algorithm 4 will solve the optimization problem or its likelihood equation Q for ni¼1 jWi j times.

844 845 846 847 848

4. Verification by simulation

849

Because the failure rate of a bathtub curve changes in an obvious phased manner over the entire lifetime cycle, a bathtub curve is suitable to test the segmented distribution and its estimation method.35,36 In order to test the performance of segmented distribution, a new bathtub curve model is put forward in this paper. Unlike what Lai et al. 37, 38 proposed,

850

ð48Þ 3.3. Flow of parameters estimation of segmented distribution If the sample set of segmented distribution is S ¼ ft1 ; t2 ; . . . ; tk ; . . . ; tm g, according to Bayesian estimation

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

851 852 853 854 855

CJA 741 27 December 2016

No. of Pages 21

16 Table 3

H. Li et al. Lifetime sample of proximity sensors in leading-edge flap.

Sample (day)

22

35

43

45

47

69

89

112

211

213

231

254

269

271

287

304

391

Frequency

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

Fig. 10

Table 4

Empirical reliability of proximity sensors in leading-edge flap estimated by product limit method.

Results of proximity sensors in leading-edge flap when li ¼ si .

Model

Bayesian estimation

n

MLE

Q

ICSD

Q

ICSD

Linear failure rate model

1 2 3

0.385113 0.391301 0.050524

998.322910 1495.454200 983.191595

0.461216 0.283457 0.603541

742.080213 1234.126147 1156.456345

Segmented Weibull distribution

1 2 3

3.013787 3.401046 1.392337

834.545652 1206.334961 929.041437

30.792925 23.915894 15.149311

329.938043 733.151436 554.069265

Table 5

Distribution parameters of proximity sensors in leading-edge flap when li ¼ si .

Method

Model

n

[a1, a2, a3]

[b1, b2, b3]

[s1, s2, s3]

[l1, l2, l3]

Bayesian estimation

Linear failure rate model

1 2 3 1 2 3

0.0003 0.0456, 0.0023 0.0432, 0.0830, 0.0118 0.0309 0.0431, 0.1000 0.0445, 0.0887, 0.1056

0.0015 0.0006, 0.0025 0.0006, 0.0004, 0.0049 1.7219 1.3331, 1.2538 1.2972, 1.0384, 1.3957

22 0, 89 0, 69, 211 0 0, 112 0, 69, 211

21.78 0, 88.11 0, 68.31, 208.89 0 0, 110.88 0, 68.31, 208.89

1 2 3 1 2 3

0.0004 0.0470, 0.0046 0.0455, 0.0664, 0.0095 0.0206 0.0287, 0.0667 0.0297, 0.0591, 0.0704

0.001 0.0003, 0.0022 0.0004, 0.0003, 0.0033 1.3776 1.0665, 1.0031 1.0378, 0.8307, 1.1165

22 22, 112 22, 69, 211 0 0, 112 0, 69, 211

10.89 21.78, 92.57 21.50, 36.58, 107.05 0 0, 55.02 0, 33.87, 104.43

Segmented Weibull distribution

MLE

Linear failure rate model

Segmented Weibull distribution

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

CJA 741 27 December 2016

No. of Pages 21

Study on segmented distribution Table 6

17

Results of proximity sensors in leading-edge flap when li ¼ 0.

Model

n

Bayesian Estimation

MLE

Q

ICSD

LnL

Q

ICSD

LnL

Linear failure rate model

1 2 3

0.023444 0.004479 0.000207

524.299062 826.590032 804.948647

259.441481 238.972833 250.450112

0.651190 0.493083 0.183909

680.265902 667.489418 765.032881

337.299737 337.507106 363.828470

Segmented Weibull distribution

1 2 3

0.736003 0.078336 0.068887

524.709035 818.894616 714.654662

259.646467 235.125125 237.598933

10.156623 6.117920 1.393693

509.523560 646.610284 842.849242

251.928567 245.101964 299.779892

Table 7

Distribution parameters of proximity sensors in leading-edge flap when li ¼ 0.

Method

Model

n

a1

Bayesian estimation

Linear failure rate model

1 2 3 1 2 3

0.032900 0.058192 0.049372 0.035840 0.043426 0.049048

1 2 3 1 2 3

0.029909 0.053672 0.043542 0.032581 0.039478 0.046568

Segmented Weibull distribution MLE

Linear failure rate model

Segmented Weibull distribution

a2

a3

b1

0.634507 0.171940

1.161458

0.006345 0.045332

0.006345

0.483327 0.018532

1.161420

0.005768 0.037705

0.007923

b2

0.000658 0.000299 0.000467 1.532463 1.327502 1.177229 0.000598 0.000257 0.000455 1.393148 1.206820 1.027791

b3

0.003251 0.001363

0.005208

4.797477 1.306900

4.797477

0.002298 0.000510

0.004735

4.361343 1.243400

3.675273

s1 0 0 0 0 0 0 0 0 0 0 0 0

s2

s3

231 211

271

231 45

231

213 112

271

231 43

211 885

856 857 858 859

the failure rate of the new bathtub curve model in this paper is not two additive Weibull distributions, but a sum of one increasing exponential function and another decreasing exponential function.

t

hðtÞ ¼ c1 ed1 t þ c2 ed2 > 0;

c1 > 0; c2 > 0; d1 > 0; d2

c1 d1 d2 > 1; t 2 ½0; 1Þ c2

RðtÞ ¼ e

4.1. A novel bathtub curve

861

A recent bathtub curve model proposed by Jiang 39 is proven to achieve a good effect, however, it is finite support. Learning from the thought of Jiang, a novel bathtub curve model is proposed in this paper, whose definition domain is ½0; þ1Þ, so as to verify the segmented distribution herein. If the failure rate decreases, h0 ðtÞ < 0, hðtÞ > 0 are obviously tenable. Therefore, supposition h0 ðtÞ ¼ d1 hðtÞ, d1 > 0 could be considered reasonable, which is a simplest differential equation. So when the failure rate decreases, the failure rate function hdec ðtÞ can be expressed as hdec ðtÞ ¼ c1 ed1 t . When the failure rate increases, then there exist h0 ðtÞ > 0, hðtÞ > 0, so assuming that hðtÞ ¼ d2 h0 ðtÞ, d2 > 0 is reasonable. After the differential equation is solved, we obtain the increasing

887 888

t c1 d t c e 1 c2 d2 ed2 d1 þc2 d2 d1 1

860

ð49Þ ð50Þ

890 891

862 863 864 865 866 867 868 869 870 871 872 873

t

874 875 876 877 878 879 880 881 882 883 884

failure rate function hinc ðtÞ ¼ c2 ed2 . If add hinc ðtÞ to hdec ðtÞ, then the failure rate function hðtÞ of the bathtub curve model is obtained as Eq. (49). Obviously, 2 ln c2 Þ hðtÞ has only one stationary point t ¼ d2 ðln cd11dd12dþ1 in the whole domain ½0; þ1Þ. hðtÞ is monotonously decreasing in the interval ½0; t Þ, but it is monotonously increasing in the interval ½t ; 1Þ. Meanwhile, it is easy to know that the reliability function and probability density function of the novel bathtub curve model can be expressed as Eqs. (50) and (51), respectively. Some curves of hðtÞ, fðtÞ, RðtÞ of the bathtub model are shown in Fig. 9

d1 t

fðtÞ ¼ ðc1 e

t d2

t c1 d t c e 1 c2 d2 ed2 d1 þc2 d2 d1 1

þ c2 e Þe

ð51Þ

4.2. Verifying segmented distribution under circumstance of bathtub curve Supposing that the random number r  Uð0; 1Þ is known, then the solution of equation

893

894 895

896 897

t

lnð1  rÞ ¼ dc11 ed1 x  c2 d2 ed2  dc11 þ c2 d2 , namely x, is a random number obeying the distribution in Eq. (50). Every 80 random numbers that obey the bathtub curve model are taken as a group of sample. Let li ¼ 0 and cpenalty is calculated by Eq. (42), and the simulation sample is estimated by the segmented model. Then randomly generate 100 groups of sample set, and the average results of estimation of 100 times are listed in Table 2. Q, ICSD, and ln L should be worked out according to Eqs. (32), (39), and (43) after each estimation, and the averages of Bayesian estimation and MLE from the estimation results of 100 times are denoted as EðQÞ, EðICSDÞ, Eðln LÞ, respectively. In addition, the results of Bayesian estimation are taken as the initial values of MLE, and directly solve the optimization problem in Eq. (44), rather than its likelihood equation. The experiment shows that the effect of segmented distribution is satisfying. (1) No matter Bayesian estimation or MLE,

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914

CJA 741 27 December 2016

No. of Pages 21

18

H. Li et al.

Fig. 11

Comparison of fitting effect of proximity sensors in leading-edge flap when li ¼ si .

Fig. 12

Comparison of fitting effect of proximity sensors in leading-edge flap when li ¼ 0.

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

CJA 741 27 December 2016

Study on segmented distribution

Fig. 13 The best segmentation fitting of proximity sensors in leading-edge flap.

936

the effect of segmented distribution is much better than that of non-segmented distribution. When the failure rate changes sharply, it is hard to fit with the Bayesian estimation under the circumstance of one and same distribution, because samples are too disorderly and unsystematic to meet the fitting prerequisites. Sometimes, even if the distribution can be estimated, the approximation error is very great. (2) It’s not like that the larger segmented number is the better the estimation will be. The number of segments in essence is determined by the nature of sample itself, and it is also closely related with the distribution model. When the distribution model is different, the optimal segment number n may also be different. If si is inconsistent with the characteristics of the sample itself, even though the segment number n is large, the estimation error is serious, too. (3) The parameters si in segmented distribution are critical and sensitive. Even though si slightly change, the results are obviously different. Such as being the case, the subtle nuance of the segmentation pattern can lead to very different results. (4) In measuring the pros and cons of segmented distribution and its estimation, Q is a more accurate index than ICSD when n is not the same. Only when the number of segments n is the same, can ICSD be comparable.

937

5. Application instance of segmented distribution

915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935

938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953

The leading-edge flap in aircraft is an airfoil, moveable device installed on the front part of the wing. It can deflect downward or (and) slide backward (forward), which is mainly utilized to increase lifting force during flight. On flaps/slats, 30 proximity sensors are fixed, among which there are eight proximity sensors in leading-edge flap. A flap/slat electronics unit (FSEU) can monitor the impedance of the proximity sensors, which is related with the positions of the sensors. Once the leadingedge flap moves, the target will move together and the impedance will change too. In this way, the proximity sensors can measure the position of the target that moves with the control surface. The proximity sensors provides the position data of the leading-edge flap to the FSEU. According to these data, the FSEU in turn controls the indicator panel of the leadingedge device and the lamps inside the cockpit, indicating the position. The FSEU also provides these data to the flight data

No. of Pages 21

19 acquisition system, a computer in a stall management yaw damper, and the proximity switch electronic unit. Therefore, the proximity sensors in leading-edge flap are very important for the manipulation of aircraft. Considering the importance of the proximity sensors in leading-edge flap to the aircraft safety, accurate evaluation of their reliability is very necessary. Here taking a B737-700 fleet of 25 aircraft as an example, the failure data of the proximity sensors in leading-edge flap from 2001 to 2012 have been collected from an airlines. After preprocessing, the lifetime samples of the proximity sensors in leading-edge flap are listed in Table 3, where the unit of the lifetime data is day. Firstly, the lifetime samples in Table 3 are estimated by the product limit method, and the estimation results, namely, failure rate and reliability are shown in Fig. 10. In Fig. 10, it can be seen that the failure rate function hðtÞ of the proximity sensors in leading-edge flap obviously is in a phase-wise manner: in the early period of their operation, the failure rate constantly keeps at a relatively low level; in the later period, the failure rate rises sharply, which is similar to the failure model B recommended by RCM from Fig. 1. Therefore, only segmented distribution can accurately be utilized to evaluate their reliability. Firstly, the distribution of the proximity sensors in the leading-edge flap is inferred by Bayesian estimation, namely Algorithm 3. Next, the parameters in segmented distribution are estimated by MLE according to Algorithm 4. Then under the condition that cpenalty is calculated with Eq. (42) and li ¼ 0, one-, two-, and three-segment linear failure rate models and Weibull models are respectively utilized to fit the reliability of the proximity sensors in leading-edge flap. The results of Bayesian estimation and MLE are respectively listed in Tables 4 and 5 when li ¼ si . In addition, estimation results are tabulated in Tables 6 and 7 when li ¼ 0. Comparison with various distributions and methods is shown in Figs. 11 and 12, after summarizing estimation results. In order to determine the most appropriate distribution of the proximity sensors in leading-edge flap, we must find out the optimal model from all the estimation results of MLE and Bayesian estimation. In this paper, the priority sequence is sorted respectively according to Q, ICSD obtained from each model. Then add the two sequence numbers of each model together. The model whose summation of sequence number is minimum is taken as the optimal model. Thus, the optimal distribution of the proximity sensors in leading-edge flap is a two-segment linear failure rate model estimated by Bayesian estimation, i.e., Eq. (17). The estimation results are tabulated in Table 4. From Table 4, we can know that the parameters of the optimal distribution are a ¼ ½0:058192; 0:634507, b ¼ ½0:000299; 0:003251, s ¼ ½0; 231, and l ¼ ½0; 0. Fig. 13 shows the best fitting of the failure rate hðtÞ of the proximity sensors in leading-edge flap. If all samples are not divided, the best distribution is a one-segment linear failure model, whose parameters are a ¼ 0:057239, b ¼ 0:001525. From Fig. 13, the error of non-segment fitting is very great, and even it is not a valid distribution. By comparison, segment distribution has an obvious advantage. In addition, the following can be concluded from this example. (1) Usually, no matter what distribution samples follow, they are measured from perspective of the whole lifetime. However, after the samples are divided, the local property of

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014

CJA 741 27 December 2016

No. of Pages 21

20

H. Li et al.

1024

the samples plays a more important role. The segmented distribution emphasizes the local characteristics of the samples, thus the expression of segmented distribution is often not the same as that of non-segmented distribution. For instance, even if the samples follow Weibull distribution, they may no longer follow segmented Weibull distribution after segmentation. (2) Although the reliability function of segmented distribution is continuous, the segmented failure rate function hðtÞ may be often discontinuous whether it is a linear failure rate mode or segmented Weibull distribution.

1025

6. Conclusion and discussion

1015 1016 1017 1018 1019 1020 1021 1022 1023

1026 1027 1028 1029 1030 1031 1032 1033

1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060

In this paper, two n-segment segmented reliability models are proposed, and then how to estimate parameters of these segmented reliability models by Bayesian estimation and MLE is discussed. Whatever in the simulation verification or in the application examples, the segmented distribution and its estimation methods have achieved good effects. Nevertheless, there are still some problems related to the segmented distribution, which deserve our further study. (1) The unknown parameters contained in the segmented distribution are relatively too many, which makes the application of the segmented distribution become more difficult. (2) Although the estimation methods for the segmented distribution are discussed in this paper, these methods are still not accurate enough in estimating li , si , so in future research, better estimation methods should be studied, especially some methods to accurately estimate li , si . (3) The segmented function hðtÞ discussed in this article is usually discontinuous, which may be inconsistent with real circumstances. Although the reliability function is continuous, errors may still occur, so continuous segmented distribution should be discussed in the future. (4) Only when there exist change-points in failure rate, the segmented distribution can take its advantage. Thus there is a problem of detecting change-points before estimation by the segmented distribution. In this paper, merely MLE is preliminarily discussed to test changepoints, which is not sufficient. In the future, methods to detect change-points should deserve a careful study, especially Bayesian method and likelihood ratio. (5) From the experiments in this paper, it’s concluded that precision of the Gibbs sampling method is not sufficiently high. In future research, it is necessary that the Gibbs sampling method should be modified to improve its accuracy.

1061

1062

Acknowledgements

1063

This work was supported by the National Natural Science Foundation of China (Nos. 60672164, 60939003, 61079013, 60879001, 90000871), the Special Project about Humanities and Social Sciences in Ministry of Education of China (No. 16JDGC008), National Natural Science Funds, and Civil Aviation Mutual Funds (Nos. U1533128 and U1233114).

1064 1065 1066 1067 1068

References

1069

1. Zhao XL, Qian T, Mei G, Kwan C, Zane R, Walsh C, et al. Active health monitoring of an aircraft wing with an embedded piezoelectric sensor/actuator network: II. Wireless approaches. Smart Mater Struct 2007;16(4):1218–25. 2. Wu B, Tian Z, Chen M. Condition-based maintenance optimization using neural network-based health condition prediction. Qual Reliab Eng Int 2012;29(8):1151–63. 3. Bae SJ, Mun BM, Kim KY. Change-point detection in failure intensity: a case study with repairable artillery systems. Comp Indust Eng 2013;64(1):11–8. 4. Toms JD, Lesperance ML. Piecewise regression: a tool for identifying ecological thresholds. Ecology 2003;84(8):2034–41. 5. Atashgar K. Identification of the change point: an overview. Int J Adv Manuf Technol 2013;64(9–12):1663–83. 6. Bae SJ, Yuan T, Ning SL, Kuo W. A Bayesian approach to modeling two-phase degradation using change-point regression. Reliab Eng Syst Safety 2015;134:66–74. 7. Kvam PH. A change-point analysis for modeling incomplete burnin for light displays. IIE Trans 2006;38:489–98. 8. Li YX, Qian LF, Zhang W. Estimation in a change-point hazard regression model with long-term survivors. Statist Probab Lett 2013;83(7):1683–91. 9. He P, Kong G, Su Z. Estimating the survival functions for rightcensored and interval-censored data with piecewise constant hazard functions. Contemp Clin Trials 2013;35(2):122–7. 10. Uhm D, Huffer FW, Park C. Additive risk model using piecewise constant hazard function. Commun Statist—Simul Comput 2011;40(9):1458–77. 11. Liu ZH, Qian LF. Changepoint estimation in a segmented linear regression via empirical likelihood. Commun Statist - Simul Comput 2009;39(1):85–100. 12. Goodman MS, Li Y, Tiwari RC. Detecting multiple change points in piecewise constant hazard functions. J Appl Statist 2011;38 (11):2523–32. 13. Suresh RP. A test for constant hazard against a change-point alternative. Commun Statist - Theory Meth 2012;41(9):1583–9. 14. Nosek K, Szkutnik Z. Change-point detection in a shape-restricted regression model. Statist: A J Theoret Appl Statist 2013;48 (3):641–56. 15. Chou K, Tang K. Burn-in time and estimation of change-point with Weibull-Exponential mixture distribution. Dec Sci 1992;23 (4):973–90. 16. Achcar JA, Loibel S. Constant hazard function models with a change point: a Bayesian analysis using Markov chain Monte Carlo methods. Biomet J 1998;40(5):543–55. 17. Boukai B. Bayes sequential procedure for estimation and for determination of burn-in time in a hazard rate model with an unknown change point parameter. Seq Anal 1987;6(1):37–53. 18. Patra K, Dey DK. A general class of change point and change curve modeling for life time data. Ann Inst Statist Math 2002;54 (3):517–30. 19. Igba J, Alemzadeh K, Anyanwu-Ebo I, Gibbons P, Friis J. A systems approach towards reliability-centred maintenance (RCM) of wind turbines. Proc Comp Sci 2013;16(1):814–23. 20. Dehghanian P, Fotuhi-Firuzabad M, Aminifar F, Billinton R. A comprehensive scheme for reliability centered maintenance in power distribution systems—Part I: methodology. IEEE Trans Power Deliv 2013;28(2):761–70. 21. Li YX, Qian LF. Likelihood ratio test for a piecewise continuous Weibull model with an unknown change point. J Math Anal Appl 2014;412(1):498–504. 22. Pandya M, Jani PN. Bayesian estimation of change point in inverse Weibull sequence. Commun Statist-Theory Meth 2006;35 (12):2223–37.

1070

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133

CJA 741 27 December 2016

Study on segmented distribution 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164

23. Zhang T, Xie M. On the upper truncated Weibull distribution and its reliability implications. Reliab Eng Syst Safety 2011;96 (1):194–200. 24. Sarhan AM, Apaloo J. Exponentiated modified Weibull extension distribution. Reliab Eng Syst Safety 2013;112(4):137–44. 25. Maboudou-Tchao Edgard M, Hawkins Douglas M. Detection of multiple change-points in multivariate data. J Appl Statist 2013;40 (9):1979–95. 26. Hu Y, Zhao K, Lian H. Bayesian quantile regression for partially linear additive models. Statist Comput 2015;25(3):651–68. 27. Alhamzawi R, Yu K. Conjugate priors and variable selection for Bayesian quantile regression. Comput Statist Data Anal 2013;64 (4):209–19. 28. Reich BJ, Smith LB. Bayesian quantile regression for censored data. Biometrics 2013;69(3):651–60. 29. Nosek K. Schwarz information criterion based tests for a changepoint in regression models. Statist Papers 2010;51(4):915–29. 30. Pan JM, Chen JH. U-statistic based modified information criterion for change point problems. Commun Statist - Theory Meth 2008;37(17):2687–712. 31. Chen JH, Pan JM. Information criterion and change point problem for regular models. Ind J Statist 2006;68(2):252–82. 32. Shen G, Ghosh JK. Developing a new BIC for detecting changepoints. J Statist Plan Inf 2011;141(4):1436–47. 33. Ninomiya Y. Information criterion for Gaussian change-point model. Statist Probab Lett 2005;72(3):237–47. 34. Hannart A, Naveau P. An improved Bayesian information criterion for multiple change-point models. Technomet A J Statist Phys Chem Eng Sci 2012;54(3):256–68. 35. El-Gohary A, Alshamrani A, Al-Otaibi AN. The generalized Gompertz distribution. Appl Math Model 2013;37(1–2):13–24.

No. of Pages 21

21 36. Mazucheli J, Coelho-Barros EA, Achcar JA. Inferences for the change-point of the exponentiated Weibull hazard function. REVSTAT 2012;10(3):309–22. 37. Xie M, Lai CD. Reliability analysis using an additive Weibull model with bathtub-shaped failure rate function. Reliab Eng Syst Safety 1996;52(1):87–93. 38. Bebbington M, Lai CD, Wellington M, Zitikis R. The discrete additive Weibull distribution: a bathtub-shaped hazard for discontinuous failure data. Reliab Eng Syst Safety 2012;106 (5):37–44. 39. Jiang R. A new bathtub curve model with a finite support. Reliab Eng Syst Safety 2013;119:44–51. Huaiyuan Li is a lecturer in Jinling Institute of Technology. He graduated and obtained a doctor degree in engineering from Nanjing University of Aeronautics and Astronautics in 2016. He received the Master of Engineering in computer application and technology from Beifang University of Nationalities in 2008. His current research interests include reliability statistics, artificial intelligence and software engineering.

1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185

Hongfu Zuo is currently a professor and PhD student supervisor in Nanjing University of Aeronautics and Astronautics. Presently, his research interests include reliability engineering, vehicle application engineering, man-machine maintainability and so on

1186

Juan Xu is a lecturer in Nanjing University Of Aeronautics And Astronautics. She is a Ph.D. student in Nanjing University of Aeronautics and Astronautics. She interest in reliability, safety and maintainability of equipment, supportability engineering.

1191

Please cite this article in press as: Li H et al. Study on segmented distribution for reliability evaluation, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j. cja.2016.12.008

1187 1188 1189 1190 1192 1193 1194 1195