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
    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):
        num = math.lgamma(alpha_a+i) + math.lgamma(beta_a+beta_b) + math.lgamma(1+i+beta_b) + math.lgamma(alpha_a+beta_a)
        den = math.log(beta_b+i) + math.lgamma(alpha_a+i+beta_a+beta_b) + math.lgamma(1+i) + math.lgamma(beta_b) + math.lgamma(alpha_a) + math.lgamma(beta_a)
        total += math.exp(num - den)
    return total

$$LANGUAGE plpythonu;