Hierarchical bernoulli model with bernoulli priors

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?

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.