Peculiar Issue with PyMC3.Potential

Hello, I have only recently got familiar with PyMC3. I have been trying to use the potential function to impose constraints on a distribution. I was trying to run the following piece of code:

*import pymc3 as pm*
*import numpy as np*

*with pm.Model() as model:*
*    2D_binary_vectors = pm.Bernoulli('2D_binary_vectors', p = 0.7, shape=2)*
*    exclude00 = pm.Potential('exclude00', pm.math.switch(pm.math.sum(2D_binary_vectors) >= 1 , 0, -np.inf))*
*    trace = pm.sample(20000)*

*trace [2D_binary_vectors]*

I get output like this:

**array([[1, 1],**
**       [1, 1],**
**       [1, 1],**
**       ...,**
**       [0, 1],**
**       [0, 1],**
**       [1, 1]], dtype=int64)**

Which is expected! There are no [0,0] which is what I expect the potential would do. However, if I decrease the value of the parameter in the Bernoulli prior to p = 0.5, the same code gives me nothing but [0,0]. So, for the following code:

**with pm.Model() as model:**
**    2D_binary_vectors = pm.Bernoulli('2D_binary_vectors', p = 0.5, shape=2)**
**    exclude00 = pm.Potential('exclude00', pm.math.switch(pm.math.sum(2D_binary_vectors) >= 1 , 0, -np.inf))**
**    trace = pm.sample(20000)**

**trace [2D_binary_vectors]**

I get:

**array([[0, 0],**
**       [0, 0],**
**       [0, 0],**
**       ...,**
**       [0, 0],**
**       [0, 0],**
**       [0, 0]], dtype=int64)**

I get nothing but [0,0], which was the only thing I was trying to exclude. I don’t understand what’s going on here.

You mention “if I decrease the value of the parameter in the Bernoulli prior to p = 0.5” but you’ve got p = 0.5 in both examples. However, in the second you’ve got np.inf whereas in the first you have -np.inf. I couldn’t run either example because I don’t have the ‘dorm1’ variable, but maybe check if the sign is making a difference.

1 Like

I edited the post to correct the mistakes. The code I tried to run looked like this

This worked for me:

with pm.Model() as model:
    _binary_vectors = pm.Bernoulli('_binary_vectors', p = 0.3, shap
      e=2)
    exclude00 = pm.Potential('exclude00', pm.math.switch(
        pm.math.sum(_binary_vectors) <= 0 , (1,1), _binary_vectors))

pred = pm.sample_prior_predictive(500, model=model)

I think sample is just taking from the original distribution of the Bernoulli distribution, not the distribution adjusted by the Potential - possibly because there aren’t any observed values to update the prior with, but I’m not sure (I’m new to pymc3 too). Note that with what I’ve got (0, 0) is being replaced with (1, 1), which would affect the probability of drawing it - you’d have to do something a bit different to keep the distribution you want

Note that forward sampling distributions (i.e sample_prior/posterior_predictive) ignore Potentials – only sample takes them into account; see this issue, for which help is welcome by the way :slight_smile:
This means that to have the constraints in your forward samples, you need to implement the contraints by hand

Thank you for the reply. I think the problem is I don’t understand the use of the switch function.

Thank you for the reply. Could you elaborate on what you mean by implementing by hand?

Switch takes the format if cond then ift else iff, as you probably know.

So what I thought was happening - which may well be incorrect, @AlexAndorra your input would be most welcome - is that the Potential would examine each pair of outputs from the Bernoulli distribution, and then the switch if their sum was 0 (cond), replace the pair with (1, 1) (ift), otherwise it would output the pair from the Bernoulli (iff).

After reading the link Alex provided, it seems Potential should influence logp of the Bernoulli. I tried

with model:
    trace = pm.sample(2000, tune=1000)
    az.summary(trace)

and got
image
The mean is quite a bit lower than 0.3, so obviously there is some effect from the Potential. However,

trace[_binary_vectors]

gives
image

And I was about to write ‘whereas pred had no [0,0] values’, but I just tried it again to make sure and it does. My apologies - not sure what I did yesterday that made me think it didn’t have them

EDIT:
pred is a dictionary with values from _binary_vectors and exclude00. So pred[‘exclude00’] looks like
image
which should be what you want

1 Like