Probability of system failure, as it requires two adjacent component to fail

Question

I am stuck with a problem in probability. It says that the system consists of 16 components distributed in corners and intersections as the following:

  s---s---s---s 
  |   |   |   | 
  s---s---s---s 
  |   |   |   | 
  s---s---s---s 
  |   |   |   | 
  s---s---s---s 

Failure of two adjacent components will lead to system failure. The probability of failure for each component is 0.0005.

I tried different ways but doesn't work, and I was thinking to do it as a simulation but I didn't have any idea

Answer 1

It's easy to see that the number of pairs of two adjacent components is 8*3=24 (each large edge has 3 pairs and there are 8 edges). Because the probability of failure is so small, we can ignore the case where three or more fail simultaneously. For failure, two have to fail simultaneously and the two have to be adjacent. The probability becomes -

$$p = {16 \choose 2}(0.0005)^2 (1-0.0005)^{14} \frac{24}{{16 \choose 2}}$$

This is a lower bound since we didn't consider 3 failing, 4 failing and so on. But those probabilities are extremely low.

If you want to include 3 failing as well your estimate becomes (24*14 because two adjacent need to fail and choose any of the remaining 14 for the third component that failed):

$$p = {16 \choose 2}(0.0005)^2 (1-0.0005)^{14} \frac{24}{{16 \choose 2}} + {16 \choose 3}(0.0005)^3 (1-0.0005)^{13} \frac{24 \times 14}{{16 \choose 3}}$$

4 onwards becomes complicated with negligible addition to the answer.

Answer 2

Another answer that uses simulation. I will provide the python code. Note that if you number the components from top-left to bottom right, two components are adjacent if the difference in the numbers they got is either 1 or 4. Using this fact:

import numpy as np
sys_fails = 0
for n in range(10000):
  failed = []
  for comp_no in range(16):
    if np.random.uniform() < .0005:
        failed.append(comp_no)
  sys_fails += did_sys_fail(failed)
prob_system_failed = sys_fails/10000


def did_sys_fail(failed):
  if len(failed)>1:
     for i1 in range(len(failed)):
       for j1 in range(i1+1,len(failed)):
          if abs(failed[i1]-failed[j1])==1 or abs(failed[i1]-failed[j1])==4:
            return 1
  return 0

Because the probability of failure is so low, you will have high variance in the result. One way to overcome this is to use importance sampling where you increase the failure probability (to say, 50%), but then adjust each simulation.