# Understanding factor potentials in PyMC

I'm trying to understand factor potentials from the PyMC documentation, but need some help on the implementation piece--or it may turn out that I am misunderstanding how potentials work altogether.

Imagine that we are building a Poisson switchpoint model as specified in PyMC's documentation tutorial. Now we want to introduce a factor potential exactly as described in the documentation, such that the difference between early and late Poisson means is less than 1.

When I tried to implement the potential as in the code below, I see that values outside the allowed constraint appear in the posterior distribution. Why does the posterior distribution of potentialCheck contain these values? Obviously, I am doing something wrong...

from pymc import DiscreteUniform, Exponential, deterministic, Poisson, Uniform, potential
import numpy as np

disasters_array =   \
np.array([ 4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0,
1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])

switchpoint = DiscreteUniform('switchpoint', lower=0, upper=110, doc='Switchpoint[year]')

early_mean = Exponential('early_mean', beta=1.)
late_mean = Exponential('late_mean', beta=1.)

@deterministic(plot=False)
def rate(s=switchpoint, e=early_mean, l=late_mean):
''' Concatenate Poisson means '''
out = empty(len(disasters_array))
out[:s] = e
out[s:] = l
return out

@potential(plot=True)
def examplePotential(em=early_mean, lm=late_mean ):
if abs(em-lm) < 1:
return e
else:
return 1.

@deterministic(plot=True)
def potentialCheck(e=early_mean, l=late_mean):
'''Replicate the Potential to plot and check if the constraint has held.'''
return abs(e-l)

disasters = Poisson('disasters', mu=rate, value=disasters_array, observed=True)

disaster_model = [switchpoint, early_mean, late_mean, examplePotential, potentialCheck]

from pymc import MCMC
M = MCMC(disaster_model)
M.sample(iter=10000, burn=1000, thin=10)

from pymc.Matplot import plot as mcplot

mcplot(M)


## 3 Answers

I agree with Chris's answer. If you make the following change to your code, it should work:

@potential(plot=True)
def examplePotential(em=early_mean, lm=late_mean ):
if abs(em-lm) < 1:
return 0. ####### No penalty applied if conditions are met
else:
return -np.inf ##### Infinite potential applied if not met


You will also have to change the initial value of early_mean and late_mean so that they satisfy the potential conditions. I just gave a value of 1 to both these quantities.

• Thanks for the answer! Exactly what I needed to get it working. – degenerate hessian Apr 23 '14 at 21:15

I'm following the same tutorial during these days so I'm far from being an expert but I saw a probable typo in your code:

@potential(plot=True)
def examplePotential(em=early_mean, lm=late_mean ):
if abs(em-lm) < 1:
return e  <<<<<<< return em !!!!!!!!
else:
return 1.


With this change it works fine for me.

• There may be a typo somewhere, but this was intentional. The logic is that I'm trying to slip either a 1 or a zero into the factor product based on the indicator value. Since potential should return log probabilities, I am using the constant e for 1 and 1 for 0. – degenerate hessian Apr 15 '14 at 14:47

The factor potential should return a log-probability, not a parameter value. Remember, it is just an arbitrary log-probability term -- it has no value itself.

For example, if you wanted to constrain a particular value to be greater than zero, you could have a factor potential that returned -inf for all values of the parameter that were non-positive.