Bernoulli with p=0.5

I’m trying to define a prior over adjacency matrices for directed unweighted graphs. Each graph is represented as a binary matrix; I tried to sample matrices as Bernoulli RV with shape (N,N).

Using p different than 0.5 gives sensible results, for example:

with pm.Model() as m:
    a = pm.Bernoulli('a',p=0.49,shape=(4,4))
    trace = pm.sample()
trace['a']

gives:

array([[[0, 1, 1, 1],
        [0, 1, 1, 0],
        [0, 0, 0, 0],
        [0, 0, 1, 0]],

       [[1, 0, 0, 0],
        [1, 0, 0, 1],
        [1, 1, 1, 1],
        [1, 1, 0, 1]],

       [[0, 1, 1, 1],
        [0, 1, 1, 0],
        [0, 0, 0, 0],
        [0, 0, 1, 0]],

       ..., 
       [[0, 0, 1, 0],
        [0, 1, 0, 0],
        [1, 1, 1, 1],
        [1, 0, 1, 1]],

       [[1, 1, 0, 1],
        [1, 0, 1, 1],
        [0, 0, 0, 0],
        [0, 1, 0, 0]],

       [[0, 0, 1, 0],
        [0, 1, 0, 0],
        [1, 1, 1, 1],
        [1, 0, 1, 1]]])

However, by setting p=0.5 the sampling “gets stuck” and only gives two unique samples:

with pm.Model() as m:
    a = pm.Bernoulli('a',p=0.5,shape=(4,4))
    trace = pm.sample()
trace['a']

yields:

array([[[1, 1, 1, 1],
    [1, 1, 1, 1],
    [1, 1, 1, 1],
    [1, 1, 1, 1]],

   [[0, 0, 0, 0],
    [0, 0, 0, 0],
    [0, 0, 0, 0],
    [0, 0, 0, 0]],

   [[1, 1, 1, 1],
    [1, 1, 1, 1],
    [1, 1, 1, 1],
    [1, 1, 1, 1]],

   ..., 
   [[0, 0, 0, 0],
    [0, 0, 0, 0],
    [0, 0, 0, 0],
    [0, 0, 0, 0]],

   [[1, 1, 1, 1],
    [1, 1, 1, 1],
    [1, 1, 1, 1],
    [1, 1, 1, 1]],

   [[0, 0, 0, 0],
    [0, 0, 0, 0],
    [0, 0, 0, 0],
    [0, 0, 0, 0]]])

Is this a bug?

Are you on master? This is recently fixed.

No, I installed pymc using conda (-forge).

If you update to master this is fixed. Otherwise, you need to wait for the next release (we are scheduling it).