Expected Loss Calculation 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
 A: 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)
mean(pmax(theta_a-theta_b, 0))

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

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...
