stable_prob_b_beat_a

create function 
prob_b_beat_a
(total_a int, conv_a int,  total_b int, conv_b int)

--INPUTS:
-- total_a = total number that saw A
-- conv_a = converters for version A
-- total_b = total number that saw B
-- conv_b = converters for version B

RETURNS float
--OUTPUT: 
--Returns the probability that variation B
--will beat variation A in the long run 
STABLE 
AS $$

    import math as math
    from scipy.special import betaln
    alpha_a = 1 + conv_a
    beta_a = 1 + total_a - conv_a
    alpha_b = 1 + conv_b
    beta_b = 1 + total_b - conv_b

    total = 0.0
    for i in range(alpha_b):
        total += math.exp(betaln(alpha_a+i, beta_b+beta_a) 
            - math.log(beta_b+i) - betaln(1+i, beta_b) - betaln(alpha_a, beta_a))
    return total

$$LANGUAGE plpythonu;