Sum of two Binomial Distributions

Of course things should work with one data point (running without an error does not mean that you’ll get the answer). The problem is that you can’t use deterministic as the likelihood. I don’t have sufficient knowledge on your exact usecase to decide the suitable likelihood for that.

Assume that suitable likelihood is Normal distribution.

n = 100
y = 50
fn = .2
fp = .1
with pm.Model() as food_model:
    positives = pm.DiscreteUniform("positives", lower=0, upper=n)
    true_positives = pm.Binomial("true_positives", n=positives, p=(1 - fn))
    false_positives = pm.Binomial("false_positives", n=y - positives, p=fp)
    model = pm.Normal('model', true_positives + false_positives, observed=y)

with food_model:
    # Draw samples
    trace = pm.sample(1000, njobs=2)
    # Plot parameter
    pm.forestplot(trace, varnames=['positives'])

Now this will work. It is because deterministic can not be a likelihood of a probabilistic model.

When I said the probabilistic model, I requested the equations (or conditional relationship) between the each random variable and how they should be mapped to likelihood.

Lets see if someone else may have a better insight on this problem