Solved – the probability P(X > Y) given X ~ Be(a1, b1), and Y ~ Be(a2, b2), and X and Y are independent

$Pr(X > Y) = int_0^1 frac{x^{a_1 – 1}(1 – x)^{b_1 – 1}}{text{Be}(a_1,b_1)} int_0^xfrac{y^{a_2 – 1}(1 – y)^{b_2 – 1}}{text{Be}(a_2,b_2)} dy dx$

$Pr(X > Y) = frac{1}{text{Be}(a_1,b_1)} int_0^1 x^{a_1 – 1}(1 – x)^{b_1 – 1}I_x(a_2, b_2) dx$

where $I_x(a, b)$ is the regularized incomplete beta function. If $a$ and $b$ are integers then

$I_x(a,b) = sum_{j=a}^{a+b-1} {(a+b-1)! over j!(a+b-1-j)!} x^j (1-x)^{a+b-1-j}.$

Substitute in, do some simple algebra, and the integral will have a closed form solution as a finite sum of beta functions.

If $a_2$ and $b_2$ aren't integers but $a_1$ and $b_1$ are, then calculate $Pr(X > Y) = 1 – Pr(Y > X)$. If neither case holds, you're pooched for an analytical solution but you can always do the integral numerically, either deterministically or by Monte Carlo.