Evaluating Best-Arm Probabilities for Two Beta Posteriors

From Two Posterior Distributions to a Probability

Suppose the uncertain response rates for two treatments are represented by independent posterior random variables $X$ and $Y$. The quantity of interest is the probability that the first treatment’s rate is larger: $P(X>Y)$. This is a probability over the posterior distributions, not a statement that the first treatment must have the larger true response rate. For continuous Beta distributions, $P(X=Y)=0$, so the two strict orderings exhaust the possibilities.

For positive shape parameters, write $X\sim\operatorname{Beta}(a,b)$ and $Y\sim\operatorname{Beta}(c,d)$. Their densities on $0<x<1$ are proportional to $x^{a-1}(1-x)^{b-1}$ and $y^{c-1}(1-y)^{d-1}$, respectively. Here $B(u,v)$ is the Beta function, defined by $B(u,v)=\int_0^1 t^{u-1}(1-t)^{v-1}\,dt$. Independence is essential: it makes the joint density the product of the two marginal densities.

Definition: For independent $X\sim\operatorname{Beta}(a,b)$ and $Y\sim\operatorname{Beta}(c,d)$, the two-arm best-arm probability for the first arm is $P(X>Y)$. If $F_Y(x)=P(Y\leq x)$, then $$ P(X>Y)=\int_0^1 f_X(x)F_Y(x)\,dx =\frac{1}{B(a,b)}\int_0^1 x^{a-1}(1-x)^{b-1}F_Y(x)\,dx. $$

To obtain the integral, condition on the value $X=x$. Given that value, the event $X>Y$ is $Y<x$, whose probability is $F_Y(x)$. Averaging that conditional probability over the density of $X$ gives the formula. Equivalently, integrating the independent joint density over $0<y<x<1$ gives

P(X>Y)= \frac{1}{B(a,b)B(c,d)} \int_0^1\int_0^x x^{a-1}(1-x)^{b-1}y^{c-1}(1-y)^{d-1}\,dy\,dx.

The assumptions used here are simply that both shape parameters for each distribution are positive, the variables have the stated Beta distributions, and the draws are independent. In a posterior comparison, independence is appropriate when the two arm-specific posterior distributions factor given the data and model. If the posterior draws are dependent, multiplying their marginal densities is not justified; the probability must instead be calculated from their joint posterior distribution.

A Finite Sum for Integer Shapes of the Second Distribution

When $c$ and $d$ are positive integers, the cumulative distribution function of $Y$ has a finite binomial-tail expression. Set $n=c+d-1$. Then $F_Y(x)=\sum_{j=c}^{n}\binom{n}{j}x^j(1-x)^{n-j}$. Substituting this expression into the integral and integrating each term gives a finite sum.

Result: For $a,b>0$ and positive integers $c,d$, $$ P(X>Y)= \sum_{j=c}^{c+d-1} \binom{c+d-1}{j} \frac{B(a+j,\ b+c+d-1-j)}{B(a,b)}. $$ The finite-sum identity requires integer $c,d$, not integer $a,b$. If all four shape parameters $a,b,c,d$ are positive integers, every Beta function in the sum can be evaluated using factorials. If $a$ or $b$ is noninteger, the sum still applies, but its Beta-function values require the Gamma function or another appropriate evaluation method.

A short justification of the cumulative-distribution identity also verifies its indices. Differentiate the proposed sum. With $n=c+d-1$, use $\frac{d}{dx}\{\binom{n}{j}x^j(1-x)^{n-j}\} =n\{\binom{n-1}{j-1}x^{j-1}(1-x)^{n-j}-\binom{n-1}{j}x^j(1-x)^{n-j-1}\}$. Summing from $j=c$ to $n$ cancels consecutive terms, leaving

n\binom{n-1}{c-1}x^{c-1}(1-x)^{n-c} = \frac{x^{c-1}(1-x)^{d-1}}{B(c,d)}.

The right side is the Beta$(c,d)$ density, since $1/B(c,d)=n!/((c-1)!(d-1)!)$. The proposed sum is zero at $x=0$, so integrating its derivative from zero to $x$ proves that it equals $F_Y(x)$. Finally, each term in the probability integral is evaluated by the defining Beta integral: $\int_0^1 x^{a+j-1}(1-x)^{b+n-j-1}\,dx=B(a+j,b+n-j)$. This proves the finite-sum result.

Worked Example: A Simple Exact Comparison

For a hypothetical comparison, let $X\sim\operatorname{Beta}(2,1)$ and $Y\sim\operatorname{Beta}(1,2)$, independently. All shapes are positive integers, so the finite sum applies and factorial evaluations are valid. Here $c=1$, $d=2$, and $c+d-1=2$.

The two terms are

P(X>Y) =2\frac{B(3,2)}{B(2,1)} +\frac{B(4,1)}{B(2,1)} =2\frac{1/12}{1/2}+\frac{1/4}{1/2} =\frac{1}{3}+\frac{1}{2} =\frac{5}{6}\approx 0.8333.

Thus, under these specified independent distributions, the probability that the draw for the first arm exceeds the draw for the second is $5/6$. As a direct check, $f_X(x)=2x$ and $F_Y(x)=2x-x^2$; integrating $2x(2x-x^2)$ from zero to one gives $4/3-1/2=5/6$, the same result.

Worked Example: A Longer Integer-Shape Sum

Now let independent hypothetical posterior draws be $X\sim\operatorname{Beta}(3,2)$ and $Y\sim\operatorname{Beta}(2,3)$. All four shape parameters are positive integers. In the sum, $c=2$, $d=3$, and $c+d-1=4$, so the terms correspond to $j=2,3,4$.

Using $B(3,2)=1/12$, the calculation is

\begin{aligned} P(X>Y) &=6\frac{B(5,4)}{B(3,2)} +4\frac{B(6,3)}{B(3,2)} +\frac{B(7,2)}{B(3,2)}\\ &=6\frac{1/280}{1/12} +4\frac{1/168}{1/12} +\frac{1/56}{1/12}\\ &=\frac{9}{35}+\frac{2}{7}+\frac{3}{14} =\frac{53}{70}\approx 0.7571. \end{aligned}

For example, $B(5,4)=4!3!/8!=1/280$, $B(6,3)=5!2!/8!=1/168$, and $B(7,2)=6!1!/8!=1/56$. The result is a probability between zero and one. It need not be $1/2$ merely because the two shape pairs reverse each other: the distributions are not identical, and symmetry of the pair of arm draws does not follow from reversing the shape parameters.

Noninteger Shapes and Appropriate Evaluation

The integral representation is valid for all positive real $a,b,c,d$. The particular binomial finite sum above is useful when $c,d$ are positive integers, but factorial simplification is a separate issue. The identity $B(r,s)=(r-1)!(s-1)!/(r+s-1)!$ applies only when both arguments are positive integers. In general, $B(r,s)=\Gamma(r)\Gamma(s)/\Gamma(r+s)$, where $\Gamma$ is the Gamma function.

Worked Example: A Fractional Shape Parameter

Let $X\sim\operatorname{Beta}(1/2,1)$ and let $Y\sim\operatorname{Beta}(1,1)$, independently. The second distribution is uniform, so $c=d=1$ and its cumulative distribution function is $F_Y(x)=x$. The finite sum has one term, but the first shape parameter of $X$ is not an integer; use Beta-function values rather than a factorial formula.

P(X>Y) =\frac{B(3/2,1)}{B(1/2,1)} =\frac{2/3}{2} =\frac{1}{3}\approx 0.3333.

The values follow from $B(r,1)=1/r$, or equivalently from the Gamma-function definition. Thus the finite-sum calculation remains valid, while the integer factorial shortcut does not apply to $B(1/2,1)$ or $B(3/2,1)$. As another check, the uniform draw is below $X=x$ with probability $x$, so $P(X>Y)=E[X]=1/3$.

Checks, Alternatives, and Limitations

A reliable calculation should be checked against the defining probability and the conditions of the formula. Useful checks include:

Bounds: the computed value must lie in $[0,1]$. A value outside this range signals an algebra, coding, or numerical error.
Reversed ordering: with independent continuous $X,Y$, $P(X>Y)+P(Y>X)=1$, because ties have probability zero. A separately calculated reverse probability should complement the first.
Identical distributions: if $X$ and $Y$ are independent and identically distributed continuous variables, exchanging their labels leaves the comparison unchanged, so each ordering has probability $1/2$.
Indices and normalization: for the finite sum, verify that the index runs from $c$ through $c+d-1$, and that the Beta denominator is $B(a,b)$. A missing coefficient or reversed CDF changes the result.

The finite sum is exact in mathematical form and is often convenient for integer shapes. Direct numerical quadrature of the integral is a general alternative when the shapes are noninteger or the finite sum is inconvenient. Simulation from the two distributions is another way to approximate the probability by the fraction of paired draws with $X>Y$; unlike the exact integral or finite sum, that approximation has simulation variability. These approaches target the same probability but differ in computational cost and numerical or sampling error.

The finite-sum expression should not be used with noninteger $c$ or $d$ as though its binomial-tail derivation still held. Nor should a factorial formula be applied just because the original distributions look like Beta distributions: factorial evaluation of every term requires all four shape parameters to be positive integers. For noninteger shapes, retain the integral or use the finite sum only when its required $c,d$ are integers, evaluating each Beta function through Gamma functions or a trusted numerical routine.

A further interpretive limitation is that $P(X>Y)$ is a posterior probability conditional on the chosen model, prior, and observed data. It is not the probability that a null hypothesis is true, and it is not by itself a frequentist error rate or a guarantee about a future patient’s outcome. Its meaning depends on the response-rate model and on the independence assumption used for the comparison.

Common Pitfall

Common mistake: “All the Beta-function terms can be evaluated with factorials because the sum is a finite sum.” Finiteness of the sum does not make its Beta-function arguments integers. The factorial expression applies to a term $B(r,s)$ only when both $r$ and $s$ are positive integers. The finite-sum identity needs positive integer $c,d$; factorial evaluation of every term additionally needs positive integer $a,b$. Otherwise, evaluate the relevant Beta functions with the Gamma function or another appropriate method.

Another frequent error is to integrate $f_X(x)\{1-F_Y(x)\}$ when the target is $P(X>Y)$. The event $X>Y$, conditional on $X=x$, requires $Y<x$, so the integrand uses $F_Y(x)$. The complementary CDF instead corresponds to $P(Y>x\mid X=x)$ and therefore reverses the comparison.

Check Your Understanding

Answer each question, then compare with the solutions below.

For independent $X$ and $Y$, which cumulative distribution function appears in the integral for $P(X>Y)$, and why?
Suppose $a=1.5$, $b=2$, $c=2$, and $d=3$. Does the integer-shape finite sum apply? Can every Beta term be evaluated with factorials?
If $X$ and $Y$ are independent and identically distributed continuous variables, what are $P(X>Y)$ and $P(Y>X)$? State the reasoning.
For independent $X\sim\operatorname{Beta}(1,1)$ and $Y\sim\operatorname{Beta}(1,1)$, evaluate $P(X>Y)$ using the integral.

Solutions. 1. Use $F_Y(x)$, since conditional on $X=x$, the event $X>Y$ is $Y<x$. 2. Yes: $c=2,d=3$ are positive integers, so the finite sum applies. No: $a=1.5$ makes some Beta-function arguments noninteger, so factorial evaluation of every term is not valid; use Gamma functions or numerical Beta-function evaluation. 3. Each probability is $1/2$. Identical distributions make the two orderings equally likely, and continuity gives zero probability of a tie; hence their probabilities sum to one and are equal. 4. Both densities are uniform, so $f_X(x)=1$ and $F_Y(x)=x$; therefore $P(X>Y)=\int_0^1x\,dx=1/2$.