Probability P(Theta1 > Theta2) for independent beta-distributed random variables

compare_binomial.m

function p = compare_binomial(x1, n1, x2, n2)
% Compute P(Beta(x1+1,n1-x1+1) > Beta(x2+1,n2-x2+1)).

% Compute p(x1, n1, x2, x2).
p = 1;
for k = 0:x2
    p = p * (x1 + 1 + k) / (n1 + 2 + k);
end

% Iterate Bayes updates with each of n2 - x2 failures.
q = p * (n1 - x1 + 1) * (x2 + 1) / (n1 + x2 + 3);
for k = 1:(n2 - x2)
    p = p + q / k;
    q = q * (n1 - x1 + 1 + k) * (x2 + 1 + k) / (k * (n1 + x2 + 3 + k));
end