Given two independent variables with Beta distribution, $X \sim \text{Be}(a_1, b_1)$ and $Y \sim \text{Be}(a_2, b_2)$, how do you find the probability that the value of X is greater than the value of Y for a given observation?
Does this probability have a name that I'm just blanking on?
$\Pr(X > Y) = \int_0^1 \frac{x^{a_1 - 1}(1 - x)^{b_1 - 1}}{\text{Be}(a_1,b_1)} \int_0^x\frac{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.
