Using PyMC to find the probability of the difference between two Binomial experiments

bernoulli_difference.py

%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
import pymc as pm
import numpy as np

p = 0.5

lambda_1 = pm.Binomial("lambda_1", 1000, p)
lambda_2 = pm.Binomial("lambda_2", 1000, p)

@pm.deterministic
def lambda_diff(lambda_1=lambda_1, lambda_2=lambda_2):
    return lambda_1-lambda_2

model = pm.Model([lambda_1, lambda_2, lambda_diff])
mcmc = pm.MCMC(model)
mcmc.sample(200000, 100000, 1)

tau_samples = mcmc.trace('lambda_diff')[:]
plt.hist(tau_samples, histtype='stepfilled', bins=30, alpha=0.85,
         label="posterior of $\lambda_1$", color="#A60628", normed=True) 
plt.show()
print 
print "The probability of having A > (B+25) is ~ %0.2f%%" % (100*float(len(tau_samples[tau_samples < 25]))/len(tau_samples))