How do you compute the P(x>y) for a joint density function in R?

Question

I'm trying to understand the Bayesian AB testing process more thoroughly. If I have two tests such that the posteriors are: $$x\sim Beta(\alpha_1, \beta_1)$$ $$y\sim Beta(\alpha_2, \beta_2)$$

Where $x$ and $y$ are posterior distributions. Then the joint pdf is of the form:

$f(x,y) = \frac{\gamma(\alpha_1+\beta_1)}{\gamma(\alpha_1)\gamma(\beta_1)}x^{\alpha_1-1}(1-x)^{\beta_1-1}\frac{\gamma(\alpha_2+\beta_2)}{\gamma(\alpha_2)\gamma(\beta_2)}y^{\alpha_2-1}(1-y)^{\beta_2-1}\quad 0\leq x \leq 1,\quad 0\leq y \leq 1$

I want to know $P(x>y)$, how do u solve the equation

$$ P(x>y) = \int_{0}^{1}\int_{y}^{1}\frac{\gamma(\alpha_1+\beta_1)}{\gamma(\alpha_1)\gamma(\beta_1)}x^{\alpha_1-1}(1-x)^{\beta_1-1}\frac{\gamma(\alpha_2+\beta_2)}{\gamma(\alpha_2)\gamma(\beta_2)}y^{\alpha_2-1}(1-y)^{\beta_2-1}dxdy$$

in R? I have the joint function coded like so

joint_dist <- function(x, y, aX, bX, aY, bY) {
  (gamma(aX + bX)/gamma(aX)*gamma(bX)*x^(aX - 1)*(1-x)^(bX-1))*(gamma(aY + 
   bY)/gamma(aY)*gamma(bY)*y^(aY-1)*(1-y)^(bY-1))
  }

How do I use the integrate function to integrate from y to 1 and then 0 to 1? Also I'm aware of Evan Miller's summation derivation, but I'm under the impression that he uses a $Beta(1,1)$ prior for both control and variation. I'd like to be able to choose a different prior for each. Maybe an uniformative prior for the variation, but a more informed prior for the control. Maybe I'm wrong here.

Answer 1

How about using MonteCarlo integration?

#You can put your value here. 
aX=3 ;bX=7; aY=3; bY=7

#MC Integration
n_sim=1000000   # number of simulation
u<-runif(n_sim) # 
C<-(1-pbeta(u,aX,bX))*dbeta(u,aY,bY)
m<-mean(C)  #integration result
s<-sd(C)    # standard deviation of the result.

You can figure out confidence interval of the calculation using the standard deviation. For example, the 95% CI would be (m-s/sqrt(n_sim)*1.96,m+s*1.96/sqrt(n_sim)).

When I set X and Y have same distribution, I got 50% as a result.

There could be way better way to calculate and maybe someone would post it as another answer.

I am not a native English speaker so any revision of this answer would be appreciated. Thank you in advance!

edit1 I would add some more detail.

Let's say

$1-F_{\alpha_1\beta_1}(y)=\int_{y}^{1}\frac{\gamma(\alpha_1+\beta_1)}{\gamma(\alpha_1)\gamma(\beta_1)}x^{\alpha_1-1}(1-x)^{\beta_1-1}dx$

and

$f_{\alpha_2\beta_2}(y)=\frac{\gamma(\alpha_2+\beta_2)}{\gamma(\alpha_2)\gamma(\beta_2)}y^{\alpha_2-1}(1-y)^{\beta_2-1}$

then, $P(x>y) =\int_{0}^{1} (1-F_{\alpha_1\beta_1}(y))f_{\alpha_2\beta_2}(y) dy$.

And this expression is actually getting expectation $E[(1-F_{\alpha_1\beta_1}(y))f_{\alpha_2\beta_2}(y)]$ when y has a uniform distribution.

So what I did is essentially simulating this expectation. I took 1000000 (n_sim) samples from the uniform distribution and get $(1-F_{\alpha_1\beta_1}(y))f_{\alpha_2\beta_2}(y)$ for every each sample (which is C) and calculated the mean and standard deviation. This mean is unbiased estimator of the integration. Thanks to the central limit theorem, the mean of C has a asymptotic normal distribution and this is the reason we can get confidence interval as well.

For addition, I guess you can interpret the integration as $E[(1-F_{\alpha_1\beta_1}(y))]$ when y has a beta distribution with $\alpha_2$, $\beta_2$.

y<-rbeta(n_sim,aY,bY)
C2<-(1-pbeta(y,aX,bX))
mean(C2)

This would give the same mean. But in this case standard deviation seems different. I haven't figured out why it is so I have to do some more research about this difference.

However, the method getting the confidence interval using the uniform distribution has the reference.