I'm trying to find the moments (or the pdf but I'm less confident there's a closed form) of $X /(X + Y)$ where X and Y are two independent Beta distributions. There's a paper from PhamGia 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.
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{1v}v)$ transforms $(X,V)$ to $(X,Y)$. The Jacobian determinant of this transformation is $x/v^2$. Also, \begin{equation*} \Big(x,x\frac{1v}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{1v}v. \end{equation*} The joint pdf of $(X,Y)$ is given by \begin{equation*} f_{X,Y}(x,y)=Cx^{a1}(1x)^{b1}y^{c1}(1y)^{d1}\,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+c1}(1x)^{b1}v^{c1}(1v)^{c1}(1rx)^{d1} \\ &\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^{c1}(1v)^{c1}J(v)I\{0 < v<1\}, \tag{1} \end{align*} where \begin{equation*} J(v):=\int_0^s x^{a+c1}(1x)^{b1}(1rx)^{d1}\,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^{\lambda1}(1x)^{\mu1}(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(1d,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^{ac}\int_0^1 t^{a+c1}(1t/r)^{b1}(1t)^{d1}\,dx \\ &=r^{ac}B(a+c,d)\,_2F_1(1b,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).

$\begingroup$ The second case has a complex part after the $\cot$? $\endgroup$ – Matt F. Nov 19 at 21:01

$\begingroup$ @MattF. : Thank you for the alert. I have redone the entire answer. (The problem was that I previously used a somewhat incorrect expression from Maximilian Janisch's answer.) $\endgroup$ – Iosif Pinelis Nov 20 at 4:30

$\begingroup$ @IosifPinelis Sorry for causing that extra work $\endgroup$ – Maximilian Janisch Nov 20 at 8:32

$\begingroup$ @MaximilianJanisch : Nothing to be sorry about. I had to be more careful. $\endgroup$ – Iosif Pinelis Nov 20 at 13:00

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}{1b}\}.$$
(Here I have the convention that $\frac{b}{b1}=\infty$ for $b=1$.)
Note that $F$ is a diffeomorphism with inverse
$$F^{1}(a,b)=\left(a, \frac{a\cdot(1b)}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$):
$$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(1b)}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^{\gamma3} (1b)^{\delta1}\cdot\int_0^{\min(1,\frac{b}{1b})} a^\alpha(1a)^{\beta1}\cdot\left(1+\frac{a\cdot(b1)}b\right)^{\delta1}\,\mathrm da. \end{split}

$\begingroup$ thank you very much for your thorough answer. I came to the same result as you do, but have trouble reducing this expression further (as I feel it could be done, but maybe I'm mistaken) and/or finding/approximating moments. $\endgroup$ – Kerighan Nov 19 at 16:53

$\begingroup$ This would be easier to read replacing $]0,1[$ with $I$, removing most of the $\cdot$'s, replacing $\operatorname{Bet}$ with the more standard $\mathrm{B}$, and replacing $const.$ with the more standard $c$. $\endgroup$ – Matt F. Nov 19 at 17:15

$\begingroup$ Hi @Kerighan I agree that this doesn't (fully) answer your question, but I think I'll leave my answer here as it might be useful to other users... As for approximating moments: I suppose that you could approximate them using numerical integration for fixed $\alpha,\beta,\gamma,\delta$. I don't have a good idea if you want an answer for every $\alpha,\beta,\gamma,\delta$ $\endgroup$ – Maximilian Janisch Nov 19 at 17:16

1$\begingroup$ @MattF. I implemented some of your suggestions $\endgroup$ – Maximilian Janisch Nov 19 at 17:21

1$\begingroup$ @MaximilianJanisch : I think your indicator for the joint pdf of $(X,\frac X{X+Y})$ is incorrect, which leads to having $1+\frac{a\cdot(b1)}b<0$ under your final integral when $b<1/2$ and $a$ is close to $1$. $\endgroup$ – Iosif Pinelis Nov 20 at 4:35
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(1b,a+c;a+c+d;\frac{t}{1t}\right)}{(t1) t B(a,b) \Gamma (ac+1)
B(c,d)} & 0<t\leq \frac{1}{2} \\
\frac{\pi (1)^{1a} \Gamma (b) \left(\frac{1}{t1}\right)^{1c} t^{c1} (\cot (\pi
(a+c))+i) \, _2\tilde{F}_1\left(a+c,1d;a+b+c;\frac{1}{t}1\right)}{B(a,b) \Gamma
(ac+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]

$\begingroup$ The second case is wrong, or at least requires clarification, because of the complex part after $\cot$. (I pointed this out to @IosifPinelos when he gave the same answer.) Sometimes Mathematica makes mistakes. $\endgroup$ – Matt F. Nov 20 at 11:37

$\begingroup$ @Matt F. However, this is a realvalued expression. See the addition in my answer. $\endgroup$ – user64494 Nov 20 at 11:59

$\begingroup$ It should be realvalued, but pdf /. {a > 1.1, b > 1.1, c > 1.1, d > 1.1, t > 0.7} gives 1.04867 + 1.74393*10^16 I instead. There's probably an equivalent expression whose factors are all realvalued, and I might be able to find it, but I have never found a simple way to get Mathematica to figure that out for me. $\endgroup$ – Matt F. Nov 20 at 16:48

$\begingroup$ @Matt F: 1.74393*10^16 I is a roundoff error when calculating a complicated expression with complexes. Don't hesitate to ask for further explanation in need. $\endgroup$ – user64494 Nov 20 at 17:28

1$\begingroup$ @MattF. : "There's probably an equivalent expression whose factors are all realvalued"  Such expressions are obtained "manually" in my answer. $\endgroup$ – Iosif Pinelis Nov 20 at 21:13