Hierarchical bernoulli model with bernoulli priors

Questions
1

I’ve defined a simple model using the following:

basic_model = pm.Model()

with pm.Model() as basic_model:
    p = pm.Bernoulli('p', p=0.5)
    out = pm.Bernoulli('out', p=0.9 * p, observed=1)
    data = pm.sample(1000)

which I would expect to observe p as always 1. However, when I sample from it, p is always 0. I can’t figure out why this would be. Any ideas? I would have assumed it was an issue with p returning an int dtype and multiplication failing, but this doesn’t seem to be the case. The default sampler used is BinaryGibbsMetropolis, which might be problematic?

2

I’m not sure what the root cause would be since I got stuck while debugging, but I will note that the problem appears to be the likelihood from the observed variable. The output of basic_model.observed_RVs[0].logp({'p':0}) is -np.inf which will lead the BinaryGibbsMetropolis proposals for p=0 to always get rejected.