# Property of two independent Beta distribution

I have been working with beta-bernoulli posteriors recently. Is it true that if $$X,Y$$ are independent rvs with $$X \sim Beta(a_1+1,b_1+1)$$ and $$Y \sim Beta(a_2+1,b_2+1)$$ then $$\mathbb{P}(X>Y)>0.5$$ iff $$\frac{a_1}{a_1+b_1}>\frac{a_2}{a_2+b_2}$$? $$a_1,a_2,b_1,b_2$$ are assumed to be whole numbers.

This makes intuitive sense as the estimated bernoulli parameter is greater in $$X$$ than in $$Y$$. I have observed this for all parameters within the range $$0 \leq a_1,b_1,a_2,b_2\leq 6$$. However, it does not seem to immediately follow from the definitions.

I have tried writing $$\mathbb{P}(X>Y)=\frac{\int_{0}^{1}\int_{0}^{x}y^{a_2}(1-y)^{b_2}x^{a_1}(1-x)^{b_1}dy dx}{B(a_1+1,b_1+1)B(a_2+1,b_2+1)}$$, but the numerator doesn't seem that easy to work with. I don't see a similar question on stackexchange. Can this be proved or can a counter example be produced for this?

Update1: It seems the result is false if $$a_1,b_1,a_2,b_2$$ are just positive real numbers as initially posed. See Henrys answer below. Can a counterexample be produced for whole numbers?

Update2: Henry has produced a whole number counterexample. This answers the question in the negative conclusively.

• This may be relevant stats.stackexchange.com/questions/436039/… Jan 12, 2023 at 14:42
• Yes, it is relevant but the given pdfs for difference of beta rvs seem daunting to integrate and make any conclusion about. Jan 12, 2023 at 14:48

It does not seem to be a correct conjecture.

It seems your condition is that the mode for $$X$$ is greater than the mode for $$Y$$. Since in non-symmetric Beta distributions, the mode is not equal to the mean, it should be possible to find a counterexample.

One way is to take a non-symmetric case where $$\frac{a_1}{a_1+b_1}=\frac{a_2}{a_2+b_2}$$, find which tends to be greater more often, and then adjust the parameters slightly to get a counter-example which works.

So, using R, we might look at

set.seed(2023)
XgreaterY <- function(a1,b1,a2,b2,cases){
mean(rbeta(cases, a1+1, b1+1) > rbeta(cases, a2+1, b2+1))
}
XgreaterY(3, 1, 30, 10, 10^7)
# 0.390946
XgreaterY(3, 1, 29, 11, 10^7)
# 0.4385228


suggesting that $$a_1=3, b_1=1, a_2=29, b_2=11$$ provides a counterexample, since $$\frac{a_1}{a_1+b_1} = 0.75 > 0.725 = \frac{a_2}{a_2+b_2}$$ but it seems $$\mathbb{P}(X>Y) < 0.44 < 0.5$$

• Could this hold if $a_1,b_1,a_2,b_2$ are all whole numbers? I had tried it only for whole numbers. Jan 12, 2023 at 14:57
• @SushantVijayan No: try $a_1=3,b_1=1,a_2=29,b_2=11$ or XgreaterY(3, 1, 29, 11, 10^7) to get a probability less than $0.44$ Jan 12, 2023 at 15:02
• @SushantVijayan I doubt it, and a similar approach should work, perhaps reversing the adjustment. The (mode) example in the answer compares $\text{Beta}(4,2)$ and $\text{Beta}(30,12)$. A (mean) example might compare $\text{Beta}(301,149)$ and $\text{Beta}(4,2)$ Jan 12, 2023 at 15:22