- Email: [email protected]
Bayesian inference for Common cause failure rate based on causal inference with missing data
Bayesian inference for Common cause failure rate based on causal inference with missing data
Journal Pre-proof
Bayesian inference for Common cause failure rate based on causal inference with missing data H.D. Nguyen, E. Gouno PII: DOI: Reference:
S0951-8320(19)30879-8 https://doi.org/10.1016/j.ress.2019.106789 RESS 106789
To appear in:
Reliability Engineering and System Safety
Received date: Revised date: Accepted date:
5 July 2019 14 November 2019 28 December 2019
Please cite this article as: H.D. Nguyen, E. Gouno, Bayesian inference for Common cause failure rate based on causal inference with missing data, Reliability Engineering and System Safety (2019), doi: https://doi.org/10.1016/j.ress.2019.106789
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.
Highlights • Introducing a new methodology using MCMC to estimate parameters of a common cause failure model. • Providing an alternative to deal with missing data in contingency tables. • Performing causality-based inference for common cause of failure. • Constructing strategies to elicit hyperparameters for a Bayesian approach. • Proposing a stochastic version of the Inverse Bayes Formulae algorithm.
1
Bayesian inference for Common cause failure rate based on causal inference with missing data Nguyen, H. D.a,b , Gouno, E.c, a Institute
of Research and Development, Duy Tan University, Danang 550000, Vietnam of Mathematics, Faculty of Information Technology, Vietnam National University of Agriculture, Hanoi, Viet Nam c University of South Brittany, LMBA, Campus de Tohannic, 56017 Vannes, France
b Department
Abstract This paper proposes a methodology to handle the causality to make inference on common cause failure (CCF) in a situation of missing data. The data are collected in the form of contingency table but the available information are only the numbers of CCF of different orders and the numbers of failure due to a given cause, i.e. the margins of the contingency table. The frequencies in each cell are unknown; we are in a situation of missing data. Assuming a Poisson model for the count, we suggest a Bayesian approach and we use the inverse Bayes formula (IBF) combined with a Metropolis-Hastings algorithm to make inference on the parameters. The performance of the resulting algorithm is evaluated through simulations. A comparison is made by analogy with results obtained from the recently proposed α-composition method. Keywords: Common-cause failure, Inverse Bayesian formula, Contingency table, Missing data, Causal inference.
1. Introduction Probability risk assessment (PRA) studies are performed to identify and understand vulnerabilities of systems. They were initially introduced in the nuclear industry. PRA is an analytical tool that allows to track potential system failures 5
and to determine the likelihood and consequences of their occurrences. This tool is now applied in many aeras and the literature on the subject is important. The NASA published a 83 pages bibliography in 2000 [1] and a procedures guide for ∗ Corresponding
author Email address: [email protected] (Gouno, E. )
Preprint submitted to Reliability Engineering & System Safety
January 3, 2020
manager and practitioners in 2011 [2]. Many books related to applications of PRA on reliability engineering problems have been published [3, 4, 5]. Kelly 10
and Smith consider a Bayesian approach of PRA [6]. A survey of application of PRA for instrumentation and control systems in nuclear power plants is given by Lu and Jiang [7]. Safie et al. [8] discuss the role of reliability and PRA at NASA. Pointing out that traditional PRAs fail to capture the dynamic behaviour of complex systems, Nielsen et al. [9] consider a dynamic PRA (DPRA)
15
and propose an optimization methods to deal with a combinatorial explosion of scenarios. More recently, Georges-William et al.[10] provide a probabilistic risk assessment of station blackout in nuclear power plants. They use an approach based on a hybrid of Monte Carlo methods, multistate modeling, and network theory as an alternative to the use of static fault tree which turns out
20
to deteriorate in applicability with increasing system size and complexity. In PRA, a common cause failure refers to “the failure or unavailable state of more than one component during the mission time and due to the same shared cause” [11]. In industry where a high reliability is required, these events are crucial challenges. The study of CCF events then plays a very important
25
role in system reliability assessments. A number of models has been developed to analyze CCF events such as the Binomial failure rate model [12, 13, 14], the α-factor model [15, 16, 17], the Beta factor model, and the multiple Greek letter model [15]. The CCF is also a major concern in several recent studies, for example, see [18, 19, 20]. However, in general, many previous studies
30
deal essentially with the magnitude of CCF events, i.e. the number of components failing simultaneously provoking the failure of the system, but they do not consider the cause of those events. This leads to some limitations for their use as discussed in [21]. Causality-based inference needs to be investigated and the U.S Nuclear Regulatory Commission has endeavored to develop
35
a database for this problem [22]. The Commission report [23] gives somes examples of possible CCF causes. For example, collecting the information from the emergency diesel generators in a set of nuclear plants, the report points out that amongst the 41 CCF events recorded, 15 events are caused by Design/Construction/Installation/Manufacture proximate cause group; 9 events
40
are qualified Internal to Component proximate cause group; 9 events are the
3
consequence of cause group related to Operational/Human Error, 5 events are the result of cause groups related to External Environment and 3 events are caused by Other proximate cause group. The report [24] provides a way to achieve further understanding of the occurrence of CCF events by presenting 45
the process of data collection and grouping the ranking of proximate failure causes. Rasmuson et al. [25] quantify the conditional probability of CCF considering the asymmetry within a common-cause component group, given independent failures or failures with common-cause potential. Using Bayesian networks,
50
Kelly et al. [26] are the first to propose the preliminary framework of a causalitybased inference, providing cause-specific quantitative insights into likely causes of failures. O’Connor and Mosleh [21] develop a general cause based methodology for analysis of common cause and dependent failures in system risk and reliability assessments. Johnston [27] presents a structured procedure to han-
55
dle dependent failure and their causes. A CCF model updating is presented in [28] where two updating processes are proposed considering the difficulty to identify failure causes immediately. Sakurahara et al. [29] utilize simulation models of physical failure mechanisms to capture underlying causalities and to generate simulation-based data for the CCF probability estimation. Zheng et al.
60
[30] suggest the α-decomposition method. In this decomposition, the authors consider a set of CCF data involving the information about not only the magnitude of external causes but also the cause occurrence frequencies of failures. This data forms a contingency table: the rows of the table record the number of components involved in each CCF event while its columns record the cause of
65
the event. As a result, the value in each cell of the table represents the number of failures due to a specific cause. In many situations it is not possible to observe complete contingency tables. We have missing data. Missing data are common in many data-driven approaches and the litterature on their statistical analysis grows up since the
70
early 1970s. Little and Rubin [31] offer a complete survey of methodologies for handling missing data problem. In particular, they provide a comprehensive study of methods for treating contingency table and they describe imputation methods in details. Tan et al. [32] present Bayesian solutions to missing data
4
problems using methods relying on EM-type algorithms, data augmentation and 75
noniterative computation. In the problem addressed here, one can not observe exactly the number of observations in each cell. Instead, only the values in the margins of the table are observed. Incomplete contingency table arises in many areas and a number of methods has been proposed in the literature. Ng et al. introduce new families of Dirichlet distributions [33, 34, 35, 36] and they use the
80
inverse Bayes formula (IBF) [37] to obtain closed-form solutions for incomplete categorical data Bayesian analysis. Several applications of the IBF can be found in [32]. The method presented in this paper is a contribution to probabilistic risk assessment where the causality of failure is taken into account. One of the key
85
points of the method is to enable the analysis of crude CCF data. That is to say that time and cost can be saved in the sense that no special procedures for collecting or correcting (completing) the data are required. We consider a scenario of CCF data mentioned in the paper of Zheng et al. [30]. The CCF data used in their paper are for illustration only and not from any real
90
database. However, one can easily imagine practical situations where for some reasons, number of components involved in the failure and cause of the failure are recorded independently. Therefore the CCF data do not contain the counting number of CCF events classified by magnitude and possible cause. We suggest a Bayesian approach to estimate the cell probabilities in this context of missing
95
data. An exact non-iterative sampling procedure proposed by Tian et al. [38] is developed. This procedure relies on the IBF. To overcome some practical difficulties specific to our problem, we use a Metropolis-Hasting (MH) algorithm to generate from conditional probability expressed with the IBF. Numerical results show the efficiency of the method. We show that the α-decomposition
100
parameters estimates are close to estimates obtained with our proposed method. The paper is organized as follows. In section 2, we present two situations of data and the corresponding distributions. The Bayesian inference for parameters of the model in the case of incomplete data using the IBF method is given in section 3. Section 4 presents a stochastic version of the IBF sampling. A
105
simulation study is investigated in section 5. Section 6 shows an application of the proposed method to the dataset provided in [30]. Finally, some concluding
5
remarks are given in section 7. 2. Data and distributions Let us consider a m × k contingency table with cell Nij . In the contex of 110
CCF, considering a m-component system observed in a time window of length T (unit of time) and k possible causes of failure C1 , . . . , Ck , Nij is a random variable corresponding to the number of CCF of order i from cause j occurring in [0, T ].
The matrix Nij , i = 1, . . . , m, j = 1, . . . , k that we denote N, represents
complete data. Each cell represents a number of failures occurring in a given time window. This number depending on time is a counting process. A natural way to model this number is to use a homogeneous Poisson process that is to say a Poisson process with a constant intensity. The consequence is then that in a given time window, the number of events follows a Poisson distribution. Then we suppose that the cell counts are independent homogeneous Poisson processes with respective intensity λij . The expected number of CCF of order i from cause j is therefore λij T . When Nij is observed, estimates for the λij ’s are easily obtained. Suppose now that the observed data are the total margins of the contingency table that is to say: Ni,• =
k X
Nij ,
j=1
the sum of the terms on the row i, corresponding to the total number of CCF events of order i, i = 1, ...., m, and N•j =
m X
Nij ,
i=1
the sum of the terms on the column j, corresponding to the total number of 115
components failed from cause j, j = 1, ...., k. In this situation, the cell counts Nij are not directly observed. We have missing data; the observed data are N1• , . . . , Nm• , N•1 , . . . , N•k , that we denote N• . We refer to N• as incomplete data and to N as complete data. Data displayed by Zheng et al. in table 1, p.24 in [30], are of this form. To
120
each line (system) of this table corresponds a contingency table where only the margins are observed. 6
Our purpose is to estimate the m × k matrix of terms λij , denoted by λ , in this situation of missing data. The distribution of N, the complete data, is obtained as a product of Poisson 125
distributions:
Y (λij T )nij
P (N = n | λ ) =
nij !
i,j
e−λij T .
(1)
Thus, the distribution of the incomplete data can be expressed as: P (N• = n• | λ ) =
X
n∈Dn•
P (N = n | λ ),
(2)
Pm where Dn• is the set of all the m×k matrices n = nij such that i=1 nij = n•j Pk and j=1 nij = ni• . In the following we aim at developing a Bayesian approach to estimate λ .
130
3. Bayesian inference: the IBF method Ng [37] introduces a method to compute posterior distributions based on the inverse Bayes formula: the IBF method. We are going to apply this method in our specific context of CCF. When the size of Dn• is not too large, the method leads to exact solutions. When it is not the case, we suggest a MH algorithm. λ), we want to compute the posterior distriAssuming a prior distribution π(λ λ | N• = n• ). Applying the Bayes theorem and the formula (2), we bution π(λ have: λ | N• = n• ) = π(λ
P
n∈Dn•
λ) P (N = n | λ ) π(λ
P (N• = n• )
.
But λ) = π(λ λ | N = n)P (N = n) P (N = n | λ ) π(λ and P (N = n)/P (N• = n• ) = P (N = n | N• = n• ). Therefore λ | N• = n• ) = π(λ
X
n∈Dn•
λ | N = n) P (N = n | N• = n• ), π(λ
(3)
We assume that the terms of λ are independent and each λij has a Gamma distribution with parameters (αij , βij ) which is a conjugate prior for the Poisson
7
distribution. The prior distribution of λ is then a product of Gamma distribu135
tions and the posterior distribution given the complete data is a product of Gamma distributions with parameters (nij + αij , T + βij ). That is, for n ∈ Dn• , Y
λ | N = n) ∝ π(λ
n +αij −1 −λij (T +βij )
λijij
e
.
(4)
i,j
We need now to compute P (N = n | N• = n• ). Let λ ∗ be an arbitrary value of λ. We have (inverse Bayes formulas): P (N = n) = and P (N• = n• ) = 140
Therefore
λ∗ ) P (N = n | λ ∗ ) π(λ ∗ λ | N = n) π(λ λ∗ ) P (N• = n• | λ ∗ ) π(λ . λ∗ | N• = n• ) π(λ
λ∗ | N• = n• ) P (N = n | λ ∗ ) π(λ P (N = n) = . λ∗ | N = n) P (N• = n• | λ ∗ ) P (N• = n• ) π(λ
(5)
Doing the summation on both sides with respect to n ∈ Dn• , we deduce: X P (N = n | λ ∗ ) P (N• = n• | λ ∗ ) = . λ∗ | N = n) λ∗ | N• = n• ) π(λ π(λ
n∈Dn•
And replacing in (5), we obtain λ∗ | N = n) P (N = n | λ ∗ )/π(λ . ∗ λ∗ | N = u) u∈Dn P (N = u | λ )/π(λ
P (N = n | N• = n• ) = P
•
Using expressions (1) and (4), we have: P (N = n | λ∗ ) Y T nij Γ(nij + αij ) ∝ = q(n) π(λ∗ | N = n) nij ! (T + βij )nij +αij i,j
(6)
n o Let us denote Dn• = n1 , . . . , nK and p` = P (N = n` | N• = n• ), that is q(n` ) , ` = 1, . . . , K. p` = PK j=1 q(nj )
(7)
Relying on this notation, we can express the posterior marginal distribution of the terms of λ as: π(λij | N• = n• ) =
K X `=1
8
p` π(λij | N = n` ),
(8)
145
which is a mixture of Gamma distribution. ˜ ij of λij are the expectation Under a quadratic loss, the Bayes estimators λ of the posterior distribution. Therefore, ˜ ij = λ
K X `=1
p`
hn
+ αij i . T + βij
`,ij
(9)
The following example illustrates the method when the size of Dn• is not too large, thus p = (p1 , . . . , p` ) can be easily computed. 150
Consider a set of CCF data from a 2-component system observed over a time window of length T = 10 unit of time. Suppose that the potential causes are divided into three main groups, say “cause 1”, “cause 2” and “cause 3” and that we have the margins, or observed data n• = {(3, 5, 2), (4, 6)}. The space Dn• can be 0 2 Dn• = 3 3 1 3 2 2 3 1 0 4
155
described as: 0 4 0 1 1 2 1 2 1 0 3 1 2 , , , , , 3 1 2 2 4 0 2 3 1 3 2 1 0 0 2 0 2 2 1 1 2 2 0 3 0 1 , , , , , 2 1 5 0 1 4 1 1 3 2 0 5 1 0 . 2
Let us consider a noninformative distribution for each λij , i.e., λij ∼ G (0.5, 0).
The probabilies p` = P (N = n` | N• = n• ), n` ∈ Dn• , are calculated from (7) and displayed in Table 1. `
p`
`
p`
1
0.10608204
7
0.13366337
2
0.07072136
8
0.04950495
3
0.12376238
9
0.06364922
4
0.07425743
10
0.14851485
5
0.04243281
11
0.12376238
6
0.06364922
Table 1: The conditional distribution P (N = n | N• = n• )
9
Applying (9), we obtain the following Bayes estimate for λ 0.19908 0.20700 0.14391 . λ˜ = 0.20092 0.39300 0.15608 4. A stochastic version of the IBF sampling When the space Dn• is too large the procedure described above is not pos160
sible since the computation of p relies on the size of Dn• . Sampling strategies are then considered. As mentioned previously, the incomplete-data marginal posterior distribution is a mixture of gamma distribution. Therefore a simple scheme to sample from this distribution is 1. choose a component of the mixture with probability p; 2. generate the λij ’s with the complete-marginal posterior
165
distributions (8). When Dn• can be identified this procedure is the exact IBF sampling suggested by Tian et al. [38]. The key in our situation is therefore to draw components of the mixture with the probability p, the conditional predictive distribution P (N = n | N• = n• ). The components of this probability has a cumbersome denominator. The MH
170
algorithm is a well-known sampling algorithm to overcome this kind of difficulty. The algorithm belongs to the family of MCMC. The principle is to build a Markov chain which stationary distribution is the distribution of interest. The algorithm is described as follows. Let Q be a specified transition matrix from one element of the space to another. Starting from n(0) ∈ Dn• , after p
175
cycles of the following scheme, n(p) is available, do: 1. Generate a candidate n(∗) ∼ Q(. | n(p) ). 2. Compute
q(n(∗) ) ρ = min 1, . q(n(p) )
3. Set n(p+1) = n(∗) with probability ρ; else, set n(p+1) = n(p) . After a number of iterations of this scheme, the sequence obtained is a sample of the conditional predictive distribution. 180
Many choices for Q, the generating-candidate matrix, are possible. In our case, the following scheme is applied to generate the candidate:
10
1. Choose randomly a 2 × 2 submatrix of n(p) with probability 1/ a11 a12 A= a21 a22
k 2
m 2
, say
2. Draw x from a uniform distribution in (max{0, a•1 −a2• }, . . . , min{a1• , a•1 }), where a1• = a11 + a12 and a•1 = a11 + a21 . 3. Replace the corresponding matrix A in n(p) by the matrix B defined as x a1• − x B= a•1 − x a2• − a•1 + x where a2• = a21 + a22 .
185
That means, we generate a new candidate by randomly changing a 2 × 2 submatrix of the current sample. The principle in step 2 and step 3 is to make sure that sum of rows and sum of columns of these sub-matrices are the same. The construction of the generating-candidate matrix is specific to our situation but it has been inspired by classical approaches that consider random walks.
190
The figure 1 compares the exact probabilities p` in Table 1 and the approximated probabilities obtained with the MH algorithm. One can see that the approximation is accurate; the algorithm is efficient. The suggested algorithm is an exact IBF sampling where Dn• is randomly covered using a MH algorithm. To summarize, the algorithm that we call the
195
stochastic IBF sampling or MH IBF sampling is described as follows: 1. Generate M
i.i.d. samples n(1) , . . . , n(M ) from the MH algorithm de-
scribed above, λ | N = n(`) ), the complete-data posterior distribution 2. Generate λ (`) ∼ π(λ defined in (4). 200
λ(`) }`=1,...,M are i.i.d. samples from the posterior π(λ λ | N• = n• ). Then {λ λ | N• = Considering the data from the previous example, sampling from π(λ n• ) using the stochastic IBF sampler and computing the empirical mean, we obtain the following estimate for λ : 0.1978159 0.2069460 λ˜ = 0.2018113 0.3933634 11
0.1431974 , 0.1561672
0.30 0.20 0.15 0.10 0.00
0.05
Probability
0.25
Theoretical distribution Equilibrium distribution
1
2
3
4
5
6
7
8
9
10
11
Elements
Figure 1: Comparison of the theoretical distribution P (N | N• = n• ) and the equilib-
rium distribution of the MH samples for the simple case
which is close to the estimate obtained with the exact distribution p` , ` = 205
1 . . . , K. We observed that in practice, the sampler is not sensitive to the choice of the burn-in size; we have roughly a steady acceptance ratio for different combinations of sample size and burn-in size. Therefore we consider the value M = 100000 and a burn-in size 50000 for this study.
5. A simulation study 210
We evaluate the performance of the proposed method through simulations. The simulation procedure is performed as follows: 1. Assign an input value for λ = (λij ) and a time window T . 2. Simulate a number in each cell of the contingency table from a Poisson distribution with parameter λij T and compute the margins.
215
3. Use the stochastic IBF sampler to estimate λ based on the two margins computed in 2. (1)
4. Repeat steps 2-3 M times to obtain the estimates λ˜
12
, . . . , λ˜
(M )
.
Guess
n11
n12
n13
n21
n22
n23
n31
n32
n33
Total
Input
6
5
1
4
3
1
4
1
2
27
Situation 1
7
4
1
3
4
2
4
3
1
29
Situation 2
4
3
1
3
2
1
2
1
1
17
Situation 3
9
7
3
6
4
2
5
3
1
40
Table 2: The average number of CCF events in each cell of contingency table M 1 X ˜ (i) ˜ 5. Compute the mean λ = λ and mean square error (MSE). M i=1
For purpose of illustration, we consider a 3 × 3 contingency table. In our
experiment, setting T = 120 and an average number of CCF events in each cell of the contingency table, the matrix input λ is deduced from λij T , an approximation of the chosen average numbers (first row in Table 2). The parameter λ used for simulation 0.0500 λ = 0.0333 0.0333
is therefore: 0.0083 0.0250 0.0083 . 0.0083 0.0167 0.0416
To specify values for prior hyperparameters, we use the strategy suggested in 220
[14]. We consider three situations of guess at the average numbers of CCF events in the time window of length T presented in Table 2. The first situation corresponds to a guess at the total number of failures close to the input number. The second and the third one correspond to situations where the guess at the total number underestimates (resp. overestimates) the input number. Then,
225
the prior hyperparameters for each λij is defined by solving a simple system of two equations matching guesses and uncertainty with mean and variance of the Gamma distribution:
α /β = ηij , ij ij
α /β 2 = ρ2 η 2 , ij ij ij ij
(10)
where ηij represents the guess at λij and ρij is the corresponding uncertainty of the guess. The hyperparameters for each situation of the guesses and the 230
corresponding assumed uncertainty are given in Table 3. They are determined as follows. We consider two levels of uncertainty: low and high. We choose ρij from a uniform distribution in [0.05, 0.3] for the low level and [0.7, 0.95] for 13
Situation 1
Situation 2
Situation 3
ρ
Low
High
Low
High
Low
High
Type of prior
A
B
C
D
E
F
(α11 , β11 )
34.59
593.07
1.30
22.41
23.09
692.70
1.99
59.98
15.18
202.52
1.36
18.16
(α12 , β12 )
51.12
2045.17
1.46
58.47
33.27
1331.10
1.43
57.50
32.18
643.70
1.83
36.65
(α13 , β13 )
17.54
526.46
1.66
49.89
24.96
1497.81
1.77
106.33
21.74
521.89
1.58
37.99
(α21 , β21 )
12.15
364.56
1.18
35.32
33.79
1351.95
1.35
54.33
11.06
189.62
1.17
20.18
(α22 , β22 )
47.81
1434.53
1.12
33.62
34.70
2082.51
1.24
74.56
12.12
363.80
1.22
36.89
(α23 , β23 )
15.16
606.50
1.52
60.89
18.18
2182.50
1.14
136.99
20.16
806.68
1.53
61.43
(α31 , β31 )
89.46
10736.35
1.88
226.11
60.09
7211.36
1.34
160.85
10.64
425.81
1.15
46.23
(α32 , β32 )
12.45
747.06
1.41
84.72
32.23
3868.41
1.13
136.52
13.71
823.03
1.30
78.10
(α33 , β33 )
18.97
2276.96
1.75
210.21
21.65
2598.92
1.20
144.76
14.12
1694.50
1.32
158.26
Table 3: Hyperparameters of prior based on the guess in Table 2.
the high level of uncertainty. The type A (resp. B) corresponds to the guess in situation 1 with low (resp. high) level of uncertainty. The type C (resp. D) 235
corresponds to the guess in situation 2 with low (resp. high) level of uncertainty. The type E (resp. F) corresponds to the guess in situation 3 with low (resp. high) level uncertainty. The type G is the noninformative case, that is to say gamma distributions with parameters (0.5, 0). The corresponding estimates are presented in Table 4. In general, when
240
the guess is close to the input, the small value of ρ leads to better estimates. Otherwise, when the guess is far from input, the large value of ρ give better results. Moreover, the noninformative prior gives acceptable results with small MSE.
6. Application 245
Zheng et al. [30] simulate CCF events and corresponding causes for 16 systems with 3 components and consider the results displayed in Table 5. They propose a method: the α- decomposition method to make inference from these data. It relies on an extension of the classical α-factor model. The classical α-factor model [15] for a m-component system considers m parameters (α1 , . . . , αm ) to
250
characterize the CCF’s. αi refers to the probability that exactly i components are involved in a CCF event given that at least one component fails.
14
λ11
λ12
λ13
λ21
λ22
λ23
λ31
λ32
λ33
Input
Type A
Type B
Type C
Type D
Type E
Type F
Type G
0.05
0.05887
0.06104
0.03604
0.04546
0.06528
0.05336
0.0612
(8.93e-5)
(4.74e-4)
(2.02e-4)
(1.80e-4)
(2.78e-4)
(3.31e-4)
(6.91e-4)
0.03236
0.03038
0.02603
0.03424
0.04973
0.03807
0.0353
(9.86e-5)
(2.56e-4)
(2.45e-4)
(1.94e-4)
(9.34e-5)
(2.07e-4)
(2.89e-4)
0.00834
0.00858
0.00836
0.00912
0.02301
0.01802
0.0160
(3.44e-8)
(2.45e-6)
(1.03e-8)
(5.40e-6)
(2.18e-4)
(1.35e-4)
(1.01e-4)
0.02495
0.02525
0.02577
0.03113
0.04761
0.03835
0.03622
(7.06e-5)
(1.31e-4)
(5.84e-5)
(1.08e-4)
(20.8e-4)
(1.50e-4)
(2.62e-4)
0.03282
0.02728
0.01699
0.020275
0.03209
0.02394
0.02781
(6.20e-5)
(1.01e-4)
(6.43e-5)
(5.98e-5)
(3.43e-5)
(5.39e-5)
(1.14e-4)
0.01653
0.01607
0.00834
0.00903
0.01592
0.01330
0.01482
(6.87e-5)
(8.85e-5)
(3.09e-8)
(6.45e-6)
(5.82e-5)
(3.67e-5)
(7.31e-5)
0.03291
0.03085
0.01768
0.02418
0.04010
0.03421
0.03155
(6.76e-6)
(1.04e-4)
(2.45e-4)
(1.30e-6)
(5.24e-5)
(1.21e-4)
(2.00e-4)
0.02399
0.02172
0.00866
0.01104
0.02418
0.02134
0.02432
(2.47e-4)
(2.08e-4)
(2.32e-8)
(1.47e-5)
(2.52e-4)
(1.93e-4)
(3.25e-4)
0.00833
0.00832
0.00854
0.01058
0.00822
0.00775
0.01481
(6.94e-5)
(7.11e-5)
(6.61e-5)
(5.01e-5)
(7.14e-5)
(8.15e-5)
(3.52e-5)
0.0416
0.0083
0.0333
0.025
0.0083
0.0333
0.0083
0.0167
Table 4: Means of Bayesian estimates and MSE for various types of prior, 100000
simulations.
15
Common causes occurence
CCF event
Cause 1
Cause 2
Cause 3
Total
Order 1
Order 2
Order 3
System 1
32
28
67
127
113
11
3
System 2
17
78
11
106
98
7
1
System 3
18
19
50
87
73
9
5
System 4
29
6
31
66
53
5
8
System 5
7
33
10
50
45
4
1
System 6
15
9
17
41
33
3
5
System 7
12
15
7
34
32
2
0
System 8
2
22
7
31
29
2
0
System 9
7
4
11
22
20
2
0
System 10
10
8
3
21
20
2
0
System 11
3
6
10
19
16
2
1
System 12
7
3
6
16
14
1
1
System 13
3
5
7
15
13
1
1
System 14
5
3
7
15
12
1
2
System 15
4
5
2
11
9
1
1
System 16
1
6
2
9
7
1
1
Table 5: A hypothetical database of CCF events provided in [30]
In the α-decomposition model, the parameters αi are decomposed as: αi =
k X
C
αi j rj , i = 1, . . . , m,
(11)
j=1
C
where αi j is the conditional probability that, given a CCF occurs from cause Cj , exactly i components are involved and rj is the probability that, given a 255
CCF occurs, it is from cause Cj . C
The decomposed parameters αi j , i = 1, . . . , m; j = 1, . . . , k are the new parameters of interest. Zheng et al. [30] present a hierarchical model and develop a Bayesian methodology to estimate the decomposed α-factors. C
The parameters αi j which are probabilities, can be roughly approximated 260
by the ratio between the number of CCF of order i from cause j and the total number of CCF of any order from cause j. These numbers can be respectively P3 approximated by λij T and λ•j T , where λ•j = i=1 λij . Therefore:
C
αi j ≈
λij λij T = , i = 1, . . . , m, j = 1, . . . , k, λ•j T λ•j
(12)
This relationship allows us to compare the results obtained with MH-IBF method 16
265
and the α-decomposition method. Since the systems are of the same nature with Poisson counting process, the data considered to apply the MH-IBF method are: n1• = 587, n2• = 53, n3• = 530, n•1 = 172, n•2 = 250, n•3 = 248. That is to say, we aggregate the data from Table 5. With noninformative priors, assuming that the observation time
270
is T = 100 (unit of time), the results characterizing the posterior distribution are displayed in Table 6. ˜ ij , it can be seen Considering mean of each posterior distribution, that is λ that ˜ 11 < λ ˜ 13 < λ ˜ 12 , λ ˜ 22 < λ ˜ 21 < λ ˜ 23 , λ
(13)
˜ 32 < λ ˜ 33 < λ ˜ 31 . λ C
Thus, the order among λij ’s is the same as the order among αi j ’s in the α275
decomposition model. Similar to the discussion in [30], cause 2 is of least CCF risk, cause 3 is best at provoking partial CCF and cause 1 is best at provoking complete CCF. However, the impact of cause 2 and cause 3 are similar as the difference between them are not significant. Based on this result, the most hazardous causes could be recognized, which has a practical engineering meaning. Figures 2-4 show the smoothed histogram as a density estimates of the posterior densities of each λij based on M = 100000 i.i.d. samples generated using the sampling MH-IBF method (the histogram is displayed only for λ11 , λ21 , λ31
4
in each graph).
1
2
3
λ11 λ12 λ13
0
Probability density function
280
0.8
1.2
1.6
2.0
2.4
2.8
Figure 2: Posterior distributions of λ1j (j = 1, 2, 3)
17
Posterior λij
Mean
SD
0.025-quantile
Median
0.975-quantile
λ11
1.4454
0.1487
1.1625
1.4442
1.7461
λ12
2.2293
0.1718
1.9016
2.2269
2.5725
λ13
2.2095
0.1728
1.8777
2.2075
2.5541
λ21
0.1825
0.0906
0.0389
0.1735
0.3839
λ22
0.1820
0.0888
0.0397
0.1708
0.3761
λ23
0.1805
0.0902
0.0379
0.1703
0.3811
λ31
0.1069
0.0541
0.0213
0.0986
0.2272
λ32
0.1041
0.0536
0.0217
0.0973
0.2262
λ33
0.1042
0.0538
0.0211
0.0983
0.2260
CCF order 1
CCF order 2
CCF order 2
2
4
λ21 λ22 λ23
0
Probability density function
6
Table 6: Summary of posterior distributions of each λij , (i, j = 1, 2, 3).
0.0
0.1
0.2
0.3
0.4
0.5
0.6
Figure 3: Posterior distributions of λ2j (j = 1, 2, 3)
18
8 6 4 2 0
Probability density function
λ31 λ32 λ33
0.0
0.1
0.2
0.3
0.4
Figure 4: Posterior distributions of λ3j (j = 1, 2, 3)
Under a quadratic loss, the Bayes 1.4454 ˜ λ = 0.1825 0.1069 C
estimate of λ is: 2.2293 2.2095 0.1820 0.1805 . 0.1041 0.1042
The parameters αi j in α-factor model can be approximated using (12). 285
C
Table 7 compares the estimates of αi j , 0.025 and 0.975 quantiles obtained from the method suggested by Zheng et al. [30] and the aproximation obtained injecting the estimates of λij ’s in (12). Mean
Values
0.025-quantile
0.975-quantile
α-decomposition
MH-IBF
α-decomposition
MH-IBF
α-decomposition
MH-IBF
α1C1
0.813
0.833
0.603
0.713
0.961
0.932
α1C2
0.910
0.886
0.826
0.804
0.978
0.955
α1C3
0.820
0.885
0.687
0.798
0.944
0.953
α2C1
0.0988
0.1052
0.00389
0.0105
0.277
0.215
α2C2
0.0729
0.0723
0.00979
0.0147
0.153
0.147
α2C3
0.112
0.0723
0.0120
0.0162
0.221
0.202
α3C1
0.0879
0.0616
0.00456
0.0124
0.218
0.198
α3C2
0.0174
0.0414
0.00052
0.0015
0.0561
0.0894
α3C3
0.0679
0.0417
0.00540
0.0088
0.151
0.139
Cj
Table 7: Compare the estimates of αi
using the α-decomposition method and our
proposed method.
19
C
The order among αi j ’s when j is fixed (similar to (13)) is preserved. The figure 5 shows that the estimates and quantiles with both methods
290
are close and overlap with each other. In other words, with consideration of uncertainty bounds, the results of the MH-IBF method are consistent with the
0.8
1.0
results of the α-decomposition method.
0.0
0.2
0.4
0.6
α−decomposition 0.025−quantile for α−decomposition 0.975−quantile for α−decomposition MH−IBF 0.025−quantile for MH−IBF 0.975−quantile for MH−IBF
C
α1 1
C
α1 2
C
α1 3
C
C
α2 1
α2 2
Cj
Figure 5: Comparison of the parameters αi
C
α2 3
C
α3 1
C
α3 2
with their corresponding approximation
using MH-IBF method
295
The MH-IBF is an plausible alternative to the α-decomposition method. The MH-IBF is constructed considering the natural assumption that we have Poisson distributions for the counts and the parameters to be estimated are the expected numbers of CCF. The α-decomposition does not have this simple direct interpretation. Therefore in the context of PRA study, MH-IBF can be
300
more adequate.
20
C
α3 3
7. Concluding remarks This paper has developped a new method for estimating a Common Cause Failure model parameters based on missing data in the form of contingency table. We have presented a procedure to characterize the relationship between 305
potential cause and failure for system reliability. The data considered are of the form of contingency tables where only the margins are observed. The unobserved cells counts are assumed to have a Poisson distribution. Relying on the inverse Bayes formulae (IBF) combined with a Metropolis-Hasting algorithm, Bayes inference on the counting rate in each cell is obtained. This method called the
310
MH-IBF method is compared with the α-decomposition method in Zheng et al. [30]. It is shown that the two methods lead to similar results. However, the model behind the MH-IBF method has the advantage to be more transparent than those implied in the α-decomposition model, where a Bayesian hierarchical model has been applied with OpenBUGS. The approach suggested in this paper
315
is then an applicable contribution to establish efficient PRA models. In the case of nuclear plants safety where usually the samples of observed failures are small, our Bayesian method which consider prior information could be an appropriate strategy for a probabilistic risk analysis. Further works, e.g. choices of hyperparameters, construction credibility region are to be investigated.
320
Acknowledgment The authors would like to thank the anonymous reviewers and the associate editor for their valuable comments and suggestions that significantly improve the paper. References
325
[1] NASA, Probabilistic Risk Assessment, 2nd Edition, NASA/SP-2000–6112, 2000. [2] NASA, Probabilistic Risk Assessment Procedures Guide for NASA Managers and Practitioners, 2nd Edition, NASA/SP-2011–3441, 2011. [3] R. R. Fullwood, Probabilistic risk assessment in the nuclear power industry:
330
fundamentals and appllications, Pergamon Press, Toronto, 1988. 21
[4] M. Stewart, Probabilistic risk assessment of engineering systems, Chapman & Hall, London, UK, 1997. [5] H. Kumamoto, Satisfying safety goals by Probabilistic Risk Assessment, Springer series in reliability enginering, Springer, 2007. 335
[6] D. Kelly, C. Smith, Bayesian inference for probabilistic risk assessment, Springer series in reliability enginering, Springer, 2013. [7] L. Lu, J. Jiang, Probabilistic safety assessment for instrumentation and control systems in nuclear power plants: an overview, Journal of Nuclear Science and Technology 41 (3) (2004) 323–330.
340
[8] F. Safie, R. Stutts, Z. Huang, Reliability and probabilistic risk assessment – how they play together, 2015 annual reliability and maintainability, Symposium (RAMS) 67 (2) (2015) 1–5. [9] J. Nielsen, A. Tokuhiro, R. Hiromoto, J. Khatry, Optimization method to branch-and-bound large SBO state spaces under dynamic probabilistic risk
345
assessment via use of LENDIT scales and S2R2 sets, journal of nuclear science and technology 51 (10) (2014) 1212–1230. [10] H. Georges-William, M. Lee, E. Patelli, Probabilistic risk assessment of station blackouts in nuclear power plants, IEEE Transactions on Reliability 67 (2) (2018) 494–512.
350
[11] NUREG/CR-5485, Guidelines on modeling common-cause failures in probabilistic risk assessment, Tech. rep., U.S. Nuclear Regulatory Commission. (1998). [12] C. Atwood, The binomial failure rate common-cause model, Technometrics 28 (2) (1986) 139–148.
355
[13] C. Atwood, D. Kelly, The binomial failure rate common-cause model with winbugs, Reliability Engineering and System Safety 46 (2009) 990–999. [14] H. Nguyen, E. Gouno, Maximum likelihood and bayesian inference for common-cause of failure model, Reliability Engineering & System Safety 182 (2019) 56–62. 22
360
[15] NUREG/CR-4780, Procedures for treating common cause failures in safety and reliability studies, Tech. rep., US NUclear Regulatory Commission (1988). [16] A. A. O’Connor, A. Mosleh, Extending the alpha factor model for cause based treament of common cause failure events in pra and event assesment,
365
Proceedings of probabilistic fafety assessment and management conference, Honolulu, HI. [17] V. Hassija, C. S. Kumar, K. Velusamy, A pragmatic approach to estimate alpha factors for common cause failure analysis, Annals of Nuclear Energy 63 (2017) 317–325.
370
[18] H. Yu, J. Yang, J. Lin, Y. Zhao, Reliability evaluation of non-repairable phased-mission common bus systems with common cause failures, Computers & Industrial Engineering 11 (2017) 445–457. [19] M. Jinhua, L. Y. Feng, P. Weiwen, H. Zhong, Reliability analysis of complex multi-state system with common cause failure based on evidential networks,
375
Reliability Engineering & System Safety 174 (2018) 71–81. [20] J. Qin, Z. Li, Reliability and sensitivity analysis method for a multistate system with common cause failure,
Complexity,
DOI:
10.1155/2019/6535726. [21] A. O’Connor, A. Mosleh, A general cause based methodology for anaysis 380
of common cause and dependent failures in system risk and reliability assessments, Reliability Engineering and System Safety 145 (2016) 341–350. [22] U.S. Nuclear Regulatory Commission, CCF Parameter Estimations, 2012 Update (2013). [23] NUREG/CR-6819, Common-cause failure event insights, U.S. Nuclear Reg-
385
ulatory Commission 1-4. [24] NUREG/CR-6268, Common-cause failure database and analysis system: Event data collection, classification, and coding, Tech. rep., U.S. Nuclear Regulatory Commission. (2007).
23
[25] D. Rasmuson, D. Kelly, Common-cause failure analysis in event assessment, 390
Journal of Risk and Reliability 222 (4) (2008) 521–532. [26] D. Kelly, S. Shen, G. DeMoss, K. Coyne, D. Marksberry, Common-cause failure treatment in event assessment: basis for a proposed new model, Proceedings of Probabilistic Safety Assessment and Management (PSAM) 10.
395
[27] B. D. Johnston, A structured procedure for dependent failure analysis (DFA), Reliability Engineering 19 (1987) 125–136. [28] M. Zhang, Z. Zhang, A. Mosleh, S. Chen, Common cause failure model updating for risk monitoring in nuclear power plants based on alpha factor model, Reliability Engineering 231 (2017) 209–220.
400
[29] T. Sakurahara, G. Schumock, S. Reihani, E. Kee, Z. Mohaghegh, Simulation informed probabilistic methodology for common cause failure analysis, Reliability Engineering & System Safety 185 (2019) 84–99. [30] X. Zheng, A. Yamaguchi, T. Takata, α-decomposition method: A new approach to the analysis of common cause failure, Reliability Engineering
405
and System Safety 116 (2013) 20–27. [31] R. Little, D. Rubin, Statistical analysis with missing data, 2nd Edition, Wiley series in probability and statistics, Wiley-interscience, 2002. [32] M. Tan, G. Tian, K. Ng, Bayesian missing data problems: EM, data augmentation and noniterative computation, 1st Edition, Biostatistics Series,
410
Chapman & Hall, 2010. [33] K. Ng, M. Tang, G. Tian, M. Tan, The nested Dirichlet distribution and incomplete categorical data analysis, Statistical Sinica 19 (2009) 251–271. [34] K. Ng, M. Tang, M. Tan, G. Tian, Grouped Dirichlet distribution: A new tool for incomplete categorical data analysis, Jounal of multivariate analysis
415
99 (2008) 490–509. [35] G. Tian, M. Tang, K. Yuen, K. Ng, Futher properties and new applications of the nested Dirichlet distribution, Computational statistics and data analysis 54 (2) (2010) 394–405. 24
[36] G. Tian, K. Ng, Z. Geng, On improved EM algorithm and confidence in420
terval constrution for incomplete r × c tables, Computational statistics and data analysis 51 (2007) 2919–2933. [37] K. Ng, Inversion of Bayes formula: explicit formulae for unconditional pdf, In Advance in the Theory and Practice in Statistics, Wiley, New York. 19 (1997) 571–584.
425
[38] G. Tian, M. Tan, K. Ng, An exact non-iterative sampling procedure for discrete missing data problems, Statistical Neerlandica 61 (2) (2007) 232– 242.
25
Declaration of interests
☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
Journal Pre-proof
Bayesian inference for Common cause failure rate based on causal inference with missing data H.D. Nguyen, E. Gouno PII: DOI: Reference:
S0951-8320(19)30879-8 https://doi.org/10.1016/j.ress.2019.106789 RESS 106789
To appear in:
Reliability Engineering and System Safety
Received date: Revised date: Accepted date:
5 July 2019 14 November 2019 28 December 2019
Please cite this article as: H.D. Nguyen, E. Gouno, Bayesian inference for Common cause failure rate based on causal inference with missing data, Reliability Engineering and System Safety (2019), doi: https://doi.org/10.1016/j.ress.2019.106789
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.
Highlights • Introducing a new methodology using MCMC to estimate parameters of a common cause failure model. • Providing an alternative to deal with missing data in contingency tables. • Performing causality-based inference for common cause of failure. • Constructing strategies to elicit hyperparameters for a Bayesian approach. • Proposing a stochastic version of the Inverse Bayes Formulae algorithm.
1
Bayesian inference for Common cause failure rate based on causal inference with missing data Nguyen, H. D.a,b , Gouno, E.c, a Institute
of Research and Development, Duy Tan University, Danang 550000, Vietnam of Mathematics, Faculty of Information Technology, Vietnam National University of Agriculture, Hanoi, Viet Nam c University of South Brittany, LMBA, Campus de Tohannic, 56017 Vannes, France
b Department
Abstract This paper proposes a methodology to handle the causality to make inference on common cause failure (CCF) in a situation of missing data. The data are collected in the form of contingency table but the available information are only the numbers of CCF of different orders and the numbers of failure due to a given cause, i.e. the margins of the contingency table. The frequencies in each cell are unknown; we are in a situation of missing data. Assuming a Poisson model for the count, we suggest a Bayesian approach and we use the inverse Bayes formula (IBF) combined with a Metropolis-Hastings algorithm to make inference on the parameters. The performance of the resulting algorithm is evaluated through simulations. A comparison is made by analogy with results obtained from the recently proposed α-composition method. Keywords: Common-cause failure, Inverse Bayesian formula, Contingency table, Missing data, Causal inference.
1. Introduction Probability risk assessment (PRA) studies are performed to identify and understand vulnerabilities of systems. They were initially introduced in the nuclear industry. PRA is an analytical tool that allows to track potential system failures 5
and to determine the likelihood and consequences of their occurrences. This tool is now applied in many aeras and the literature on the subject is important. The NASA published a 83 pages bibliography in 2000 [1] and a procedures guide for ∗ Corresponding
author Email address: [email protected] (Gouno, E. )
Preprint submitted to Reliability Engineering & System Safety
January 3, 2020
manager and practitioners in 2011 [2]. Many books related to applications of PRA on reliability engineering problems have been published [3, 4, 5]. Kelly 10
and Smith consider a Bayesian approach of PRA [6]. A survey of application of PRA for instrumentation and control systems in nuclear power plants is given by Lu and Jiang [7]. Safie et al. [8] discuss the role of reliability and PRA at NASA. Pointing out that traditional PRAs fail to capture the dynamic behaviour of complex systems, Nielsen et al. [9] consider a dynamic PRA (DPRA)
15
and propose an optimization methods to deal with a combinatorial explosion of scenarios. More recently, Georges-William et al.[10] provide a probabilistic risk assessment of station blackout in nuclear power plants. They use an approach based on a hybrid of Monte Carlo methods, multistate modeling, and network theory as an alternative to the use of static fault tree which turns out
20
to deteriorate in applicability with increasing system size and complexity. In PRA, a common cause failure refers to “the failure or unavailable state of more than one component during the mission time and due to the same shared cause” [11]. In industry where a high reliability is required, these events are crucial challenges. The study of CCF events then plays a very important
25
role in system reliability assessments. A number of models has been developed to analyze CCF events such as the Binomial failure rate model [12, 13, 14], the α-factor model [15, 16, 17], the Beta factor model, and the multiple Greek letter model [15]. The CCF is also a major concern in several recent studies, for example, see [18, 19, 20]. However, in general, many previous studies
30
deal essentially with the magnitude of CCF events, i.e. the number of components failing simultaneously provoking the failure of the system, but they do not consider the cause of those events. This leads to some limitations for their use as discussed in [21]. Causality-based inference needs to be investigated and the U.S Nuclear Regulatory Commission has endeavored to develop
35
a database for this problem [22]. The Commission report [23] gives somes examples of possible CCF causes. For example, collecting the information from the emergency diesel generators in a set of nuclear plants, the report points out that amongst the 41 CCF events recorded, 15 events are caused by Design/Construction/Installation/Manufacture proximate cause group; 9 events
40
are qualified Internal to Component proximate cause group; 9 events are the
3
consequence of cause group related to Operational/Human Error, 5 events are the result of cause groups related to External Environment and 3 events are caused by Other proximate cause group. The report [24] provides a way to achieve further understanding of the occurrence of CCF events by presenting 45
the process of data collection and grouping the ranking of proximate failure causes. Rasmuson et al. [25] quantify the conditional probability of CCF considering the asymmetry within a common-cause component group, given independent failures or failures with common-cause potential. Using Bayesian networks,
50
Kelly et al. [26] are the first to propose the preliminary framework of a causalitybased inference, providing cause-specific quantitative insights into likely causes of failures. O’Connor and Mosleh [21] develop a general cause based methodology for analysis of common cause and dependent failures in system risk and reliability assessments. Johnston [27] presents a structured procedure to han-
55
dle dependent failure and their causes. A CCF model updating is presented in [28] where two updating processes are proposed considering the difficulty to identify failure causes immediately. Sakurahara et al. [29] utilize simulation models of physical failure mechanisms to capture underlying causalities and to generate simulation-based data for the CCF probability estimation. Zheng et al.
60
[30] suggest the α-decomposition method. In this decomposition, the authors consider a set of CCF data involving the information about not only the magnitude of external causes but also the cause occurrence frequencies of failures. This data forms a contingency table: the rows of the table record the number of components involved in each CCF event while its columns record the cause of
65
the event. As a result, the value in each cell of the table represents the number of failures due to a specific cause. In many situations it is not possible to observe complete contingency tables. We have missing data. Missing data are common in many data-driven approaches and the litterature on their statistical analysis grows up since the
70
early 1970s. Little and Rubin [31] offer a complete survey of methodologies for handling missing data problem. In particular, they provide a comprehensive study of methods for treating contingency table and they describe imputation methods in details. Tan et al. [32] present Bayesian solutions to missing data
4
problems using methods relying on EM-type algorithms, data augmentation and 75
noniterative computation. In the problem addressed here, one can not observe exactly the number of observations in each cell. Instead, only the values in the margins of the table are observed. Incomplete contingency table arises in many areas and a number of methods has been proposed in the literature. Ng et al. introduce new families of Dirichlet distributions [33, 34, 35, 36] and they use the
80
inverse Bayes formula (IBF) [37] to obtain closed-form solutions for incomplete categorical data Bayesian analysis. Several applications of the IBF can be found in [32]. The method presented in this paper is a contribution to probabilistic risk assessment where the causality of failure is taken into account. One of the key
85
points of the method is to enable the analysis of crude CCF data. That is to say that time and cost can be saved in the sense that no special procedures for collecting or correcting (completing) the data are required. We consider a scenario of CCF data mentioned in the paper of Zheng et al. [30]. The CCF data used in their paper are for illustration only and not from any real
90
database. However, one can easily imagine practical situations where for some reasons, number of components involved in the failure and cause of the failure are recorded independently. Therefore the CCF data do not contain the counting number of CCF events classified by magnitude and possible cause. We suggest a Bayesian approach to estimate the cell probabilities in this context of missing
95
data. An exact non-iterative sampling procedure proposed by Tian et al. [38] is developed. This procedure relies on the IBF. To overcome some practical difficulties specific to our problem, we use a Metropolis-Hasting (MH) algorithm to generate from conditional probability expressed with the IBF. Numerical results show the efficiency of the method. We show that the α-decomposition
100
parameters estimates are close to estimates obtained with our proposed method. The paper is organized as follows. In section 2, we present two situations of data and the corresponding distributions. The Bayesian inference for parameters of the model in the case of incomplete data using the IBF method is given in section 3. Section 4 presents a stochastic version of the IBF sampling. A
105
simulation study is investigated in section 5. Section 6 shows an application of the proposed method to the dataset provided in [30]. Finally, some concluding
5
remarks are given in section 7. 2. Data and distributions Let us consider a m × k contingency table with cell Nij . In the contex of 110
CCF, considering a m-component system observed in a time window of length T (unit of time) and k possible causes of failure C1 , . . . , Ck , Nij is a random variable corresponding to the number of CCF of order i from cause j occurring in [0, T ].
The matrix Nij , i = 1, . . . , m, j = 1, . . . , k that we denote N, represents
complete data. Each cell represents a number of failures occurring in a given time window. This number depending on time is a counting process. A natural way to model this number is to use a homogeneous Poisson process that is to say a Poisson process with a constant intensity. The consequence is then that in a given time window, the number of events follows a Poisson distribution. Then we suppose that the cell counts are independent homogeneous Poisson processes with respective intensity λij . The expected number of CCF of order i from cause j is therefore λij T . When Nij is observed, estimates for the λij ’s are easily obtained. Suppose now that the observed data are the total margins of the contingency table that is to say: Ni,• =
k X
Nij ,
j=1
the sum of the terms on the row i, corresponding to the total number of CCF events of order i, i = 1, ...., m, and N•j =
m X
Nij ,
i=1
the sum of the terms on the column j, corresponding to the total number of 115
components failed from cause j, j = 1, ...., k. In this situation, the cell counts Nij are not directly observed. We have missing data; the observed data are N1• , . . . , Nm• , N•1 , . . . , N•k , that we denote N• . We refer to N• as incomplete data and to N as complete data. Data displayed by Zheng et al. in table 1, p.24 in [30], are of this form. To
120
each line (system) of this table corresponds a contingency table where only the margins are observed. 6
Our purpose is to estimate the m × k matrix of terms λij , denoted by λ , in this situation of missing data. The distribution of N, the complete data, is obtained as a product of Poisson 125
distributions:
Y (λij T )nij
P (N = n | λ ) =
nij !
i,j
e−λij T .
(1)
Thus, the distribution of the incomplete data can be expressed as: P (N• = n• | λ ) =
X
n∈Dn•
P (N = n | λ ),
(2)
Pm where Dn• is the set of all the m×k matrices n = nij such that i=1 nij = n•j Pk and j=1 nij = ni• . In the following we aim at developing a Bayesian approach to estimate λ .
130
3. Bayesian inference: the IBF method Ng [37] introduces a method to compute posterior distributions based on the inverse Bayes formula: the IBF method. We are going to apply this method in our specific context of CCF. When the size of Dn• is not too large, the method leads to exact solutions. When it is not the case, we suggest a MH algorithm. λ), we want to compute the posterior distriAssuming a prior distribution π(λ λ | N• = n• ). Applying the Bayes theorem and the formula (2), we bution π(λ have: λ | N• = n• ) = π(λ
P
n∈Dn•
λ) P (N = n | λ ) π(λ
P (N• = n• )
.
But λ) = π(λ λ | N = n)P (N = n) P (N = n | λ ) π(λ and P (N = n)/P (N• = n• ) = P (N = n | N• = n• ). Therefore λ | N• = n• ) = π(λ
X
n∈Dn•
λ | N = n) P (N = n | N• = n• ), π(λ
(3)
We assume that the terms of λ are independent and each λij has a Gamma distribution with parameters (αij , βij ) which is a conjugate prior for the Poisson
7
distribution. The prior distribution of λ is then a product of Gamma distribu135
tions and the posterior distribution given the complete data is a product of Gamma distributions with parameters (nij + αij , T + βij ). That is, for n ∈ Dn• , Y
λ | N = n) ∝ π(λ
n +αij −1 −λij (T +βij )
λijij
e
.
(4)
i,j
We need now to compute P (N = n | N• = n• ). Let λ ∗ be an arbitrary value of λ. We have (inverse Bayes formulas): P (N = n) = and P (N• = n• ) = 140
Therefore
λ∗ ) P (N = n | λ ∗ ) π(λ ∗ λ | N = n) π(λ λ∗ ) P (N• = n• | λ ∗ ) π(λ . λ∗ | N• = n• ) π(λ
λ∗ | N• = n• ) P (N = n | λ ∗ ) π(λ P (N = n) = . λ∗ | N = n) P (N• = n• | λ ∗ ) P (N• = n• ) π(λ
(5)
Doing the summation on both sides with respect to n ∈ Dn• , we deduce: X P (N = n | λ ∗ ) P (N• = n• | λ ∗ ) = . λ∗ | N = n) λ∗ | N• = n• ) π(λ π(λ
n∈Dn•
And replacing in (5), we obtain λ∗ | N = n) P (N = n | λ ∗ )/π(λ . ∗ λ∗ | N = u) u∈Dn P (N = u | λ )/π(λ
P (N = n | N• = n• ) = P
•
Using expressions (1) and (4), we have: P (N = n | λ∗ ) Y T nij Γ(nij + αij ) ∝ = q(n) π(λ∗ | N = n) nij ! (T + βij )nij +αij i,j
(6)
n o Let us denote Dn• = n1 , . . . , nK and p` = P (N = n` | N• = n• ), that is q(n` ) , ` = 1, . . . , K. p` = PK j=1 q(nj )
(7)
Relying on this notation, we can express the posterior marginal distribution of the terms of λ as: π(λij | N• = n• ) =
K X `=1
8
p` π(λij | N = n` ),
(8)
145
which is a mixture of Gamma distribution. ˜ ij of λij are the expectation Under a quadratic loss, the Bayes estimators λ of the posterior distribution. Therefore, ˜ ij = λ
K X `=1
p`
hn
+ αij i . T + βij
`,ij
(9)
The following example illustrates the method when the size of Dn• is not too large, thus p = (p1 , . . . , p` ) can be easily computed. 150
Consider a set of CCF data from a 2-component system observed over a time window of length T = 10 unit of time. Suppose that the potential causes are divided into three main groups, say “cause 1”, “cause 2” and “cause 3” and that we have the margins, or observed data n• = {(3, 5, 2), (4, 6)}. The space Dn• can be 0 2 Dn• = 3 3 1 3 2 2 3 1 0 4
155
described as: 0 4 0 1 1 2 1 2 1 0 3 1 2 , , , , , 3 1 2 2 4 0 2 3 1 3 2 1 0 0 2 0 2 2 1 1 2 2 0 3 0 1 , , , , , 2 1 5 0 1 4 1 1 3 2 0 5 1 0 . 2
Let us consider a noninformative distribution for each λij , i.e., λij ∼ G (0.5, 0).
The probabilies p` = P (N = n` | N• = n• ), n` ∈ Dn• , are calculated from (7) and displayed in Table 1. `
p`
`
p`
1
0.10608204
7
0.13366337
2
0.07072136
8
0.04950495
3
0.12376238
9
0.06364922
4
0.07425743
10
0.14851485
5
0.04243281
11
0.12376238
6
0.06364922
Table 1: The conditional distribution P (N = n | N• = n• )
9
Applying (9), we obtain the following Bayes estimate for λ 0.19908 0.20700 0.14391 . λ˜ = 0.20092 0.39300 0.15608 4. A stochastic version of the IBF sampling When the space Dn• is too large the procedure described above is not pos160
sible since the computation of p relies on the size of Dn• . Sampling strategies are then considered. As mentioned previously, the incomplete-data marginal posterior distribution is a mixture of gamma distribution. Therefore a simple scheme to sample from this distribution is 1. choose a component of the mixture with probability p; 2. generate the λij ’s with the complete-marginal posterior
165
distributions (8). When Dn• can be identified this procedure is the exact IBF sampling suggested by Tian et al. [38]. The key in our situation is therefore to draw components of the mixture with the probability p, the conditional predictive distribution P (N = n | N• = n• ). The components of this probability has a cumbersome denominator. The MH
170
algorithm is a well-known sampling algorithm to overcome this kind of difficulty. The algorithm belongs to the family of MCMC. The principle is to build a Markov chain which stationary distribution is the distribution of interest. The algorithm is described as follows. Let Q be a specified transition matrix from one element of the space to another. Starting from n(0) ∈ Dn• , after p
175
cycles of the following scheme, n(p) is available, do: 1. Generate a candidate n(∗) ∼ Q(. | n(p) ). 2. Compute
q(n(∗) ) ρ = min 1, . q(n(p) )
3. Set n(p+1) = n(∗) with probability ρ; else, set n(p+1) = n(p) . After a number of iterations of this scheme, the sequence obtained is a sample of the conditional predictive distribution. 180
Many choices for Q, the generating-candidate matrix, are possible. In our case, the following scheme is applied to generate the candidate:
10
1. Choose randomly a 2 × 2 submatrix of n(p) with probability 1/ a11 a12 A= a21 a22
k 2
m 2
, say
2. Draw x from a uniform distribution in (max{0, a•1 −a2• }, . . . , min{a1• , a•1 }), where a1• = a11 + a12 and a•1 = a11 + a21 . 3. Replace the corresponding matrix A in n(p) by the matrix B defined as x a1• − x B= a•1 − x a2• − a•1 + x where a2• = a21 + a22 .
185
That means, we generate a new candidate by randomly changing a 2 × 2 submatrix of the current sample. The principle in step 2 and step 3 is to make sure that sum of rows and sum of columns of these sub-matrices are the same. The construction of the generating-candidate matrix is specific to our situation but it has been inspired by classical approaches that consider random walks.
190
The figure 1 compares the exact probabilities p` in Table 1 and the approximated probabilities obtained with the MH algorithm. One can see that the approximation is accurate; the algorithm is efficient. The suggested algorithm is an exact IBF sampling where Dn• is randomly covered using a MH algorithm. To summarize, the algorithm that we call the
195
stochastic IBF sampling or MH IBF sampling is described as follows: 1. Generate M
i.i.d. samples n(1) , . . . , n(M ) from the MH algorithm de-
scribed above, λ | N = n(`) ), the complete-data posterior distribution 2. Generate λ (`) ∼ π(λ defined in (4). 200
λ(`) }`=1,...,M are i.i.d. samples from the posterior π(λ λ | N• = n• ). Then {λ λ | N• = Considering the data from the previous example, sampling from π(λ n• ) using the stochastic IBF sampler and computing the empirical mean, we obtain the following estimate for λ : 0.1978159 0.2069460 λ˜ = 0.2018113 0.3933634 11
0.1431974 , 0.1561672
0.30 0.20 0.15 0.10 0.00
0.05
Probability
0.25
Theoretical distribution Equilibrium distribution
1
2
3
4
5
6
7
8
9
10
11
Elements
Figure 1: Comparison of the theoretical distribution P (N | N• = n• ) and the equilib-
rium distribution of the MH samples for the simple case
which is close to the estimate obtained with the exact distribution p` , ` = 205
1 . . . , K. We observed that in practice, the sampler is not sensitive to the choice of the burn-in size; we have roughly a steady acceptance ratio for different combinations of sample size and burn-in size. Therefore we consider the value M = 100000 and a burn-in size 50000 for this study.
5. A simulation study 210
We evaluate the performance of the proposed method through simulations. The simulation procedure is performed as follows: 1. Assign an input value for λ = (λij ) and a time window T . 2. Simulate a number in each cell of the contingency table from a Poisson distribution with parameter λij T and compute the margins.
215
3. Use the stochastic IBF sampler to estimate λ based on the two margins computed in 2. (1)
4. Repeat steps 2-3 M times to obtain the estimates λ˜
12
, . . . , λ˜
(M )
.
Guess
n11
n12
n13
n21
n22
n23
n31
n32
n33
Total
Input
6
5
1
4
3
1
4
1
2
27
Situation 1
7
4
1
3
4
2
4
3
1
29
Situation 2
4
3
1
3
2
1
2
1
1
17
Situation 3
9
7
3
6
4
2
5
3
1
40
Table 2: The average number of CCF events in each cell of contingency table M 1 X ˜ (i) ˜ 5. Compute the mean λ = λ and mean square error (MSE). M i=1
For purpose of illustration, we consider a 3 × 3 contingency table. In our
experiment, setting T = 120 and an average number of CCF events in each cell of the contingency table, the matrix input λ is deduced from λij T , an approximation of the chosen average numbers (first row in Table 2). The parameter λ used for simulation 0.0500 λ = 0.0333 0.0333
is therefore: 0.0083 0.0250 0.0083 . 0.0083 0.0167 0.0416
To specify values for prior hyperparameters, we use the strategy suggested in 220
[14]. We consider three situations of guess at the average numbers of CCF events in the time window of length T presented in Table 2. The first situation corresponds to a guess at the total number of failures close to the input number. The second and the third one correspond to situations where the guess at the total number underestimates (resp. overestimates) the input number. Then,
225
the prior hyperparameters for each λij is defined by solving a simple system of two equations matching guesses and uncertainty with mean and variance of the Gamma distribution:
α /β = ηij , ij ij
α /β 2 = ρ2 η 2 , ij ij ij ij
(10)
where ηij represents the guess at λij and ρij is the corresponding uncertainty of the guess. The hyperparameters for each situation of the guesses and the 230
corresponding assumed uncertainty are given in Table 3. They are determined as follows. We consider two levels of uncertainty: low and high. We choose ρij from a uniform distribution in [0.05, 0.3] for the low level and [0.7, 0.95] for 13
Situation 1
Situation 2
Situation 3
ρ
Low
High
Low
High
Low
High
Type of prior
A
B
C
D
E
F
(α11 , β11 )
34.59
593.07
1.30
22.41
23.09
692.70
1.99
59.98
15.18
202.52
1.36
18.16
(α12 , β12 )
51.12
2045.17
1.46
58.47
33.27
1331.10
1.43
57.50
32.18
643.70
1.83
36.65
(α13 , β13 )
17.54
526.46
1.66
49.89
24.96
1497.81
1.77
106.33
21.74
521.89
1.58
37.99
(α21 , β21 )
12.15
364.56
1.18
35.32
33.79
1351.95
1.35
54.33
11.06
189.62
1.17
20.18
(α22 , β22 )
47.81
1434.53
1.12
33.62
34.70
2082.51
1.24
74.56
12.12
363.80
1.22
36.89
(α23 , β23 )
15.16
606.50
1.52
60.89
18.18
2182.50
1.14
136.99
20.16
806.68
1.53
61.43
(α31 , β31 )
89.46
10736.35
1.88
226.11
60.09
7211.36
1.34
160.85
10.64
425.81
1.15
46.23
(α32 , β32 )
12.45
747.06
1.41
84.72
32.23
3868.41
1.13
136.52
13.71
823.03
1.30
78.10
(α33 , β33 )
18.97
2276.96
1.75
210.21
21.65
2598.92
1.20
144.76
14.12
1694.50
1.32
158.26
Table 3: Hyperparameters of prior based on the guess in Table 2.
the high level of uncertainty. The type A (resp. B) corresponds to the guess in situation 1 with low (resp. high) level of uncertainty. The type C (resp. D) 235
corresponds to the guess in situation 2 with low (resp. high) level of uncertainty. The type E (resp. F) corresponds to the guess in situation 3 with low (resp. high) level uncertainty. The type G is the noninformative case, that is to say gamma distributions with parameters (0.5, 0). The corresponding estimates are presented in Table 4. In general, when
240
the guess is close to the input, the small value of ρ leads to better estimates. Otherwise, when the guess is far from input, the large value of ρ give better results. Moreover, the noninformative prior gives acceptable results with small MSE.
6. Application 245
Zheng et al. [30] simulate CCF events and corresponding causes for 16 systems with 3 components and consider the results displayed in Table 5. They propose a method: the α- decomposition method to make inference from these data. It relies on an extension of the classical α-factor model. The classical α-factor model [15] for a m-component system considers m parameters (α1 , . . . , αm ) to
250
characterize the CCF’s. αi refers to the probability that exactly i components are involved in a CCF event given that at least one component fails.
14
λ11
λ12
λ13
λ21
λ22
λ23
λ31
λ32
λ33
Input
Type A
Type B
Type C
Type D
Type E
Type F
Type G
0.05
0.05887
0.06104
0.03604
0.04546
0.06528
0.05336
0.0612
(8.93e-5)
(4.74e-4)
(2.02e-4)
(1.80e-4)
(2.78e-4)
(3.31e-4)
(6.91e-4)
0.03236
0.03038
0.02603
0.03424
0.04973
0.03807
0.0353
(9.86e-5)
(2.56e-4)
(2.45e-4)
(1.94e-4)
(9.34e-5)
(2.07e-4)
(2.89e-4)
0.00834
0.00858
0.00836
0.00912
0.02301
0.01802
0.0160
(3.44e-8)
(2.45e-6)
(1.03e-8)
(5.40e-6)
(2.18e-4)
(1.35e-4)
(1.01e-4)
0.02495
0.02525
0.02577
0.03113
0.04761
0.03835
0.03622
(7.06e-5)
(1.31e-4)
(5.84e-5)
(1.08e-4)
(20.8e-4)
(1.50e-4)
(2.62e-4)
0.03282
0.02728
0.01699
0.020275
0.03209
0.02394
0.02781
(6.20e-5)
(1.01e-4)
(6.43e-5)
(5.98e-5)
(3.43e-5)
(5.39e-5)
(1.14e-4)
0.01653
0.01607
0.00834
0.00903
0.01592
0.01330
0.01482
(6.87e-5)
(8.85e-5)
(3.09e-8)
(6.45e-6)
(5.82e-5)
(3.67e-5)
(7.31e-5)
0.03291
0.03085
0.01768
0.02418
0.04010
0.03421
0.03155
(6.76e-6)
(1.04e-4)
(2.45e-4)
(1.30e-6)
(5.24e-5)
(1.21e-4)
(2.00e-4)
0.02399
0.02172
0.00866
0.01104
0.02418
0.02134
0.02432
(2.47e-4)
(2.08e-4)
(2.32e-8)
(1.47e-5)
(2.52e-4)
(1.93e-4)
(3.25e-4)
0.00833
0.00832
0.00854
0.01058
0.00822
0.00775
0.01481
(6.94e-5)
(7.11e-5)
(6.61e-5)
(5.01e-5)
(7.14e-5)
(8.15e-5)
(3.52e-5)
0.0416
0.0083
0.0333
0.025
0.0083
0.0333
0.0083
0.0167
Table 4: Means of Bayesian estimates and MSE for various types of prior, 100000
simulations.
15
Common causes occurence
CCF event
Cause 1
Cause 2
Cause 3
Total
Order 1
Order 2
Order 3
System 1
32
28
67
127
113
11
3
System 2
17
78
11
106
98
7
1
System 3
18
19
50
87
73
9
5
System 4
29
6
31
66
53
5
8
System 5
7
33
10
50
45
4
1
System 6
15
9
17
41
33
3
5
System 7
12
15
7
34
32
2
0
System 8
2
22
7
31
29
2
0
System 9
7
4
11
22
20
2
0
System 10
10
8
3
21
20
2
0
System 11
3
6
10
19
16
2
1
System 12
7
3
6
16
14
1
1
System 13
3
5
7
15
13
1
1
System 14
5
3
7
15
12
1
2
System 15
4
5
2
11
9
1
1
System 16
1
6
2
9
7
1
1
Table 5: A hypothetical database of CCF events provided in [30]
In the α-decomposition model, the parameters αi are decomposed as: αi =
k X
C
αi j rj , i = 1, . . . , m,
(11)
j=1
C
where αi j is the conditional probability that, given a CCF occurs from cause Cj , exactly i components are involved and rj is the probability that, given a 255
CCF occurs, it is from cause Cj . C
The decomposed parameters αi j , i = 1, . . . , m; j = 1, . . . , k are the new parameters of interest. Zheng et al. [30] present a hierarchical model and develop a Bayesian methodology to estimate the decomposed α-factors. C
The parameters αi j which are probabilities, can be roughly approximated 260
by the ratio between the number of CCF of order i from cause j and the total number of CCF of any order from cause j. These numbers can be respectively P3 approximated by λij T and λ•j T , where λ•j = i=1 λij . Therefore:
C
αi j ≈
λij λij T = , i = 1, . . . , m, j = 1, . . . , k, λ•j T λ•j
(12)
This relationship allows us to compare the results obtained with MH-IBF method 16
265
and the α-decomposition method. Since the systems are of the same nature with Poisson counting process, the data considered to apply the MH-IBF method are: n1• = 587, n2• = 53, n3• = 530, n•1 = 172, n•2 = 250, n•3 = 248. That is to say, we aggregate the data from Table 5. With noninformative priors, assuming that the observation time
270
is T = 100 (unit of time), the results characterizing the posterior distribution are displayed in Table 6. ˜ ij , it can be seen Considering mean of each posterior distribution, that is λ that ˜ 11 < λ ˜ 13 < λ ˜ 12 , λ ˜ 22 < λ ˜ 21 < λ ˜ 23 , λ
(13)
˜ 32 < λ ˜ 33 < λ ˜ 31 . λ C
Thus, the order among λij ’s is the same as the order among αi j ’s in the α275
decomposition model. Similar to the discussion in [30], cause 2 is of least CCF risk, cause 3 is best at provoking partial CCF and cause 1 is best at provoking complete CCF. However, the impact of cause 2 and cause 3 are similar as the difference between them are not significant. Based on this result, the most hazardous causes could be recognized, which has a practical engineering meaning. Figures 2-4 show the smoothed histogram as a density estimates of the posterior densities of each λij based on M = 100000 i.i.d. samples generated using the sampling MH-IBF method (the histogram is displayed only for λ11 , λ21 , λ31
4
in each graph).
1
2
3
λ11 λ12 λ13
0
Probability density function
280
0.8
1.2
1.6
2.0
2.4
2.8
Figure 2: Posterior distributions of λ1j (j = 1, 2, 3)
17
Posterior λij
Mean
SD
0.025-quantile
Median
0.975-quantile
λ11
1.4454
0.1487
1.1625
1.4442
1.7461
λ12
2.2293
0.1718
1.9016
2.2269
2.5725
λ13
2.2095
0.1728
1.8777
2.2075
2.5541
λ21
0.1825
0.0906
0.0389
0.1735
0.3839
λ22
0.1820
0.0888
0.0397
0.1708
0.3761
λ23
0.1805
0.0902
0.0379
0.1703
0.3811
λ31
0.1069
0.0541
0.0213
0.0986
0.2272
λ32
0.1041
0.0536
0.0217
0.0973
0.2262
λ33
0.1042
0.0538
0.0211
0.0983
0.2260
CCF order 1
CCF order 2
CCF order 2
2
4
λ21 λ22 λ23
0
Probability density function
6
Table 6: Summary of posterior distributions of each λij , (i, j = 1, 2, 3).
0.0
0.1
0.2
0.3
0.4
0.5
0.6
Figure 3: Posterior distributions of λ2j (j = 1, 2, 3)
18
8 6 4 2 0
Probability density function
λ31 λ32 λ33
0.0
0.1
0.2
0.3
0.4
Figure 4: Posterior distributions of λ3j (j = 1, 2, 3)
Under a quadratic loss, the Bayes 1.4454 ˜ λ = 0.1825 0.1069 C
estimate of λ is: 2.2293 2.2095 0.1820 0.1805 . 0.1041 0.1042
The parameters αi j in α-factor model can be approximated using (12). 285
C
Table 7 compares the estimates of αi j , 0.025 and 0.975 quantiles obtained from the method suggested by Zheng et al. [30] and the aproximation obtained injecting the estimates of λij ’s in (12). Mean
Values
0.025-quantile
0.975-quantile
α-decomposition
MH-IBF
α-decomposition
MH-IBF
α-decomposition
MH-IBF
α1C1
0.813
0.833
0.603
0.713
0.961
0.932
α1C2
0.910
0.886
0.826
0.804
0.978
0.955
α1C3
0.820
0.885
0.687
0.798
0.944
0.953
α2C1
0.0988
0.1052
0.00389
0.0105
0.277
0.215
α2C2
0.0729
0.0723
0.00979
0.0147
0.153
0.147
α2C3
0.112
0.0723
0.0120
0.0162
0.221
0.202
α3C1
0.0879
0.0616
0.00456
0.0124
0.218
0.198
α3C2
0.0174
0.0414
0.00052
0.0015
0.0561
0.0894
α3C3
0.0679
0.0417
0.00540
0.0088
0.151
0.139
Cj
Table 7: Compare the estimates of αi
using the α-decomposition method and our
proposed method.
19
C
The order among αi j ’s when j is fixed (similar to (13)) is preserved. The figure 5 shows that the estimates and quantiles with both methods
290
are close and overlap with each other. In other words, with consideration of uncertainty bounds, the results of the MH-IBF method are consistent with the
0.8
1.0
results of the α-decomposition method.
0.0
0.2
0.4
0.6
α−decomposition 0.025−quantile for α−decomposition 0.975−quantile for α−decomposition MH−IBF 0.025−quantile for MH−IBF 0.975−quantile for MH−IBF
C
α1 1
C
α1 2
C
α1 3
C
C
α2 1
α2 2
Cj
Figure 5: Comparison of the parameters αi
C
α2 3
C
α3 1
C
α3 2
with their corresponding approximation
using MH-IBF method
295
The MH-IBF is an plausible alternative to the α-decomposition method. The MH-IBF is constructed considering the natural assumption that we have Poisson distributions for the counts and the parameters to be estimated are the expected numbers of CCF. The α-decomposition does not have this simple direct interpretation. Therefore in the context of PRA study, MH-IBF can be
300
more adequate.
20
C
α3 3
7. Concluding remarks This paper has developped a new method for estimating a Common Cause Failure model parameters based on missing data in the form of contingency table. We have presented a procedure to characterize the relationship between 305
potential cause and failure for system reliability. The data considered are of the form of contingency tables where only the margins are observed. The unobserved cells counts are assumed to have a Poisson distribution. Relying on the inverse Bayes formulae (IBF) combined with a Metropolis-Hasting algorithm, Bayes inference on the counting rate in each cell is obtained. This method called the
310
MH-IBF method is compared with the α-decomposition method in Zheng et al. [30]. It is shown that the two methods lead to similar results. However, the model behind the MH-IBF method has the advantage to be more transparent than those implied in the α-decomposition model, where a Bayesian hierarchical model has been applied with OpenBUGS. The approach suggested in this paper
315
is then an applicable contribution to establish efficient PRA models. In the case of nuclear plants safety where usually the samples of observed failures are small, our Bayesian method which consider prior information could be an appropriate strategy for a probabilistic risk analysis. Further works, e.g. choices of hyperparameters, construction credibility region are to be investigated.
320
Acknowledgment The authors would like to thank the anonymous reviewers and the associate editor for their valuable comments and suggestions that significantly improve the paper. References
325
[1] NASA, Probabilistic Risk Assessment, 2nd Edition, NASA/SP-2000–6112, 2000. [2] NASA, Probabilistic Risk Assessment Procedures Guide for NASA Managers and Practitioners, 2nd Edition, NASA/SP-2011–3441, 2011. [3] R. R. Fullwood, Probabilistic risk assessment in the nuclear power industry:
330
fundamentals and appllications, Pergamon Press, Toronto, 1988. 21
[4] M. Stewart, Probabilistic risk assessment of engineering systems, Chapman & Hall, London, UK, 1997. [5] H. Kumamoto, Satisfying safety goals by Probabilistic Risk Assessment, Springer series in reliability enginering, Springer, 2007. 335
[6] D. Kelly, C. Smith, Bayesian inference for probabilistic risk assessment, Springer series in reliability enginering, Springer, 2013. [7] L. Lu, J. Jiang, Probabilistic safety assessment for instrumentation and control systems in nuclear power plants: an overview, Journal of Nuclear Science and Technology 41 (3) (2004) 323–330.
340
[8] F. Safie, R. Stutts, Z. Huang, Reliability and probabilistic risk assessment – how they play together, 2015 annual reliability and maintainability, Symposium (RAMS) 67 (2) (2015) 1–5. [9] J. Nielsen, A. Tokuhiro, R. Hiromoto, J. Khatry, Optimization method to branch-and-bound large SBO state spaces under dynamic probabilistic risk
345
assessment via use of LENDIT scales and S2R2 sets, journal of nuclear science and technology 51 (10) (2014) 1212–1230. [10] H. Georges-William, M. Lee, E. Patelli, Probabilistic risk assessment of station blackouts in nuclear power plants, IEEE Transactions on Reliability 67 (2) (2018) 494–512.
350
[11] NUREG/CR-5485, Guidelines on modeling common-cause failures in probabilistic risk assessment, Tech. rep., U.S. Nuclear Regulatory Commission. (1998). [12] C. Atwood, The binomial failure rate common-cause model, Technometrics 28 (2) (1986) 139–148.
355
[13] C. Atwood, D. Kelly, The binomial failure rate common-cause model with winbugs, Reliability Engineering and System Safety 46 (2009) 990–999. [14] H. Nguyen, E. Gouno, Maximum likelihood and bayesian inference for common-cause of failure model, Reliability Engineering & System Safety 182 (2019) 56–62. 22
360
[15] NUREG/CR-4780, Procedures for treating common cause failures in safety and reliability studies, Tech. rep., US NUclear Regulatory Commission (1988). [16] A. A. O’Connor, A. Mosleh, Extending the alpha factor model for cause based treament of common cause failure events in pra and event assesment,
365
Proceedings of probabilistic fafety assessment and management conference, Honolulu, HI. [17] V. Hassija, C. S. Kumar, K. Velusamy, A pragmatic approach to estimate alpha factors for common cause failure analysis, Annals of Nuclear Energy 63 (2017) 317–325.
370
[18] H. Yu, J. Yang, J. Lin, Y. Zhao, Reliability evaluation of non-repairable phased-mission common bus systems with common cause failures, Computers & Industrial Engineering 11 (2017) 445–457. [19] M. Jinhua, L. Y. Feng, P. Weiwen, H. Zhong, Reliability analysis of complex multi-state system with common cause failure based on evidential networks,
375
Reliability Engineering & System Safety 174 (2018) 71–81. [20] J. Qin, Z. Li, Reliability and sensitivity analysis method for a multistate system with common cause failure,
Complexity,
DOI:
10.1155/2019/6535726. [21] A. O’Connor, A. Mosleh, A general cause based methodology for anaysis 380
of common cause and dependent failures in system risk and reliability assessments, Reliability Engineering and System Safety 145 (2016) 341–350. [22] U.S. Nuclear Regulatory Commission, CCF Parameter Estimations, 2012 Update (2013). [23] NUREG/CR-6819, Common-cause failure event insights, U.S. Nuclear Reg-
385
ulatory Commission 1-4. [24] NUREG/CR-6268, Common-cause failure database and analysis system: Event data collection, classification, and coding, Tech. rep., U.S. Nuclear Regulatory Commission. (2007).
23
[25] D. Rasmuson, D. Kelly, Common-cause failure analysis in event assessment, 390
Journal of Risk and Reliability 222 (4) (2008) 521–532. [26] D. Kelly, S. Shen, G. DeMoss, K. Coyne, D. Marksberry, Common-cause failure treatment in event assessment: basis for a proposed new model, Proceedings of Probabilistic Safety Assessment and Management (PSAM) 10.
395
[27] B. D. Johnston, A structured procedure for dependent failure analysis (DFA), Reliability Engineering 19 (1987) 125–136. [28] M. Zhang, Z. Zhang, A. Mosleh, S. Chen, Common cause failure model updating for risk monitoring in nuclear power plants based on alpha factor model, Reliability Engineering 231 (2017) 209–220.
400
[29] T. Sakurahara, G. Schumock, S. Reihani, E. Kee, Z. Mohaghegh, Simulation informed probabilistic methodology for common cause failure analysis, Reliability Engineering & System Safety 185 (2019) 84–99. [30] X. Zheng, A. Yamaguchi, T. Takata, α-decomposition method: A new approach to the analysis of common cause failure, Reliability Engineering
405
and System Safety 116 (2013) 20–27. [31] R. Little, D. Rubin, Statistical analysis with missing data, 2nd Edition, Wiley series in probability and statistics, Wiley-interscience, 2002. [32] M. Tan, G. Tian, K. Ng, Bayesian missing data problems: EM, data augmentation and noniterative computation, 1st Edition, Biostatistics Series,
410
Chapman & Hall, 2010. [33] K. Ng, M. Tang, G. Tian, M. Tan, The nested Dirichlet distribution and incomplete categorical data analysis, Statistical Sinica 19 (2009) 251–271. [34] K. Ng, M. Tang, M. Tan, G. Tian, Grouped Dirichlet distribution: A new tool for incomplete categorical data analysis, Jounal of multivariate analysis
415
99 (2008) 490–509. [35] G. Tian, M. Tang, K. Yuen, K. Ng, Futher properties and new applications of the nested Dirichlet distribution, Computational statistics and data analysis 54 (2) (2010) 394–405. 24
[36] G. Tian, K. Ng, Z. Geng, On improved EM algorithm and confidence in420
terval constrution for incomplete r × c tables, Computational statistics and data analysis 51 (2007) 2919–2933. [37] K. Ng, Inversion of Bayes formula: explicit formulae for unconditional pdf, In Advance in the Theory and Practice in Statistics, Wiley, New York. 19 (1997) 571–584.
425
[38] G. Tian, M. Tan, K. Ng, An exact non-iterative sampling procedure for discrete missing data problems, Statistical Neerlandica 61 (2) (2007) 232– 242.
25
Declaration of interests
☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: