Bayesian Statistics: An Introduction for the Practicing Reliability Engineer

A common way to measure the reliability of a system is to determine the probability that it will pass, or survive, a stress test. This typically requires several system tests, but it may be prohibitively expensive to conduct many tests of a sophisticated system, such as an air-to-air missile¹ or medical device.²

Bayesian statistical methods can help in such a situation, since they enable inclusion of other types of data (such as computer simulation experiments or subject-matter-expert opinions). Bayesian methods may also be necessary because many modern systems do not fail during testing. With no failures, it is difficult for classical statistics to accurately quantify the probability of failure.

This article begins by briefly reviewing Bayesian methods. It then discusses how these methods can be used to learn more about the probability of a system surviving a test, and it concludes by describing a straightforward algorithm for calculating a total system’s reliability once it has been tested. This algorithm is simple, produces an exact answer, and is not mentioned in the Bayesian reliability literature. The efficiency of this algorithm is illustrated with an example.

This section introduces the basics of Bayesian principles and Bayesian statistical methodology. The most effective way to introduce this concept is to contrast it with the principles and methodology of classical statistics. The biggest difference between Bayesian and classical statistics is in how probability is defined. In classical statistics, probability is the long-run frequency of an event. So for a fixed (and unknown) parameter such as a population mean, μ,

\[
\mathbb{P}(3.66 \leq \mu \leq 4.11) =
\begin{cases}
1 & \text{if true} \\
0 & \text{if not}
\end{cases}.
\tag{1}
\]

In words, Eq. 1 states that the fixed parameter μ is either in the stated interval or it is not.

Bayesian statisticians think about probability in a different way. In Bayesian statistics, probability is the belief that a statement is true. So if one believes (based on their experience and/or the data that they have seen) that μ is within the stated interval with 95% probability, it would be fair to say that

\[
\mathbb{P}(3.66 \leq \mu \leq 4.11) = .95\text{.} \tag{2}
\]

The objective and point of Bayesian statistics is to calculate probabilities like the one in Eq. 2 and to ensure that this calculation is scientifically respected.

To calculate such a probability, a Bayesian statistician begins with a prior distribution. Assuming that the unknown parameter of interest is θ, this prior distribution is typically denoted as π(θ). The prior distribution indicates where the user believes the parameter θ to be before data are observed or collected. Assume, for example, that we purchased a coin at a magic shop. Upon the purchase, the shop owner tells us that the coin will more often turn up heads than tails. In this case, we will let θ = ℙ(H) and define the prior distribution π(θ) for all values of θ between 0 and 1. This prior will also be more heavily weighted toward values of 1 to indicate that, a priori, the coin is expected to turn up heads more frequently than tails.

Once the prior is formulated, data are collected. The distribution of the data conditioned on a value of θ is written as p(x|θ), i.e.,

\[
(X_1, X_2, \ldots, X_n) \sim p(x_1, x_2, x_3, \ldots, x_n \mid \theta) = p(\mathbf{x} \mid \theta)\text{,}
\]

where x = (x₁, x₂, . . . , x_n). The function p(x|θ) is also referred to as the likelihood of θ.

With the prior and the likelihood, the posterior distribution, typically denoted as π(θ|x), can be calculated. The posterior is calculated using Bayes’s rule.³ This calculation is shown below:

\[
\pi(\theta \mid \mathbf{x}) = \frac{p(\mathbf{x}, \theta)}{p(\mathbf{x})}
= \frac{p(\mathbf{x} \mid \theta) \pi(\theta)}{\int_{\Theta} p(\mathbf{x} \mid \theta) \pi(\theta)\, d\theta}
\propto p(\mathbf{x} \mid \theta) \pi(\theta)\text{,} \tag{3}
\]

where Θ is the set of all possible values of θ. The formula given in Eq. 3 makes sense: the posterior is proportional to the prior distribution of θ (where we thought θ was before collecting data) times the likelihood (where the data suggests θ to be).

This article focuses on how to use Bayesian methods to learn more about the probability that a system survives a test. The subsection titled The Prior and Posterior of One Subsystem specifically discusses how Bayesian methods are used to learn about the survival probability of one system, and the subsection following that one, The Prior and Posterior of the Entire System, discusses the methods necessary to learn about the survival probability of one system composed of multiple subsystems. The algorithm discussed in the latter subsection is elementary and uses no approximations when calculating its answer.

The Prior and Posterior of One Subsystem

Let us first assume that we are interested in learning about the probability that a system passes an endurance test of some sort. We will denote this probability as θ, and we will conduct n trials/tests on this system and record X, the number of times (out of the n trials) that the system passes a test.

To do a Bayesian analysis on θ, we begin by specifying a prior distribution for it. The beta distribution is often used as the prior for the probability of success in a sequence of success/failure trials.^4,5 The beta distribution is specified by two parameters and is especially convenient in cases such as this since it is a conjugate prior (i.e., it produces a posterior distribution of the same form). The beta prior takes the form

\[
\pi(\theta) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \, \theta^{\alpha - 1}(1 - \theta)^{\beta - 1} \quad 0 \leq \theta \leq 1\text{,}
\]

where Γ( · ) is the gamma function. This prior has mean

\[
\text{Prior Mean} = \frac{\alpha}{\alpha + \beta}
\]

and variance

\[
\text{Prior Var} = \frac{\alpha \beta}{(\alpha + \beta)^2 (\alpha + \beta + 1)}\text{.}
\]

The values of α and β (α, β > 0) are selected to reflect the user’s prior belief. This prior belief is often informed in a variety of ways, such as expert opinion, computer simulation, or prior experiments. A user who believes that θ is small (<0.5) would set α < β (making the prior mean < 0.5). A user who believes that θ is large (>0.5) would set α > β. The confidence in these prior beliefs is, of course, reflected in the variance of the prior. If a user wanted to set the prior mean of θ to 0.4, they could set α = 2 and β = 3, making the prior variance 0.04. If a user wanted to elevate the confidence in this statement (that the prior mean of θ is 0.4), they could adjust the values of α and β to α = 20 and β = 30, making the prior variance 0.004. And if a user knew absolutely nothing about θ, they would set α = β = 1, in which case the prior for θ is uniform over the interval (0, 1).

There is no consistent selection of α and β in reliability studies. Leoni et al.⁶ set α = 3 and β = 1 in one of their reliability studies, Burke and Harman⁷ set α = 7.2 and β = 0.8 in one of their reliability studies, and Martz et al.⁸ set α = 27.3 and β = 0.5. When priors are informed by simulation results or previous experiments, analysts sometimes set

\[
\alpha = n_{\text{pr}} \cdot \hat{\theta}^{\text{pr}} + 1
\quad \text{and} \quad
\beta = n_{\text{pr}} \cdot (1 - \hat{\theta}^{\text{pr}}) + 1\text{,}
\]

where \(\hat{\theta}^{\text{pr}}\) is a prior estimate of θ, and n_pr is some positive number that represents the confidence the analyst has in the simulation or experiment informing the prior;⁹ think of n_pr as the effective sample size that informs the prior. The greater this effective sample size that informs the prior, the more peaked the prior distribution is near the prior estimate of θ. If no confidence exists in the simulation informing the prior, then n_pr = 0 and the prior would be flat.

Figures 1 and 2 illustrate what these prior distributions look like. The prior in Figure 1 puts large probability on low values of θ and does so by setting the value of α to be significantly less than the value of β. The prior in Figure 2 puts large probability on high values of θ and does so by setting the value of β to be smaller than the value of α. Also observe that the prior is much more peaked for small values of θ in Figure 1 than it is for large values of θ in Figure 2. This is a consequence of the difference in the values between α and β. The difference is larger for the prior in Figure 1 than it is for the prior in Figure 2.

Let us now assume that we observe x successes of the system out of n tests conducted. In this case, the likelihood is binomial,

\[
p(x \mid \theta) = \binom{n}{x} \theta^x (1 - \theta)^{n - x}\text{,}
\]

making the posterior distribution of θ

\[
\pi(\theta \mid x) = \frac{p(x \mid \theta) \, \pi(\theta)}{\int_{\Theta} p(x \mid \theta) \, \pi(\theta) \, d\theta}
= \frac{\binom{n}{x} \theta^x (1 - \theta)^{n - x} \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \theta^{\alpha - 1}(1 - \theta)^{\beta - 1}}{\int_0^1 \left[ \binom{n}{x} \theta^x (1 - \theta)^{n - x} \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \theta^{\alpha - 1}(1 - \theta)^{\beta - 1} \right] d\theta}\text{.}
\]

There is a simple way to calculate the posterior distribution in Eq. 4 without having to evaluate the integral in the denominator. First observe that the expression in the denominator is not a function of θ; it is a normalizing constant independent of θ, and for this reason we can write

\[
\pi(\theta \mid x) \propto p(x \mid \theta) \, \pi(\theta)\text{.}
\]

Eliminating all multiplicative constants in p(x|θ) π(θ) that do not depend on θ, we get that

\[
\pi(\theta \mid x) = c \cdot \theta^{x + \alpha - 1} (1 - \theta)^{n - x + \beta - 1}\text{,}
\]

where c is some constant such that

\[
\int_0^1 c \cdot \theta^{x + \alpha - 1} (1 - \theta)^{n - x + \beta - 1} \, d\theta = 1\text{.}
\]

The posterior π(θ|x) takes the form of a beta distribution, making

\[
c = \frac{\Gamma(\alpha + \beta + n)}{\Gamma(\alpha + x) \, \Gamma(\beta + n - x)}\text{.}
\]

The posterior of θ is thus a beta distribution with parameters α^pst and β^pst where

\[
\begin{aligned}
\alpha^{\text{pst}} &= \alpha + x, \quad \text{and} \\
\beta^{\text{pst}} &= \beta + n - x\text{.}
\end{aligned}
\]

The plots in Figures 3 and 4 show the posteriors corresponding to the priors shown in Figures 1 and 2, respectively. In Figure 3, one success was observed after 10 trials, emphasizing even more that the value of θ is small. Observe how the posterior in this case is more peaked at small values of θ than the prior was. In Figure 4, two successes were observed in 10 trials, indicating that the probability of success was much smaller than the prior anticipated. Observe how, in this case, the peak of the posterior has significantly shifted to smaller values of θ.

The Prior and Posterior of the Entire System

Let us now put this problem in the context of one large system that is composed of several subsystems. If all the subsystems have to work for the entire system to work, how do the posterior distributions of the subsystem reliabilities inform the distribution of the total system’s reliability? And how would testing the entire system (as a whole) affect the posterior of the total system reliability? The next two subsections address these questions.

Subsystem Test Sizing

This subsection focuses on how the distributions of the subsystem survival probabilities affect the distribution of the total system’s survival probability. If there are S subsystems, and the entire system fails if any one of the subsystems fail (i.e., the subsystems work in series), then the success probability of the total system, θ_{Tot Sys}, is calculated as

\[
\begin{aligned}
\theta_{\text{Tot Sys}} 
&= \mathbb{P}(\text{Success of Subsys 1}) \times \mathbb{P}(\text{Success of Subsys 2}) \times \cdots \\
&\quad \times \mathbb{P}(\text{Success of Subsys } S) \\
&= \prod_{j=1}^{S} \theta_j,
\end{aligned}
\]

where θ_j is the success probability of the jth subsystem.

If S = 5 and the block diagram of the subsystems is as shown in Figure 5, then the system fails if subsystem 1, 4, 5, or both 2 and 3 fail. In this case, the success probability of the entire system would be calculated as

\begin{align*}
\theta_{\text{Tot Sys}} 
&= \mathbb{P}(\text{Success of Subsys 1}) 
\times \mathbb{P}(\text{Success of Subsys 2 or 3}) 
\times \mathbb{P}(\text{Success of Subsys 4}) 
\times \mathbb{P}(\text{Success of Subsys 5}) \\
&= \theta_1 \left( \theta_2 + \theta_3 - \theta_2 \theta_3 \right) \theta_4 \theta_5\text{.}
\end{align*}

The value of θ_{Tot Sys} is thus the product and/or sum of beta random variables. The resulting distribution of a random variable such as θTot Sys has been derived in a number of publications, but this distribution is very complicated and thus difficult to work with analytically.^{1,10,11,12,13} The distribution of θ_{Tot Sys} is easy to work with and understand, however, using Monte Carlo methods. Since the posterior distribution of all the components of the system take the form of a beta distribution with known parameters, assuming independence of the subsystems, we can easily simulate n_Sim values from the prior distribution of θ_{Tot Sys}. This requires simulating n_Sim S-tuples of (θ₁, θ₂, θ₃,..., θ_S). With each simulated S-tuple, we can calculate a value of θ_{Tot Sys}. The algorithm for generating n_Sim values of θ_{Tot Sys} for a system in series is given in Procedure 1; that for a system as shown in Figure 5 is given in Procedure 2.

Figures 6–12 illustrate how the posterior distributions of subsystem reliability affect the distribution of θ_{Tot Sys}. In the simulations performed, we assumed that the entire system was composed of three subsystems (S = 3) and that these subsystems worked in series. The priors of the three subsystems are shown in black in Figures 6–8. The subsystems are then tested with n₁ = 2, n₂ = 5, and n₃ = 4, where n_j is the number of times the jth subsystem is tested. The resulting posteriors are shown in red in Figures 6–8, and 10,000 draws from the resulting distribution of θ_{Tot Sys} are shown in Figure 9. The subsystems were also tested at n₁ = 11, n₂ = 14, and n₃ = 12, and the corresponding posteriors are shown in blue in Figures 6–8. Observe that these posteriors are more peaked (more informed) than the others since the sample sizes are larger. The distribution of θ_{Tot Sys} corresponding to these larger sample sizes is shown in Figure 10. Observe how the variance of this posterior is smaller than that shown in Figure 9; this is because the subsystem sample sizes are larger.

Total System Test Sizing

Let us now investigate how testing the entire system (and not just its individual components) affects the posterior distribution of θ_{Tot Sys}. Updating the prior distribution of θ_{Tot Sys} given test results on the total system is more challenging than updating the subcomponent values of θ because, in this case, the original distribution of θ_{Tot Sys} is not a beta distribution. Recall that the prior distribution of θ_{Tot Sys} was analytically challenging to work with and, as a result, was obtained using Monte Carlo methods. It is not uncommon for practitioners to approximate this prior with another (perhaps beta) distribution^{1,11,14,15,16,17} to make the posterior analysis simpler and more convenient. Others redefine the priors of the independent components entirely just so the prior of the total system’s reliability is analytically tractable. Zoh et al.,¹⁸ for example, set the priors of the components to negative log-gamma distributions, yet they still had to execute complicated Markov chain Monte Carlo methods to do posterior analysis of the total system.

These workarounds and approximations to make posterior analysis easier are not necessary. This article describes a simple and quick Monte Carlo method that updates the prior distribution of θ_{Tot Sys} after testing the entire system. This method is a simple application of the algorithm proposed by Rubin.¹⁹ Rubin observed that a sample from the posterior distribution of a parameter can be obtained by first generating values from its prior and then generating data conditioned on these sampled values. Those values of the parameter for which the generated data match the observed data follow the posterior distribution.

To apply this algorithm in our case, we begin by writing the posterior for θ_{Tot Sys} as

\[
\pi(\theta_{\text{Tot Sys}} \mid x_{\text{TS}}) \propto p(x_{\text{TS}} \mid \theta_{\text{Tot Sys}}) \, \pi(\theta_{\text{Tot Sys}})\text{,}
\]

where

\[
p(x_{\text{TS}} \mid \theta_{\text{Tot Sys}}) = 
\binom{n_{\text{TS}}}{x_{\text{TS}}} \, 
\theta_{\text{Tot Sys}}^{x_{\text{TS}}} 
(1 - \theta_{\text{Tot Sys}})^{n_{\text{TS}} - x_{\text{TS}}}\text{,} \tag{6}
\]

and n_TS and x_TS are the number of tests (and successes) of the total system. Given that \(x_{TS}^\ast\) successes have been observed from n_TS trials of the total system, we sample from the posterior \(\pi(\theta_{\text{Tot Sys}} \mid x_{TS} = x_{TS}^\ast)\) by first simulating from the prior of θ_{Tot Sys}, π(θ_{Tot Sys}). (An algorithm similar to the ones shown in Procedure 1 or 2 could be used to sample from the prior, π(θ_{Tot Sys}).) We then condition on these sampled values of θ_{Tot Sys} to generate candidate values of x_TS from the likelihood shown in Eq. 6. The simulated values of θ_{Tot Sys} for which the likelihood generates \(x_{\text{TS}} = x_{\text{TS}}^\ast\) are then considered to be an exact sample from the posterior. The details of this algorithm (assuming the subsystems work in series; minor changes to the first “for” loop of the algorithm would be necessary if the subsystems did not work in series) are given in Procedure 3.

The plots in Figures 11 and 12 show how the distribution of θ_{Tot Sys} changes when tests on the entire system are executed. The plot in Figure 11 shows how the distribution of θ_{Tot Sys} changes from the distribution in Figure 9 when four successes are observed out of four tests on the entire system. Observe that with this extra evidence of success, the distribution of θ_{Tot Sys} shifts to the right. The same story is told in Figure 12. It shows how the distribution of θ_{Tot Sys} changes from the distribution in Figure 10 when five successes are observed out of seven tests on the entire system. This distribution moves to the right and is also more peaked.

It is critical to understand the distinction between the two examples discussed above. In the first case (with the resulting distribution of θ_{Tot Sys} shown in Figure 11), the total system was tested four times and four successes were observed. Because the system works in series, a successful test of the entire system implies a successful test of each component. The posterior distribution of θ_{Tot Sys} can thus easily be calculated by simply updating the posteriors of the system’s three components and then applying Procedure 1. This is not true for the second example (with the resulting distribution of θ_{Tot Sys} shown in Figure 12). Recall that in the second example, the entire system was tested seven times, but only five successes were observed. Since it is not clear which component(s) failed (causing the failure of the entire system), the posteriors of the components cannot be updated, and Procedure 1 cannot be applied. In this case, Procedure 3 is necessary in calculating the posterior of θ_{Tot Sys}.

This example illustrates the efficiency of Rubin’s algorithm when calculating total system reliability. Consider two types of systems/fault trees, each with m different types of components. The first fault tree works in series, and in the second, every other component operates in series. Figures 13 and 14 illustrate these fault trees.

The efficiency of Rubin’s algorithm is demonstrated with a Monte Carlo study. We initially place Beta(α,β) priors on all the components with α = 999 and β = 0.5. (Such priors assume a prior reliability mean of 0.9995, which is not uncommon for highly reliable systems.) For each system/fault tree, and for a specific value of m, n_j (the number of trials for subsystem j) and x_j (the number of successful trials for subsystem j) can then be simulated. The number of trials is simulated from a Poisson(λ = 5) distribution, and the number of successful trials is simulated from a binomial distribution with parameters n_j and p_j = α/(α + β) = 999/(999 + 0.5). With these simulated test sizes, and assuming that the entire system is tested twice (with both tests being a success), the time it took to generate 10,000 draws from the posterior of the total system’s reliability can be calculated. Table 1 shows the results of doing this 100 times and calculating the average time it took to generate these posterior values for m = 10, 20, 30, 40, and 50.

It is clear from Table 1 that this algorithm efficiently calculates the posterior of the total system’s reliability.

This article reviews some of the methodologies related to Bayesian reliability. It initially focuses on success/failure data of systems and their subsystems. The article specifically addresses how the number of subsystem (or total system) tests affects the reliability of the entire system. It also presents a simple and efficient Monte Carlo method that can be employed to update the prior of a total system’s reliability when only data from the total system are available.