In [67]:

    

from scipy.special import betaln as lbeta
from scipy.stats import beta as beta_dist, norm as norm_dist
from numpy import log, exp


    

In [113]:

    

def probability_B_beats_A_evan(α_A, β_A, α_B, β_B):
    total = 0.0
    for i in range(α_B):
        total += exp(lbeta(α_A+i, β_B+β_A) - log(β_B+i) - lbeta(1+i, β_B) - lbeta(α_A, β_A))
    return total


    

In [114]:

    

def probability_B_beats_A_mc(α_A, β_A, α_B, β_B, N=1000000):
    p_A = beta_dist.rvs(α_A, β_A, size=N)
    p_B = beta_dist.rvs(α_B, β_B, size=N)
    return (p_B>p_A).mean()


    

In [115]:

    

def beta_mean(a, b):
    return a/(a+b)
def beta_var(a, b):
    return a*b/((a+b)**2*(a+b+1))


    

In [116]:

    

def probability_B_beats_A(α_A, β_A, α_B, β_B):
    mu = beta_mean(α_A, β_A) - beta_mean(α_B, β_B)
    sigma = (beta_var(α_A, β_A) + beta_var(α_B, β_B))**.5
    return norm_dist.cdf(0, mu, sigma)


    

In [121]:

    

α_A, β_A, α_B, β_B = 21, 200, 31, 200
print(probability_B_beats_A_evan(α_A, β_A, α_B, β_B))
print(probability_B_beats_A(α_A, β_A, α_B, β_B))


    

    




0.906704245199
0.905669069182

In [122]:

    

%%timeit
probability_B_beats_A_evan(α_A, β_A, α_B, β_B)


    

    




1000 loops, best of 3: 461 µs per loop

In [123]:

    

%%timeit
probability_B_beats_A(α_A, β_A, α_B, β_B)


    

    




10000 loops, best of 3: 96.9 µs per loop