Let $X,Y$ be independent random variables such that $X \sim \operatorname{Gamma}(a,b)$ and $Y\sim \operatorname{Gamma}(c,b)$. We denote $M = E_{X,Y}\left[\frac{X}{X+Y}\right]$.
We know that in the above framework $\frac{X}{X+Y}$ is beta distributed, and hence $M$ is the mean of a beta random variable.
Now I need to obtain $$M1 = E_{X,Y}\left[\frac{X^2}{X+Y}\right].$$
Is is correct to claim that $$M1 = \frac{c_{a}}{c_{a-1}}E_{X,Y}\left[\frac{X}{X+Y}\right],$$ where $X \sim \operatorname{Gamma}(a+1,b)$, and $Y \sim \operatorname{Gamma}(c,b)$, and $c_a$ and $c_{a-1}$ are normalization constants for the PDFs of $X$ with parameters $a$ and $a-1$, respectively?
Any other solutions to $M1$ are welcome.
For independent $X \sim \operatorname{Gamma}(a,b)$ and $Y \sim \operatorname{Gamma}(c,b)$ $$\tag{1}\label{1} E_{X,Y}\left[\frac{X^2}{X+Y}\right] = \int_0^\infty \int_0^\infty \frac{x^2}{x+y} \frac{b^a}{\Gamma(a)}x^{a-1}e^{-bx}\frac{b^c}{\Gamma(c)}y^{c-1}e^{-by}\,dx\,dy $$ For independent $U \sim \operatorname{Gamma}(a+1,b)$ and $V \sim \operatorname{Gamma}(c,b)$ $$\tag{2}\label{2} E_{U,V}\left[\frac{U}{U+V}\right] = \int_0^\infty \int_0^\infty \frac{x}{x+y} \frac{b^{a+1}}{\Gamma(a+1)}x^{a}e^{-bx}\frac{b^c}{\Gamma(c)}y^{c-1}e^{-by}\,dx\,dy $$ One can see that r.h.s.'s in (\ref{1}) and (\ref{2}) are differ by constant only: $$ E_{U,V}\left[\frac{U}{U+V}\right] =\frac{b}{a}E_{X,Y}\left[\frac{X^2}{X+Y}\right]. $$ Then $$ M1=E_{X,Y}\left[\frac{X^2}{X+Y}\right]=\frac{a}{b}E_{U,V}\left[\frac{U}{U+V}\right]=\frac{a}{b}\cdot\frac{a+1}{a+1+c} $$ where $U \sim \operatorname{Gamma}(a+1,b)$ and $V \sim \operatorname{Gamma}(c,b)$ are independent r.v.