Sum of two Binomial Distributions

@Nadheesh is right, the problem here is that you cannot assign observed to pm.Deterministic, in this case it is also difficult to get a close form likelihood function for the observed Positive.

We are working to implement a likelihood-free algorithm (Approximate Bayesian Computation) for these cases. Currently, this is an ongoing project by @agustinaarroyuelo, @aloctavodia and me - Hopefully, this will be selected as one of the GSOC 2018 projects.

In this relatively simple case, you can try a reject sampling implemented in numpy:

n = 100
y = 50
fn = .2
fp = .1

# reject sampling
Nsample = 1000
yhat = 0
TrueP = []
for i in range(Nsample):
    # only accept if the simulated data is the same as the observed
    while yhat != y:
        # simulate from latent true postive, can be replaced with informative prior
        trueP = np.random.randint(0, n)
        # simulate observation
        yhat = np.random.binomial(trueP, 1 - fn) +\
               np.random.binomial(n - trueP, fp)
    TrueP.append(trueP)
    yhat = 0
TrueP = np.asarray(TrueP)

# display result
plt.hist(TrueP, np.arange(n));

image

2 Likes