09-Probabilities#

When quantifying the risk from a specific failure mode, it can be important to understand several event likelihoods or probabilities such as:

  • probability of the event occurring in a given time period (or given number of cycles)

  • probability of a given energy level being exceeded (such as an energy level that would cause a failure) during that event

  • probability that given the energy level is exceeded, that the failure mode or failure outcome will occur

  • probability that an individual or group of individuals will be exposed to the event

  • etc.

Perhaps a simple example is the probability of breaking a hip after age 65. You’d include the probability of falling down as one impetus and perhaps a bike crash as another and then for each of those you’d have the probability of a break given the event occurred. We are part of such calculations as we pay insurance premiums. Actuarial science is the discipline that deals with the quantification of risk using probabilities.

In the discussion here, we will focus on quantitatively estimating event probabilities, or the probability an event will occur in a given time period.

Event Probability from Poisson Distribution#

There are many statistical distributions including the normal distribution, chi squared, exponential, etc. that can be used to model event probabilities. One of the most common distributions used in reliability engineering is the Poisson distribution. The Poisson distribution is often used to model the number of events that occur in a fixed interval of time or space. The probability of observing \(k\) events in a given time period is given by the Poisson probability mass function:

\[ P(k) = \frac{\lambda^k e^{-\lambda}}{k!} \]

where \(\lambda\) is the rate of occurrences.

In the treatment here, we will create a \(\lambda'\) prime that is the rate per unit time and we will also only care that k or the number of events is greater than or equal to 1. This is because we are interested in the probability of at least one event occurring in a given time period. As such, the following is true of the probability estimate:

\[\begin{split} P = 1 - e^{-\lambda' t} \\ P = 1 - e^{-\mu t} \end{split}\]

where the second equation is how it is given in Crowl and Louvar’s Chemical Process Safety 4th edition.

Multiple Event Probability Estimates#

For clarity, the example shown above has a mean time between failures (1/\(\mu\)) of 6000 seconds, 100 minutes, or 1.67 hours. The probability of a failure occurring in 1 second is 1.67e-4 which is very close to \(\mu\) = 1.67e-4. However for longer times like 1 month, the probability of failure is very close to 1 but \(\mu\) is 432 failures expected per month. Don’t confuse \(\mu\) or the failure rate with P, the probability of failure in a given time period.

import numpy as np
P = 1 - np.exp(-1/6000*1)
print(f'The probability of failure in 1 second is {P:.8f}, which is approximately equal to mu = {1/6000:.8f}')

The probability of failure in 1 second is 0.00016665, which is approximately equal to mu = 0.00016667

mu = 1/100*60*24*30 # rate per month
P = 1 - np.exp(-mu*1)
print(f'The probability of failure in 1 month is {P:.3f}, which is not equal to mu = {mu:.0f}')

The probability of failure in 1 month is 1.000, which is not equal to mu = 432

Important

Please be very familiar with Examples 12-1 and 12-2 in Crowl and Louvar’s Chemical Process Safety 4th edition. Those examples are great in helping you work with reliability, probability, failure rates and mean time between failures for and and or systems.

Monte-Carlo#

Monte-Carlo methods can be used when there are more complex scenarios ini terms of interactions between events, dependencies, etc. In this case, the probability of an event occurring is estimated by simulating the scenario many times and counting the number of times the event occurs. The probability is then estimated as the number of times the event occurs divided by the total number of simulations.

Example Monte-Carlo Simulation

#import needed packages
import numpy as np
import matplotlib.pyplot as plt

# Number of random numbers to generate
N = int(1e4)

MTBF_A = 1 # mean time between failures, unit time for Event A
muA = 1/MTBF_A # average events per unit time for A
1-np.exp(-muA) # probability of at least one event in unit time (poisson process)

0.6321205588285577

muA*np.exp(-muA) # probability of exactly one event in unit time (poisson process)

0.36787944117144233

unif = np.random.uniform(0, 1, N) # generate N random numbers between 0 and 1
# 1 if event, 0 if no event, where the if statement is the Monte Carlo approach
eventA = [1 if x < 1-np.exp(-muA)  else 0 for x in unif]

Each entry in that list is an “event” over a duration of unit time and it can either be a failure or a success.

pA = sum(eventA)/len(eventA) # fraction of time with at least one event, probability of event in time unit
pA

0.6253

#determine the average rate of events per unit time
-np.log(1-pA)

0.9816295731824953

# now let's do the same thing for a Poisson process with a second 'event'
unif = np.random.uniform(0,1,N) #generate new random numbers
MTBF_B = 2 # mean time between failures
muB = 1/MTBF_B # average events per unit time
eventB = [ 1 if x < 1-np.exp(-muB) else 0 for x in unif ]

# now consider eventA and eventB as two independent Poisson processes in an OR gate
eventAorB = [min(1,sum([a,b])) for a,b in zip(eventA, eventB)]

pAorB = sum(eventAorB)/len(eventAorB) # fraction of time with at least one event
pAorB

0.7728

#combination of the two events for the average rate of events per unit time (mu)
-np.log(1-pAorB)

1.481924592135141

# now consider eventA and eventB as two independent Poisson processes in an AND gate
eventAandB = [a*b for a,b in zip(eventA, eventB)]
pAandB = sum(eventAandB)/len(eventAandB) # fraction of time with at least one event
pAandB

0.2489

#Calculation of the probability of A and B occurring at the same time (product of the two probabilities)
pab = (1-np.exp(-muB))*(1-np.exp(-muA))
pab

0.2487200592643541

#combination of the two events for the average rate of events per unit time (mu) from simulation
-np.log(1-pAandB)

0.286216480290171

#combination of the two events for the average rate of events per unit time (mu) from product of probabilities
-np.log(1-pab)

0.28597693937029134