Let $U\sim \mathcal U(0, 1)$ be a random variable uniformly distributed over the interval $[0, 1]$. Let $X_1, X_2\sim \Gamma(a, b)$ be two iid random variables with a Gamma distribution. Now it is straightforward to calculate the mean and variance of the random linear combination $$Y = U\cdot X_1 + (1-U)\cdot X_2,$$ which is simply $\mathbb E(Y)=\mathbb E(X_1)$ and $var(Y)=\frac23 var(X_1)$. However, it would be very interesting for me to know the precise distribution of $Y$ but I didn't find any result for this particular problem. Does anyone know what the distribution is of $Y$?

  • $\begingroup$ This is not a random mixture. A mixture of two identical distributions is the distribution itself, whence making a random choice of the weights will not change the distribution. This is a random linear combination or, if you like, a random convex combination of the two values. So, which question do you want answered: the one about a random mixture or the one you have written mathematically? $\endgroup$
    – whuber
    Commented Jul 18, 2022 at 14:49
  • 1
    $\begingroup$ Thanks for clarifying the notation - I honestly did not know how to label this combination. I would like to know the answer to the mathematically formulated problem. $\endgroup$
    – jfiedler
    Commented Jul 19, 2022 at 8:32
  • $\begingroup$ to clarify @whuber's point a mixture distribution would be something like $p_Y(y) = w*p_{X_1}(y) + (1-w)*p_{X_2}(y)$ where $p_Y, p_{X_1}, p_{X_2}$ are density functions and as whuber points out if $p_{X_1}=p_{X_2}$ then $p_Y=p_{X_1}=p_{X_2}$ for any value of $w$. $\endgroup$
    – bdeonovic
    Commented Jul 19, 2022 at 16:08

1 Answer 1


A direct attack (via integration) on computing the density looks intractable.

Instead, we may more easily compute the characteristic function of $Y$ (when the scale factor $b=1,$ which we may assume without any loss of generality simply by changing the units in which we express $Y$) as

$$\begin{aligned} \phi_Y(t;a) &= E\left[e^{itY}\right] = E\left[e^{it(UX_1+(1-U)X_2)}\right] = E_U\left[E_{X_1,X_2}\left[e^{it(UX_1+(1-U)X_2)}\right]\right]\\ &= E_U\left[ E_{X_1} \left[e^{it UX_1}\right]\, E_{X_2} \left[e^{it (1-U)X_2}\right]\right]\\ &= E_U\left[(1 - i U t)^{-a}\, (1 - i(1-U)t)^{-a}\right]\\ &= \int_0^1 \left((1 - iut)(1 - i(1-u)t)\right)^{-a}\,\mathrm{d}u\\ &= \frac{_2F_1\left(1, 2-2a, 2-a; \frac{i}{2i+t}\right)\ {-}\ _2F_1\left(1, 2-2a, 2-a; \frac{i+t}{2i + t}\right)}{(a-1)t(2i + t)(1 - it)^{a-1}}. \end{aligned}$$

The hypergeo library for R will efficiently compute these hypergeometric functions, as shown in these graphs of $\phi_Y$ for a range of $a.$

Figure 1

The solid graph plots the real part of $\phi_Y$ and the dotted graph plots its imaginary part.

(There is a technical problem: the foregoing formula is not defined for positive integral values of $a.$ You can get close enough by computing the average of $\phi_Y(t;a\pm2\epsilon)$ when $a$ is within $\epsilon$ of an integer and $\epsilon$ is small. It doesn't have to be terribly small for this to be highly accurate; I used $\epsilon = 10^{-5}.$ Smaller values will cause too much loss of precision in double precision floating point calculations.)

That answers the question. However, often it's more useful to have a density or distribution function. Either can be obtained through a Fourier transform. The density, for instance, is

$$f_Y(y;a) = \frac{1}{2\pi}\int_{\mathbb R} e^{- i y t} \phi_Y(t;a)\,\mathrm{d}t.$$

Numerical integration is possible, albeit a little slow for small values of $a.$ For $a=3/2$ it took half a minute to generate 101 points on the curve below, which matches the histogram of 100,000 values of $Y$ generated directly from its definition.

Figure 2

For $a=6$ the time dropped to four seconds; for $a=7.5$ to under half a second; etc. The timing is a rather irregular function of $a$ (odd multiples of $1/2$ tend to be very quick to compute) so take some care if you need many evaluations of $f_Y.$

This is the full R code.

# The characteristic function of Y.
phi. <- function(t, a) {
  (hypergeo(1, 2 - 2*a, 2 - a, 1.i / (2.i + t)) -
     hypergeo(1, 2 - 2*a, 2 - a, (1.i + t) / (2.i + t))) /
    (t * (2.i + t) * (1 - 1.i * t)^(a-1) * (1 - a))
phi <- function(t, a, eps=1e-5) {
  ifelse(t == 0, 1, 
         if(isTRUE(abs(a - round(a)) > eps)) phi.(t, a) else {
           (phi.(t, a - 2*eps) + phi.(t, a + 2*eps)) / 2
# Test `phi` by plotting it.
for (s in c(bquote(1/4), bquote(3/2), bquote(10))) {
  a <- eval(s)
  plot(c(-3,3), c(-1,1), type="n", 
       ylab="", yaxp=c(-1,1,2),
       xlab=expression(italic(t)), cex.lab = 1.5,
       main=bquote(italic(a)==.(s)), cex.main = 1.5)
  curve(Im(phi(x, a)), add=TRUE, lty = 3, n=201)
  curve(Re(phi(x, a)), add=TRUE, n=201)
# The density of Y.
f <- Vectorize(function(x, a, ...) {
  h <- function(t) Re(exp(-1.i * t * x) * phi(t, a)) / (2 * pi)
  integrate(h, -Inf, Inf, ...)$value
}, "x")
# Test with direct simulation of Y.
a <- 5/2 # Watch out for small or integral values -- integration can take time
b <- 1
n <- 1e5
X.1 <- rgamma(n, a, b)
X.2 <- rgamma(n, a, b)
U <- runif(n)
Y <- U * X.1 + (1 - U) * X.2

hist(Y, freq=FALSE, breaks=100)
  # Many subdivisions are needed for small values of `a`.
  curve(f(x, a, subdivisions=5e3, rel.tol=1e-3, abs.tol=1e-4), n=101, lwd=2, add = TRUE)
  • 1
    $\begingroup$ Ah, many thanks for this very comprehensive answer! This helps a lot and also shows me that there are most likely no simple closed form solutions - maybe only for some special choices of Gamma-distribution parameters. Again, many thanks! $\endgroup$
    – jfiedler
    Commented Jul 20, 2022 at 7:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.