This article has been accepted for inclusion in a future issue of this journal.

Content is final as presented, with the exception of pagination.

IEEE TRANSACTIONS ON RELIABILITY 1

A Machine Learning Algorithm

for Reliability Analysis
Maria Luz Gámiz , Fernando Jesus Navas-Gómez, and Rocío Raya-Miranda

Abstract—In this article, we build a statistical model able to pre- be determined when either the minimal paths or the minimal
dict the reliability of the system based on a dataset. Our objective cuts are known.
is double. On the one hand, we aim at constructing a function
Several procedures have been developed to generalize the con-
that classifies the system in one of the two categories (operative
or failed) based on the knowledge of components states. On the cept of binary coherent structure to a multistate configuration,
other hand, we present a statistical test to decide the order of and then, the structure function has been specified by means
importance of components in terms of the effect each one has of a finite set of boundary points (see, for example, [29] and
on the system performance. We present a supervised algorithm references therein).
involving isotonic smooth logistic regression and cross-validation Aven [3] justifies the introduction of multistate models by
techniques. Our method is completely data-driven not lying in
any parametric assumptions. The method is illustrated through necessity in some areas of application in engineering, such
an extensive simulation study. as gas/oil production and transport systems, where a binary
approach would give a poor representation of the real world.
Index Terms—Cross-validation, importance reliability
measures, isotonic smoothing, logistic regression (LR),
This author investigates the problem of calculating certain mea-
quasi-likelihood. sures of the performance of a monotone system with multiple
states and presents comparative studies on the accuracy of his
calculations by a Monte Carlo simulation study.
I. INTRODUCTION Baxter [7] introduced continuous models for system reliabil-
ity and since then a wide variety of performance measures have
MATHEMATICAL representation of the logic of a system
A is one of the main objectives pursued in reliability analysis.
Since system performance assessment can be a complex problem
been defined and calculated to be valid for binary, multistate and
continuous systems (see [10] for a detailed summary).
For continuous systems, if the structure function cannot be de-
in practice, even for simple structures, it seems reasonable that
termined based on qualitative characteristics (e.g., series and/or
a skillful procedure for modeling the relationship between the
parallel structures) or by boundary point analysis, approximation
state of the system and its components can help efficiently in
methods are required. To this end, Lisnianski and Levitin [23]
assessing the reliability of complex systems. Thus, the system
explore an approximation based on the technique of the uni-
structure function, which represents the link function between
versal generating function. The method consists of a discrete
the state of the system and all its components, has been a
approximation of the system performance with a continuous set
major topic in the field of system reliability. It is defined as
of states using a finite multistate system and the purpose is to
a mapping that relates each possible combination of states of all
construct upper and lower bounds for the reliability measures of
components to a system state.
the continuous system.
Most reliability models have traditionally focused on a binary
In view of the great difficulty inherent in analytical perfor-
formulation of system behavior, that is, the models allow only
mance evaluation (even more for continuous systems), new
two levels of performance for a system and its components:
approaches to empirical methods have emerged in the last two
perfect operation (1) and complete failure (0). The total number
decades. Brunelle and Kapur [10] proposed a multivariate in-
of combinations of components states (vector states) in this case
terpolation method by which, starting from a dataset, they can
is 2p , being p the size of the system. The structure function can
construct the structure function of the continuous system to
represent the value of the state of the system given the states
of its components.
Manuscript received December 16, 2019; revised May 23, 2020; accepted
July 18, 2020. This work was supported by the Spanish Ministry of Economy
Another approach is presented in [19] where a new method-
and Competitiveness under Grant RTI2018-099723-B-I00 and Grant EDU2016- ology that includes regression methods to approximate the
79322-R. Associate Editor: K. Goseva-Popstojanova. (Corresponding author: structure function of a system based on empirical information
Maria Luz Gámiz.)
The authors are with the Department of Statistics and Operational Re-
about the states of a sample of systems and its corresponding
search, University of Granada, 18071 Granada, Spain (e-mail: mgamiz@[Link]; components states is proposed. The main idea is the use of
fernandonavas2222@[Link]; rraya@[Link]). nonparametric monotone regression techniques, given that the
Color versions of one or more of the figures in this article are available online
at [Link]
nature of the problem requires the regression function must be
Digital Object Identifier 10.1109/TR.2020.3011653 monotone in each explanatory variable. The adjustment method

0018-9529 © 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
See [Link] for more information.

Authorized licensed use limited to: University of New South Wales. Downloaded on August 09,2020 at 08:58:37 UTC from IEEE Xplore. Restrictions apply.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

2 IEEE TRANSACTIONS ON RELIABILITY

is developed taking into account two steps, the first one in- and U is a random variable for which we assume that E[U ] = 0.
volves the use of local regression techniques for constructing Note that the function ψ is not the structure function of the
the surface that best fits the data, and then it is considered a nu- system as it is understood in the usual literature.
merical method of isotonization of the response. The algorithm Unlike [19], we are not interested in quantifying the level
of pool adjacent violators, which was introduced to solve the of performance of the system. Rather we aim at predicting
one-dimensional case, and subsequently has been generalized the probability that the system provides an acceptable level
for larger dimension problems (see [12]). of performance as well as evaluating the effect of each single
In this context of using expert data to construct the structure component on the system behavior.
function of a multistate system, the recently published paper by We use a machine learning procedure to build the model able
Zaitseva and Levashenko [34] proposes the use of fuzzy decision to predict the reliability of the system on the basis of a dataset.
trees. These authors have into account that data collected from The final model is determined by a tuning parameter that is
the real world is uncertain. This uncertainty may be due to selected by a validation technique. We use leave-one-out cross-
inaccuracy or error in measurements or, in a more general validation where the learning algorithm is applied once for each
sense, to some random factors that are inherent to the physical sample item, using the rest of the sample as a training set and
behavior of the system however not controlled by the analyst. using the selected item as a single-item test set. Recently, as can
In their approach, Zaitseva and Levashenko [34] consider that be seen in [1] and the references therein have demonstrated the
uncertainty in the data is not random in nature and then it cannot potential of a machine learning approach for industrial reliability
be indicated in a quantitative form by probability theory and analysis.
they propose fuzzy logic to define the structure function. Within a probabilistic framework, LR is the most widely
Motivated by the fact that there may be uncertainty affecting used method for classification problems. However, many other
the system behavior, meaning by this, for example, lack of algorithms are very popular in machine learning such as support
knowledge or limited time that does not allow a detailed analysis, vector machine (SVM), random forest (RF), neural network
Coolen and Coolen-Maturi [13] propose the generalization of (NN), among others. One of the arguments against LR is that it
the traditional deterministic conception of the structure function cannot handle nonlinearities in the data. In such a case, SVM, for
to a more flexible formulation of the concept as a predictive example, is recommended. However, in this article, we consider
(imprecise) probability. LR from a nonparametric viewpoint, so we do not impose a linear
We believe that taking into account uncertainty is crucial model for the underlying logit function. Although classification
when describing system structure models. In this article, we performance is almost identical in both cases LR and SVM,
reconsider the viewpoint of [19] and assume that the system local LR can provide class probabilities whereas SVM is a
performance level is measured with the uncertainty that is in- deterministic classifier. This is an argument in favor of LR since
deed of random nature, and we construct the structure function our objective in this article is not limited to assign elements to
using as input variables the states of the components that are one in a set of categories, rather we want to quantify the strength
modeled as random variables ranging in an interval. An extra of the relationship between the outcome and every single feature
term of random error is added to the model. As in [19], we use in the model. In [28], Musa presents a comparative study on
regression techniques to evaluate the probability that the system performance between SVM and LR. The conclusion of the study
is functioning. is that LR has higher interpretability, whereas SVM is considered
Our objective is double. On the one hand, considering the to be a black-box predictor. It neither makes its prediction
context of structural reliability, we want to construct a limit implicit nor gives incite in the rules governing its prediction,
state function that allows us to divide the space of states of which is not the case in LR. Therefore, in the case of considering
components into the failure and safe regions. To this end, we classification only, each of them can be used while when the
propose the use of classification techniques based on logistic interpretation is necessary such LR should be used. On the other
regression (LR). hand, LR can be seen as a special and simple case of NN, but
On the other hand, we want to rank the components of the in general, NNs work as a black-box that makes it difficult
system according to their effect on the system performance to interpret, so it is not easy to identify important predictors
without any knowledge of the system design. That is, the only and relative rank, moreover they require more computational
information required is the current state of each component as time [6]. Decision trees and their generalization as RFs are also
well as the state of the system itself, but no information about intuitive and easy to implement but do not provide a probability
the system design is necessary. for the outcome as a direct result. Generally speaking (see [33]),
Without any loss of generality, we represent the state of for classification purposes none of the above-mentioned meth-
the ith component of the system by a random variable Xi ods is clearly superior, it rather depends on the problem at
(i = 1, 2, . . . , p), which is directly observable and that takes hand. However, when further inferences about the model are
values on the interval [0, 1]. Given a particular configuration needed, as is the case of this article, LR is the recommended
of the components, the system state Y ∗ is a random variable that method.
fits the following model (see [19]): The rest of this article is organized as follows. Section II
proposes a logistic model to estimate the reliability of the
Y ∗ = ψ (X1 , X2 , . . . , Xp ) + U
system. In Section III, the problem is solved using nonparametric
where ψ is a function usually assumed to be linear, that statistics techniques. In Section IV model is evaluated by means
is, ψ(X1 , X2 , . . . , Xp ) = β0 + β1 X1 + β2 X2 + · · · + βm Xp ; of sensitivity analysis. In Section V, a version of the Birnbaum

Authorized licensed use limited to: University of New South Wales. Downloaded on August 09,2020 at 08:58:37 UTC from IEEE Xplore. Restrictions apply.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

GÁMIZ et al.: MACHINE LEARNING ALGORITHM FOR RELIABILITY ANALYSIS 3

importance measure is given. Section VI is devoted to numerical the states of all components in the system, p is the size of the
applications. In Section VII, a brief discussion on some aspects system and Xi ∈ [0, 1] for all i = 1, 2, . . . , p.
about system complexity is presented. Finally, Section VIII The reliability of the system is the probability that the system
concludes this article. is operative. We consider the reliability of the system as a
function of the state vector, i.e.,
II. MONOTONE LR FOR RELIABILITY FITTING R (x1 , . . . , xp ) = E [Y |x1 , . . ., xp ]
Our first interest is in estimating the probability that a structure
eβ0 +β1 x1 +···+βp xp
performs as required under specified service conditions. To = . (1)
estimate this probability we propose the following procedure. 1 + eβ0 +β1 x1 +···+βp xp
The goal is estimating the vector of regression parameters β =
A. LR Model Based on a Latent Variable (β0 , β1 , . . ., βp ) from sampling information about the system
and its components.
Let X = (X1 , X2 , . . . , Xp )t be the input vector of variables For a given dataset {(xi ; yi ) ∈ [0, 1]p × [0, 1]; i = 1, . . ., n},
that gives the states of all components. The level of performance we estimate β by maximum-likelihood and obtain an estimate of
of the system is determined by the function the system reliability according to (1). Alternatively, the model
Y ∗ = ψ (X) + U can be expressed as
R(x)
where U is a symmetric random variable with E[U ] = 0 and g(x) = log = β0 + β1 x 1 + · · · + β p x p .
1 − R(x)
Var(U ) = σ 2 . Let y0 be a threshold value such that when Y ∗ >
y0 the system is performing its function satisfactorily; otherwise, Maximum-likelihood methods lead to the following optimiza-
i.e., when Y ∗ ≤ y0 , the system is failing. In this situation, we tion problem:
define Y equal to 1 when Y ∗ > y0 , and equal to 0 for Y ∗ ≤  = arg max {log L(β)}
β
y0 . When only information related to variable Y is registered β
whereas Y ∗ is not directly observable, Y ∗ can be interpreted as  n


a latent variable. = arg max Yi log Ri + (1 − Yi ) log(1 − Ri )
The probability that a particular configuration X = x leads β i=1
to a good behavior in the system, i.e., Y = 1, can be expressed
as +with Ri = R(xi ). To be consistent with the coherence require-
ments of the system, we need the reliability to be a nonde-
P (Y = 1|X = x) = P (Y ∗ > y0 |X = x) creasing function in each of its arguments. So ∂R(x) ∂xj > 0, for
= P (ψ(X) + U > y0 |X = x) j = 1, . . ., p. Note that the sign of the jth partial derivative is
determined by the sign of βj , for all j = 1, 2, . . . , p. Besides,
= P (U > y0 − ψ(x)) the proper extrema condition must be held. This means that
= 1 − F (y0 − ψ(x)) R(0, 0, . . . , 0) = 0 and R(1, 1, . . . , 1) = 1. It is then needed
to introduce some order restrictions in the above optimization
where F is the distribution function corresponding to U . problem. Therefore, the problem is stated as follows.
As a first approach, we take ψ a linear function of the input
Solve
variables, i.e.,
   = arg max {log L(β)}
β
y0 − ψ(x) β
P (Y = 1|X = x) = P U0 >
σ s.t.
= 1 − F0 (β0 + β1 x1 + · · · + βm xp ) 1) βj > 0 (j = 1, . . . , p)
2) β0 < 0
= 1 − F0 (βx ) 3) β0 + β1 + · · · + βp > 0.
where we denote β = (β0 , β1 , . . . , βm ), and F0 is the distribu- The two last restrictions are needed to meet the proper extreme
tion function of the standardized random variable U , which is conditions.
assumed symmetric around 0. Perhaps the most popular option
eu
is F0 (u) = Flogis (u) = 1+e u , with u ∈ R, the standard logistic III. LOCAL LOGISTIC MODEL FOR RELIABILITY
distribution function. Note that an identification problem is
Usually, parametric regression models entail functional as-
involved since the parameters y0 , and σ cannot be estimated
sumptions that can be restrictive in empirical applications. In
separately from vector β.
Section II, we assume that the logit transformation of R is a
linear function depending on x. However, in many practical
B. LR Model Under Coherence Restrictions occasions, this functional relationship may not be appropriate
According to the previous section, let Y be the observable leading to wrong decisions, and then it is more convenient not
random variable specifying by the value of 1 that the system to speculate with any particular parametric model to fit the data.
is operative and by 0 that the system is out of work; X = Fully nonparametric regression models allow for more flexi-
(X1 , . . ., Xp )t is a p-dimensional random vector that specifies bility although they are rarely used for the estimation in binary

Authorized licensed use limited to: University of New South Wales. Downloaded on August 09,2020 at 08:58:37 UTC from IEEE Xplore. Restrictions apply.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

4 IEEE TRANSACTIONS ON RELIABILITY

response applications. One major problem of nonparametric where xt = (x1 , x2 , . . . , xp ) and with g a smooth function in the
regression in higher dimensions is the so-called curse of di- sense of derivability. We can use the first-order Taylor expansion
mensionality referring to the existence of sparse data in the of g around a point x0 , and write
neighborhood of an estimation point. One of the consequences
g(x) = g (x0 ) + (x − x0 )t g  (x0 ) + o ( x − x0 )
of this problem is the high variance of estimators. In the case
of binary response regression, the curse of dimensionality does where o(·) is a function that satisfies that o(h)/h −→ 0, when
not disappear, however, the finite-sample variance problems are h → 0. The approximation above is valid for x in a neighbor-
ameliorated due to the boundedness of the response (see [18]). hood of x0 conveniently chosen, that is for all x such that
In our context of reliability systems, there is an extra issue to x − x0 < h, with · the Euclidean norm and h > 0 suf-
∂g
have into account, that is, the fitted model must meet the coher- ficiently small. We denote g  = ( ∂x , ∂g , . . . , ∂x
1 ∂x2
∂g t
p
) the gra-
ence order restrictions. Nonparametric regression produces esti- ∂g(x0 )
dient vector of g. With notation β0 = g(x0 ), y βj = ∂xj ,
mators with good asymptotic properties such as consistency and
normal distribution that are reached at a reasonable convergence j = 1, . . . , p, the logit can be written as
rate. However, these estimators are not necessarily monotone. g (x) ≈ β0 + β1 (x1 − x01 ) + · · · + βp (xp − x0p ) . (3)
To get monotonicity of the final estimator, we propose a two
stages sequential procedure: first isotonize the data on {X, Y } With (2) and (3), we can see the nonparametric model as a locally
and then smooth the resulting isotonic response using a local log-linear model, and the parameters can be estimated locally
quasi-likelihood logit method. by maximum-likelihood techniques.
At this second stage, we do not work with the original data,
instead, we have the following dataset {(xi , Yi ); i = 1, . . . , n},
A. First Step: Isotonization
where {Y1 , . . . , Yn } are the responses isotonized with the pro-
There are good reasons to first isotonize and then smooth cedure of the previous section, and that are obtained by con-
instead of proceeding in the opposite order. They are explained veniently averaging some of the adjacent original data Yi .
in the work of [27]. On the one hand, as explained by the author, Then, Y |X=x is not a binary random variable and the like-
frequently isotonic regression produces too many “flat spots,” lihood method cannot be applied as it is stated in previous
which is an argument against having a good estimate of the sections.
slope of the response, an objective pursued in this article to However, we can assume that E[Y |X=x ] ≈ E[Y |X=x ] and
estimate the importance measures defined in Section V. On the Var[Y |X=x ] = c Var[Y |X=x ] = c V (R(x)), with V a known
other hand, the usual kernel functions used in the smoothing step function. In our case V (u) = u · (1 − u). So, although we can-
lead to nondecreasing smooth regression functions, moreover, not specify a full likelihood function based on the isotonized
the convergence rate of isotonic estimators can be improved data, we can specify the relationship between the mean and the
according to the study in [27]. variance and we can solve the estimation problem considering
As explained in [24], the isotonization problem can be seen a quasi-likelihood approach (see [16]). The quasi-likelihood
as a quadratic programming problem. The interest is in fitting a function satisfies
scatterplot from the model yi = f (xi ) + i , i = 1, . . . , n, where ∂ y−u
no parametric form is assumed for f although it is known Q(u, y) = .
∂u c · V (u)
that it is nondecreasing. Defining zi = f (xi ), the problem is to
find z ∈ Rn minimizing ||y − z|| subject to Az ≥ 0, for matrix The above score function satisfies similar properties as the
A of dimension (n − 1) × n and with elements A(i, i) = −1, likelihood score function that are as follows.
∂
A(i, i + 1) = 1, and A(i, j) = 0, for j = i, i + 1. The set C = 1) E[ ∂u Q(u, y)] = 0.
∂
{z ∈ Rn : Az ≥ 0} is a polyhedral convex cone in Rn . This 2) Var[ ∂u Q(u, y)] = c V1(u) .
2
problem was first solved by using the pool adjacent violators ∂ 1
3) −E[ ∂u 2 Q(u, y)] = c V (u) .
algorithm (PAVA) algorithm proposed by Brunk [11]. In this ar- In this context, we want to find the value of β that maximizes
ticle, we will use the hinge algorithm presented in [24] and [26], the local quasi-likelihood function
which is an algorithm for cone projection simple and intuitive
n

and faster than PAVA [26]. The result is a new dataset {(xi , yi )},
0 (β) = wh,0 (xi ) Q R0 (xi ) , Yi
that satisfies the required order restriction: xi ≺ xj ⇒ yi ≤ yj ,
i=1
for all i, j ∈ {1, 2, . . . n}, with ≺ denoting the usual partial order
β t (x −x )
in the p-dimensional space Rp . We take this dataset as the input where R0 (xi ) = e βt (x i 0
; β = (β0 , β1 , . . . , βm )t ; and
1+e i −x0 )
for the local smoothing algorithm. {wh,0 (xi ), i = 1, . . . , n} determine the weights of the obser-
vations around the estimation point x0 . This is to ensure that the
B. Second Step: Local Quasi-Likelihood Logit linear approximation of g given in (3) is considered only into
The random variable Y |X=x follows a binomial distribution the window where it is valid. We define
B(1, R(x)). From a nonparametric point of view, we assume wh,0 (xi ) = Kh ( xi − x0 )
the following:
with Kh (·) = (1/h)K(·/h), being h the bandwidth parameter
R(x) that controls the amount of smoothing, i.e., it determines the
log = g(x) (2)
1 − R(x) size of the window around the point x0 where the local linear

Authorized licensed use limited to: University of New South Wales. Downloaded on August 09,2020 at 08:58:37 UTC from IEEE Xplore. Restrictions apply.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

GÁMIZ et al.: MACHINE LEARNING ALGORITHM FOR RELIABILITY ANALYSIS 5

approximation is valid. Usually, K is a density function with indicating that there exists a threshold along which both classes
a compact support. The quasi-likelihood solution for β comes can be separated without error. An uninformed classifier that
from solving the following score equations, for j = 0, 1, . . . , p: performs no better than chance would report AUC  = 0.5.
 n
∂0
(β) = wh,0 (xi ) V. IMPORTANCE OF COMPONENTS
∂βj i=1
The main purpose of the measures of importance in reliability
Yi − R0 (xi ) ∂R0 analysis is to rank the components of the system in order to detect
(xi ) = 0.
c V (R0 (xi )) ∂βj those components that due to location or other particular char-
acteristics have a higher probability of causing system failure.
∂0
The gradient vector with elements given by ∂βj , for j = A number of measures have been introduced considering the
1, . . . , p, is written as importance a component has for the system behavior according
n
 to different points of view (see, e.g., [5]).
0 (β) = {c}−1 wh,0 (xi ) Yi − R0 (xi )
i=1 A. Birnbaum Importance Measures
· (xi − x0 ) . (4) Birnbaum in 1969 was first to introduce the concept of impor-
Solving the problem of nonlinear equations, we can obtain a tance measures. In [2], [8], and [21], the importance measures
local estimate of the function g, that is, the estimated value of are categorized into three types based on the knowledge for
β0 provides an estimate of g(x0 ) and we also obtain the estima- determining them. These three classes are as: structure impor-
tions of the corresponding partial derivatives. The asymptotic tance measures, reliability importance measures, and lifetime
properties of the local estimator are studied in the Appendix. importance measures.
Then, Birnbaum [8] defines the importance of the jth com-
IV. SENSITIVITY ANALYSIS: ROC CURVE ponent in terms of the probability that the system is in a state
in which the functioning of the jth component is critical, i.e.,
Until now we have constructed a model that estimates the the system functions if the jth component functions and fails
probability that the system fails. We want to predict the state of otherwise. This measure is probably the reference index used by
the system then the next step is to define a classification rule to practitioners to evaluate the relative importance of components
map the LR probabilities into the set {0, 1}. To this purpose, we in complex systems and in this sense it is continuously utilized
use the area under the receiver operating characteristic (ROC) in practical applications where it is revised and refined (see, for
curve (AUC) to tune the classification threshold parameter. The example, [32]). Let r = (r1 , r2 , . . . , rp ) be the vector whose ith
ROC curve is a graphical evaluation of the performance of the element is the reliability of the ith component of the system.
model as a classifier. In the case of independent components, this measure is the
In our context, the value of the AUC can be seen as an partial derivative of the system reliability function R, and in
estimate of the probability that the predicted reliability of a failed the literature, it is common to refer to the Birnbaum measure
system is not higher than the predicted reliability for a system by this derivative (see, e.g., [4] and [31]). Hence, the Birnbaum
that is operative. To approximate this probability we proceed as importance measure I B is given by
follows. First, we split the sample into two sets. On the one side,
we take the failed systems (Yi = 0) and, on the other side, we IjB (R) =
∂R
(r1 , r2 , . . . , rp ) .
consider the operative systems (Yj = 1). So let define the two ∂rj
sets Y0 = {i : Yi = 0} and Y1 = {j : Yj = 1}, with respective
sizes n0 and n1 , (n0 + n1 = n). The LR model also classifies Note that this quantity measures the change in the system
the systems of the sample in one of the two classes 0 and 1 reliability in terms of component reliabilities. However, there
according to the value of the predicted system reliability R(x). are situations when only the design of a system is known, but
This classifier is doing reasonably well if the predicted reliability no information is available about the component reliabilities. In
is higher for any system indexed in Y1 than for systems indexed this situation, Birnbaum [8] defines what is called “structural
in Y0 . Then, we define the following statistic importance.” For a binary system, it is defined as
 
 = 1  i ) < R(x j) . IjB (φ) = {φ (1j ; x) − φ (0j ; x)} .
AUC I R(x (5)
n0 · n1 x
i∈Y0 j∈Y1
This index quantifies the effect that the change of a component
 is an estimator of the probability that
It is obvious that AUC state has on the state of the system. As argued in [8], if nothing
a randomly chosen element of class 0, which is indexed in the is known about the reliabilities of the components, and for lack
set Y0 is ranked by the regression model lower than a randomly of better knowledge, it is assumed that all vectors x ∈ {0, 1}p
chosen element of class 1, indexed in the set Y1 . are equally probable and then the reliability importance reduces
 can thus be used as an assessment of the dis-
The value of AUC to the structural importance in the sense that
criminatory power of this classifier (how well it can separate two
 = 1,
classes). Obviously a perfect classifier would lead to AUC IjB (R) = E [φ (1j ; x) − φ (0j ; x)] = IjB (φ).

Authorized licensed use limited to: University of New South Wales. Downloaded on August 09,2020 at 08:58:37 UTC from IEEE Xplore. Restrictions apply.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

6 IEEE TRANSACTIONS ON RELIABILITY

B. New Version of the Birnbaum Importance Measure 

according to a factor eβj Δ . Therefore, we choose as the most
In this article, we see the system reliability as a function relevant component for the function of the system, the one
of the states of the components instead of the reliabilities. for which a fixed amount of improvement in its state-level
So, if we denote R(x|xj = a) the probability that the system produces the highest increment in the odds index. In other words,
works for a configuration of components where the state of component j is more relevant for the system performance than
the jth component is a, then we can interpret the following component k if βj > βk , for j, k = 1, 2, . . . , p.
index R(x|xj = 1) − R(x|xj = 0) in the following sense. With
a binary approach, a ∈ {0, 1}, this index measures the effect D. Inference About IjB
that the event “failure of the jth component” has on the system
The conclusion from Section V-C is that the importance
reliability. In a continuous system, the state of the jth component
of a component with respect to the good performance of the
is xj = a ∈ [0, 1], for j = 1, . . . , p and R(x) = P r(Y = 1|x)
system can be established according to the magnitude of the
the probability that the system is operative when the vector state
corresponding partial derivative of the logit function. So the
takes value x = (x1 , . . ., xp ) ∈ [0, 1]p , and we can generalize to
ranking of components is made just by the decreasing order of
the following R(x|xj = a + Δ) − R(x|xj = a), which quan-
the βj coefficients, for j = 1, 2, . . . , p. In the parametric case of
tifies the change in system performance caused by an increment
(1), the β coefficients are constant, and then the ranking is valid
of magnitude Δ in the state of component j. So, we focus, on
for all range of x. Based on the logistic model, it is interesting to
the one side, on the reliability of the system, and the states of
carry out hypothesis testing to statistically decide whether two
components, on the other, to define the following.
components have the same effect for the system performance or
Definition. A version of the Birnbaum measure IjB (R, x)
not. In other words, let j = k, we can use the data and the LR
Let S be a coherent system of size p and let us denote
model to decide between H0 : βj = βk , against H1 : βj = βk .
R(x1 , . . ., xp ) = P (Y = 1|x1 , . . ., xp ). The local importance
The usual Wald test, score test, or likelihood-ratio test can be
Birnbaum measure of the jth component is defined as
used to solve the problem.
∂R (x) In the case of a local-logit model, we fit the curve locally
IjB (R, x) = . (6)
∂xj and then, the local-slopes, i.e., the βj coefficients for j > 1, are
defined to depend on the estimation point. Thus, the importance
This measure allows us to rank the components of the system
of a particular component is estimated locally. We want to test
according to the magnitude of the effect that a gradual change in
the null hypothesis
the state of the components has on the reliability of the system.
H0 : βj (x) = βk (x)
C. Local Estimation of IjB (R, x)
for a particular system configuration x = (x1 , x2 , . . . , xp ).
We can use the LR model just defined to estimate the reliability
Local-likelihood estimates are asymptotically normally dis-
of the system and then an estimate of the above importance
tributed (see, for example, [15]). To solve the problem, we
measure (6) can be derived accordingly. First, we obtain  with Ejk vector
consider the following statistic Tjk = Ejk β,

∂R eβx of 0 s, except for positions j and k, which are, respectively,
(x) =  2 j
β 1 and −1. We can use the results of Section V III to derive
∂xj (1 + eβx )
the distribution of the statistics Tjk . Since we take h → 0, we
with x = (x1 , . . ., xp ). Since Y is a random variable with bino- can assume that the bias is small enough to be ignored, and the
mial distribution B(1, R(x)), we can equally write  t .
variance of the statistics is Var(Tjk ) = Ejk Var(β)E jk

∂R
IjN B (R, x) = (x) = V (Y )βj . VI. NUMERICAL APPLICATION
∂xj
In this section, we summarize the steps of the algorithm
Using the LR model, for a given x, the system components can be
and present some numerical results. Let us consider a dataset
ranked according to the value of the corresponding β-coefficient,
{(xi , yi ); i = 1, 2, . . . , n} as previously denoted.
or, equivalently according to the value of the odds − ratio,

that is, ORj = eβj , which is a central index provided by the
statistical software most used by the practitioner to implement A. Algorithm: Local-Logit Fit and Components Comparison
LR fitting. Note that in our context, OR is always greater than Step 1: Isotonization (Section III-A)
1. Isotonize the data matrix according to the partial
The ratio P r{Y = 1|x}/P r{Y = 0|x} measures the odds order of components states. This routine implements
(chance) for the system to be operative against to be failed. the hinge algorithm for cone projection to minimize
The bigger the odds index is, the better the components are y − θ 2 over a cone C of the form {θ : Aθ ≥ 0} as
configured for the system behavior. If we concentrate on com- explained in Section III-A. In particular we use part
ponent j, whose state is at level xj and, controlling for the rest of the routines coneproj and sparsebasis that can
of components, let us now consider that the state is improved be found in [25].
to value xj + Δ, then the odds index is expected to increase Step 2: Smoothing (Section III-B)

Authorized licensed use limited to: University of New South Wales. Downloaded on August 09,2020 at 08:58:37 UTC from IEEE Xplore. Restrictions apply.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

GÁMIZ et al.: MACHINE LEARNING ALGORITHM FOR RELIABILITY ANALYSIS 7

We construct a local logistic model to fit the iso-

tonized dataset obtained in Step 1. We built a routine
to carry out a quasi-likelihood estimation procedure.
To solve the nonlinear estimating equations we use
the functions implemented in the CRAN package
nleqslv. As a result, for a given level of components Fig. 1. Series–parallel system.
performance (location), we estimate the reliability of
the system as well as the partial derivatives at the
given location. detailed in Appendix. For a given location, i.e., a par-
Step 3: Model selection ticular level of components state values, we develop
The final model is determined by a tuning param- the estimator of the model and evaluate its accuracy in
eter (the bandwidth parameter) that is selected by terms of bias and variance as explained in Appendix.
a validation technique. We use the leave-one-out The ranking of components is established locally as
cross-validation. We let the model learn from data follows.
as follows. We select one element of the sample and Let us consider a system of p components.
apply the estimation algorithm considering the rest of 1) Define a grid of points in the interval [0, 1],
the sample as the training set, the selected element is that is XM = {x0 , x1 , . . . , xM }. Construct an
treated as a single unit test set. More specifically, for M × p matrix XM , whose jth column is the
[i]

β
[i]
vector (x0 , x1 , . . . , xM )t , and denote xi· its ith
each i ∈ {1, 2, . . . , n}, define Ri,h (xi ) = e [i] ,
h,0


β row. Note that the ith row of matrix XM has
1+e h,0
where βh,0 = g[i] (xi ) is obtained as explained in all its elements equal, indicating that we want to
[i]

Section III-B for a bandwidth h and using all the evaluate at a configuration of the system where all
data except xi . Define the cross-validation score as components are at the same level of performance.
2) For each i = 1, 2, . . . , M , fit the local-logit
n
 2 model following the steps explained previously.
Yi − Ri,h (xi )
[i]
Q(h) = n−1 Let us denote β(x  i· ), the corresponding estimator
i=1
at a given state vector xi· = (xi , xi , . . . , xi ).
then, the cross-validated bandwidth is defined as 3) Choose two components in the system. In other
hCV = arg minQ(h). words, choose two different elements from set
h
Step 4: Model evaluation and limit state function (Sec- {1, 2, . . . , p}. Let us denote them j and k, with
tion IV) j = k. Define the statistic Tjk (xi· ) = βj (xi· ) −
To evaluate the model we develop sensitivity anal- βk (xi· ). The normalized statistic Zjk (xi· ) =
jk (xi· )
Tjk (xi· )−B
ysis based on the ROC curve. A measure of the follows a normal distribution
 jk (xi· )
Σ
goodness-of-fit is given by the area under the curve
with mean βj (xi· ) − βk (xi· ), where we consider
that is approximated by the expression (5).
(9) and (10) to obtain bias and variance estima-
To define the limit state function we first find the
tions.
p0 value of probability (p0 ∈ (0, 1)) that provides
4) For each i = 1, 2, . . . , M , we consider a hypoth-
the point of the ROC curve closest to the point of
esis testing problem
coordinates (0, 1), that is, we estimate the optimal
cut probability for correct classification. For a given H0 : βj (xi· ) = βk (xi· )
x ∈ [0, 1]p define the limit state function as f (x) =
H1 : βj (xi· ) > βk (xi· ).

I[R(x) > p0 ].
Step 5: Ranking of components (Section V) 5) Under the null hypothesis, Zjk (xi· ) follows a
We develop a series of hypothesis testing aimed to standard Normal distribution. Positive and large
compare locally the Birnbaum importance measures values of these statistics will lead to a rejecting
for each pair of components of the systems. In short, decision of H0 in favor of the alternative.
for a given state vector x0 ∈ [0, 1]p , we are interested
in testing the following hypothesis: B. Simulation Results

H0 : IjN B (x0 ) = IkN B (x0 ) We consider three practical examples.

1) Example 1. Series–Parallel System of Size 3: Our first ex-
for all 1 < j < k < p. This question can be equiv- ample considers a series–parallel combination system of p = 3
alently formulated in terms of the coefficients of components, as displayed in Fig. 1. Component 1 is in series with
the model, i.e., β(x0 ). A comparison test for the a subsystem made of components 2 and 3, which are placed in
parameters of the model can be developed according parallel. In these conditions, the structure function of the system
to Section V − D. We obtain the distribution of the is derived from ψ(X1 , X2 , X3 ) = min{X1 , max{X2 , X3 }},
test statistic from the asymptotic normal distribution where Xj denotes the state of component jth, j = 1, 2, 3. We

Authorized licensed use limited to: University of New South Wales. Downloaded on August 09,2020 at 08:58:37 UTC from IEEE Xplore. Restrictions apply.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

8 IEEE TRANSACTIONS ON RELIABILITY

TABLE I
SERIES–PARALLEL SYSTEM

TABLE II
SUMMARY STATISTICS FROM M = 1000 SIMULATED SAMPLES FROM THE
SERIES–PARALLEL SYSTEM

Fig. 2. Bridge system.

TABLE III
BRIDGE SYSTEM
simulate a sample of n = 100 such systems so the data consist of
a matrix with p + 1 columns. The first p columns report the states
of the components as the value of independent random variables
with distribution U(0, 1). The state of the system is simulated
considering a latent variable that is not directly observed and TABLE IV
is assumed to be Y ∗ ∼ N (ψ(x1 , x2 , x3 ), σ), for a particular SUMMARY STATISTICS FROM M = 1000 SIMULATED SAMPLES FROM THE
configuration of the state vector (x1 , x2 , x3 ). We fix σ = 0.2 BRIDGE SYSTEM
and finally, the information about the state of the system is
simulated from a Binomial distribution with event probability
given by π(x1 , x2 , x3 ) = P {Y ∗ > y0 }, with y0 = 0.5.
We have run a simulation study with a total of M = 1000
repetitions of the experiment. Table I lists the results of model
fitting considering both methods presented in previous sections: As expected, the results provided by the parametric model
parametric and local-logistic approach. The values presented in and displayed in Table II suggest that component 1 has a higher
the table are the average along the M samples of the corre- effect on the system behavior than components 2 and 3, which
sponding indexes at the table header. Specifically, we consider are similar.
the following. 2) Example 2. Bridge System: In this case, we consider the
1) True positive rate (TPR): The systems correctly classified structure represented in Fig. 2. The size of the system is p = 5.
as operative (i.e., Sensitivity). The structure function of the system is related to the function
2) True negative rate (TNR): The systems correctly classified
as nonoperative (i.e., Specificity). ψ(X1 , X2 , X3 , X4 , X5 ) = max{min{X1 , X4 },
3) Model predictive capability (MPC): The systems correctly min{X2 , X5 }, min{X1 , X3 , X5 }, min{X2 , X3 , X4 }}
classified.
4) AUC. where Xj denotes the state of component jth, j = 1, 2, 3, 4, 5.
From Table I, we see that the two approaches (nonparametric We simulate a sample of n = 100 such systems so the data
and parametric model) fit the data with acceptable accuracy, as consist of a matrix with p + 1 = 6 columns. The first p columns
the AUC is largely above 80%. The local model provides more report the states of the components as the value of independent
reliable predictions since the MPC increases up to 3% for all random variables with distribution U(0, 1). The state of the
cases. system is simulated considering a latent variable that is not
Table II provides the average and standard deviation of the directly observed and is assumed to be Y ∼ N (0.5, 0.2). Finally,
estimated β coefficients obtained along the M = 1000 samples. what we observe at the system level is the value of a Binomial
Denote similar to the former case.
M
We have simulated a total of M = 1000 samples. Table III lists

βj
(m) the results of model fitting considering both methods: parametric
Avj = (1/M )
m=1
and local-logistic approach. The values presented in the table are
the average along the M samples of the corresponding indexes
and at the table header.
 M
1/2 In this case, we notice that the results from the parametric fit
 2
βj
(m) are not so good as in the previous model. Although the measure
SDj = (1/M ) − Avj
m=1 of the goodness-of-fit given by the AUC index is acceptable,
the predictability capacity decreases for all cases (Y = 0 and
for j = 1, 2, 3, and being βj , the estimated coefficient for the
(m)
Y = 1) to 75%. The local approach also provides results slightly
jth component based on the mth sample, for j = 1, 2, 3 and worse results as compared with the previous example, however
m = 1, 2, . . . , M . we still get a good predictability power in this case.

Authorized licensed use limited to: University of New South Wales. Downloaded on August 09,2020 at 08:58:37 UTC from IEEE Xplore. Restrictions apply.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

GÁMIZ et al.: MACHINE LEARNING ALGORITHM FOR RELIABILITY ANALYSIS 9

TABLE V
NON SERIAL/PARALLEL SYSTEM

Fig. 3. Comparison test for components 1 and 3 at four performance levels in

the bridge system.

Fig. 4. Non series–parallel system.

Similar to the previous example, Table IV provides the av-

erage and standard deviation of the estimated β coefficients
obtained along the M = 1000 samples. Denote Avj and SDj ,
and βj , similar to the previous example, for the jth component
(m)

based on the mth sample, for j = 1, . . . , 5 and m = 1, 2, . . . , M .

The results provided by the parametric model and displayed
in Table IV suggest that component 3 is the least important for
the system performance, according to the Birnbaum criterion. Fig. 5. Non series–parallel system: Comparison test for components
importance.
Using the local-logit model, we have considered the following
particular testing problem:
state vector with all its elements equal to x0 ∈ [0, 1]. Denote x0
H0 : I1N B (x) = I3N B (x) as the p-dimensional vector whose elements are all equal to x0 ,
with p = 6. Focusing, for example, on components C and F , we
H1 : I1N B (x) = I3N B (x) set the following testing problem:
for x = x · 1 and taking different levels of components states, H0 : IFN B (x0 ) = ICN B (x0 )
specifically we take x = 0.6, 0.7, 0.8, 0.9. For each sample,
we have computed the value of the test statistic Z1,3 as ex- H1 : IFN B (x0 ) > ICN B (x0 ).
plained above. We have computed the associated p-value for
Table V presents a summary of values obtained for the test
each case, so we have a total of M = 1000 p-values. As usual
statistics ZF C along a simulation study based on a total of
this quantity is computed as p−value(m) = 2 Pr{Z > Z1,3 },
M = 1000 repetitions of the experiment. We have simulated
for m = 1, 2, . . . , M . We reject the null hypothesis when the
in this case samples of size n = 200.
p−value(m) < 0.05. Fig. 3 presents the boxplots obtained from
At a significance level of α = 0.05, we have to compare the
the p-values compared with the threshold value of 0.05. The
value of the statistics ZF C with the limit value 1.65. Fig. 5 gives
conclusion is that we cannot reject that components 1 and 3 are
a graphical summary of the simulation study where it is high-
equally important for the system performance at the significance
lighted by the red line the lower limit of the acceptance region.
level of α = 0.05.
Notice that there is more evidence in favor of the alternative
3) Example 3. A Non Series–Parallel System [20]: Our third
hypothesis as the level of performance increases. This means
case study is focused on an example taken from the book [20]
that unit F is less relevant for the system behavior at the lowest
where it is used to illustrate an algorithm to compute the relia-
levels.
bility of non serial/parallel systems. The system is represented
in Fig. 4 where six modules and five nodes are considered.
VII. HIGH-DIMENSIONAL SYSTEMS
The components of the system are the six modules named
{A, B, C, D, E, F }, since it is considered that the system works The complexity of systems encountered in the real world is
properly once the communication from node #1 to node #5 is usually higher than the models considered in this article, in the
established. sense of the number of single components in the system. In this
Let us consider that all components of the system are per- section, we first present a short simulation study to assess the
forming at the same level. In other words, we concentrate on a feasibility of the algorithm when system complexity increases.

Authorized licensed use limited to: University of New South Wales. Downloaded on August 09,2020 at 08:58:37 UTC from IEEE Xplore. Restrictions apply.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

10 IEEE TRANSACTIONS ON RELIABILITY

TABLE VI
COMPUTATIONAL COMPLEXITY OF THE ALGORITHM MEASURED AS THE
AVERAGED RUNNING TIME ALONG R = 100 EXECUTIONS OF THE
ALGORITHM AS A FUNCTION OF SYSTEM SIZE (p)

Fig. 6. Airplane cockpit system.

the reliability of each module from the state of its compo-

nents. For each j = 1, 2, . . . , J, we get R j (xj1 , . . . , xjp ), with
j
Then, we outline some proposals to be addressed in the case of p1 + p2 + · · · + pJ = p, the total size of the system, where J is
high-dimensional systems. We do not present a full development the total number of subsystems considered.
since it is beyond the scope of this article. Now, we can write the system reliability as a function
RS (R1 , R2 , . . . , RJ ), and use the algorithm as above to estimate
A. Study of Feasibility as a Function of System Complexity it. Then, for component i in subsystem j, the importance measure
S ∂Rj
To check the feasibility of our algorithm with the increasing as defined in Section V-B is calculated as IjN B = ∂R ∂Rj ∂xi .
complexity of the system, we have run a short simulation study To illustrate, in this section, we analyze the structure for the
aimed to measure the computational cost involved as the number information system in an airplane cockpit. It was first presented
of units in the system increases. We have considered two scenar- in a report from NASA [30]. Cockpit instrumentation is the
ios. In the first case (Case 1), we start with the series/parallel sys- required minimum instrumentation for general aviation aircraft
tem of Fig. 1, that is, p = 3. We focus on the testing problem H0 : flying for satisfying conditions of instrument flight rules as
I1NB (0.5) = I2NB (0.5) against H1 : I1NB (0.5) = I2NB (0.5). To defined in Federal Aviation. A reliability block diagram of the
draw samples we take the same specifications as in Example cockpit information system is shown in [22, Appendix]. For
1 (see Section VI-B), except for n = 250 and M = 100. For simplicity, we represent in Fig. 6 a sketch of the diagram without
each sample, we derive the p-value associated to the alternative specifically labeling the different units in the system. For a
hypothesis. We are interested in particular, in execution time, detailed description of the system, see [30] or [22]. The only
that is the difference between time at the start and the end of required information is the state of single units, the state of the
one single repetition. We obtain a total of M = 100 running system (0 or 1), and the state of each macrocomponent (0 or 1).
times and consider the mean value. Then, we build a system by The number of units in the system is p = 30, but the algorithm
adding a single unit disposed in series with the three-component is run for J = 4 different subsystems of respective sizes equal
series/parallel system. We consider the increased system, p = 4, to p1 = 4, p2 = 9, p3 = 6 and p4 = 10. It is worth to mention
and formulate exactly the same testing problem. We carry out that the internal configuration of the modules is not required in
the same simulation study and obtain the same summary of the procedure proposed.
execution times. Then, a new component is added in series to
the ensemble to build a system of size p = 5. We repeat all the C. Lack of Knowledge on System Structure
process until p = 10. The averaged execution times are reported
Unless some expertise reveals some rational of the system de-
in Table VI. The simulations have been run using a 3.60 GHz
sign that allows some strategies for model-fit as explained in the
Intel Core i5-8600 K.
previous section, the estimation problem has too many features
In the second scenario (Case 2), we also start from the se-
to deal with, the dimensions of the space grow exponentially and
ries/parallel system but additional units are located sequentially
the available data becomes sparse. This is the so-called curse of
(one each time) in parallel with the subsystem formed with units
dimensionality.
2 and 3 of Fig. 1. The testing problem is the same as in Case
In this case we assume that the unique information is the
1. The results obtained using a 2.00 GHz Intel Core
states of components and the system itself. Although we can
i3-6006 U are presented in the second row of Table VI. It
feed our model using all features at hand the result will have
can be noticed that the running time increases linearly with p.
so many terms that it becomes hard to interpret. Not to mention
the elevated computational cost that it will involve. Dimension
B. Example With Some Expert Knowledge
reduction is suggested but it can be difficult to decide which
In some cases, some information on the system design is features are truly worth keeping in our model. We need an
available making it possible to consider the modular decom- algorithm able to transform the large set of inputs into an ideal
position approach. A new coherent structure can be formed set of inputs. When the number of features is very high it is
where subsystems of the original system are treated as single reasonable to assume some correlations between the system
units and the system reliability is computed as a function (to be units can arise in the latent structure underlying the data and then
estimated) of the performance of these macrounits, which are principal component analysis (PCA) might be an appropriate
in turn derived from the information about the units that form tool. As a result of PCA, new independent inputs (or variables)
them. Treating each module separately. More specifically, with from a linear combination of existing variables. In summary, we
the only information about the state of a module (0 = failure are adding a new step to the algorithm.
and 1 = operative) as well as the state of the units inside the For a given dataset {X, Y }n×p , we have the following
module, we use the algorithm proposed to estimate (separately) conditions.

Authorized licensed use limited to: University of New South Wales. Downloaded on August 09,2020 at 08:58:37 UTC from IEEE Xplore. Restrictions apply.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

GÁMIZ et al.: MACHINE LEARNING ALGORITHM FOR RELIABILITY ANALYSIS 11

1) Isotonize the response to obtain {X, Y }n×p .  then, after a convenient Taylor expansion, we get
0 is β
2) Apply an algorithm for PCA to transform the data to a
reduced set {Z, Y }n×p , with p < p, and Z = XΓ, being  =  (β) +  (β)(β
0 = 0 (β)  − β) (7)
0 0
Γ the PCA coefficient-matrix of dimension p × p .
where we denote 0 the Hessian matrix, whose (j, k)-element is
3) Fit a local-logistic model based on the reduced dataset. ∂2
 denote the estimated vector of coefficients.
Let b ∂βj ∂βk , for j, k = 0, 1, . . . , p, i.e.,

The system reliability estimation, at a given x0 , is n



0 (β) = {c}−1 wh,0
0 (x) = eβ(x−x0 ) i=1
R 
1 + eβ(x−x 0)
(xi ) R0 (xi ) (1 − R0 (xi )) (xi − x0 ) (xi − x0 )t .
=b
with β  Γ. For convenience, we introduce the following notation: X, is the
design matrix with dimensions n × (p + 1)
VIII. CONCLUSION ⎛ ⎞
1 x11 − x01 . . . x1p − x0p
In this article, we built a statistical model able to predict the ⎜ 1 x21 − x01 . . . x2p − x0p ⎟
⎜ ⎟
reliability of the system based on a dataset. We used a supervised X=⎜. .. .. .. ⎟
⎝ .. . . . ⎠
learning algorithm based on local LR and the model was trained
considering cross-validation techniques. The method was en- 1 xn1 − x01 . . . xnp − x0p
tirely data-driven, not relying on any parametric assumptions. the vector of transformed responses
Once the model was validated through sensitivity analysis, a
function that classifies the system in one of two categories,  0 = (Y1 − R0 (x1 ), . . . , Yn − R0 (xn ))t
Y
operative or failed, was defined based on the knowledge of
components states. Furthermore, we presented a statistical test the diagonal matrix Wn whose ith element is Wn (i, i) =
to decide the order of importance of components in terms of the wh,0 (xi ), for i = 1, . . . , n; and the diagonal matrix Vn whose
effect each one has on the system performance. The importance ith element is Vn (i, i) = R0 (xi )(1 − R0 (xi )), for i = 1, . . . , n.
of a component was defined by the effect an increment in the Then, the above-mentioned equations can be written in the
state of the component has on the system reliability. To illustrate matrix form as
the method, we carried out an extensive simulation study. We
0 (β) = Xt Wn Y0 and 0 (β) = Xt Wn Vn X.
considered three practical scenarios with different complexity,
that is, the size of the systems considered is, respectively, Then, from (7), we have
p = 3, 5, and 6. For the simplest structure (the series–parallel
combination with p = 3), the parametric approach gives a sat-  − β = − ( (β))−1  (β) .
β (8)
0 0
isfactory solution, as can be deduced from the results provided
by sensitivity analysis. However, for the bridge system (p = 5), We use this expression to approximate the bias and the vari-
we see that the accuracy given by the parametric model is not so ance of the estimator.
good and so the conclusions from the model can be misleading.
As stated in [21], importance measures were used in various Bias:
fields to evaluate the relative importance of various objects such According to (4) and (8), to approximate an expression for
as components in a system. The absolute values of importance the bias we have to handle with the following:
measures may not be as important as their relative rankings. In
t
their book, the authors highlighted the necessity of addressing eg(xi ) eβ (xi −x0 )
E [Yi ] − R0 (xi ) = − t
importance ranking of components to provide the most efficient 1 + eg(xi ) 1 + eβ (xi −x0 )
way to identify weaknesses in the system structure, which can
be decisive to minimize the damage caused by an unpredictable since xi − x0 < h, then g(x) = β t (xi − x0 ) + o(h). Now
ex
accident. This can be of primary importance for some industries we can use a Taylor expansion for the function f (x) = 1+e x

like power plants, nuclear plants, or electric power industry (see, getting that
for example, [9], [14], or [21]). These fields of application can t
eβ (xi −x0 )
benefit from the contribution of this article. E [Yi ] − R0 (xi ) =   · o(h)
1 + eβ
t
(xi −x0 ) 2

APPENDIX = Vn (i, i) · o(h)

ASYMPTOTIC PROPERTIES OF THE LOCAL ESTIMATOR
where o(h)/h → 0, as h → 0. Finally, we have
We follow similar arguments as in [17] to approximate the
 
bias and variance of the estimator of the logit and its partial  = − Xt Wn Vn X −1 Xt Wn Vn 1n o(h)
Bias(β) (9)
derivatives. Let us denote β  the vector the estimates of g and its
first partial derivatives evaluated at x0 . Since the maximizer of with 1n = (1, 1, . . . , 1)t .

Authorized licensed use limited to: University of New South Wales. Downloaded on August 09,2020 at 08:58:37 UTC from IEEE Xplore. Restrictions apply.
 This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.

12 IEEE TRANSACTIONS ON RELIABILITY

V ariance: [12] O. Burdakov, A. Grimvall, and M. Hussian, “A generalised PAV algorithm

for monotonic regression in several variables,” in COMPSTAT: Proceed-
Again using expression (8), we have ings in Computational Statistics, J. Antoch, Ed. New York, NY, USA:
Springer, 2004, pp. 761–767.
 = ( (β))
Var(β)
−1
Var [0 (β)] (0 (β))
−1
[13] F. P. A. Coolen and T. Coolen-Maturi, “The structure function for system
0
reliability as predictive (imprecise) probability,” Rel. Eng. Syst. Saf.,
we obtain vol. 154, pp. 180–187, 2016.
[14] J. F. Espiritu, D. W. Coit, and U. Prakash, “Component criticality impor-
Var [0 (β)] = (Wn X)t Var (Y0 ) Wn X tance measures for the power industry,” Elect. Power Syst. Res., vol. 77,
pp. 407–420, 2007.
= Xt Wn Vn Wn X. [15] J. Fan and I. Gijbels, Local Polynomial Modelling and Its Applications.
London, U.K.: Chapman & Hall, 1996.
Then, we have [16] J. Fan, N. E. Heckman, and M. P. Wand, “Local polynomial kernel
regression for generalized linear models and quasi-likelihood functions,”

Var(β) J. Amer. Statist. Assoc., vol. 90, no. 129, pp. 141–150, 1995.
[17] J. Fan, M. Farmen, and I. Gijbels, “Local maximum likelihood estimation
 −1 t  −1
= Xt Wn Vn X X W n V n W n X Xt W n V n X . and inference,” J. Roy. Statist. Soc. B, vol. 60, no. 3, pp. 591–608, 1998.
[18] M. Frolich, “Non-parametric regression for binary dependent variables,”
(10) Econometrics J., vol. 9, no. 3, pp. 511–540, 2006.
[19] M. L. Gamiz and M. D. Martinez-Miranda, “Regression analysis of the
ACKNOWLEDGMENT structure function for reliability evaluation of continuous-state system,”
Rel. Eng. Syst. Saf., vol. 95, pp. 134–142, 2010.
The authors are grateful for constructive comments from three [20] I. Koren and C. M. Krisna, Fault-Tolerant Systems. San Mateo, CA, USA:
Morgan Kaufman, 2007.
anonymous reviewers and the associate editor. [21] W. Kuo and X. Zhu, Importance Measures in Reliability, Risk and Opti-
mization: Principles and Applications. Hoboken, NJ, USA: Wiley, 2012.
REFERENCES [22] M. Lenz and J. Rhodin, Reliability Calculations for Complex Systems.
Linkoping, Sweden: Linkoping Univ., 2011.
[1] E. F. Alsina, M. Chica, K. Trawinski, and A. Regattieri, “On the use of [23] I. Lisnianski and G. Levitin, Multi-State System Reliability: Assessment,
machine learning methods to predict component reliability from data- Optimization and Applications. Singapore: World Scientific, 2001.
driven industrial case studies,” Int. J. Adv. Manuf. Technol., vol. 94, [24] M. C. Meyer, “Inference for multiple isotonic regression,” 2010. [Online].
pp. 2419–2433, 2018. Available: [Link]/research/TechnicalReports/2010/2010
[2] K. P. Amrutkar and K. K. Kamalja, “An overview of various importance [Link]
measures of reliability system,” Int. J. Math., Eng. Manage. Sci., vol. 2, [25] M. C. Meyer, “Isotonicmultiple.R,” 2012. [Online]. Available: https:
no. 3, pp. 150–171, 2017. //[Link]/ meyer/isotonicmultiple.R
[3] T. Aven, “On performance measures for multistate monotone systems,” [26] M. C. Meyer, “Semi-parametric additive constrained regression,” J. Non-
Rel. Eng. Syst. Saf., vol. 41, no. 2, pp. 59–266, 1993. parametric Statist., vol. 25, no. 3, pp. 715–730, 2013.
[4] T. Aven and U. Jensen, Stochastic Models in Reliability. New York, NY, [27] H. Mukerjee, “Monotone nonparametric regression,” Ann. Statist., vol. 16,
USA: Springer, 2013. no. 2, pp. 741–750, 1988.
[5] T. Aven and T. E. Nockland, “On the use of uncertainty importance [28] A. B. Musa, “Comparative study on classification performance between
measures in reliability and risk analysis,” Rel. Eng. Syst. Saf., vol. 95, support vector machine and logistic regression,” Int. J. Mach. Learn.
no. 2, pp. 127–133, 2010. Cybern., vol. 4, pp. 13–24, 2013.
[6] T. Ayer, J. Chhatwal, O. Alagoz, E. K. Charles, R. W. Woods, and [29] B. Natvig, Multi-State Systems Reliability Theory and Applications. Hobo-
E. S. Burnside, “Comparison of logistic regression and artificial neural ken, NJ, USA: Wiley, 2011.
network models in breast cancer risk estimation,” RadioGraphics, vol. 30, [30] D. Pettit and A. Turnbull, General Aviation Aircraft Reliability Study.
pp. 13–22, 2010. Hampton, VA, USA: NASA/CR-2001-210647, 2001.
[7] L. A. Baxter, “Continuum structures I,” J. Appl. Probab., vol. 21, [31] M. Rausand and A. Hoyland, System Reliability Theory: Models, Statisti-
pp. 802–815, 1984. cal Methods, and Applications, 2nd ed. Hoboken, NJ, USA: Wiley, 2004.
[8] Z. W. Birnbaum, “On the importance of different components in a mul- [32] S. Si, M. Liu, Z. Jiang, T. Jin, and Z. Cai, “System reliability allocation
ticomponent system,” in Multivariate Analysis-II, P. R. Krishnaiah, Ed. and optimization based on generalized Birnbaum importance measure,”
New York, NY, USA: Academic, 1969, pp. 581–592. IEEE Trans. Rel., vol. 68, no. 3, pp. 831–843, Sep. 2019.
[9] S. Bisanovic, M. Samardzic, and D. Aganovic, “Application of component [33] W. Yoo, B. A. Ference, M. L. Cote, and A. Schwartz, “A comparison of
criticality importance measures in design scheme of power plants,” Int. J. logistic regression, logic regression, classification tree, and random forests
Elect. Comput. Eng., vol. 6, no. 1, pp. 63–70, 2016. to identify effective gene-gene and gene-environmental interactions,” Int.
[10] R. D. Brunelle and K. C. Kapur, “Continuous-state system-reliability: An J. Appl. Sci. Technol., vol. 2, no. 7, pp. 268–294, 2012.
interpolation approach,” IEEE Trans. Rel., vol. 47, no. 2, pp. 181–187, [34] E. Zaitseva and V. Levashenko, “Construction of a reliability structure
Jun. 1998. function based on uncertain data,” IEEE Trans. Rel., vol. 65, no. 4,
[11] H. D. Brunk, “Maximum likelihood estimates of monotone parameters,” pp. 1710–1723, Dec. 2016.
Ann. Math. Statist., vol. 26, pp. 607–616, 1955.

Authorized licensed use limited to: University of New South Wales. Downloaded on August 09,2020 at 08:58:37 UTC from IEEE Xplore. Restrictions apply.