Some problems with estimating total size for Binomial distribution

Dear all,

I have a quite simple problem: if I know p for Binomial distribution, the observed counts C, then I want to predict the total size. First, if I try to find MAP (or MLE), then it simple does not run at all:

with pm.Model() as model:
    i0 = pm.HalfFlat('i0', testval = 100.0, shape=6)
    C = pm.Data('C', np.array([1, 2, 3, 4, 6, 8]))
    pm.Binomial('like', i0, .3, observed = C)
    sol = pm.find_MAP()

results in:

logp = -147.2, ||grad|| = 0: 100%|██████████| 2/2 [00:00<00:00, 940.85it/s]
{'i0_log__': array([4.60517019, 4.60517019, 4.60517019, 4.60517019, 4.60517019,
        4.60517019]), 'i0': array([100., 100., 100., 100., 100., 100.])}

The sampling using NUTS

    sol = pm.sample()

is slow and results in divergences.

Could you point me out where I am wrong? Because I’ve tried to re-write it as negative binomial, etc., but results were the same. If I need to estimate p from using the total size, then I have no problems. But I wonder why I cannot do the opposite and try to estimate the total size by knowing the value of p.

Thank you in advance for your time and checking my question.