Bayesian Nonparametric System Reliability using Sets of Priors

An imprecise Bayesian nonparametric approach to system reliability with multiple types of components is developed. This allows modelling partial or imperfect prior knowledge on component failure distributions in a flexible way through bounds on the functioning probability. Given component level test data these bounds are propagated to bounds on the posterior predictive distribution for the functioning probability of a new system containing components exchangeable with those used in testing. The method further enables identification of prior-data conflict at the system level based on component level test data. New results on first-order stochastic dominance for the Beta-Binomial distribution make the technique computationally tractable. Our methodological contributions can be immediately used in applications by reliability practitioners as we provide easy to use software tools.


Introduction
System reliability analysis is concerned with estimating the lifetime T sys of complex systems.Usually, the goal is to determine the system reliability function R sys (t) = P (T sys > t) based on the lifetime distributions of system components.
A critique of the methodological approach to a reliability analysis may often encompass a few common concerns.First, in a parametric setting, there may be no particularly strong reason to believe that the small part of component model space covered by a particular probability distribution necessarily contains the 'true' component lifetime distribution.Further, Bayesian methods may be invoked in order to incorporate expert opinion or other knowledge which falls outside the specific testing data under consideration.The classic concern here is in whether one can truly express ones beliefs with the exactness a prior distribution requires.Finally, it would be valuable in application to have a means of identifying when the prior choice is having a strong effect and when not.Any method hoping to address these concerns must do so whilst enabling realistic system models (with heterogeneous component types) and remain computationally tractable.
Herein, we make steps toward addressing these concerns by developing a nonparametric method which utilises imprecise probability [3,20] to model more vague or imperfect prior beliefs using upper and lower probabilities.This overcomes the concern about component lifetimes being outside a particular parametric family, uses a more flexible prior modelling framework and leads to an easy method of detecting conflicts between prior assumptions and observed failure times in test data.In the general context of Bayesian methods, this phenomenon is known as prior-data conflict, see, e.g., [13] or [6].
Furthermore, the method is based on the survival signature [8], a recent development which naturally accommodates heterogeneous component types laid out in any arbitrary manner.By separating the (time-invariant) system structure from the time-dependent failure probabilities of components, it allows straightforward and efficient computation of R sys (t).
Our imprecise probability approach provides bounds for R sys (t) by computing, for each t in an arbitrarily fine grid of time points T , the posterior predictive probability interval for the event T sys > t.Assuming the number of functioning components for each type and time t as binomially distributed, the intervals are derived from an imprecise Bayesian model using sets of conjugate Beta priors which allow to specify weak or partial prior information in an intuitive way.The width of the resulting posterior predictive probability intervals reflects the precision of the corresponding probability statements: a short range indicates that the system functioning probability can be quantified quite precisely, while a large range will indicate that our (probabilistic) knowledge is indeterminate.In particular, prior-data conflict leads to more cautious probability statements: When there is not enough data to overrule the prior, it is unclear whether to put more trust to prior assumptions or to the observations, and posterior inferences clearly reflect this state of uncertainty by larger ranges.
While the use of imprecise probability methods can often lead to tractability issues, new results on first-order stochastic dominance for the Beta-Binomial distribution keep the need for numerical optimization in our model to a minimum.
In Section 2 we review the survival signature and in Section 3 we review the nonparametric approach to Bayesian reliability analysis upon which our work builds [2].Section 4 details the reparameterisation of that approach which enables the natural formulation of the system reliability bounds, leading to nice closed form results in some later theory.Section 5 lays the ground work to incorporate imprecise probability, culminating in the main results and contributions of this work, detailed in Section 6. Section 7 provides details on the software contributions of this work and shows two worked examples demonstrating the practicality and usefulness of the method.

Survival Signature
In the mathematical theory of reliability, the main focus is on the functioning of a system given the functioning, or not, of its components and the structure of the system.The mathematical concept which is central to this theory is the structure function [4].For a system with m components, let state vector x = (x 1 , x 2 , . . ., x m ) ∈ {0, 1} m , with x i = 1 if the ith component functions and x i = 0 if not.The labelling of the components is arbitrary but must be fixed to define x.The structure function φ : {0, 1} m → {0, 1}, defined for all possible x, takes the value 1 if the system functions and 0 if the system does not function for state vector x.Most practical systems are coherent, which means that φ(x) is non-decreasing in any of the components of x, so system functioning cannot be improved by worse performance of one or more of its components.The assumption of coherent systems is also convenient from the perspective of uncertainty quantification for system reliability.It is further logical to assume that φ(0) = 0 and φ(1) = 1, so the system fails if all its components fail and it functions if all its components function.
For larger systems, working with the full structure function may be complicated, and one may particularly only need a summary of the structure function in case the system has exchangeable components of one or more types.We use the term 'exchangeable components' to indicate that the failure times of the components in the system are exchangeable [12].Coolen and Coolen-Maturi [8] introduced such a summary, called the survival signature, to facilitate reliability analyses for systems with multiple types of components.In case of just a single type of components, the survival signature is closely related to the system signature [17], which is well-established and the topic of many research papers during the last decade.However, generalization of the signature to systems with multiple types of components is extremely complicated (as it involves ordering order statistics of different distributions), so much so that it cannot be applied to most practical systems.In addition to the possible use for such systems, where the benefit only occurs if there are multiple components of the same types, the survival signature is arguably also easier to interpret than the signature.
Consider a system with K ≥ 1 types of components, with m k components of type k ∈ {1, . . ., K} and K k=1 m k = m.Assume that the random failure times of components of the same type are exchangeable [12].Due to the arbitrary ordering of the components in the state vector, components of the same type can be grouped together, leading to a state vector that can be written as x = (x 1 , x 2 , . . ., x K ), with x k = (x k 1 , x k 2 , . . ., x k m k ) the sub-vector representing the states of the components of type k.
The survival signature for such a system, denoted by Φ(l 1 , . . ., l K ), with l k = 0, 1, . . ., m k for k = 1, . . ., K, is defined as the probability for the event that the system functions given that precisely l k of its m k components of type k function, for each k ∈ {1, . . ., K} [8].Essentially, this creates a Kdimensional partition for the event T sys > t, such that R sys (t) = P (T sys > t) can be calculated using the law of total probability: where C k t ∈ {0, 1, . . ., m k } denotes the random number of components of type k functioning at time t.
For calculating the survival signature based on the structure function, observe that there are denote the set of these state vectors for components of type k and let S l 1 ,...,l K denote the set of all state vectors for the whole system for which m k i=1 x k i = l k , k = 1, . . ., K. Due to the exchangeability assumption for the failure times of the m k components of type k, all the state vectors x k ∈ S k l k are equally likely to occur, hence [8] It should be emphasized that when using the survival signature, there are no restrictions on dependence of the failure times of components of different types, as the probability P ( K k=1 {C k t = l k }) can take any form of dependence into account, for example one can include common-cause failures quite straightforwardly into this approach [9].However, there is a substantial simplification if one can assume that the failure times of components of different types are independent, and even more so if one can assume that the failure times of components of type k are conditionally independent and identically distributed with CDF F k (t).With these assumptions, we get We will employ both assumptions in this paper, leading to C k t having a Beta-Binomial distribution, giving us a closed form expression for P (C k t = l k ) for all t, k, and l k .
The main advantage of the survival signature, in line with this property of the signature for systems with a single type of components [17], is that the information about the system structure is fully separated from the information about functioning of the components, which simplifies related statistical inference as well as considerations of optimal system design.In particular for study of system reliability over time, with the structure of the system, and hence the survival signature, not changing, this separation also enables relatively straightforward statistical inferences.
There are several relatively straightforward generalizations of the use of the survival signature.The probabilities for the numbers of functioning components can be generalized to lower and upper probabilities, as e.g.done by Coolen et al. [10] within the nonparametric predictive inference (NPI) framework of statistics [7], where lower and upper probabilities for the events C k = l k are inferred from test data on components of the same types as those in the system.This is an approach that is also followed in the current paper, but with the use of generalized Bayesian inference instead of NPI.Like Coolen et al. [10], we will utilize the monotonicity of the survival signature for coherent systems to simplify computations.

Nonparametric Bayesian Approach for Component Reliability
Let us denote the random failure time of component number i of type k by T k i , i = 1, . . ., m k .The failure time distribution can be written in terms of the cdf F k (t) = P (T k i ≤ t), or in terms of the reliability function , also known as the survival function.For a nonparametric description of R k (t), we consider a set of time points t, t ∈ T = {t 1 , . . ., t max }.
At each time point t, the operational state of a single component of type k is Bernoulli distributed (functioning: 1, failed: 0) with parameter p k t , so that That is, We can also express this failure time distribution through the probability mass function (pmf) and discrete hazard function, The time grid T can be chosen to be appropriately dense for the application at hand, with the natural extension between grid points by taking R k (•) to be the right continuous step function induced by the grid values, R k (t) = p k t j , t ∈ [t j , t j+1 ), or by taking p k t j and p k t j+1 as upper and lower bounds for R k (t), t ∈ [t j , t j+1 ).
The independence assumption for components of the same type immediately implies that the number of functioning components of type k in the system is binomially distributed, t 's can, in theory, be directly chosen to arbitrarily closely approximate any valid lifetime pdf on [0, ∞), for example matching a bathtub curve for the corresponding hazard rate h k (t j ).Naturally, p k t j ≥ p k t j+1 should hold (assuming no repair).However, such direct specification is non-trivial, neglects any inherent uncertainty in the particular choice, and cannot be easily combined with test data.To account for the uncertainty, one can express knowledge about p k t through a prior distribution.A convenient and natural choice is p k t ∼ Beta(α k t , β k t ), particularly because in a Bayesian inferential setting this is the conjugate prior which leads to a Beta posterior.
Let the lifetime test data collected on component k be t k = (t k 1 , . . ., t k n k ).At each fixed time t ∈ T , this corresponds to an observation from the Binomial model described above, s k t = n k i=1 I(t k i > t).The posterior is then . The combination of a Binomial observation model with a Beta prior is often called Beta-Binomial model.
The posterior predictive distribution for the number of components surviving at time t in a new system, based on the lifetime data and the prior information, is then given by a so-called Beta-Binomial distribution, That is, we have , where B(•, •) is the Beta function.This posterior predictive distribution is essentially a Binomial distribution where the functioning probability p k t takes values in [0, 1] with the probability given by the posterior on p k t .Consequently, in order to calculate the system reliability according to (1), for each component type k we need lifetime data t k , and have to choose 2 × |T | parameters to specify the prior distribution for the discrete survival function of T k i .In total, values for 2 × |T | × K parameters must be chosen.

Reparametrisation of the Beta Distribution
The parametrisation of the Beta distribution used above is common, and allows α k t and β k t to be interpreted as hypothetical numbers of functioning and failed components of type k at time t, respectively.However, when we generalise to sets of priors in the sequel, it is useful to consider a different parametrisation.
For clarity of presentation we will temporarily drop the super-and subscript k and t indices for component type and time.Instead of α and β, we consider the parameters n (0) ∈ [0, ∞) and y (0) ∈ [0, 1], where or equivalently, α = n (0) y (0) and β = n (0) (1 − y (0) ).The upper index (0) is used to identify these as prior parameter values, in contrast to their posterior values n (n) and y (n) obtained after observing n failure times (see below).n (0) and y (0) are sometimes called canonical parameters, identified from rewriting the density in canonical form; see for example [5, pp. 202  From the properties of the Beta distribution, it follows that y (0) = E[p] is the prior expectation for the functioning probability p, and that larger n (0) values lead to greater concentration of probability measure around y (0) , since Var(p) = y (0) (1−y (0) ) n (0) +1 .Consequently, n (0) represents the prior strength and moreover can be directly interpreted as a pseudocount, as will become clear.Indeed, consider the posterior given that s out of n components function: by conjugacy p | s is Beta distributed with updated parameters Thus, after observing that s out of n components function (at time t), the posterior mean y (n) for p is a weighted average of the prior mean y (0) and s/n (the fraction of functioning components in the data), with weights proportional to n (0) and n, respectively.Therefore n (0) takes on the same role for the prior mean y (0) as the sample size n does for the observed mean s/n, leading to the notion of it being a pseudocount.
Reintroducing time and component type indices, the posterior predictive Beta-Binomial probability mass function (pmf) can be written in terms of the updated parameters as with the corresponding cumulative mass function (cmf) given by The parameterisation in terms of prior mean and prior strength (or pseudocount) makes clear that in this conjugate setting, learning from data corresponds to averaging between prior and data.This form is attractive not only because it enhances the interpretability of the model and prior specification, but crucially it also makes clear what should be a serious concern in any Bayesian analysis: when observed data differ greatly from what is expressed in the prior, this conflict is simply averaged out and is not reflected in the posterior or posterior predictive distributions.
As a simple example, imagine that we expect p k t to be about 0.75 for a certain k and t, so we choose y (0) k,t = 0.75, and that we value this choice of mean functioning probability with n (0) k,t = 8, i.e., equivalently to having seen 8 observations with a mean 0.75.If we observe n k = 16 components of type k in the test data and s k t = 12 function at time t, then s k t /n k = 0.75 as we expect, so that the updated parameters are n k,t = 0.75.However, in contrast, unexpectedly observing that no component functions at time t instead leads to parameters n k,t = 0.25.The prior and the posteriors based on these two scenarios are depicted in the left panels of Figure 1, along with their corresponding predictive Beta-binomial pmf and cmf for the case m k = 5 (right panels).
Due to symmetry, we see that both posteriors have the same variance, although arising from two fundamentally different scenarios.Posterior 1 is based on data exactly according to prior expectations; the increase in confidence on p k t ≈ 0.75 is reflected in a more concentrated posterior density, and the posterior predictive is changed only slightly.However, it may be cause for concern to see the same degree of confidence in Posterior 2, which is based on data that is in sharp conflict with prior expectations.Posterior 2 places most probability weight around 0.25, averaging between prior expectation and data, with the same variance as Posterior 1. Accordingly, rather than conveying the conflict between observed and expected functioning probabilities with increased variance, Posterior 2 instead gives a false sense of certainty.
To enable diagnosis of when this undesirable behaviour occurs, we propose to use an imprecise probability approach based on sets of Beta priors, described in the following section.

Sets of Beta Priors
As was shown by Walter and Augustin [22], we can have both tractability and meaningful reaction to prior-data conflict by using sets of priors M (0) k,t produced by parameter sets Π (0) k,t according to Bayes' Rule.This element-by-element updating can be rigorously justified as ensuring coherence [20, §2.5], and was termed "Generalized Bayes' Rule" by Walley [20,§6.4].Due to conjugacy, k,t according to (4), leading to the set of updated parameters Examples for parameter sets Π p t k q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 k,t = 0.75, and posteriors based on n k t = 16 observations with s k t = 12 (Posterior 1) and s k t = 0 (Posterior 2), respectively.Data for Posterior 1 confirm prior assumptions, while data for Posterior 2 are in conflict with the prior.However, this conflict is averaged out, and Posterior 1 and Posterior 2 have the same spread, both in the posterior pdf/cdf and the posterior predictive pmf/cmf, such that Posterior 2 gives a false sense of certainty despite the massive conflict between prior and data.
k,t has the 'spotlight' shape (left); in case of priordata conflict ( k,t has the 'banana' shape (right), leading to a large degree of imprecision in the y k,t , thus reflecting increased uncertainty about the functioning probability p k t due to the conflict between prior assumptions and observed data.
desirable model properties and ease of elicitation (see Walter [21, pp. 123f] or Troffaes et al. [19]).For each component type k and time point t, one need only specify the four parameters n k,t (so in total 4 × |T | parameters are needed to define the set of prior distributions on the survival function of each component).
A desirable inference property arising from this setup is that the posterior parameter set Π (n) k,t is not rectangular in the way that the prior parameter set is.Indeed, the shape of Π (n) k,t depends on the presence or absence of priordata conflict, which is naturally operationalised as k,t ]: that is, prior-data conflict is defined to occur when, at time t, the observed fraction of functioning components is outside its a priori expected range.
First, in the absence of prior-data conflict, Π k,t shrinks in the y k,t dimension; how much it shrinks depending on n , leading to the so-called spotlight shape depicted in Figure 2   k,t interval compared to the no conflict case, reflecting the extra uncertainty due to prior-data conflict.In other words, the posterior sets make more cautious probability statements about p k t , as desired in this scenario.Based on these shapes and (4), it is possible to deduce the following expressions for the lower and upper bounds of y Note that the lower bound for y k,t ], both the lower and the upper bounds for y k,t , corresponding to the spotlight shape.However, when k,t ], the banana shape indicates that one of the bounds for y The different locations and sizes of Π (n) k,t in the conflict versus no conflict case are then, in turn, also reflected in the corresponding sets of Beta cdfs and Beta-Binomial cmfs.As an example, those corresponding to the parameter sets in Figure 2 are depicted in Figure 3.
In the no conflict case (Posterior 1, top row), the reduction of the y k,t leads to a much smaller set of Beta and Beta-Binomial distributions.For example, the range of predictive probabilities that two out of a set of five components of type k function at time t has changed from [0.10, 0.28] a priori to [0.11, 0.14] a posteriori.This reflects the gain in precision due to test data in accordance with prior assumptions.
In contrast, for the prior-data conflict case (Posterior 2, bottom row), the wide y k,t leads to a set of Beta and Beta-Binomial distributions that is much larger than in the no conflict case.Here, the range of posterior predictive probabilities that two out of a set of five components of type k function at time t is now [0.86, 1.00] a posteriori, i.e., less precise than in the no conflict case.Using sets of Beta priors, the resulting set of posterior predictive Beta-Binomial distributions reflects the precision of prior information, the amount of data, and prior-data conflict.p t k cdf q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 l k cmf Item2 q q q q q q q q Prior Posterior k,t (or their posterior counterparts) as solid lines.The top row depicts the set of prior cdfs/cmfs and the set of posterior cdfs/cmfs for the case where data confirm prior assumptions (see left panel of Figure 2); the bottom row depicts the (identical) set of prior cdfs/cmfs and the set of posterior cdfs/cmfs in case of prior-data conflict (see right panel of Figure 2).The set of posterior cdfs and cmfs is much larger in case of prior-data conflict: uncertainty due to this conflict is reflected through increased imprecision.
Furthermore, with sets of Beta priors it is also possible to express prior ignorance by letting y (0) k,t → 0 and y (0) k,t → 1 for some or all t ∈ T .(Note that it is not advisable to choose y (0) k,t = 0 and y (0) k,t = 1, as this can lead to improper posterior predictive distributions.For example, at any t < min(t k ), we would have y (n) k,t = 1, leading to one argument of the Beta function in the denominator of (5) being zero.)These limits for y (0) k,t imply we are only prepared to give trivial bounds for the functioning probability and do not wish to commit to any specific knowledge about p k t a priori.This provides a more natural choice of 'noninformative' prior over [0, 1] than the usual choice of a Beta prior with k,t = 0.5).Such a prior for all t ∈ T actually reflects a belief that the component reliability function is on average 1/2 for all t, which is not an expression of ignorance, but rather a very specific (and arguably peculiar) prior belief.
In a near-noninformative setting, the choice of n (0) k,t is not relevant, because (8) implies both lower and upper bound for y k,t < 1 can be chosen such that k,t ] for all t ∈ min(t k ), max(t k ) .Naturally, one cannot have prior-data conflict in cases of near prior ignorance.

Sets of System Reliability Functions
The elements reviewed and extended above culminate hereinafter in the primary contribution of the current work, providing a framework in which the nonparametric Bayesian system reliability approach developed in [2] is extended to sets of system reliability functions by incorporating the sets of priors approach of Walter and Augustin [22].This allows for partial or vague specification of prior component reliability functions, and enables diagnosis of prior-data conflict which is consequential at the system level.

Computation of bounds
To obtain the lower and upper bound for the system reliability function R sys (t), we now need to minimise and maximise Equation (1) over Π (0) 1,t , . . ., Π (0) K,t for each t, where the posterior predictive probabilities for C k t are given by the Beta-Binomial pmf (5).We therefore have and similarly maximising for R sys (•).Note that Φ(•) is non-decreasing in each of its arguments l 1 , . . ., l K , thus if there is first-order stochastic ordering on k,t , s k t ) for each k, then this ordering can be used to determine the elements of Π (0) k,t which minimise and maximise the overall system reliability function without resorting to computationally expensive exhaustive searches or numerical optimisation.
We therefore start by providing the following result, where indices are suppressed for readability.We use ≥ st to denote first-order stochastic dominance.
Theorem 1.Let β y denote the Beta-Binomial distribution with probability mass function parameterised as: with n, m, s, and N fixed and unknown.Then β y ≥ st β y ∀ y > y with y, y ∈ (0, 1).
The proof is provided in Appendix A, p.28.Consequently, for each component, the posterior predictive Beta Binomial distributions with larger prior functioning probability stochastically dominate those with smaller prior functioning probability, providing rigorous proof which accords with intuition.Applying this result to the sets of system reliability functions, together with the monotonicity in the survival signature, means that R sys (•) is attained when y (0) k,t = y (0) k,t and R sys (•) is attained when y    k,t is more subtle, because stochastic dominance is not guaranteed at a single value.The following Theorem provides simple sufficient conditions under which an upper or lower limit has firstorder stochastic dominance and has virtually no computational overhead to test.
Theorem 2. Let β n denote the Beta-Binomial distribution with probability mass function parameterised as: with y, m, s, and N fixed and unknown.Then, The proof is provided in Appendix A, p.29.If s N +m−1 < y < s+m−1 N +m−1 then Theorem 2 cannot determine stochastic dominance.The following Lemma which is slightly more computationally costly, but still much faster than an exhaustive search, may be able to determine first-order stochastic dominance in such situations.
The proof is provided in Appendix A, p.31.In the cases where neither Theorem 2 or Lemma 3 apply, the entire posterior system reliability function must be optimised to find the minima/maxima.In practice, in the examples to be presented in the sequel, Theorem 2 and Lemma 3 do provide guarantees of first-order stochastic dominance for the vast majority of time points, t, substantially lowering the computational costs of performing the minimisation/maximisation involved in finding the sets of system reliability functions compared to either numerical optimisation or an exhaustive grid search (which would get exponentially slower in the number of different components). Thus, where The result for R sys (•) is completely analogous.It is interesting to note that if m k = 1 the bounds are sharp on stochastic dominance.In particular, when indicates the lower bound is not in conflict with the observed data, whilst y (0) k,t > s k t n k is in conflict since the observed empirical probability of functioning at time t is below the prior lower bound.Consequently, note that n (0) k,t is used only when the prior comes into conflict with the data.Since n (0) k,t controls the prior certainty, this accords with the intuition that the least certain prior bound is invoked when in a conflict setting and the more certain prior bound used when the data agrees.

Prior parameter choice
In the following, we will give some guidelines on how to choose the parameter sets Π k,tmax which define the set of prior discrete reliability functions for components of type k.We advocate that this is much easier in terms of n (0) and y (0) than it would be in terms of α and β.
As mentioned in Section 3, the functioning probabilities p k t must satisfy p k t j ≥ p k t j+1 .This naturally translates to conditions on the prior for p k t , so that for example y should hold.Because s k t /n k is decreasing in t, the weighted average property of the update step in Equation ( 4) for y k,t ensures that y . In situations where one has a high degree of certainty about the functioning probability for low t, but less certainty about what happens for larger t, then one can let y (0) k,t drop to (almost) 0, but clearly y , and k,t j ] (meaning there is no prior-data conflict), then should there be no observed failures in [t j , t j+1 ], so that s k t j+1 /n k = s k t j /n k , then and Again, this follows from (4), the weighted average property.It is possible to construct similar examples with regard to the lower bound n (0) k,t j .Therefore, we advise taking the same n (0) k,t bounds for all t as far as possible.If they do change, it must be very gradual and we recommend diagnosing any problems as above.
Generally, the interpretation as pseudocount or prior strength should guide the choice of bounds for n k,t ] will be dominated by the location of s k t /n k .Furthermore, the length of [y (n) k,t , y k,t values.Specifically, in a no-conflict situation, when n k,t ].In contrast, high values for n (0) k,t will lead to slower learning and wider y   k,t with help of the half-width rule as described above.
As mentioned in Section 5, it is not advisable to choose y (0) k,t = 0 and y (0) k,t = 1.For any t ∈ min(t k ), max(t k ) , this can lead to improper posterior predictive distributions.However, it is possible to choose values close to 0 and 1, respectively, and due to the linear update step (4) for y  k,t = 0.9999.Likewise, our nonparametric method does not cause unintuitive tail behaviour as some parametric methods do; there is no problem, for example, with assigning y (n) k,t near-zero for large t if prior knowledge suggests so.While it is possible to set the bounds y (0) k,t and y (0) k,t for each t ∈ T individually, in practice this will be often too time-consuming when T forms a dense grid.Switching to a coarser time grid will waste information from data, as then failure times in the test data are rounded up to the next t ∈ T .
In the examples here we elicit bounds for a subset of T and fill up the time grid with the least committal bounds, i.e., taking y (0) k,t equal to last (in the time sequence) elicited y (0) k,t , and likewise y (0) k,t equal to next (in time sequence) elicited y (0) k,t .A possible elicitation procedure in this vein could be to start with eliciting y (0) k,t bounds for a few 'central' time points t, filling up the grid as described above accordingly, and then to further refine the obtained bounds as deemed necessary by the expert.

Software
The methods of this paper have been implemented in the R [16] package ReliabilityTheory [1], providing an easy to use interface for reliability practitioners.The primary function, which computes the upper and lower posterior predictive system survival probabilities as in (10), is named nonParBayesSystemInferencePriorSets().The user specifies the times at which to evaluate the bounds, the survival signature (Φ(•)), the component test data (t 1 , . . ., t K ), and the prior parameter set for each component type and time (Π Table 1: Survival signature for the bridge system from Figure 4, omitting all rows with T3 = 0, since Φ = 0 for these. the CPU has multiple cores and making automatic use of the theoretical results in Section 6 where applicable, performing exhaustive search in the few cases they are not.Note that computation of the system signature itself can be simplified by expressing the structure of the system as an undirected graph using the computeSystemSurvivalSignature() function in the same package, leaving only data and prior to be handled.These publicly available functions have been used in computing all the following examples for reproducibility.See Appendix B for further details of how to use this software.

Examples 7.2.1. Toy example
As a toy example, consider a 'bridge' type system layout with three types of components T1, T2 and T3, as depicted in Figure 4.The survival signature for this system is given in Table 1.All rows with T3 = 0 have been omitted; without T3, the system cannot function, thus Φ = 0.For component types T1 and T2, we consider a near-noninformative set of prior reliability functions.For components of type T3, we consider an informative set of prior reliability functions as given in Table 2.This set could result from eliciting prior functioning probabilities at times 0, 1, 2, 3, 4, 5 only, and filling up the rest.These prior assumptions, together with sets of posterior reliability functions resulting from three different scenarios for test data for  component type T3, are illustrated in Figures 5, 6 and 7; test data for components of type T1 and T2 are invariably taken as t 1 = (2.2,2.4, 2.6, 2.8) and t 2 = (3.2,3.4, 3.6, 3.8), respectively.
In Figure 5, test data for component type T3 is t 3 = (0.5, 1.5, 2.5, 3.5), and so in line with expectations.The posterior set of reliability functions for each component type and the whole system is considerably smaller compared to the prior set (due to the low prior strength intervals [n t ] = [1,4]) and so giving more precise reliability statements.We see that posterior lower and upper functioning probabilities drop at those times t when there is a failure time in the test data, or a drop in the prior functioning probability bounds.Note that the lower bound for the prior system reliability function is zero due to the prior lower bound of zero for T1; for the system to function, at least two components of type T1 must function.
In Figure 6, test data of component type T3 is t 3 = (0.6, 0.7, 0.8, 0.9), and so earlier than expected.Compared to Figure 5, posterior functioning intervals for T3 are wider between t = 1 and t = 3.5, reflecting additional imprecision due to prior-data conflict.For t > 1, it is clearly visible how y 3,t = 4 and n 3 = 4), while y (n) 3,t is one-fifth of y (0) 3,t (weights n 3,t = 1 and n 3 = 4).Note that the posterior system functioning probability is constant for t ∈ [1,2] because in that interval the prior functioning probability is constant and there are no observed failures.
In Figure 7, test data of component type T3 is t 3 = (4.1,4.2, 4.3, 4.4), and so observed failures are later than expected.Here we see that for t ∈ [2,4], posterior functioning bounds for T3 are even wider than prior functioning bounds.The width turns back to being half the prior width only after the four failures.The imprecision carries over to the system bounds, where we see wider bounds as compared to the other two scenarios especially between t = 2 and t = 4.In particular, also note that at the system level posterior bounds are a subset of prior bounds after t = 2.6, although prior-data conflict for the individual component type T3 extends well beyond t = 4.This demonstrates the power of this technique to identify prior-data conflict which is actually consequential at the system level, not just the component levelin other words, for mission times t > 2.6, we can diagnose that the prior-data conflict need not be of elevated concern for this system viewed as a whole.Nevertheless, the posterior system reliability bounds are wider than in the no-conflict case for t ∈ [1, 4.4], signalising the general need for caution in this scenario.
bounds for M are based on a Weibull cdf with shape 2.5 and scales 6 and 8 for the lower and upper bound, respectively.The prior bounds for P can be seen as the least committal bounds derived from an expert statement of y P,t ∈ [0.5, 0.65] for t = 5 only.For H, near-noninformative prior functioning probability bounds have been selected; with the upper bound for P being approximately one for t ≤ 5 as well, the prior upper system reliability bound for t ≤ 5 is close to one, too, since the system can function on H and one of P1 -P4 alone.Note that the posterior functioning probability interval for M is wide not only due to the limited number of observations, but also because n (0) M,t = 8 and the prior-data conflict reaction.
Posterior functioning probability bounds for the complete system are much more precise than the prior system bounds, reflecting the information gained from component test data.The posterior system bounds can be also seen to reflect location and precision of the component bounds; for example, the system bounds drop drastically between t = 2.5 and t = 3.5 mainly due to the drop of the bounds for P at that time.
It is also interesting to note that the prior-data conflict which is consequential at the system level occurs over roughly the same range in t as there is prior-data conflict for component type P. Indeed, this occurs despite there being prior-data conflict in both M and C over much larger ranges, giving valuable insight into which prior requires further expert attention -thus the technique avoids wasted time addressing prior-data conflict in components which may not be relevant when propagated to the uncertainty in the whole system.

Conclusions
In this paper we have contributed an imprecise Bayesian nonparametric approach to system reliability with multiple types of components.The approach allows modelling partial or imperfect prior knowledge on component failure distributions in a flexible way through bounds on the functioning probability for a grid of time points, and combines this information with test data in an imprecise Bayesian framework.Component-wise predictions on the number of functioning components are then combined to bounds for the system survival probability by means of the survival signature.New results on first-order stochastic dominance for the Beta-Binomial distribution enable closed-form solutions for these bounds in most cases and avoid exponential growth in the complexity of computing the estimate as the number of components grows.The widths of the resulting system reliability bounds reflect the amount of test data, the precision of prior knowledge, and crucially provide an easily used method to identify whether these two information sources are in conflict in a way which is of consequence to the whole system reliability estimate.
These methodological contributions can be immediately used in applications by reliability practitioners as we provide easy to use software tools.
An important next step is to extend the model to include right-censored observations which are common in the reliability setting.In particular, this allows to use component failure observations from a running system to calculate its remaining useful life.We see two potential approaches.First, to obtain lower and upper system reliability bounds one can assume that a component either fails immediately after censoring or continues to function during the entire time horizon.This minimal assumption will be simple to implement but will lead to high imprecision.Alternatively, one can assume exchangeability with other surviving components at the moment of censoring.This approach will be more complex to accomodate but will lead to less imprecision.Indeed, this assumption lies at the core of the Kaplan-Meier estimator [14], and has already been adopted by Coolen and Yan [11] in an imprecise probability context.
Upscaling the survival signature to large real-world systems and networks, consisting of thousands of components, is a major challenge.However, even for such systems the fact that one only needs to derive the survival signature once for a system is an advantage, and also the monotonicity of the survival signature for coherent systems is very useful if one can only derive it partially.
The survival signature and its use for uncertainty quantification for system reliability can be generalized quite straightforwardly, mainly due to the simplicity of this concept.For example, one may generalize the system structure function from a binary function to a probability, to reflect uncertainty about system functioning for known states of its components, with a further generalization to imprecise probabilities possible. Thus, However, unlike the case for the y parameter in Theorem 1, neither β n nor β n can be guaranteed to dominate for all possible values for the other parameters, so that necessary conditions for monotonicity (either increasing or decreasing) must be established.We require, After extensive routine algebra, this can be conveniently expressed as This limit is hardest to satisfy for l = m − 1 since n − n > 0 (note l = m since we are evaluating L(l + 1)/L(l), so m − 1 is the maximal value l can take).Thus, for monotonicity to hold for all l, we require Since n − n > 0 by definition, we have a monotonically increasing likelihood ratio only when y By a similar argument, the likelihood ratio is only monotonically decreasing when y(N + m − 1) − s < 0.
Thus, Proof of Lemma 3, p16.(A.1) and (A.2) are sufficient but not necessary conditions.Using theory in [15] we can sharpen these conditions to provide first-order stochastic dominance conditions for a larger range of parameter values.Proposition 2.1, p.399 of [15] proves that half-monotone likelihood ratio ordering -i.e.monotonicity of L(l + 1)/L(l) -together with left and right tail conditions on L(•) imply first order stochastic dominance.

. , m}
Thus we can conclude that L(•) is half monotone decreasing.
Tail conditions on L(•) It is not difficult to derive the same loose bounds as in Theorem 2 using the tail conditions.However, it is also easy to see that these are sufficient but not necessary.Sharpening these bounds in terms of the other parameter values involves seemingly intractable algebra, so we leave the tail condition as the alternative slightly more costly numerical check when the conditions of Theorem 2 are not satisfied.Evaluation of L(•) at two values is still orders of magnitude less costly than reevaluation of R sys (•) or R sys (•).

B. Software details
Functions which make it easy to use the methods of this paper have been added to the R package ReliabilityTheory [1].There are two functions of particular note: computeSystemSurvivalSignature and, implementing the result from Appendix A above, nonParBayesSystemInferencePriorSets.

B.1. Computing the survival signature
The function computeSystemSurvivalSignature allows easy computation of the survival signature if the system is expressed as an undirected graph with 'start' and 'terminal' nodes (which are not considered components for survival signature computation).The system is considered to work if there is a path from the start to the terminal node passing only through functioning components.
Graph representations of systems are most simply defined by using the graph.formulafunction.The 'start' node should be denoted s and the 'terminal' node should be denoted t and intermediate nodes (representing actual components) should be numbered and connected by edges denoted by -, where the numbering denotes physically distinct components.Component numbers can be repeated to include multiple links.For example, to build a simple three component series system: sys <-graph.formula(s-1 -2 -3 -t) and to build a three component parallel system: sys <-graph.formula(s-1 -t, s -2 -t, s -3 -t) There is an additional shorthand which indicates a link exists to a list of multiple components separated by the : operator, so that the parallel system can be also be expressed more compactly by: sys <-graph.formula(s-1:2:3 -t) Therefore, the simple bridge system of Figure 4 can be constructed with: sys <-graph.formula(s-1 -2 -3 -t, s -4 -5 -3 -t, 1:4 -6 -2:5) The grid of times, at.times, is specified as simply a vector of time points.The test.data argument is a list of component type names (as the tag) and corresponding lifetime data (as the value), for example a toy sized dataset for each component would be expressed as: test.data=list("T1"=c(0.19, 0.73, 1.87, 1.17), "T2"=c(0.22,0.27, 0.63, 1.80, 1.25, 1.95), "T3"=c(1.33,0.65, 1.59)) Finally, there are multiple options for specifying the prior parameter sets.Each of the nLower, nUpper, yLower and yUpper arguments can be specified as: • a single value for a homogeneous prior across time and components.e.g.nLower=2 =⇒ n (0) • a vector of values of length |T | (length(at.times)), for a time inhomogeneous prior which is identical across component types.
• a data frame of size 1 × K, where each column is named the same as in the survival.signatureand test.dataarguments, for a time homogeneous prior which varies across component types.
• a data frame of size |T | × K, where each column is named the same as in the survival.signatureand test.dataarguments, for a time inhomogeneous prior which varies across component types.
With these arguments supplied, nonParBayesSystemInferencePriorSets will then compute the posterior sets automatically in parallel across the cores of a multicore CPU and return a list with two objects, named lower and upper, containing respectively the lower and upper bound for the system reliability function R sys (t).
as in(7) are depicted in Figure 2.Such rectangular prior parameter sets Π (0) k,t have been shown to balance

Figure 1 :
Figure 1: Beta densities (top left) and cdfs (bottom left), with the corresponding Betabinomial predictive probability mass functions (top right) and cumulative mass functions (bottom right), for a prior with n (0) k,t = 8, y

Figure 3 :
Figure 3: Sets of Beta pdfs (left) and Beta-Binomial cmfs (right, for m k = 5) corresponding to the prior and posterior parameter sets in Figure 2. The sets are depicted as shaded areas, with the distributions corresponding to the four corners of the prior parameter set Π are obtained with n . In particular, y (0) k,t > 0 and y for all possible n values.The analogous result for n

Lemma 3 .
Let β n denote the Beta-Binomial distribution as in Theorem 2. Define L n,n (l) := p(l | y, n, m, s, N ) p(l | y, n, m, s, N ) should not increase.It is inadvisable to express certainty in the expected functioning probabilities with n bounds that vary substantially over the range of t.With (strongly) differing n bounds, monotonicity of the y bounds cannot be guaranteed.For example, if y

(
; low values for n as compared to the test sample size n k give low weight to the prior expected functioning probability intervals [y (0) k,t , y (0) k,t ], and the location of posterior intervals [y (n) k,t , y (n) k,t intervals, which means more cautious posterior inferences.The difference betweeen n determines the strength of the prior-data conflict sensitivity; as is clear from Figure2and (8), the wider the n interval, the wider [y (n) k,t , y (n) k,t ] in case of conflict.So it seems useful to choose n , posterior inferences are not overly sensitive to whether y

Figure 4 :
Figure 4: Reliability block diagram for a 'bridge' system with three component types.

Figure 5 :Figure 6 :Figure 7 :
Figure 5: Prior and posterior sets of reliability functions for the 'bridge' system and its three component types, with failure times as expected for component type Test data failure times are denoted with tick marks near the time axis.

Figure 8 :
Figure 8: Prior and posterior sets of reliability functions for a simplified automotive brake system with layout as depicted in the lower right panel.

Table 2 :
Lower and upper prior functioning probability bounds for component type T3 in the 'bridge' system example.