We can now implement this algorithm in Python. Let's start by defining the sample methods using the NumPy function np.random.binomial(1, p), which draws a sample from a Bernoulli distribution with probability p:
import numpy as npdef X1_sample(p=0.35): return np.random.binomial(1, p)def X2_sample(p=0.65): return np.random.binomial(1, p)def X3_sample(x1, x2, p1=0.75, p2=0.4): if x1 == 1 and x2 == 1: return np.random.binomial(1, p1) else: return np.random.binomial(1, p2) def X4_sample(x3, p1=0.65, p2=0.5): if x3 == 1: return np.random.binomial(1, p1) else: return np.random.binomial(1, p2)
At this point, we can implement the main cycle. As the variables are Boolean, the total number of probabilities is 16, so we set ...