4
$\begingroup$

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)
$\endgroup$
1
$\begingroup$

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.

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

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – degenerate hessian Apr 15 '14 at 14:47
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.