Let $\epsilon_\pm \ge 1$ be real numbers. Consider a following random variable: \begin{equation} {\mathcal R}:= \frac{1-Z}{Z} \cdot \Xi \quad (i) \end{equation} where $Z \in (0,1)$ is a random variable with a density $\rho_Z(z) = z^{\epsilon_+-1} (1-z)^{\epsilon_--1} /B(\epsilon_-,\epsilon_+)$ and $\Xi$ is a uniform random variable , i.e. $\Xi = U(0,1)$. Both variables $Z$ and $\Xi$ are independent and $B(\cdot,\cdot)$ is the beta function.
We have shown that the probability density of the variable $R$ is given as follows: \begin{equation} \rho_R(x) = \frac{x^{-1-\epsilon_+}}{(1+\epsilon_+) B(\epsilon_-,\epsilon_+)} \cdot _2F_1\left(\epsilon_++1,\epsilon_-+\epsilon_+;\epsilon_++2;-\frac{1}{x}\right) \cdot 1_{x \ge 0} \quad (ii) \end{equation}
Now, the natural thing is to check the normalization of the pdf above. If we now use the functional identity http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/17/02/09/0002/ and the series expansion of the hypergeometric function and integrate term by term then after some manipulations we arrive at a following identity:
\begin{eqnarray} B_{\frac{1}{2}}(\epsilon_--1,\epsilon_++1)-B_{\frac{1}{2}}(\epsilon_-,\epsilon_+)-B_{\frac{1}{2}}(\epsilon_+,\epsilon_-)+B_{\frac{1}{2}}(\epsilon_++1,\epsilon_--1) = -\frac{(\epsilon_--\epsilon_+-1) \Gamma (\epsilon_--1) \Gamma (\epsilon_+)}{\Gamma (\epsilon_-+\epsilon_+)} \quad (iii) \end{eqnarray} where $B_z(\cdot,\cdot)$ is the incomplete beta function.
In[566]:= {em, ep} = RandomReal[{1, 10}, 2, WorkingPrecision -> 50];
( (ep - em + 1) Gamma[em - 1] Gamma[ep])/
Gamma[em +
ep] + (NIntegrate[(t^(ep - 1) - t^(em - 2)) (1 -
t) (1 + t)^(-em - ep), {t, 0, 1}, WorkingPrecision -> 30])
((-1 + em - ep) Gamma[-1 + em] Gamma[ep])/
Gamma[em + ep] - (-Beta[1/2, em - 1, ep + 1] + Beta[1/2, em, ep] +
Beta[1/2, ep, em] - Beta[1/2, ep + 1, em - 1])
Out[567]= 0.*10^-32
Out[568]= 0.*10^-50
Now, I have two questions. The first one is simple, i.e.how do we prove the identity $(iii)$ otherwise?
The second question is related to the moments of the distribution of ${\mathcal R}$. Take some $m \ge 0$ and real. Then from the definition $(i)$ we clearly have: \begin{equation} E\left[ {\mathcal R}^m \right] = \frac{B(\epsilon_-+m,\epsilon_+-m)}{B(\epsilon_-,\epsilon_+)} \cdot \frac{1}{m+1} \end{equation}
Can we actually prove the same result by using the closed form expression $(ii)$ for the pdf of ${\mathcal R}$ ?
Here is the answer to the first question.
By using the identity http://functions.wolfram.com/GammaBetaErf/Beta3/17/02/03/0001/ we can write: \begin{eqnarray} lhs &=& B(\epsilon_--1,\epsilon_++1) - B(\epsilon_-,\epsilon_+) \\ &=&\frac{\Gamma(\epsilon_--1) \Gamma(\epsilon_++1)}{\Gamma(\epsilon_-+\epsilon_+)} - \frac{\Gamma(\epsilon_-) \Gamma(\epsilon_+)}{\Gamma(\epsilon_-+\epsilon_+)}\\ &=& \left(\epsilon_+-\epsilon_-+1\right)\frac{\Gamma(\epsilon_--1) \Gamma(\epsilon_+)}{\Gamma(\epsilon_-+\epsilon_+)} = rhs \end{eqnarray}
Now we note the following identity: \begin{equation} F_{2,1} \left[ 1,b,c, z \right] = (c-1) z^{1-c} \cdot (1-z)^{-b-1+c} \cdot \left( B(b-c+1,c-1) - B_{1-z}(b-c+1,c-1) \right) \quad (I) \end{equation} for $b,c \ge 1$ and $-1 < z < 1$.
Here is the answer to the second question. \begin{eqnarray} &&E\left[ {\mathcal R}^m \right]= \int\limits_0^\infty x^m \cdot \rho_{\mathcal R}(x) dx \\ && =\int\limits_0^1 x^m \left( \frac{\epsilon_+}{\epsilon_--1} + \frac{x^{-1-\epsilon_-}}{(1+\epsilon_+)B(\epsilon_+,\epsilon_-)} F_{2,1} \left[-1+\epsilon_-,\epsilon_-+\epsilon_+,\epsilon_-,-x \right]\right) dx + \\ && \int\limits_1^\infty x^m \left( \frac{x^{-1-\epsilon_+}}{(1+\epsilon_+) B(\epsilon_+,\epsilon_-)} F_{2,1} \left[ 1+\epsilon_+,\epsilon_-+\epsilon_+,2+\epsilon_+, -\frac{1}{x}\right] \right) dx \\ &&= \frac{\epsilon_+}{\epsilon_--1} \cdot\frac{1}{m+1}+\\ &&\frac{1}{(1-\epsilon_-)B(\epsilon_-,\epsilon_+)} \frac{1}{m+1} \left( \, _2F_1(\epsilon_--1,\epsilon_-+\epsilon_+;\epsilon_-;-1)-\frac{(\epsilon_--1) \, _2F_1(\epsilon_-+\epsilon_+,\epsilon_-+m;\epsilon_-+m+1;-1)}{\epsilon_-+m} \right)+\\ &&\frac{1}{(1+\epsilon_+) B(\epsilon_+,\epsilon_-)} \frac{1}{m+1} \left( \frac{(\epsilon_++1) \, _2F_1(\epsilon_-+\epsilon_+,\epsilon_+-m;\epsilon_+-m+1;-1)}{\epsilon_+-m}-\, _2F_1(\epsilon_++1,\epsilon_-+\epsilon_+;\epsilon_++2;-1) \right) \\ &&= \frac{1}{m+1} \cdot \frac{B(\epsilon_-+m,\epsilon_+-m)}{B(\epsilon_-,\epsilon_+)} \end{eqnarray} In the second line we used http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric2F1/17/02/09/0002/ and in the third line we integrated term by term the power series expansions. In the last line we firstly used the Pfaff transformation https://en.wikipedia.org/wiki/Hypergeometric_function#Transformation_formulas to reduce the hypergeometric functions to values at one half and then we used the identity $(I)$ to express the later quantities (i.e. the hypergeometric functions at one half) through both beta and incomplete beta functions. Finally we used http://functions.wolfram.com/GammaBetaErf/Beta3/17/02/03/0001/ to simplify the result.