How to calculate the PDF of the 'difference' between two Beta distributions? I start with two Beta distributions:
$$\mathrm{Beta_A}(p; \alpha_A, \beta_A) = \frac{p^{\alpha_A-1}\,(1-p)^{\beta_A-1}}{\mathrm{B}(\alpha_A, \beta_A)}$$
$$\mathrm{Beta_B}(p; \alpha_B, \beta_B) = \frac{p^{\alpha_B-1}\,(1-p)^{\beta_B-1}}{{\mathrm{B}(\alpha_B, \beta_B)}}$$
where in the context of Bernoulli trials, $\alpha$ can be interpreted as $1 + \mathrm{successes}$ and $\beta$ can be interpreted as $1 + \mathrm{fails}$. $\mathrm{B}$ is the Beta function.
I then define the 'difference' between $\mathrm{Beta_A}$ and $\mathrm{Beta_B}$ as:
$$F(x; \alpha_A, \beta_A, \alpha_B, \beta_B) = \mathrm{Beta_A}(p) - \mathrm{Beta_B}(p)$$
Questions


*

*what is the PDF of $F(x)$?

*what family of probability density distributions does $F(x)$ belong to?


Example and illustration
For example for $\alpha_A=3, \beta_A=9$ (2 successes from 8 Bernoulli trials) and $\alpha_A=1, \beta_A=5$ (0 successes from 4 Bernoulli trials) the distribution of values that $p$ can take is:

If I then take $n$ random values $X_A \sim \mathrm{Beta_A}$ and $X_B \sim \mathrm{Beta_B}$, and find the differences between each $i^\mathrm{th}$ element, $X_{A,i} - X_{B,i}$, and plot these $n$ differences in a histogram, I am essentially sampling $F(x)$ - the underlying distribution of $\mathrm{Beta_A} - \mathrm{Beta_B}$ which can only be defined for $x \in [-1,+1]$.
With $n = 5 \times 10^7$ random samples and bin widths of $\Delta x = 0.004$, $F(x)$ takes the following form:

What is the PDF of $F(x)$?

Notes


*

*more verbose version of question

*Kullback-Leibler divergence only gives a scalar value of difference measure

 A: I know this is a bit of an old question but for what it's worth there is an established closed-form solution to this problem, found by Pham-Gia, Turkkan, and Eng in 1993. It's a piecewise solution that relies on the Appell F1 hypergeometric function. Given
$$
\begin{align}
\theta_1 &= \text{beta}(\alpha_1, \beta_1) \\
\theta_2 &= \text{beta}(\alpha_2, \beta_2) \\
\theta_d &= \theta_0 - \theta_1
\end{align}
$$
Then the probability of the difference of $\theta_d$ is piecewise over $\theta_d$. I've re-written it here, in case you can't access the paper. I'm using $\cdot$ to indicate multiplication when I need to break the equation over several lines.
Given $A = \text{Beta}(\alpha_1, \beta_1)\text{Beta}(\alpha_2, \beta_2)$.
For $-1 \leq \theta_d < 0$:
$$
\begin{align}
p(\theta_d) = &\text{Beta}(\alpha_2, \beta_1)\theta_d^{\beta_1 + \beta_2 - 1}(1 - \theta_d)^{\alpha_2 + \beta_1 - 1} \\
&\cdot F_1(\beta_1, \alpha_1 + \beta_1 + \alpha_2 + \beta_2 - 2, 1 - \alpha_1; \beta_1 + \alpha_2; 1 - \theta_d, 1 - \theta_d^2) \\
&/ A
\end{align}
$$
For $0 < \theta_d \leq 1$:
$$
\begin{align}
p(\theta_d) = &\text{Beta}(\alpha_1, \beta_2)\theta_d^{\beta_1 + \beta_2 - 1}(1 - \theta_d)^{\alpha_1 + \beta_2 - 1} \\
&\cdot F_1(\beta_2, 1 - \alpha_2, \alpha_1 + \beta_1 + \alpha_2 + \beta_2 - 2; \alpha_1 + \beta_2; 1 - \theta_d^2, 1 + \theta_d) \\
&/ A
\end{align}
$$
And for $\theta_d = 0$, $\alpha_1 + \alpha_2 > 1$, $\beta_1 + \beta_2 > 1$:
$$
p(\theta_d) = \text{Beta}(\alpha_1 + \alpha_2 - 1, \beta_1 + \beta_2 -1) / A
$$
As for the family of this distribution, I'm honestly not sure. Somebody else may be able to jump in there.
