Here’s a neat example from chapter 7 of Risk Assessment and Decision Analysis with Bayesian Networks, by Fenton & Neil to check the numbers are coming out as expected.
import pymc3 as pm
import numpy as np
import theano
import theano.tensor as tt
names = ['martin_oversleeps', 'train_strike', 'martin_late', 'norman_late']
# P(martin_late|train_strike, martin_oversleeps) -------------------------
martin_late_lookup = theano.shared(np.asarray([
[[.7, .3], [.4, .6]],
[[.4, .6], [.2, .8]]]))
def p_martin_late(train_strike, martin_oversleeps):
return martin_late_lookup[train_strike, martin_oversleeps]
# P(norman_late|train_strike) --------------------------------------------
p_norman_late_lookup = theano.shared(np.asarray([[0.9, .1], [.2, .8]]))
def p_norman_late(train_strike):
return p_norman_late_lookup[train_strike]
# Build model ------------------------------------------------------------
with pm.Model() as m:
martin_oversleeps = pm.Categorical('martin_oversleeps', [.6, .4])
train_strike = pm.Categorical('train_strike', [.9, .1])
norman_late = pm.Categorical('norman_late',
p_norman_late(train_strike))
martin_late = pm.Categorical('martin_late',
p_martin_late(train_strike, martin_oversleeps))
prior = pm.sample_prior_predictive(500_000)
Then added some helper functions
def marginal_distributions(trace, names):
"""Samples from the prior.
trace = samples from prior
names = list of strings of node names
"""
for name in names:
print(f"P({name}) = {trace[name].mean()*100:.2f}%")
def conditional_probabilities(trace, names, evidence):
"""Output the probabilities of all the variables, conditioned on the evidence.
trace = samples from prior
names = list of strings of node names
evidence = {'norman_late': 1, 'martin_late': 1} (for example)
"""
n_samples = len(list(trace.values())[0])
for name in names:
# Calculate Boolean subset of samples confiorming to the evidence
subset = np.full(n_samples, True, dtype=bool)
for key, val in evidence.items():
subset = (subset) & (trace[key] == val)
# Calculate mean of variable of interest consistent with the evidence
prob = trace[name][subset == 1].mean()
# Compose string
string = f"P({name} | "
for key, value in evidence.items():
string += key + f" = {value}"
string += ", "
string += f") = {prob*100:.2f}%"
print(string)
So now you can do things like
marginal_distributions(prior, names)
and get
P(martin_oversleeps) = 40.00%
P(train_strike) = 10.07%
P(martin_late) = 44.58%
P(norman_late) = 17.01%
or condition upon evidence provided in the form of a dictionary
conditional_probabilities(prior, names, {'norman_late': 1, 'martin_late': 1})
and get
P(martin_oversleeps | norman_late = 1, martin_late = 1, ) = 51.51%
P(train_strike | norman_late = 1, martin_late = 1, ) = 59.41%
P(martin_late | norman_late = 1, martin_late = 1, ) = 100.00%
P(norman_late | norman_late = 1, martin_late = 1, ) = 100.00%
Pretty cool. Although I’ve not checked it’s robust, eg for nodes with >2 category levels