Bayesian Networks in Python

I am implementing two bayesian networks in this tutorial, one model for the Monty Hall problem and one model for an alarm problem. A bayesian network (BN) is a knowledge base with probabilistic information, it can be used for decision making in uncertain environments.

Bayesian networks is a systematic representation of conditional independence relationships, these networks can be used to capture uncertain knowledge in an natural way. Bayesian networks applies probability theory to worlds with objects and relationships. Conditional independence relationships among variables reduces the number of probabilities that needs to be specified in order to represent a full joint distribution. A full joint distribution can answer any question but it will become very large as the number of variables increases.

A bayesian network is created as a directed acyclic graph (DAG) with nodes, edges and conditional probabilities. Conditional probabilities is calculated with Bayes theorem, calculations is based on joint probability distributions that we create when we build the network. Nodes represents variables (Alarm, Burglary) and edges represents the links (connections) between nodes.

We can ask questions to a bayesian network and get answers with estimated probabilities for events. It is possible to use different methods for inference, some is exact and slow while others is approximate and fast. The library that I use have the following inference algorithms: Causal Inference, Variable Elimination, Belief Propagation, MPLP and Dynamic Bayesian Network Inference.

Libraries

I am using pgmpy, networkx and pylab in this tutorial. I had some problems when installing pgmpy as it requires torch, the installation of torch failed. I installed torch to Python 3.7 with: pip install https://download.pytorch.org/whl/cpu/torch-1.1.0-cp37-cp37m-win_amd64.whl.

Monty Hall Problem

This problem is about a contest in which a contestant can select 1 of 3 doors, it is a price behind one of the doors. The host of the show (Monty) opens a empty door after the contestant has selected a door and asks the contestant if he want to switch to the other door. The question is if it is best to stick with the selected door or switch to the other door. It is best to switch to the other door because it is a higher probability that the price is behind that door.

# Import libraries
import pgmpy.models
import pgmpy.inference
import networkx as nx
import pylab as plt

# Create a bayesian network
model = pgmpy.models.BayesianModel([('Guest', 'Monty'), 
                                    ('Price', 'Monty')])

# Define conditional probability distributions (CPD)

# Probability of guest selecting door 0, 1 and 2
cpd_guest = pgmpy.factors.discrete.TabularCPD('Guest', 3, [[0.33, 0.33, 0.33]])

# Probability that the price is behind door 0, 1 and 2
cpd_price = pgmpy.factors.discrete.TabularCPD('Price', 3, [[0.33, 0.33, 0.33]])

# Probability that Monty selects a door (0, 1, 2), when we know which door the guest has selected and we know were the price is
cpd_monty = pgmpy.factors.discrete.TabularCPD('Monty', 3, [[0, 0, 0, 0, 0.5, 1, 0, 1, 0.5], 
                                                           [0.5, 0, 1, 0, 0, 0, 1, 0, 0.5], 
                                                           [0.5, 1, 0, 1, 0.5, 0, 0, 0, 0]], 
                                              evidence=['Guest', 'Price'], 
                                              evidence_card=[3, 3])

# Add CPDs to the network structure
model.add_cpds(cpd_guest, cpd_price, cpd_monty)

# Check if the model is valid, throw an exception otherwise
model.check_model()

# Print probability distributions
print('Probability distribution, P(Guest)')
print(cpd_guest)
print()
print('Probability distribution, P(Price)')
print(cpd_price)
print()
print('Joint probability distribution, P(Monty | Guest, Price)')
print(cpd_monty)
print()

# Plot the model
nx.draw(model, with_labels=True)
plt.savefig('C:\\DATA\\Python-data\\bayesian-networks\\monty-hall.png')
plt.close()

# Perform variable elimination for inference
# Variable elimination (VE) is a an exact inference algorithm in bayesian networks
infer = pgmpy.inference.VariableElimination(model)

# Calculate probabilites for doors including price, the guest has selected door 0 and Monty has selected door 2
posterior_probability = infer.query(['Price'], evidence={'Guest': 0, 'Monty': 2})

# Print posterior probability
print('Posterior probability, Guest(0) and Monty(2)')
print(posterior_probability)
print()

Probability distribution, P(Guest)
+----------+------+
| Guest(0) | 0.33 |
+----------+------+
| Guest(1) | 0.33 |
+----------+------+
| Guest(2) | 0.33 |
+----------+------+

Probability distribution, P(Price)
+----------+------+
| Price(0) | 0.33 |
+----------+------+
| Price(1) | 0.33 |
+----------+------+
| Price(2) | 0.33 |
+----------+------+

Joint probability distribution, P(Monty | Guest, Price)
+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+
| Guest    | Guest(0) | Guest(0) | Guest(0) | Guest(1) | Guest(1) | Guest(1) | Guest(2) | Guest(2) | Guest(2) |
+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+
| Price    | Price(0) | Price(1) | Price(2) | Price(0) | Price(1) | Price(2) | Price(0) | Price(1) | Price(2) |
+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+
| Monty(0) | 0.0      | 0.0      | 0.0      | 0.0      | 0.5      | 1.0      | 0.0      | 1.0      | 0.5      |
+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+
| Monty(1) | 0.5      | 0.0      | 1.0      | 0.0      | 0.0      | 0.0      | 1.0      | 0.0      | 0.5      |
+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+
| Monty(2) | 0.5      | 1.0      | 0.0      | 1.0      | 0.5      | 0.0      | 0.0      | 0.0      | 0.0      |
+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+

Finding Elimination Order: : : 0it [00:00, ?it/s]
0it [00:00, ?it/s]
Posterior probability, Guest(0) and Monty(2)
+----------+--------------+
| Price    |   phi(Price) |
+==========+==============+
| Price(0) |       0.3333 |
+----------+--------------+
| Price(1) |       0.6667 |
+----------+--------------+
| Price(2) |       0.0000 |
+----------+--------------+

Alarm Problem

A person has installed a new alarm system that can be triggered by a burglary or an earthquake. This person also have two neighbors (John and Mary) that are asked to make a call if they hear the alarm. This problem is modeled in a bayesian network with probabilities attached to each edge. Alarm has burglary and earthquake as parents, JohnCalls has Alarm as parent and MaryCalls has Alarm as parent. We can ask the network: what is the probability for a burglary if both John and Mary calls.

# Import libraries
import pgmpy.models
import pgmpy.inference
import networkx as nx
import pylab as plt

# Create a bayesian network 
model = pgmpy.models.BayesianModel([('Burglary', 'Alarm'), 
                                    ('Earthquake', 'Alarm'),
                                    ('Alarm', 'JohnCalls'), 
                                    ('Alarm', 'MaryCalls')])

# Define conditional probability distributions (CPD)

# Probability of burglary (True, False)
cpd_burglary = pgmpy.factors.discrete.TabularCPD('Burglary', 2, [[0.001, 0.999]])

# Probability of earthquake (True, False)
cpd_earthquake = pgmpy.factors.discrete.TabularCPD('Earthquake', 2, [[0.002, 0.998]])

# Probability of alarm going of (True, False) given a burglary and/or earthquake
cpd_alarm = pgmpy.factors.discrete.TabularCPD('Alarm', 2, [[0.95, 0.94, 0.29, 0.001], 
                                                           [0.05, 0.06, 0.71, 0.999]], 
                                              evidence=['Burglary', 'Earthquake'], 
                                              evidence_card=[2, 2])

# Probability that John calls (True, False) given that the alarm has sounded
cpd_john = pgmpy.factors.discrete.TabularCPD('JohnCalls', 2, [[0.90, 0.05], 
                                                           [0.10, 0.95]], 
                                              evidence=['Alarm'], 
                                              evidence_card=[2])

# Probability that Mary calls (True, False) given that the alarm has sounded
cpd_mary = pgmpy.factors.discrete.TabularCPD('MaryCalls', 2, [[0.70, 0.01], 
                                                           [0.30, 0.99]], 
                                              evidence=['Alarm'], 
                                              evidence_card=[2])

# Add CPDs to the network structure
model.add_cpds(cpd_burglary, cpd_earthquake, cpd_alarm, cpd_john, cpd_mary)

# Check if the model is valid, throw an exception otherwise
model.check_model()

# Print probability distributions
print('Probability distribution, P(Burglary)')
print(cpd_burglary)
print()
print('Probability distribution, P(Earthquake)')
print(cpd_earthquake)
print()
print('Joint probability distribution, P(Alarm | Burglary, Earthquake)')
print(cpd_alarm)
print()
print('Joint probability distribution, P(JohnCalls | Alarm)')
print(cpd_john)
print()
print('Joint probability distribution, P(MaryCalls | Alarm)')
print(cpd_mary)
print()

# Plot the model
nx.draw(model, with_labels=True)
plt.savefig('C:\\DATA\\Python-data\\bayesian-networks\\alarm.png')
plt.close()

# Perform variable elimination for inference
# Variable elimination (VE) is a an exact inference algorithm in bayesian networks
infer = pgmpy.inference.VariableElimination(model)

# Calculate the probability of a burglary if John and Mary calls (0: True, 1: False)
posterior_probability = infer.query(['Burglary'], evidence={'JohnCalls': 0, 'MaryCalls': 0})

# Print posterior probability
print('Posterior probability of Burglary if JohnCalls(True) and MaryCalls(True)')
print(posterior_probability)
print()

# Calculate the probability of alarm starting if there is a burglary and an earthquake (0: True, 1: False)
posterior_probability = infer.query(['Alarm'], evidence={'Burglary': 0, 'Earthquake': 0})

# Print posterior probability
print('Posterior probability of Alarm sounding if Burglary(True) and Earthquake(True)')
print(posterior_probability)
print()

Probability distribution, P(Burglary)
+-------------+-------+
| Burglary(0) | 0.001 |
+-------------+-------+
| Burglary(1) | 0.999 |
+-------------+-------+

Probability distribution, P(Earthquake)
+---------------+-------+
| Earthquake(0) | 0.002 |
+---------------+-------+
| Earthquake(1) | 0.998 |
+---------------+-------+

Joint probability distribution, P(Alarm | Burglary, Earthquake)
+------------+---------------+---------------+---------------+---------------+
| Burglary   | Burglary(0)   | Burglary(0)   | Burglary(1)   | Burglary(1)   |
+------------+---------------+---------------+---------------+---------------+
| Earthquake | Earthquake(0) | Earthquake(1) | Earthquake(0) | Earthquake(1) |
+------------+---------------+---------------+---------------+---------------+
| Alarm(0)   | 0.95          | 0.94          | 0.29          | 0.001         |
+------------+---------------+---------------+---------------+---------------+
| Alarm(1)   | 0.05          | 0.06          | 0.71          | 0.999         |
+------------+---------------+---------------+---------------+---------------+

Joint probability distribution, P(JohnCalls | Alarm)
+--------------+----------+----------+
| Alarm        | Alarm(0) | Alarm(1) |
+--------------+----------+----------+
| JohnCalls(0) | 0.9      | 0.05     |
+--------------+----------+----------+
| JohnCalls(1) | 0.1      | 0.95     |
+--------------+----------+----------+

Joint probability distribution, P(MaryCalls | Alarm)
+--------------+----------+----------+
| Alarm        | Alarm(0) | Alarm(1) |
+--------------+----------+----------+
| MaryCalls(0) | 0.7      | 0.01     |
+--------------+----------+----------+
| MaryCalls(1) | 0.3      | 0.99     |
+--------------+----------+----------+

Finding Elimination Order: : 100%|█████████████| 2/2 [00:00<00:00, 1002.94it/s]
Eliminating: Alarm: 100%|███████████████████████| 2/2 [00:00<00:00, 401.10it/s]
Posterior probability of Burglary if JohnCalls(True) and MaryCalls(True)
+-------------+-----------------+
| Burglary    |   phi(Burglary) |
+=============+=================+
| Burglary(0) |          0.2842 |
+-------------+-----------------+
| Burglary(1) |          0.7158 |
+-------------+-----------------+

Finding Elimination Order: : 100%|█████████████| 2/2 [00:00<00:00, 2008.77it/s]
Eliminating: MaryCalls: 100%|███████████████████| 2/2 [00:00<00:00, 501.08it/s]
Posterior probability of Alarm sounding if Burglary(True) and Earthquake(True)
+----------+--------------+
| Alarm    |   phi(Alarm) |
+==========+==============+
| Alarm(0) |       0.9500 |
+----------+--------------+
| Alarm(1) |       0.0500 |
+----------+--------------+