# Statistical moments of X / (X + Y) when X and Y are two independent Beta distributions

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 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.

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 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, then $$0, $$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, 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).

• The second case has a complex part after the $\cot$? – Matt F. Nov 19 at 21:01
• @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.) – Iosif Pinelis Nov 20 at 4:30
• @IosifPinelis Sorry for causing that extra work – Maximilian Janisch Nov 20 at 8:32
• @MaximilianJanisch : Nothing to be sorry about. I had to be more careful. – Iosif Pinelis Nov 20 at 13:00
• thank you for this great answer. – Kerighan Nov 20 at 17:19

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$$):

$$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}$$

• 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. – Kerighan Nov 19 at 16:53
• 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$. – Matt F. Nov 19 at 17:15
• 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$ – Maximilian Janisch Nov 19 at 17:16
• @MattF. I implemented some of your suggestions – Maximilian Janisch Nov 19 at 17:21
• @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(b-1)}b<0$ under your final integral when $b<1/2$ and $a$ is close to $1$. – 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(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 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] • 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. – Matt F. Nov 20 at 11:37
• @Matt F. However, this is a real-valued expression. See the addition in my answer. – user64494 Nov 20 at 11:59
• It should be real-valued, 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 real-valued, 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. – Matt F. Nov 20 at 16:48
• @Matt F: 1.74393*10^-16 I is a round-off error when calculating a complicated expression with complexes. Don't hesitate to ask for further explanation in need. – user64494 Nov 20 at 17:28
• @MattF. : "There's probably an equivalent expression whose factors are all real-valued" -- Such expressions are obtained "manually" in my answer. – Iosif Pinelis Nov 20 at 21:13