I’ve put a lot of effort into writing software for evaluating random inequality probabilities with beta distributions because such inequalities come up quite often in application. For example, beta inequalities are at the heart of the Thall-Simon method for monitoring single-arm trials and adaptively randomized trials with binary endpoints.
It’s not easy to evaluate P(X > Y) accurately and efficiently when X and Y are independent random variables. I’ve seen several attempts that were either inaccurate or slow, including a few attempts on my part. Efficiency is important because this calculation is often in the inner loop of a simulation study. Part of the difficulty is that the calculation depends on four parameters and no single algorithm will work well for all parameter combinations.
Let g(a, b, c, d) equal P(X > Y) where X ~ beta(a, b) and Y ~ beta(c, d). Then the function g has several symmetries.
- g(a, b, c, d) = 1 – g(c, d, a, b)
- g(a, b, c, d) = g(d, c, b, a)
- g(a, b, c, d) = g(d, b, c, a)
The first two relations were published by W. R. Thompson in 1933, but as far as I know the third relation first appeared in this technical report in 2003.
For special values of the parameters, the function g(a, b, c, d) can be computed in closed form. Some of these special cases are when
- one of the four parameters is an integer
- a + b + c + d = 1
- a + b = c + d = 1.
The function g(a, b, c, d) also satisfies several recurrence relations that make it possible to bootstrap the latter two special cases into more results. Define the beta function B(a, b) as Γ(a, b)/(Γ(a) Γ(b)) and define h(a, b, c, d) as B(a+c, b+d)/( B(a, b) B(c, d) ). Then the following recurrence relations hold.
- g(a+1, b, c, d) = g(a, b, c, d) + h(a, b, c, d)/a
- g(a, b+1, c, d) = g(a, b, c, d) – h(a, b, c, d)/b
- g(a, b, c+1, d) = g(a, b, c, d) – h(a, b, c, d)/c
- g(a, b, c, d+1) = g(a, b, c, d) + h(a, b, c, d)/d
For more information about beta inequalities, see these papers: