Resampling Approach to the Estimation of Reliability Systems

# Resampling Approach to the Estimation of Reliability Systems

Maxim Fioshin, Helen Fioshina
Institute of Transport Machines and Technologies,
Riga Technical University,
Lomonosova Str. 1, LV1019, Riga
Latvia
Maksims.Fiosins@rtu.lv, Jelena.Fiosina@rtu.lv
###### Abstract

The article is devoted to the resampling approach application to the reliability problems. This approach to reliability problems was first proposed by Ivnitsky (1967). Resampling is intensive statistical computer method, which is non-parametrical, that uses initial samples data in different combinations to simulate the process many times and get finally the estimator of the characteristics of interest. At the present paper simple resampling, hierarchical resampling, the case of one sample for several variables, the case of partially known distributions, analysis of degradation flow, analysis of degradation-renewal process, construction of confidence intervals are described. All those resampling application cases can be applied successfully to solve the reliability problems as an alternative to classical methods.111This work has been partly supported by the European Social Fund within the National Programme ”Support for the carrying out doctoral study program’s and post-doctoral researches” project ”Support for the development of doctoral studies at Riga Technical University”

Resampling Approach to the Estimation of
Reliability Systems

Maxim Fioshin, Helen Fioshina Institute of Transport Machines and Technologies, Riga Technical University, Lomonosova Str. 1, LV1019, Riga Latvia Maksims.Fiosins@rtu.lv, Jelena.Fiosina@rtu.lv

## Introduction

Computers play an important role in the development of modern statistical science. Many statistical methods developed in the last time require a big amount of computations. Such methods are usually called intensive statistical methods. Some authors speak about separate discipline called ”computational statistics” [Gentle_02].

In statistical problems of reliability we often have difficulties with application of classical methods. In the case of complex systems, small samples, unknown distributions of system characteristics classical methods work not so good. In this case, methods of computational statistics can be used as an alternative.

The field of computational statistics at present time includes a big number of methods. Jackknife method was suggested by M. Quenouille in 1949 [Quen_49]. It uses the estimator, which is the combination of the estimator’s obtained using all data and the estimators, obtained using only part of the same data. In 1967 V. Ivnitsky [Ivn_67] suggested to use resampling for the estimation of the systems’ reliability by simulation. B. Efron [Efron_79_2] suggested bootstrap method, which is the generalization of the jackknife method. Resampling approach was investigated from 1995 under prof. A. Andronov supervision. During this investigation, simple and hierarchical resampling [Andr_95] and their implementations in reliability theory, queuing theory [AndrF_99_2, Afan_02], stochastic processes [Andr_00_2, Afan_05_1, Andr_06], optimization tasks [Andr_00] and construction of confidence intervals [Andr_02, AndrF_04] were considered. The present paper is devoted to the description of the main results connected with the reliability problems.

The classical ”plug-in” approach consists in the choosing of the forms of unknown distributions and estimating their unknown parameters on the base of available sample populations. Then the estimated distributions are used in formulas for calculation of unknown parameters instead of real values or they are used for generation of the pseudo-random numbers, which are used in simulation instead of latter variables. Here the estimation of probability distributions, espetially in the case of small initial samples, leads to making mistakes in choosing the form of a distribution and estimating its parameters. It can lead to the bias of estimated parameters and results of the simulation.

Resampling is an alternative approach. We do not estimate unknown distributions, but extract values from initial sample populations at random and use them as values of input variables during simulation. This method does not require any preliminary information about the type of the distributions, so it is non-parametrical. It uses the available information in different combinations, which require a big number of computations, but allows to get good estimators in situations, where classical estimators are not good.

The first section of the present article is devoted to the main principles of the resampling approach. The resampling algorithm and the principals of the resampling estimator calculation are described. The variance of the estimator is supposed to be criterion of obtained estimator efficiency. In the case of biased estimators it is better to take the mean squared error as the criterion of efficiency, which is calculated on the base of variance and bias. The variance is calculated on the base of resampling pairs (-pair, -pair, -pair). The resampling pair shows the common elements structure of two different resamples, so the type of the pair depends on the problem. In this section the resampling is applied to the estimation of the reliability function of logical schema. Different variants of this method are described: simple resampling, hierarchical resampling, resampling of one sample for several distributions and the case of partially known distributions.

In the second part the resampling approach implementation to stochastic processes analysis is described. Two problems are considered there. The first problem is analysis of the degradation process with the accumulation of damages. The second problem is the comparison of two renewal processes implementing to the analysis of degradation-renewal process analysis.

The last part of the article is devoted to resampling application to the interval estimation of logical system and the algorithm of true coverage probability calculation of the constructed confidence interval.

Finally the concluding remarks are made. All the parts contain only the brief description of the considered problems and the main results. For more detailed description please refer to a bibliography, which is presented at the end of the article.

## 1 Resampling estimation of reliability function of logical schemes

### 1.1 Simple resampling

Let we have a known function of real arguments. This function can represent some characteristics of -element logical system with working times of elements (for example, the total working time of the system, the indicator that the system works at the given time etc.).

Let we have a vector of independent r.v. X=(,,…,), each component with unknown cdf ; only sample populations H=(,,…,), are available for each , . Let our parameter of interest is the expectation of the function , which argument is vector X:

 Θ=Eϕ(X). (1)

We have to estimate on the base of given sample populations H.

The classical nonparametric plug-in approach is the following. We calculate the empirical cdf , =1,2,…, on the base of samples H, take formulas for calculation of (which include ) and use instead of . This gives formula for plug-in estimator of . Usually this estimator has good properties, but in some cases (the same samples for the same variables, small samples, etc.) it can have big variance or be biased.

The resampling approach is nonparametric approach to systems simulation and estimation [Andr_00]. It is a simple alternative to classical methods. We perform the following iterative procedure. On the -th step we extract at random one element from each population , . Let be the indices of extracted elements in the corresponding populations H, so . The estimator of is calculated using vector as argument of the function :

 Θ∗q=Eϕ(X∗q). (2)

The procedure is repeated times, and as the estimator of an average of all is taken:

 Θ∗=1rr∑q=1Θ∗q. (3)

The obtained estimator (3) is unbiased:

 EΘ∗=Θ. (4)

So, in order to provide proper application of resampling method, we should know other properties of . In most cases estimators’ variance can be taken as efficiency criterion. Although, in some situations the estimator is biased and then bias and mean squared error also should be calculated.

The variance of the resampling estimator is:

 VarΘ∗=1rVarΘ∗q+r−1rCov(Θ∗q,Θ∗q′)=1rμ2+r−1rμ11−μ2,q≠q′, (5)

where , , .

In this formula, the variance , the second moment and the expectation depend only on the properties of the function and r.v. X. Only the covariance and the mixed moment depend on the applied resampling procedure. So our goal is the calculation of the variance (5) depending on the used resampling procedure.

In order to calculate the covariance or the mixed moment , we introduce the notation of the -pair.

Let . We say that two vectors and produce the -pair, , when , if and , if . So, the -pair shows indices of the arguments which have the same elements from the initial sample populations in two different realizations and .

For example, let we have function and the resample vectors j and j. In this case, they will form -pair.

Let us denote the event ”-pair is happened”. Let = be the probability of this event. We can calculate the covariance or the mixed given .

Let , be conditional covariance given -pair. Let , be conditional mixed moment given -pair. Then can be calculated as follows:

 Cov(Θ∗q,Θ∗q′)=∑\boldmathω∈MP{\boldmathω}C(\boldmathω). (6)

The mixed moment can be calculated as follows:

 μ11=∑\boldmathω∈MP{\boldmathω}μ11(\boldmathω). (7)

Now let us show how to obtain the probabilities . As the arguments of the function are independent and the probability to extract the same element from twice on the -th realization and on the -th realization is , the probability can be calculated as follows:

 P{\boldmathω}=∏i∈\boldmathω1ni∏i∉\boldmathω(1−1ni). (8)

The construction of all -pairs is simple combinatorial problem. By constructing all the -pairs we can calculate the variance (5), using formulas (6)-(7).

As an example, let us consider a system ”2 of 3” [Fio_03]. It consists on 3 elements, and it works if at least 2 of 3 elements work. The problem is to estimate the probability that at the time moment the system works.

The function of interest is the following indicator function :

 ϕt(x1,x2,x3)={1if at least 2 % elements of {x1,x2,x3} are greater than t,0otherwise.

Our aim is to estimate the probability that system works at time moment as the expectation of the function:

In our case the following 8 -pairs are possible:

The probabilities of the -pairs and the mixed moments can be easily calculated. This allows us to use the resampling approach and to calculate the variance of obtained estimator.

### 1.2 Hierarchical resampling

Often the simulated system is complex, but structured. It can be split onto subsystems, which can be simulated separately. Then, the results of subsystems simulation can be used to simulate the whole system. In this case we can use hierarchical resampling.

Simple resampling can also be used for analysis of hierarchical systems. We can extract, like in previous section, values for input data from initial samples, calculate the value of function of interest without paying attention to its hierarchical structure and get the simple resampling estimator. But hierarchical resampling has some advantages in the described situation [Andr_95] and allows: parallel calculations, sample sizes optimization, more clear analysis of the efficiency.

Let the function has hierarchical structure, i.e. it can be split into subfunctions . The results of subfunctions are used as arguments of the functions on the higher layer.

It is convenient to represent function by the calculation tree. The root of this tree has index and corresponds to the function . Nodes with indices correspond to input variables . The rest nodes are intermediate ones, which correspond to intermediate functions .

Let us denote the set of nodes, from which arcs go to the node , and - the set of initial variables, such that the node depends upon them. Note that for vertices of the same level the sets do not cross and the sets also do not cross.

Now, the function can be calculated by so-called ”wave” algorithm [Andr_95]. We create samples sequentially. The sample , is calculated by extracting values from corresponding samples and using them as arguments of the function . It is clear that

 Xv,q=ϕ(X∗qi|i∈Iv),q=1,2,…,nv. (9)

Finally, the estimator of can be calculated as an average of the sample elements at the root of the calculation tree:

 Θ∗=1nknk∑i=1Xk,i. (10)

Now our purpose is to calculate variance . It can be calculated by formulas (5)-(7), but using slightly generalized -pair notation.

Let us consider the value , calculated by formula (9). It was calculated using only one value from each sample , . Let us define these values indices in the initial samples by vector , . Some elements of vectors and can be equal. So, we can use almost the same definition of the -pair.

Let M={1,2,…,}. We will say that two vectors and produce the -pair, , when , if and , if .

Let be the event ”the -pair takes place at the node ”. Let be the probability of this event. The values of is calculated recurrently. If the node is on zero level , then can be calculated easily:

 Pv{\boldmathω}={0if \boldmathω≠∅,1otherwise. (11)

Now let us consider a node on another level. Let us consider a sample , from which the node depends. The probability that some element from is chosen twice: for and for is equal to .

The event can happen if in each node one of two events occurs:

• We selected different elements from for and (with the probability ) and selected elements produced - pair (with the probability );

• We selected the same element from for and (with the probability ). In this case all elements from will be extracted twice. The event can happen in this case only if all elements of belong to : .

Now let us denote :

 δi,\boldmathω={1ifIi0⊂\boldmathω,0otherwise. (12)

Now we can write formula for the calculation of :

 Pv{\boldmathω}=∏i∈Iv((1−1/ni)Pi{\boldmathω}+(1/ni)δi,\boldmathω). (13)

Our goal is to calculate for all . This will allow us to use formula (6) or (7) for calculation the value or . Then we will use (5) for calculation of the estimator variance.

Let us consider an example. Let we have a system that consists of 6 elements. The 1-st and the 2-nd elements are connected in parallel, the 3-d and the 4-th elements are connected sequentially, the 5-th and the 6-th elements are connected in parallel, but the 6-th element is switched on when the 5-th element fails (cold reservation). Our purpose is to estimate the probability that the system will work at the time moment .

Our function of interest can be denoted as follows:

 ϕt(x1,x2,x3,x4,x5,x6)={1,if min(max(x1,x2),min(x3,x4),x5+x6)

The function can be represented by using the calculation tree. The tree has 6 leaves, three elements on the first level correspond to the functions , and and the element on the second level (root of this tree) corresponds to the function .

In this case we have 4 -pairs in each node on the first level and 8 -pairs in the root of the tree. This allows us using the resampling approach and calculating the variance of obtained estimator.

### 1.3 The case of the only one sample for several r.v.

In this section we will show how to deal with in the case, when only one sample is available for several logical elements. This situation can often appear when elements are considered to have identical characteristics.

Let we have a function , but some arguments are considered to be statistically identical. We have only one sample for all identical arguments. Let the numeration of arguments corresponds to the order of sample numbers. So, let we have samples and the sample is used for arguments (and we put ). Let be a number of arguments, for which the sample is used, .

The resampling approach here is the following: on the -th realization we extract one value for each argument from the corresponding population. The values from one sample are extracted without replacement. The extracted values are used as arguments of the function , and the resampling estimator for -th realization is calculated by formula (2). Finally, the resampling estimation is calculated by formula (3). Note that obtained estimator is unbiased. Our goal is to show how to calculate its variance.

The variance is calculated by formula (5). In order to calculate the covariance we generalize the notation of -pair, introducing the -pair. The idea lying behind this definition is the following: in the previous cases, where -pair was used, the same element from initial sample in two different realizations and was possible only for the argument . In the present case the same element in two different realizations and is possible for different arguments, namely any of arguments can have the same element, because it is extracted from for all of them. The -pair shows what arguments have the same elements extracted.

Let be the indices of elements extracted on the -th step. Note that are indices in the sample and they are different, .

Let . We say that two vectors and produce the -pair when for all if and , then , otherwise .

For example, if we have a function and two samples and , the first sample is for arguments (this means ) and the second sample is for arguments (this means ). Now let we have resample vectors j and j. In this case, they will produce -pair.

Note that the -pair is the specific case of the -pair. -pair becomes -pair when only two variants are possible: or . Also, zeros in the -pair are not stored.

The variance of obtained estimator can be calculated by formula (5). In order to calculate the covariance or the mixed moment , we can use formula (6) or (7), where we use the probability of -pair, conditional covariance given -pair and the conditional mixed moment given -pair.

Now let us show how to calculate the probability of the -pair. The arguments can be split into independent blocks. The block corresponds to arguments which belong to the sample . Let be a number of non-zero elements at the block . Then the probability can be calculated using the hypergeometrical distribution

 P(\boldmathβ)=k∏i=1(miαi)(ni−mimi−αi)(nimi), (15)

Note that in the case when , it is impossible for two samples not to have common elements. For this case, let us put , if .

The calculation of depends on the function structure. Note that for the calculation of the probability we need only the information about a number of non-equal elements in each block. In general, the information about the indices of equal elements is used for a calculation of .

Now let us consider a specific case of the above described situation - the case when some elements are not simply equivalent, but their influence to the system work is equivalent [Fio_00, Fio_02]. It can happen in the situation, when the function is commutative by the block arguments, i.e. changing the order of arguments inside the block does not change the function result. Note that a reliability function, which includes identical elements, can often be commutative by these arguments, because can include sum, min, max etc. of these arguments.

Here, we do not need to store for each element of , what element of contains the same element of corresponding sample. We need only to know how much the same elements were selected for each block from the corresponding sample. So, instead of -element -pair we introduce the -element -pair.

Formally, we say that two vectors and produce the -pair if for all , . This means that -pair stores only the number of common elements inside each block. The probability can be calculated by formula (15).

For example, if we have a function and two samples and , the first sample is for arguments (this means ) and the second sample is for arguments (this means ). Now let we have resample vectors and . In this case, they will form -pair.

### 1.4 The case of partially known distributions

In the previous sections we supposed that distributions of all r.v. X={} are unknown, but only sample populations H are available for each variable. In many practical situations the distributions of some r.v. are known, but the distributions of the rest r.v. are unknown, and only sample populations H are available. The problem is, how to use the available information about r.v. distribution in the most efficient way [AndrF_99_3].

Let be a known function of independent r.v. X={} with unknown distributions, Z={} with known distributions: . The problem consists in estimation of the expectation

The idea of the estimation of is the following: we use the conditional expectation of given :

 g(X)=E(ϕ(X,Z)|X). (16)

It is clear that

Two situations are possible for the conditional expectation (: either has known functional form for any X or the functional form of is unknown.

In the first case we can operate as in usual simple resampling, but we estimate instead of . In each realization we create a resample and use it to estimate : . Then, the estimator can be obtained by formula (3).

Now we consider the case when the conditional expectation (16) is unknown. In this case we are able to estimate , because the distributions of random variables are known. On the -th step we generate mutually independent realizations of vector Z, and estimate in the following way:

 Θ∗q=1NN∑i=1ϕ(X∗q,Z∗q,i). (17)

Then the estimator can be obtained by formula (3).

Now let us see how hierarchical resampling can be applied in the case of partially known distributions. Let the function be calculated by a calculation tree. The leaves of the tree correspond to the variables . Note that the intermediate functions depend on the values of child nodes and the variables .

Instead of the conditional expectation , we must know the conditional distribution function of each given X: . As in the case of the simple resampling, two variants are possible here: either the distribution function can be calculated in each node or it is unknown.

Let us consider the first variant when during a sequential calculation of our function we are able to find conditional distributions of its subfunctions. Let is known distribution function of . Then can be calculated in usual way:

 Fv,X(y)=∫∫…∫ϕv(X,Z)≤y∏jdGj(zj). (18)

As the result, we have the following procedure. Let us consider the first (not initial) level. We choose values from samples , . It allows us to calculate a realization of conditional distribution function by formula (18). We repeat this procedure times and get a sample population , which elements are the distribution functions , .

Now let a vertex be an intermediate one. We have to extract the distribution functions instead of simple values. Let , be the set of the distribution functions, extracted on the -th step from the child samples. Then we calculate the distribution function

 Fqv,Y(y)=∫∫…∫ϕv(Y,Z)≤y∏iF∗qi(yi)∏jdGj(zj) (19)

and use it as -th element of the sample .

Finally for the root of the calculation tree we calculate the estimator by the following formula:

 Θ∗=1rr∑q=1∞∫0ydFqk,Y(y), (20)

where .

Now we consider the second variant when the conditional distributions of the available subfunctions are unknown. For each vertex and realization number , we have a sequence , of independent realizations of subfunction with the same distribution function . Therefore, the -th element of the sample population is the vector that ”represents” unknown distribution .

On the -th step we extract a vector from each sample , , forming resample , . Then for each random variable we generate (by the random number generator) a vector of its independent realizations, forming vector of realizations: in accordance to the distribution function . Further we calculate values and form a vector , which becomes -th element of the sample .

When the root of the calculation tree will be reached, we are able to estimate by analogy with (17):

 Θ∗=1rNr∑q=1N∑ξ=1Xqk,ξ, (21)

where .

Note that previously has denoted as the sample of function values. Now we have more general case: denotes either the set of the conditional distributions or the set of vectors which represent these distributions.

Let us consider the same system as in the Section 1.2, but with partially known distributions. Let us denote variables from that problem as . Let the distributions of random variables , and are unknown, but the distributions of , and are known. In our present definitions, , , ; , , . Our characteristic of interest is the expectation of function

 ϕt(x1,x2,x3,z1,z2,z3)={1ifmin{max{x1,z1},min{x2,z2},x3+z3}

Here we have the case when the conditional expectation (16) is known. It can be calculated easily:

 gt(x1,x2,x3)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩1,ifx2t,x1t,x1>t,x3t,x1t,¯¯¯¯¯¯G2(t),ifx2>t,x1>t,x3>t.

Now we are able to use the resampling approach to estimate . In this case, we have 8 -pairs, which are subsets of M={1,2,3}. We can calculate the values and .

## 2 Resampling estimation of stochastic process parameters

### 2.1 The failure model with accumulation of damages

Above we consider statistical models where time factor is absent. On the contrary stochastic processes have a dynamic character. Here the efficiency investigation of resampling estimator encounters great difficulties [Andr_00_2]. Let’s consider some statistical models, which are implemented to different reliability problems.

The model with accumulation of damages was considered in [Afan_02] and [Andr_06]. It is based on the failure model and its modifications, presented in [Andr_Gerc_72, Andr_94, Ger_00]. The model supposes two types of failures - initial and terminal failures. The initial failures (damages) appear according to a homogeneous Poisson process with rate . Each initial failure degenerates into a terminal failure after a random time . So if -th initial failure appears at time then the corresponding terminal failure appears at the time instant . Terminal failure and the corresponding initial failure are eliminated instantly. We assume that are i.i.d. r.v’s, independent on with cdf . We take interest in the number of initial failures which did not degenerate into terminal failures at time (further initial failures) and the number of terminal failures observed up to time . Let and be the corresponding expectations, , be the corresponding probability distributions, .

It is well known that and are mutually independent r.v’s, by that:

 EXt=λ∫t0(1−F(x))dx, EYt=λ∫t0F(x)dx, (22)
 PXt(i)=1i!(EXt)iexp(−EXt),i=0,1,…. (23)

The probability is calculated analogously by formula (23) where is replaced by . In future all formulas will be obtained for , but for all of them can be calculated analogously.

The rate and the cdf are unknown. Two samples are given: the sample =( of the intervals between initial failures and the sample =( of degeneration times. We need to estimate ,, and using samples and .

In order to estimate values (22) and (23) we use the resampling approach. On the -th realization we extract (without replacement) elements from the sample obtaining the -th resample , where . Then we calculate the instants of initial failures , ,  . Then we extract values from the sample (suppose that ), obtaining the -th resample , where . This allows us to calculate the instants of terminal failures and keep in mind the number of failures of each type up to time . Further all extracted values are returned into the initial samples and the described procedure is reiterated times.

Let be the indicator function of the event: ”The -th initial failure occurred, but did not degenerate into a terminal failure up to the time ”:

 ζ∗qj(t)={1 if τ∗qj≤t<τ∗qj+B∗qj,0 otherwise. (24)

Then the number of initial failures at time for the -th realization is

The resampling-estimator of can be obtained from formula (3), where .

Now we need to calculate the resampling-estimators of the probabilities . Let be the indicator function of the event . The resampling-estimators of the probabilities can be determined by formula (3) taking .

Let us calculate the expectations of the resampling-estimators (note that they are biased). Obviously , . These expectations can be calculated taking into account the following reasoning: 1) The probability of the event that the number of initial failures occurred during time is equal to can be found by a Poisson distribution: ; 2) It is known, that if the number is fixed then the moments of initial failures are independent and uniformly distributed on (0, ); 3) If is the time moment of an initial failure appearance then in the time instant with probability it remains initial; 4) The probability that at time the considered failure remains initial is: .

Therefore, the expectation of the resampling-estimator is calculated as follows:

 E(E∗Xt)=p1nA∑j=1jdt(j)+p1nA∞∑j=nA+1dt(j). (25)

We also can find the expectation of the estimator :

 EP∗Xt(i)=nA∑j=idt(j)(ji)pi1(1−p1)j−i+(nAi)pi1(1−p1)l−i∞∑nA+1dt(j), i=1,2,…. (26)

Let us illustrate the idea of the calculation of the variance of . It can be obtained from formula (5), where .

For that purpose we need to calculate the covariance for two different realizations and . Let us consider the second mixed moment . We have:

 EX∗qtX∗q′t=EnA∑j=1ζ∗qj(t)nA∑i=1ζ∗q′i(t)=nA∑j=1nA∑i=1Eζ∗qj(t)ζ∗q′i(t). (27)

Now we need to calculate . We have to take into account that and can be formed by the same intervals between initial failures and by the same degeneration times. Let be the number of the same intervals between initial failures on realizations and ’. Let be the random event ”the -th initial failure in the -th realization and the -th initial failure in the ’-th realization have the same degeneration time.” Then

 Eζ∗qj(t)ζ∗q′i(t)=j∑ν=0P{αA=ν}(P{RjB=i}E