Example of direct sampling

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 np

def 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 Nsamples to 5000 (smaller values are also acceptable):

N = 4
Nsamples = 5000

S = np.zeros((N, Nsamples))
Fsamples = {}

for t in range(Nsamples):
    x1 = X1_sample()
    x2 = X2_sample()
    x3 = X3_sample(x1, x2)
    x4 = X4_sample(x3)
    
    sample = (x1, x2, x3, x4)
    
    if sample in Fsamples:
        Fsamples[sample] += 1
    else:
        Fsamples[sample] = 1

When the sampling is complete, it's possible to extract the full joint probability:

samples = np.array(list(Fsamples.keys()), dtype=np.bool_)
probabilities = np.array(list(Fsamples.values()), dtype=np.float64) / Nsamples

for i in range(len(samples)):
    print('P{} = {}'.format(samples[i], probabilities[i]))

P[ True False  True  True] = 0.0286
P[ True  True False  True] = 0.024
P[ True  True  True False] = 0.06
P[False False False False] = 0.0708
P[ True False  True False] = 0.0166
P[False  True  True  True] = 0.1006
P[False False  True  True] = 0.054
...

We can also query the model. For example, we could be interested in P(X₄=True). We can do this by looking for all the elements where X₄=True, and summing up the relative probabilities:

p4t = np.argwhere(samples[:, 3]==True)
print(np.sum(probabilities[p4t]))

0.5622

This value is coherent with the definition of X₄, which is always p >= 0.5. The reader can try to change the values and repeat the simulation.