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

P(X > Y + \delta) &=& \int_\delta^1 \int_0^{x-\delta} f_X(x) \, f_Y(y)\, dy,dx \ &=& \int_\delta^1 f_X(x)\, F_Y(x-\delta) \, dx

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.

Click to learn more about Bayesian statistics consulting


One thought on “Beta inequalities in R

  1. Hi John,
    Very much appreciate to your vignettes on R.V. inequalities. They are really helpful to my research on Bayesian response adaptive randomization.
    I’ve tried an extended version solving P(X > max(Y, Z)), where X, Y, Z are beta R.V.s according to your publication(Stochastic inequality probabilities for adaptively randomized clinical trials. Cook JD, Nadarajah S).
    three_brv_max_prob_solver <- function(a1, b1, a2, b2, a3, b3) {
    ibeta <- function(x, a, b) {
    pbeta(x, a, b)
    integrand <- function(x) {
    x ^ (a1 – 1) * (1 – x) ^ (b1 – 1) * ibeta(x, a2, b2) * ibeta(x, a3, b3)
    result <- 1 / beta(a1, b1) * integrate(integrand, 0, 1, rel.tol = 1e-4)$value
    Hope it can help others if I didn't make any mistakes.

Leave a Reply

Your email address will not be published. Required fields are marked *