Sum of inverse gamma random variables with unity parameters? (challenging integral)

definite integralsexponential functiongamma distributionsmooth-functionsspecial functions

I came across a problem on this site (which I now cannot find) asking about the distribution of sums of i.i.d. inverse gamma random variables for paramters $\alpha=\beta=1$. Defining $Z=X+Y$ with $X,Y\sim\operatorname{InvGamma}(1,1)$ we can apply the usual formula to obtain the density
$$
f_Z(z)=\int_0^z(z-x)^{-2}\exp\left(-\tfrac{1}{z-x}\right)x^{-2}\exp\left(-\tfrac{1}{x}\right)\,\mathrm dx,\qquad z>0.
$$

This integral seems deceptively hard which is mainly due to the exponential terms. Now we do have the useful fact that
$$
\int(z-x)^{-2}\exp\left(-\tfrac{1}{z-x}\right)\,\mathrm dx=\exp\left(-\tfrac{1}{z-x}\right),
$$

however, I don't think integrating by parts gets us very far. Another observation is that the integrand is symmetric about $x=z/2$ so that
$$
f_Z(z)=2\int_0^{z/2}(z-x)^{-2}\exp\left(-\tfrac{1}{z-x}\right)x^{-2}\exp\left(-\tfrac{1}{x}\right)\,\mathrm dx.
$$

But again I'm not sure this accomplishes much. How can we deal with the exponential terms to put this integral in a more useful form for evaluating? Taylor series expansions are certainly out. Maybe a clever substitution? The integrand is quite interesting and reminds me of a bump function. Here it is for $z=1$:

enter image description here

Best Answer

Substitute $y=z^2/[4x(z-x)]$. Then $0<x<z/2$ maps to $1<y<\infty$, and $$f_Z(z)=\frac{8}{z^3}\int_1^\infty\left(1-\frac1y\right)^{-1/2}\exp\left(-\frac{4y}{z}\right)\,{\rm d}y=-\frac{4}{z^3}g'\left(\frac2z\right),$$ where $$g(a)=\int_1^\infty\frac{e^{-2ay}\,{\rm d}y}{\sqrt{y(y-1)}}\ \underset{y=(1+t)/2}{\phantom{\Big[}=\phantom{\Big]}}\ \int_1^\infty\frac{e^{-a(1+t)}}{\sqrt{t^2-1}}\,{\rm d}t=e^{-a}K_0(a)$$ uses the modified Bessel function $K_0$. Hence $f_Z(z)=(4/z^3)e^{-2/z}\big(K_0(2/z)+K_1(2/z)\big)$.