Statistical moments of $\frac X{X + Y}$ when $X$ and $Y$ are two independent random variables with a Beta distribution I'm trying to find the moments (or the pdf but I'm less confident there's a closed form) of $\frac X{X + Y}$ where $X$ and $Y$ are two independent random variables with a Beta distribution. There's a paper from Pham-Gia that I tried to read, and a similar (but yet different) question posted here, but they didn't help.  I find the resulting distribution to be a lot similar to the Beta distribution.
 A: Suppose that $X$ and $Y$ are independent beta random variables (r.v.'s) with parameters $(a,b)$ and $(c,d)$, respectively. Let 
\begin{equation*}
 V:=\frac X{X+Y}. \tag{0}
\end{equation*}
The transformation $(x,v)\mapsto(x,x\frac{1-v}v)$ transforms $(X,V)$ to $(X,Y)$. 
The Jacobian determinant of this transformation is $-x/v^2$. Also,
\begin{equation*}
\Big(x,x\frac{1-v}v\Big)\in(0,1)\times(0,1)
\iff
 (0 < v<1\ \&\ 0 < x < s),
\end{equation*}
where 
\begin{equation*}
 s:=1\wedge\frac1r=\min\Big(1,\frac1r\Big),\quad r:=r_v:=\frac{1-v}v. 
\end{equation*}
The joint pdf of $(X,Y)$ is given by 
\begin{equation*}
 f_{X,Y}(x,y)=Cx^{a-1}(1-x)^{b-1}y^{c-1}(1-y)^{d-1}\,I\{0<x<1\ \&\ 0<y<1\}, 
\end{equation*}
where $I$ denotes the indicator and 
\begin{equation*}
 C:=\frac1{B(a,b)B(c,d)}. 
\end{equation*}
So, joint pdf of $(X,V)$ is given by 
\begin{align*}
 f_{X,V}(x,v)&=f_{X,Y}(x,xr)x/v^2 \\ 
 &=
 Cx^{a+c-1}(1-x)^{b-1}v^{-c-1}(1-v)^{c-1}(1-rx)^{d-1} \\ 
 &\times I\{0 < v<1\ \&\ 0 < x < s\}. 
\end{align*}
So, the pdf of $V$ is given by 
\begin{align*}
 f_V(v)&=\int_{\mathbb R}f_{X,V}(x,v)\,dx \\ 
& =Cv^{-c-1}(1-v)^{c-1}J(v)I\{0 < v<1\}, \tag{1}
\end{align*}
where 
\begin{equation*}
 J(v):=\int_0^s x^{a+c-1}(1-x)^{b-1}(1-rx)^{d-1}\,dx. 
\end{equation*}
To evaluate $J(v)$, use formula 3.197.3 of Table of Integrals, Series, and Products,
Seventh Edition by Gradshteyn and Ryzhik: 
\begin{equation*}
 \int_0^1 x^{\lambda-1}(1-x)^{\mu-1}(1-\beta x)^{-\nu}\,dx
 =B(\lambda,\mu)\,_2F_1(\nu,\lambda;\lambda+\mu;\beta) \tag{*}
\end{equation*}
for $\lambda>0$, $\mu>0$, $|\beta|<1$, where $_2F_1$ is the hypergeometric function. 
If $1/2<v<1$, then $0<r<1$, $s=1$, and (*) immediately yields 
\begin{equation*}
 J(v)=B(a+c,b)\,_2F_1(1-d,a+c;a+c+b;r). \tag{2}
\end{equation*}
If $0<v<1/2$, then $r>1$, $s=1/r\in(0,1)$, and the substitution $t=rx$ together with (*) 
yields 
\begin{align*}
 J(v)&=r^{-a-c}\int_0^1 t^{a+c-1}(1-t/r)^{b-1}(1-t)^{d-1}\,dx \\ 
 &=r^{-a-c}B(a+c,d)\,_2F_1(1-b,a+c;a+c+d;1/r). \tag{3}
\end{align*}
Formulas (1), (2) (for $v\in(1/2,1)$), and (3) (for $v\in(0,1/2)$) give the pdf $f_V$ of $V$, defined by (0). 
A: Not a complete answer but too long for a comment:
Let $I=]0,1[$ be the unit interval. Let $$F:I^2\to M, (x,y)\mapsto \left(x,\frac{x}{x+y}\right),$$
where $$M=\{(a,b)\in I\times I\mid \frac{a}{a+1}\le b\}=\{(a,b)\in I\times I\mid a\le\frac{b}{1-b}\}.$$
(Here I have the convention that $\frac{b}{b-1}=\infty$ for $b=1$.)
Note that $F$ is a diffeomorphism with inverse
$$F^{-1}(a,b)=\left(a, \frac{a\cdot(1-b)}b\right).$$
So, where $\operatorname{Jac}$ denotes the Jacobian determinant,
$$|(\operatorname{Jac}F^{-1})(a,b)|=\frac{a}{b^2}.$$
Use the following result from probability Theory (here, $f_Z$ denotes the density of the random vector $Z$):
$$\bbox[15px,border:1px groove navy]{f_{F\circ(X,Y)}(a,b)=f_{(X,Y)}(F^{-1}(a,b))\cdot|(\operatorname{Jac}F^{-1})(a,b)|\cdot1_{\operatorname{Im}(F)}(a,b).}$$
In our case, $$f_{(X,\frac{X}{X+Y})}(a,b)=f_{(X,Y)}\left(a,\frac{a\cdot(1-b)}b\right)\cdot\frac{a}{b^2}\cdot 1_{M}(a,b).$$
Note that, since $X$ and $Y$ are independent, the joint density is obtained simply as a product of the individual densities of a Beta distribution. Say $X\sim \operatorname{Beta}(\alpha,\beta)$ and $Y\sim\operatorname{Beta}(\alpha,\beta)$. Let $$c=\operatorname{B}(\alpha,\beta)\cdot \operatorname{B}(\gamma, \delta),$$ where $\operatorname{B}$ is the Beta function.
Then
\begin{split}
f_{\frac{X}{X+Y}}(b)&=\int_\mathbb R f_{(X,\frac{X}{X+Y})}(a,b)\,\mathrm da\\
&=1_{I}(b)\cdot c^{-1} b^{\gamma-3} (1-b)^{\delta-1}\cdot\int_0^{\min(1,\frac{b}{1-b})} a^\alpha(1-a)^{\beta-1}\cdot\left(1+\frac{a\cdot(b-1)}b\right)^{\delta-1}\,\mathrm da.
\end{split}
A: The command of Mathematica 12.0
pdf=PDF[TransformedDistribution[x/(x + y), {x \[Distributed] BetaDistribution[a, b], 
y \[Distributed] BetaDistribution[c, d]}], t] // TeXForm

performs
$$ \begin{cases}
 -\frac{\pi  \left(\frac{1}{t}-1\right)^{-a} \Gamma (d) \csc (\pi  (a+c)) \,
   _2\tilde{F}_1\left(1-b,a+c;a+c+d;\frac{t}{1-t}\right)}{(t-1) t B(a,b) \Gamma (-a-c+1)
   B(c,d)} & 0<t\leq \frac{1}{2} \\
 \frac{\pi  (-1)^{1-a} \Gamma (b) \left(\frac{1}{t-1}\right)^{1-c} t^{-c-1} (\cot (\pi
   (a+c))+i) \, _2\tilde{F}_1\left(a+c,1-d;a+b+c;\frac{1}{t}-1\right)}{B(a,b) \Gamma
   (-a-c+1) B(c,d)} & \frac{1}{2}<t<1
\end{cases},$$ where $ _2\tilde{F}_1 $ is the regularized hypergeometric function. Unfortunately, TeXForm omitted zero value of the PDF under consideration for $t\notin [0,1]$.
Addition. In response to @Matt F comment, let us build its plot for concrete noninteger values of the parameters:
Plot[pdf /. {a -> 1.1, b -> 1.1, c -> 1.1, d -> 1.1}, {t, -2, 2}, PlotStyle -> Thick]


