How do I subtract one Beta distribution from another? The context is that I want to evaluate risk mitigation efforts. 
I can get three-point (optimistic-likely-pessimistic) estimates of each of


*

*the un-mitigated risk,  

*the cost of mitigation and

*the residual risk (after mitigation)


and convert each of these to Beta PERT which is a decent approximation for a Beta distribution. 
So the value of the mitigation is the original risk minus the residual risk, 
and I'd like to subtract the mitigation cost to get a marginal value for the mitigation. 
That involves subtracting a beta distribution from another, and I can't find a recognised means to do that. 
Can you help please?
My preferred tool is a spreadsheet (because I want others in my firm to be able to do this easily), or alternatively R. 
 A: Because the R code mentioned above does not seem to be available anymore, here is one such implementation in R for the difference in two independent beta random variables.  (This is all based on the answer by @passerby51.)
# Set parameters
a1 <- 1    # alpha1
b1 <- 1    # beta1
a2 <- 20   # alpha2
b2 <- 17   # beta2

# Define integrand
integrand <- function(y, z, a1, b1, a2, b2) {
  (z+y)^(a1-1) * (1-z-y)^(b1-1) * y^(a2-1) * (1-y)^(b2-1)
}

# Construct pdf of difference of two independent beta distributions
  pdf <- NULL
  n <- 100 # Number of equally-spaced evaluations of the pdf between -1 and +1
  z <- -1 + 2*c(0:n)/n
  for (i in 1:(n/2)) {
      pdf[i] <- integrate(integrand, -z[i], 1, z=z[i], a1=a1, b1=b1, a2=a2, b2=b2)[[1]]
  }
  for (i in (n/2+1):(n+1)) {
      pdf[i] <- integrate(integrand, 0, 1-z[i], z=z[i], a1=a1, b1=b1, a2=a2, b2=b2)[[1]]
  }
  pdf <- pdf/(beta(a1, b1)*beta(a2, b2))

# Plot pdf and compare against histogram with a large random sample
  zz <- rbeta(1000000, a1, b1) - rbeta(1000000, a2, b2)
  hist(zz, breaks=100, freq=FALSE, las=1, xlim=c(-1, 1))
  lines(z, pdf, col="red")


Note that there are closed-form solutions for the pdf but only if the parameters are all positive integers.
A: Let $X$ and $Y$ be independent variables with densities $f_X$ and $f_Y$ and let $Z = X-Y$. Then, $Z$ will have the following density
$$
f_Z(z) = \int f_X(z+y) f_Y(y) dy
$$
which is a form of convolution of the two densities. If $X$ and $Y$ have Beta distributions with parameters $(\alpha,\beta)$ and $(a,b)$, then
\begin{align*}
f_Z(z) &= \frac{1}{B(\alpha,\beta)B(a,b)} \int_0^1 (z+y)^{\alpha-1} (1-z-y)^{\beta-1} 1\{z+y \in(0,1)\} y^{a-1} (1-y)^{b-1} dy \\
&= \frac{1}{B(\alpha,\beta)B(a,b)} \int_{\max\{0,-z\}}^{\min\{1,1-z\}} (z+y)^{\alpha-1} (1-z-y)^{\beta-1} y^{a-1} (1-y)^{b-1} dy.
\end{align*}
The density will be supported on $z \in [-1,1]$. This perhaps can be numerically integrated.
The easiest way to get estimates is perhaps by Monte Carlo simulation: Generate independent draws $(X_i,Y_i), i=1,\dots,n$ and form $Z_i = X_i - Y_i$. Then we can look at the histogram of $Z_i$ or obtain kernel density estimates. Here is some R code.
