[SOLVED] How to make Conditional Probability Tables (CPTs) for Bayesian networks with pymc

How to make Conditional Probability Tables (CPTs) for Bayesian networks with pymc

I would like to build a Bayesian network of discrete (pymc.Categorical) variables that are dependent on other categorical variables. As a simplest example, suppose variables a and b are categorical and b depends on a

Here is an attempt to code it with pymc (assuming a takes one of three values and b takes one of four values). The idea being that the CPT distributions would be learned from data using pymc.

import numpy as np
import pymc as pm
aRange = 3
bRange = 4

#make variable a
a = pm.Categorical('a',pm.Dirichlet('aCPT',np.ones(aRange)/aRange))

#make a CPT table as an array of 
CPTLines = np.empty(aRange, dtype=object)
for i in range(aRange):
    CPTLines[i] = pm.Dirichlet('CPTLine%i' %i,np.ones(bRange)/bRange)

#make a deterministic node that holds the relevant CPT line (dependent on state1)
@pm.deterministic
def selectedCPTLine(CPTLines=CPTLines,a=a):
    return CPTLines[a]

#make a node for variable b 
b=pm.Categorical('b', selectedCPTLine)

model = pm.MCMC([a, b, selectedCPTLine])

If we draw this model it looks like this

However, running this code we get an error:

Probabilities in categorical_like sum to [ 0.8603345]

Apparently, pymc can take a Dirichlet variable as the parameter of a Categorical variable. When the Categorical variable gets a Dirichlet variable as its parameter, it knows to expect a k-1 vector of probabilities with the assumption that the kth probability sums the vector to 1. This breaks down, however, when the Dirichlet variable is the output of a deterministic variable, which is what I need to make a CPT.

Am I going about this the right way? How can the representation mismatch be solved? I should mention that I'm relatively new to pymc and Python.

This question is related to a previous question on making a discrete state Markov model with pymc

Solution

OK, thanks. The problem is that, while usually, PyMC will recognize a Dirichlet as the parent of a Categorical and completes the probability simplex, here your Categoricals are embedded in a Container, and the Categorical does not make the automatic adjustment needed. The following code does this for you:

import numpy as np
import pymc as pm
aRange = 3
bRange = 4

aCPT = pm.Dirichlet('aCPT', np.ones(aRange))

#make variable a
a = pm.Categorical('a', aCPT)

#make a CPT table as an array of
CPTLines = [pm.Dirichlet('CPTLine%i' %i, np.ones(bRange)) for i in range(aRange)]

#make a node for variable b
@pm.stochastic(dtype=int)
def b(value=0, CPT=CPTLines, a=a):
    return pm.categorical_like(value, p=pm.extend_dirichlet(CPT[a]))

model = pm.MCMC([a, b, CPTLines])

Hope that helps.