Welcome!
There’s something going on with the sum calculated over samples, though I am not 100% what. Sampling works:
idata = pm.sample()
idata.posterior["p"].mean()
yields:
<xarray.DataArray 'p' ()>
array(0.5000402)
And this:
# Generate the Bernoulli samples
#samples = pm.Bernoulli('samples', p=p, shape=N)
samples = pm.Binomial('samples', p=p, n=N)
yields:
{'p_interval__': array(0.),
'samples': array(50),
'p': array(0.5),
'samples_sum': array(50)}