# Beta inequalities in R

Someone asked me yesterday for R code to compute the probability P(X > Y + δ) where X and Y are independent beta random variables. I’m posting the solution here in case it benefits anyone else.

For an example of why you might want to compute this probability, see A Bayesian view of Amazon resellers.

Let X be a Beta(a, b) random variable and Y be a Beta(c, d) random variable. Denote PDFs by f and CDFs by F. Then the probability we need is

If you just need to compute this probability a few times, here is a desktop application to compute random inequalities.

But if you need to do this computation repeated inside R code, you could use the following.

```beta.ineq <- function(a, b, c, d, delta)
{
integrand <- function(x) { dbeta(x, a, b)*pbeta(x-delta, c, d) }
integrate(integrand, delta, 1, rel.tol=1e-4)\$value
}```

The code is as good or as bad as R’s `integrate` function. It’s probably accurate enough as long as none of the parameters a, b, c, or d are near zero. When one or more of these parameters is small, the integral is harder to compute numerically.

There is no error checking in the code above. A more robust version would verify that all parameters are positive and that `delta` is less than 1.

Here’s the solution to the corresponding problem for gamma random variables, provided `delta` is zero: A support one-liner.

And here is a series of blog posts on random inequalities.