Question: How to prove this result?

A classical probability result says that if G₁ and G₂ are two independent Gamma random variables with same scale parameter (let's say 1 to simplify) and shape parameters a₁ and a₂ respectively, then G_k / (G₁ + G₂ ) is a Beta random variable with parameters (a_k , a_3-k ) (k=1..2).

In the attached file it is shown that (Maple 2015) function Statistics:-PDF fails in computing the PDF of G_k.
Noting strange here if you observe that even in the extremely simple case Z = X / (X+Y), where both X and Y are independent Uniform random variable with support [0, 1), Maple 2015 already fails in computing PDF(Z).

An alternative to Statistics:-PDF is to write explicitely the double integration which defines CDF(Z) (to begin with, and later PDF(G_k)) and ask Maple to do the integrations.
This approach works for Z but requires helping Maple when X and Y are still independent Uniform random variables but with respective non instanciated supports [0, a₁) and [0, a₂).

Applying to the Gamma-Gamma case the recipies I introduced in the Uniform-Uniform case does not give any result, unless in the very particular case where the shape parameters a₁ and a₂ are (strictly) positive integers.

All the details are in X_over_(X_plus_Y).mw

Do you have any idea how to prove with Maple the probability result mentioned at the head of this question?

PS: The "classical method" to do compute PDF(Z) consists in changing the integration variables < x₁, x₂ > into < x₁ = v₁v₂, x₂ = v₂ (1-v₁) > (see for instance Stack exchange)... but even after having dome it I still cannot get the desired result.

Thanks in advance.