effective-engineering-1.py

import pymc3 as pm

da_points_solved = [47,55,39,36,27,37,31,33,39,38,20,52,36,40,10,41,42,11,20,15,26,30,11,22,25,19,51,29,17,29,39,30,54,35]
da_points = [82,68,64,45,27,48,43,46,52,46,39,59,55,59,35,65,61,18,27,37,38,42,34,27,30,37,51,36,30,42,47,47,62,57]

sw_points_solved = [8,11,6,6,12,12,7,15,13,6,17,12,3,11,6,7,6,9,10,0,0,6,17,16,19,9,15,10,14,8,24,17]
sw_points = [11,17,12,12,12,12,16,16,13,16,22,20,3,11,6,17,13,15,10,10,10,6,17,16,24,14,15,15,24,15,24,17]

with pm.Model() as model:
    p_sw = pm.Beta('completion_rate_sw', alpha=1, beta=1)
    solved_sw = pm.Binomial('solved_points_sw', n=sw_points, p=p_sw, observed=sw_points_solved)
    
    p_da = pm.Beta('completion_rate_da', alpha=1, beta=1)
    solved_da = pm.Binomial('solved_points_da', n=da_points, p=p_da, observed=da_points_solved)
    
    diff = pm.Deterministic('diff', p_sw - p_da)
    
    trace = pm.sample(10000)

ax = pm.plot_posterior(trace, varnames=['completion_rate_sw', 'completion_rate_da'], figsize=(15,5))
ax[0].set_title('$p_{SW}$')
ax[1].set_title('$p_{DA}$');
pm.plot_posterior(trace, varnames=['diff'], ref_val=0, text_size=16, figsize=(7,5));
print('P(p_sw > p_da) = %.2f' % ((trace['diff'] > 0).mean()))