Evaluate the Beta(3,1) and Beta(91,11) probabilities In here there is one diagram

And a statement:

After observing 90 positive reviews and 10 negative reviews, our posterior estimate on $θ_A$ has a $beta(91, 11)$ distribution. After observing 2 positive reviews, our posterior estimate on $θ_B$ has a $beta(3, 1)$ distribution. The probability that a sample from $θ_A$ is bigger than a sample from $θ_B$ is 0.713.

Now I am trying to understand how to get this 0.713?
What I know for Beta distribution the first moment (mean) is:
$E(X) = \Large \frac{a}{a + b}$
I can find the mean for the $\mu_{A} = \frac{91}{102} =0.892$ (red)
I can find the mean for the $\mu_{B} = \frac{3}{4} =0.75 $ (blue)

PDF:
$f(x)=\frac{1}{B(a, b)} x^{a-1}(1-x)^{b-1}, \quad x \in(0,1)$
Probability density function PDF(A):
$\left\{\begin{array}{ll} \frac{1}{B(91,11)} x^{90}(1-x)^{10} & 0<x<1 \\ 0 & \text { (otherwise) }\end{array}\right.$
Probability density function PDF(B):
$\left\{\begin{array}{ll}3 x^{2} & 0<x<1 \\ 0 & \text { (otherwise) }\end{array}\right.$

I can approximate the blue curve with the line $y=3x$, and the average $y$ value will be 1.5
The red curve I can possible approximate with Normal distribution, with the mean of 0.892, and variance 0.03056 because of
$ Beta(a, b)$ as $Normal \left(\frac{a}{a+b}, \sqrt{\frac{ab}{\left(a+b\right)^{2}\left(a+b+1\right)}}\right)$
Now the peak of the normal distribution should be $\frac{1}{\sigma\sqrt{2\pi}}$ is ~13. If I create a triangle from normal distribution with the height 13 and bottom $2*\sigma$, this means the average $y=6.5$.
Again I cannot find 0.713.
Any idea how John Cook found it?
 A: If $\theta_A$ has the PDF $\frac{101!}{90!10!} x^{90} (1-x)^{10}$ for $x \in [0,1]$, and $\theta_B$ has the PDF $3x^2$ for $x \in [0,1]$, and they're independent, then their joint PDF is
$$
     \frac{3 \cdot 101!}{90! 10!} x^{90} (1-x)^{10} y^2
$$
for $(x,y) \in [0,1]^2$. To find the probability $\theta_A > \theta_B$, integrate over the region where $x > y$:
$$
   \int_0^1 \int_0^x \frac{3 \cdot 101!}{90! 10!} x^{90} (1-x)^{10} y^2 \,dy\,dx = \frac{4991}{7004} \approx 0.713.
$$
We can also approximate by an easier integral: since $\theta_A$ has expected value $\frac{91}{102}$, and is pretty concentrated around that value, we could estimate $\Pr[\theta_A > \theta_B]$ as $\Pr[\frac{91}{102} > \theta_B] = \int_0^{91/102} 3x^2\,dx$. This gives us $(\frac{91}{102})^3 \approx 0.710$, which is pretty close.
A: The variables $\theta_A$ and $\theta_B$ are the posterior estimates for the proportion of positive reviews.  If these are independent (which is assumed to be the case), then $$\Pr[\theta_A > \theta_B] = \int_{\theta = 0}^1 \Pr[\theta_A > \theta]f_{\theta_B}(\theta) \, d\theta = \int_{\theta = 0}^1 \int_{\psi = \theta}^1 f_{\theta_A}(\psi) f_{\theta_B}(\theta) \, d\psi \, d\theta,$$ where $f_{\theta_A}$ and $f_{\theta_B}$ are the respective posterior densities, which we are told are beta distributed.  Thus
$$\Pr[\theta_A > \theta_B] = \frac{\Gamma(91+11)}{\Gamma(91)\Gamma(11)}\frac{\Gamma(3+1)}{\Gamma(3)\Gamma(1)} \int_{\theta=0}^1 \int_{\psi = \theta}^1 \psi^{90} (1-\psi)^{10} \theta^2 \, d\psi \, d\theta.$$
