What is the convolution of a normal distribution with a gamma distribution?

Question

Is there a closed form expression for the convolution of a normal distribution (ND) with a gamma distribution (GD)? There does not seem to be a direct method of solving this convolution.

Answer 1

Often, convolving something with itself gives a solution even when the more direct convolution of two different distributions has no obvious answer. To convolve a ND and a GD, I used Pearson III and convolved two Pearson III distributions with themselves after reparameterization of those Pearson III distributions to be ND and GD using Mathematica. $$\text{PDF}[\text{PearsonDistribution}[3,a,b,x,y,z],t]=\begin{array}{cc} & \begin{cases} \dfrac{\sqrt{\frac{a}{z}} e^{-\dfrac{(a t+b)^2}{2 a z}}}{\sqrt{2 \pi }} & y=0\land a z>0 \\ \dfrac{a \ e^{-\dfrac{a \left(t+\dfrac{z}{y}\right)}{y}} \left(\dfrac{a \left(t+\dfrac{z}{y}\right)}{y}\right)^{\dfrac{a z}{y^2}-\dfrac{b}{y}}}{y \Gamma \left(-\dfrac{b}{y}+\dfrac{a z}{y^2}+1\right)} & y^2>0\land a (t y+z)>0 \\ \end{cases} \\ \end{array}$$

Then ND from Pearson III is $$\text{PDF}\left[\text{PearsonDistribution}\left[3,1,-\mu ,x,0,\sigma ^2\right],t\right]=\dfrac{ e^{-\dfrac{(t-\mu )^2}{2 \sigma ^2}}}{\sigma\sqrt{2 \pi }}\,, \\ $$ And GD from Pearson III is $$\text{PDF}\left[\text{PearsonDistribution}\left[3,\beta ,1,x,1,\frac{\alpha }{\beta }\right],t-\frac{\alpha }{\beta }\right]=\begin{array}{cc} & \begin{cases} \dfrac{\beta e^{-\beta \ t} (\beta \ t)^{\alpha -1}}{\Gamma (\alpha )} & \beta t>0 \\ 0 & \text{Otherwise} \\ \end{cases} \\ \end{array}\,.$$ Then ND*GD is $$f(s)=\text{Convolve}\left[\text{PDF}\left[\text{PearsonDistribution}\left[3,\beta ,1,x,1,\frac{\alpha }{\beta }\right],t-\frac{\alpha }{\beta }\right],\text{PDF}\left[\text{PearsonDistribution}\left[3,1,-\mu ,x,0,\sigma ^2\right],t\right],t,s\right]=2^{-\frac{\alpha }{2}} \beta ^{\alpha } \sigma ^{\alpha -2} e^{-\frac{(s-\mu )^2}{2 \sigma ^2}} \left(\frac{\sigma \, _1F_1\left(\frac{\alpha }{2};\frac{1}{2};\frac{\left(\beta \sigma ^2-s+\mu \right)^2}{2 \sigma ^2}\right)}{\sqrt{2} \Gamma \left(\frac{\alpha +1}{2}\right)}+\frac{\left(-\beta \sigma ^2-\mu +s\right) \, _1F_1\left(\frac{\alpha +1}{2};\frac{3}{2};\frac{\left(\beta \sigma ^2-s+\mu \right)^2}{2 \sigma ^2}\right)}{\Gamma \left(\frac{\alpha }{2}\right)}\right)\,.$$ That is, after a lot of simplifying. Note, $_1F_1(a;b;z)$ is the confluent hypergeometric function of the first kind.

Answer 2

For integral shape parameters there is a relatively simple form.

Without any loss of generality, we may suppose the Normal distribution is a standard Normal, because any shifting of it translates directly to a shift of the convolution; and we may incorporate all scale factors in the Gamma distribution (and then effect a rescaling at the end if needed). This means we may limit the analysis to the standard Normal distribution $\Phi$ with density

$$\phi(z) = \frac{1}{\sqrt{2\pi}}\exp\left(-z^2/2\right)$$

and Gamma densities that are zero on negative values and otherwise are given by

$$\gamma_{a,\sigma}(t) = \frac{1}{\Gamma(a)\sigma^a} t^{a-1}\exp\left(-t/\sigma\right).$$

The convolution can be computed as

$$f(z; a, \sigma) = \left(\phi \star \gamma_{a,\sigma}\right)(z) = \frac{1}{\sqrt{2\pi}\Gamma(a)\sigma^a}\int_0^\infty t^{a-1}\exp\left(-t/\sigma\right)\exp\left(-(z-t)^2/2\right)\,\mathrm{d}t.$$

The instances with integral $a$ promise to yield simple expressions. Let's begin at $a=1,$ for which the calculation yields to the usual expedient of completing the square:

$$\begin{aligned} f(z;1,\sigma) &= \frac{1}{\sqrt{2\pi}\sigma} \int_0^\infty \exp\left(-t/\sigma - (z-t)^2/2\right)\,\mathrm{d}t\\ &= \frac{\exp\left(-\left(z-\frac{1}{2\sigma}\right)/\sigma\right)}{\sqrt{2\pi}\sigma}\int_0^\infty \exp\left(-(t - (1/\sigma-z))^2\right)/2)\,\mathrm{d}t\\ &= \frac{\exp\left(-\left(z-\frac{1}{2\sigma}\right)/\sigma\right)}{\sqrt{2\pi}\sigma}\int_{-\infty}^{z-1/\sigma} \exp\left(-t^2/2\right)\,\mathrm{d}t\\ &=\exp\left(-\left(z-\frac{1}{2\sigma}\right)/\sigma\right)\Phi\left(z-1/\sigma\right). \end{aligned}$$

That, at least, can be considered a "closed form." We can exploit it to find other convolution densities. Returning to the convolution integral, observe

$$\begin{aligned} f^\prime(z;a,\sigma) &= \frac{1}{\sqrt{2\pi}\Gamma(a)\sigma^a}\int_0^\infty \frac{\mathrm d}{\mathrm{d}z} t^{a-1}\exp\left(-t/\sigma\right)\exp\left(-(z-t)^2/2\right)\,\mathrm{d}t\\ &= \frac{1}{\sqrt{2\pi}\Gamma(a)\sigma^a}\int_0^\infty t^{a-1}\exp\left(-t/\sigma\right)(t-z)\exp\left(-(z-t)^2/2\right)\,\mathrm{d}t\\ &= \frac{1}{\sqrt{2\pi}\Gamma(a)\sigma^a}\left[\int_0^\infty t^{a}e^\left(-t/\sigma\right)e^\left(-(z-t)^2/2\right)\,\mathrm{d}t - z\int_0^\infty t^{a-1}e^\left(-t/\sigma\right)e^\left(-(z-t)^2/2\right)\,\mathrm{d}t\right]\\ &= \frac{\sigma\Gamma(a+1)}{\Gamma(a)} f(z;a+1,\sigma) - z f(z;a,\sigma)\\ &= \sigma a f(z;a+1,\sigma) - z f(z;a,\sigma). \end{aligned}$$

Solving for $f(z;a+1,\sigma)$ and writing $D=\mathrm{d}/\mathrm{d}z$ for the differential operator produces the fundamental recurrence

$$f(z;a+1,\sigma) = \frac{f^\prime(z;a,\sigma) + z f(z;a,\sigma)}{\sigma a} = \frac{D+z}{\sigma a} f(z;a,\sigma)$$

Repeating this $\lfloor a\rfloor - 1$ times gives

$$\begin{aligned} f(z;a,\sigma) &= \frac{D+z}{\sigma (a-1)}\frac{D+z}{\sigma (a-2)}\cdots\frac{D+z}{\sigma(a-\lfloor a\rfloor+1)} f(z;a-\lfloor a\rfloor+1,\sigma) \\ &= \frac{\Gamma(a-\lfloor a \rfloor + 1)}{\sigma^{\lfloor a\rfloor-1}\Gamma(a)} (D+z)^{\lfloor a\rfloor-1} f(z;a-\lfloor a\rfloor+1,\sigma). \end{aligned}$$

For integral $a$ it simplifes further to

$$f(z;a,\sigma) = \frac{1}{\sigma^{a-1}\Gamma(a)} (D+z)^{a-1} f(z;1,\sigma).$$

Just to see the patterns, let $\sigma=1,$ drop the "$z-1$" arguments from $\Phi$ and $\phi$ (reducing the visual clutter in the formulas) and use $\phi^\prime(z-1) = -(z-1) \phi(z-1)$ to compute

$$\begin{aligned} f(z;1,1) &= e^{1/2-z} \Phi\\ f^\prime(z;1,1) &= e^{1/2-z}\left(-\Phi + \phi\right)\\ f(z;2,1) &= \frac{1}{1}\left(f^\prime(z;1,1) + z f(z;1,1)\right) = e^{1/2-z}\left((z-1)\Phi + \phi\right)\\ f^\prime(z;2,1) &= e^{1/2-z}\left(-(z-1)\Phi-\phi + \Phi + (z-1)\phi - (z-1)\phi\right)\\ &= e^{1/2-z}\left((2-z)\Phi - \phi\right)\\ f(z;3,1) &= \frac{1}{2}\left(f^\prime(z;2,1) + z f(z;2,1)\right) \\ &= \frac{1}{2}e^{1/2-z}\left((z^2-2z+2)\Phi + (z-1)\phi\right)\\ etc. \end{aligned}$$

Generally, the form of these densities is

$$f(z;a,1) = e^{(1/(2\sigma)-z)/\sigma}\left[P(z;a,\sigma)\Phi(z-1/\sigma) + Q(z;a,\sigma)\phi(z-1/\sigma)\right]$$

where $P$ and $Q$ are polynomials in $z.$ The basic recurrence (for $\sigma=1$) translates to

$$\begin{aligned} P(z;a+1) &= \left[(z-1)P(z;a) + P^\prime(z;a)\right]/a\\ Q(z;a+1) &= \left[P(z;a) + Q^\prime(z;a)\right]/a \end{aligned}.$$

It is now immediate that

For positive integral $a,$ the density $f(z;a,\sigma)$ is $\exp(1/(2\sigma) - z)/\sigma)$ times a polynomial linear combination of $\Phi(z-1)$ and $\phi(z-1).$ The degree of $P$ is $a-1$ and the degree of $Q$ is $a-2$ (when $a\ge 2$).

Here are histograms of simulated values for $a=1,3,10$ overplotted with the graphs of $f(z;a,1)$ to demonstrate the agreement.

---

By comparing these to all the Wikipedia listings of continuous distributions supported on the real line we can establish that they are not among those lists. Most of this is obvious, with the exception perhaps of the Pearson system of distributions. But these can be ruled out since (by definition) the logarithmic derivative of a Pearson density is a rational function of its argument with at most one zero and two (complex) poles. The recurrence relation for $f(z;a,\sigma)$ readily demonstrates that no such relationship holds.