I am trying to find a closed form expression of the Gini concentration ratio (G) for a mixture of $M$ Gamma distributions. Thus, I am trying to solve:
\begin{equation} G=\dfrac{1}{2\mu}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}dx dy|x-y|P(x)P(y) \end{equation}
where
\begin{equation} P(·) = \sum_{i=1}^M\omega_iG(·;\alpha_i,\beta_i) \end{equation}
and
\begin{equation} G(x;\alpha_i,\beta_i)=\dfrac{\beta_i^{\alpha_i}}{\Gamma(\alpha_i)}x^{\alpha_i-1}e^{-\beta_ix} \end{equation}
is the Gamma PDF. $\alpha$ and $\beta$ can take any non-negative real value.
I have managed to solve most of it, but I am stuck in an integral of the form:
\begin{equation} \int_0^\infty dx G(x;\alpha_i,\beta_i) F(x;\alpha_j,\beta_j) \end{equation}
where $F(x;\alpha_j,\beta_j)$ is the CDF of $G(x;\alpha_j,\beta_j)$. I recognize that this is the expectation value of the CDF and I know how to solve the case where $i=j$. The case for $i\neq j$ is a different story altogether. Integrating by parts seems to only lead to a recursive solution, and trying to find the relation between $G_j(x;\alpha_j,\beta_j)$ and $G_i(x;\alpha_i,\beta_i)$, so that I can solve the case $i=j$, seems to worsen things.
Hence, what I would like to find is, if posible, a closed form expression for this final integral.
Given that somebody might be interested in this problem in a future, I will post my work-around solution to solve this integral, for which I have not found yet an exact analytical expression. We can rewrite our integral of interest as:
\begin{equation} \int_0^\infty dx G(x;\alpha_i,\beta_i) F(x;\alpha_j,\beta_j)e^{x}e^{-x}=\int_0^\infty dx f(x)e^{-x} \end{equation}
This is the type of integrals that can be approximated by the Laguerre-Gauss quadrature method. Hence, the approximated solution to the original integral is given by:
\begin{equation} \int_0^\infty dx G(x;\alpha_i,\beta_i) F(x;\alpha_j,\beta_j)\approx \sum_{k=1}^n w_k G(x_k;\alpha_i,\beta_i) F(x_k;\alpha_j,\beta_j)e^{x_k} \end{equation}
where $x_k$ is the $k$-th root of the Laguerre polynomial $L_n(x)$ and the weights $w_k$ are given by:
$$w_k=\dfrac{x_k}{(n+1)^2[L_{n+1}(x_k)]^2} $$
It is not a closed form expression, but it is closed enough for my purposes.