Find_MAP when optimizing a Potential

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)}