Summation of independent Gamma distributions?

Question

Let $X_k$ be independent gamma distribution: $$X_k =\frac{X}{k}$$

where $X$ is a gamma distributions on $(0, \infty)$.

How to find the distribution of the summation of $X_k$:

$$\sum_{k=0}^{+\infty} X_k $$

Can anyone provide some detailed calculation steps ?

Answer 1

Let $Y:=\sum_{k\geqslant 1}X_k$, then its easy to check that $\operatorname{Var}[Y]=\sum_{k\geqslant 1}\operatorname{Var}[X_k]=\sum_{k\geqslant 1}\frac{4}{\pi^4k^4}<\infty $ so there is a theorem that ensures that $Y$ is well-defined and almost sure finite. Indeed we can see that $\operatorname{E}[Y]=\frac2{3}$.

Also observe that $\varphi _k(t):=\operatorname{E}[e^{-tX_k}]=(1+2t/(\pi^2k^2))^{-2}$, so it follows that $$ \varphi (t):=\operatorname{E}[e^{-tY}]=\prod_{k\geqslant 1}\varphi _k(t)=\left[\prod_{k\geqslant 1}\left(1+\frac{2t}{\pi^2k^2}\right)\right]^{-2}=\frac{2t}{(\sinh(\sqrt{2t}))^2} $$ Now, the inverse Laplace transform of $\varphi$ will be the density of $Y$. I dont see a way to compute this inverse Laplace transform by hand so I tried to compute it using Wolfram Mathematica, and it gives as a result that

$$ f_Y(s)=\pi^2\sum_{n\geqslant 1}((n\pi)^2 s-3)n^2 e^{-(n\pi)^2 s/2}\tag1 $$

A little easier is to find the CDF of $Y$, as it is the inverse Laplace transform of $\frac1{t}\varphi (t)$, then we get (with the help of Wolfram Mathematica) that

$$ F_Y(s)=1-2\sum_{n\geqslant 1}e^{-(n\pi)^2 s/2}((n\pi)^2s-1)\tag2 $$

After some work I can confirm that both expressions in (1) and (2) given by Wolfram Mathematica are correct for all $s>0$. Its annoying but they doesn't hold for $s=0$. You can justify them rigorously with a bit of Fourier, real and complex analysis.

Using a numerical computation of the inverse Laplace transform I get

This is the code that I used in Julia to do the figure

using InverseLaplace, GLMakie

# Laplace transform of density of Y
f(t) = 2t/(sinh(sqrt(2t)))^2

# Inverse Laplace transform of Y
pdf1 = ILT(t -> f(t), 64)
CDF1 = ILT(t -> f(t)/t)

# Drawing
fig = Figure()
ax = Axis(fig[1, 1],
    xlabel = "Support",
    ylabel = "Probability"
)
x = range(0, 3, length=100)
y1 = pdf1.(x)
y2 = CDF1.(x)
lines!(ax,x,y1, color = :blue, label = "pdf of Y")
lines!(ax,x,y2, color = :red, label ="CDF of Y")
axislegend()
fig
save("./Julia/figure.png", fig)
````