# Inconsistencies between conditional probability calculations by hand and with pgmpy (Bayesian Graphical Models)

I am teaching myself about Bayesian graphical networks. I'm attempting to use the python package pgmpy to generate the networks in python. This seems like a great resource.

For my first test, I generated a simple network depicted below (I set the known probabilities and conditional probabilities to infer the unconditional probabilities):

Now, I entered the probabilities for $A$ and $B$, as well as the probabilities for $P(C|A,B)$ into a bayesian graphical model structure in pgmpy:

IN:
#These are based on the the Monte Hall example found at https://github.com/pgmpy/pgmpy/blob/dev/examples/Monte%20Hall%20Problem.ipynb

from pgmpy.models import BayesianModel
from pgmpy.factors import TabularCPD

# Defining the network structure
model = BayesianModel([('A', 'C'), ('B', 'C')])

# Defining the CPDs:
cpd_p = TabularCPD('A', 2, [[0.99, 0.01]])
cpd_a = TabularCPD('B', 2, [[0.9, 0.1]])
cpd_t = TabularCPD('C', 2, [[0.9, 0.5, 0.4, 0.1],
[0.1, 0.5, 0.6, 0.9]],
evidence=['A', 'B'], evidence_card=[2, 2])

# Associating the CPDs with the network structure.

# Some other methods
model.get_cpds()

OUT:
[<TabularCPD representing P(A:2) at 0x10e24cfd0>,
<TabularCPD representing P(B:2) at 0x10e24cf10>,
<TabularCPD representing P(C:2 | A:2, B:2) at 0x10df9a750>]


However, when I calculate the probabilities in pgmpy I get:

IN:
# Infering the posterior probability
from pgmpy.inference import VariableElimination

print 'P(B|A=1,C=1)'
infer = VariableElimination(model)
posterior_p = infer.query(['B'], evidence={'A': 1, 'C': 1})
print(posterior_p['B'])

print 'P(B|C=1)'
posterior_p = infer.query(['B'], evidence={'C': 1})
print(posterior_p['B'])

print 'probs'
posterior_p = infer.query(['B','C','A'])
for entry in posterior_p:
print posterior_p[entry]

OUT:
P(B|A=1,C=1)
+-----+----------+
| B   |   phi(B) |
|-----+----------|
| B_0 |   0.8333 |
| B_1 |   0.1667 |
+-----+----------+
P(B|C=1)
+-----+----------+
| B   |   phi(B) |
|-----+----------|
| B_0 |   0.6082 |
| B_1 |   0.3918 |
+-----+----------+
probs
+-----+----------+
| A   |   phi(A) |
|-----+----------|
| A_0 |   0.9900 |
| A_1 |   0.0100 |
+-----+----------+
+-----+----------+
| C   |   phi(C) |
|-----+----------|
| C_0 |   0.8461 |
| C_1 |   0.1539 |
+-----+----------+
+-----+----------+
| B   |   phi(B) |
|-----+----------|
| B_0 |   0.9000 |
| B_1 |   0.1000 |
+-----+----------+


However, when I calculate the conditional probabilities (or total probabilities for variable $C$) shown above by hand (calculated from the unconditional probabilities), I get different answers:

$$P(B=1|A=1,C=1) = \frac{P(B=1|A=1,C=1)}{P(B=0|A=1,C=1)+P(B=1|A=1,C=1)} = \frac{0.0009}{0.0054+0.0009} = 0.143$$

$$P(B=1|C=1) = \frac{P(B=1,C=1)}{P(C=1)} = \frac{0.0495+0.0009}{0.0891+0.0495+0.0054+0.0009} = 0.348$$

$$P(C=1) = 0.0891 + 0.0495 + 0.0054 + 0.0009 = 0.145$$

Obviously this is different than the probabilities output by pgmpy.

Does anyone have a clue about where I'm going wrong? Either with my "by-hand" calculations or with the coding.

THANKS!!!

Seems as thought I figured this out. I'm posting it here just in case it helps someone.

Well, it turns out that I wasn't making an error with my "by-hand" calculation, but it was indeed a problem with the pgmpy package. For some reason, it was inferring independence between $A$ and $B$ given $C$ (they should be independent without knowing C).

I updated the package to the developer edition and it gives the desired result.

• Thank you SO much for answering this here, I've been pulling my freaking HAIR OUT trying to figure out why pgmpy is giving me the wrong conditional probabilities! I updated to dev like you said, and it works now. It's really unacceptable to have such a major issue in the production version of such a widely used package. It was causing some seriously egregious errors. See here: stackoverflow.com/questions/47087068/… Nov 3 '17 at 5:37

I tried installing the new package and the bug remains. So i looked at the Bayesian network code itself and found the bug.

Its in line 530 of BayesianModel.py when finding the descendents. Because neighbors is an iterator, it is empty after visit.extend and so descendants.extend does nothing and the result is empty. I fixed it in place by casting neighbors as a list.

    Since Bayesian Networks are acyclic, this is a very simple dfs
which does not remember which nodes it has visited.
"""
descendents = []
visit = [node]
while visit:
n = visit.pop()
neighbors = self.neighbors(n)
visit.extend(neighbors)
**descendents.extend(neighbors)**
return descendents


Confirmed this works well (matches the hand calculations) in version 0.1.14. Note I had to make changes to the code, to fix both errors in how pgmpy was being used (incorrect shapes) and for Python 3. Here is the updated code.

from pgmpy.models import BayesianModel
from pgmpy.factors.discrete import TabularCPD

# Defining the network structure
model = BayesianModel([('A', 'C'), ('B', 'C')])

# Defining the CPDs:
cpd_p = TabularCPD('A', 2, [[0.99], [0.01]])
cpd_a = TabularCPD('B', 2, [[0.9], [0.1]])
cpd_t = TabularCPD('C', 2, [[0.9, 0.5, 0.4, 0.1],
[0.1, 0.5, 0.6, 0.9]],
evidence=['A', 'B'], evidence_card=[2, 2])

# Associating the CPDs with the network structure.

# Some other methods
model.get_cpds()

from pgmpy.inference import VariableElimination

print('P(B|A=1,C=1)')
infer = VariableElimination(model)
posterior_p = infer.query(['B'], evidence={'A': 1, 'C': 1})
print(posterior_p)

print('P(B|C=1)')
posterior_p = infer.query(['B'], evidence={'C': 1})
print(posterior_p)

print('probs')
posterior_p = infer.query(['B','C','A'])
print(posterior_p)