How do you solve the integral $$E(L(\theta_A,\theta_B)) = \int_0^1\int_{\theta_B}^1(\theta_A - \theta_B)f(\theta_A)f(\theta_B)d\theta_Ad\theta_B$$ where $\theta_A \sim Beta(\alpha_1, \beta_1)$ and $\theta_B \sim Beta(\alpha_2,\beta_2)$. Essentially looking at the loss function $$ L(\theta_A,\theta_B) = \begin{cases} \theta_A - \theta_B &\quad\theta_A\ge\theta_B \\ 0 &\quad\theta_A<\theta_B \end{cases} $$ If you have examples in R that'd be greatly appreciated! Thanks

  • 1
    $\begingroup$ Are you asking for a closed-form solution (not likely) or R code? $\endgroup$ – jbowman Sep 21 '18 at 1:35
  • $\begingroup$ Yeah probably not a closed form solution. Looking for a way to at least estimate it numerically whether in R or python. $\endgroup$ – jjt3 Sep 21 '18 at 15:43

Here are two methods, a nested numerical integration and a Monte Carlo integration, in R. The nested numerical integration has an outer integration that goes from 0 to 1 with respect to $\theta_B$ and an inner integration that goes from $\theta_B$ to 1 with respect to $\theta_A$; note that since the integrate function in R is vectorized, we have to un-vectorize the outer integration in order to call the inner integration (no doubt better programming of the bar function would overcome this.)

The parameters of the integrate function are the function to be integrated, the lower and upper bounds, and various parameters that get passed to the function to be integrated. integrate(bar, theta_b[j], 1, ...) means "integrate the function bar from theta_b[j] to 1 passing the list of variables in ... to bar."

bar <- function(theta_a, theta_b, a1, b1) {
   (theta_a - theta_b) * dbeta(theta_a, a1, b1)

foo <- function(theta_b, a1, b1, a2, b2) {
   retval <- rep(0, length(theta_b))
   for (j in 1:length(theta_b)) {
      retval[j] <- integrate(bar, theta_b[j], 1, theta_b=theta_b[j], a1=a1, b1=b1)$value
   retval * dbeta(theta_b, a2, b2)

a1 <- 5
b1 <- 4
a2 <- 7
b2 <- 3

# The integration itself:
integrate(foo, 0, 1, a1=a1, b1=b1, a2=a2, b2=b2)

# Here's what we get:
0.03124686 with absolute error < 3.5e-16

# Check via simulation
theta_a <- rbeta(1e7, a1, b1)
theta_b <- rbeta(1e7, a2, b2)

#  Calculate the average of max(theta_a - theta_b, 0)

# Here's what we get:
[1] 0.03125562

As you can see, the two results are the same to three decimal places, and off by only a little in the fourth.

Naturally this can return results that may be problematic with extreme parameterizations, for example:

> a1 <- 2
> b1 <- 40
> a2 <- 50
> b2 <- 3

> integrate(foo, 0, 1, a1=a1, b1=b1, a2=a2, b2=b2)
1.167158e-24 with absolute error < 4.8e-25
> # Check via simulation
> theta_a <- rbeta(1e7, a1, b1)
> theta_b <- rbeta(1e7, a2, b2)
> mean(pmax(theta_a-theta_b))
[1] 0

as even with $10^7$ simulations, none of the $\theta_A$ values were $\geq$ the corresponding $\theta_B$ value. In this case, the numerical integration result probably isn't bad, as the correct answer would appear to be "very, very close to 0", but still...

  • 1
    $\begingroup$ (+1) These integrals can be expressed in terms of generalized hypergeometric functions. For integral parameters they have rational values (the integrands are polynomial, after all). For instance, the first one evaluates to $1243/39780\approx 0.03124685771744595$ and the second one to $547/468597582252207969843412746 \approx 1.16731\times 10^{-24}.$ $\endgroup$ – whuber Sep 21 '18 at 18:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.