# Sum of truncated Gammas

I have a set of i.i.d. variables $X_i$ that are distributed according to a truncated $\text{Gamma}(\alpha,\beta)$ distribution, with support on $[0,w)$ where $w$ is a known constant. What's the distribution of $Y=\sum_{i=1}^NX_i$?

I've been trying to get the distribution for $N=2$ to generalise from there. The pdf is:

$$f(x) = \frac{\beta^\alpha}{\gamma(\alpha,w\beta)}x^{\alpha-1}e^{-x\beta}$$

where $\gamma$ is the lower incomplete gamma function and $f$ is supported on $[0,w)$. Then for $N=2$ the distribution of $Y$ is:

\begin{aligned} g(y) &= (f*f)(y) \\ &=\int_{-\infty}^{+\infty}f(x)f(y-x)\text dx \\ &=\int_{-\infty}^{+\infty}\mathbb{1}_{0\leq x < w}\frac{\beta^\alpha}{\gamma(\alpha,w\beta)}x^{\alpha-1}e^{-x\beta}\mathbb{1}_{0\leq y-x < w}\frac{\beta^\alpha}{\gamma(\alpha,w\beta)}(y-x)^{\alpha-1}e^{-(y-x)\beta}\text dx \\ &=\left(\frac{\beta^\alpha}{\gamma(\alpha,w\beta)}\right)^2 e^{-y\beta} \int_{\max(0,y-w)}^{\min(y,w)} x^{\alpha-1} (y-x)^{\alpha-1}\text dx \\ &=\frac{\beta^{2\alpha}}{\gamma(\alpha,w\beta)^2} e^{-y\beta}y^{2\alpha-1} \left[\mathbb{1}_{0\leq y < w}\text B_{\frac x y}(\alpha,\alpha)\middle|_{0}^{y} + \mathbb{1}_{w\leq y< 2w}\text B_{\frac x y}(\alpha,\alpha)\middle|_{y-w}^{w} \right] \end{aligned}

The limits of integration can be derived as described here and $\text B_{x}(\alpha,\beta)$ is the incomplete beta function.

\begin{aligned} \left[\text B_{\frac x y}(\alpha,\alpha)\right]_{0}^{y} &= \text B_{1}(\alpha,\alpha) - \text B_{0}(\alpha,\alpha) \\ &= \text B(\alpha,\alpha)\\ &= \frac{\Gamma(\alpha)^2}{\Gamma(2\alpha)} \\ \left[\text B_{\frac x y}(\alpha,\alpha)\right]_{y-w}^{w} &= \text B_{\frac w y}(\alpha,\alpha) - \text B_{1 - \frac w y}(\alpha,\alpha) \\ &= \text B(\alpha,\alpha) \left( I_{\frac w y}(\alpha,\alpha) - \text I_{1 - \frac w y}(\alpha,\alpha) \right) \\ &= \text B(\alpha,\alpha) \left(\text I_{\frac w y}(\alpha,\alpha) - \left( 1 - \text I_{\frac w y}(\alpha,\alpha) \right) \right) \\ &= \text B(\alpha,\alpha) \left(2\text I_{\frac w y}(\alpha,\alpha) - 1 \right) \\ &= \frac{\Gamma(\alpha)^2}{\Gamma(2\alpha)} \left(2\text I_{\frac w y}(\alpha,\alpha) - 1 \right) \\ \end{aligned}

where $\text I_x(\alpha,\beta)$ is the regularized incomplete beta function and $\text B(\alpha,\beta)$ is the complete beta function. Therefore:

\begin{aligned} g(y) = \frac{\Gamma(\alpha)^2}{\Gamma(2\alpha)} \frac{\beta^{2\alpha}}{\gamma(\alpha,w\beta)^2} e^{-y\beta}y^{2\alpha-1} &\left[\mathbb{1}_{0\leq y < w} + \mathbb{1}_{w\leq y< 2w}\left(2\text I_{\frac w y}(\alpha,\alpha) - 1 \right) \right] \end{aligned}

So now let $Y = X_1 + X_2$ and let $Z = Y + X_3$ (in other words, $N=3$). The distribution of $Z$ is:

\begin{aligned} h(z) &= \left(f*g\right)(z) \\ &= \int_{-\infty}^{+\infty}f(z-y)g(y)\text dy \\ &= \int_{-\infty}^{+\infty} \mathbb{1}_{0\leq z-y < w}\frac{\beta^\alpha}{\gamma(\alpha,w\beta)}(z-y)^{\alpha-1}e^{-(z-y)\beta} \frac{\beta^{2\alpha}}{\gamma(\alpha,w\beta)^2} e^{-y\beta}y^{2\alpha-1} \mathbb{1}_{0\leq y < w}\frac{\Gamma(\alpha)^2}{\Gamma(2\alpha)} \text dy\\ &+ \int_{-\infty}^{+\infty} \mathbb{1}_{0\leq z-y < w}\frac{\beta^\alpha}{\gamma(\alpha,w\beta)}(z-y)^{\alpha-1}e^{-(z-y)\beta} \frac{\beta^{2\alpha}}{\gamma(\alpha,w\beta)^2} e^{-y\beta}y^{2\alpha-1} \mathbb{1}_{w\leq y< 2w}\frac{\Gamma(\alpha)^2}{\Gamma(2\alpha)} 2\text I_{\frac w y}(\alpha,\alpha) \text dy \\ &- \int_{-\infty}^{+\infty} \mathbb{1}_{0\leq z-y < w}\frac{\beta^\alpha}{\gamma(\alpha,w\beta)}(z-y)^{\alpha-1}e^{-(z-y)\beta} \frac{\beta^{2\alpha}}{\gamma(\alpha,w\beta)^2} e^{-y\beta}y^{2\alpha-1} \mathbb{1}_{w\leq y< 2w}\frac{\Gamma(\alpha)^2}{\Gamma(2\alpha)} \text dy \\ &= \frac{\beta^{3\alpha}}{\gamma(\alpha,w\beta)^3} e^{-z\beta} \frac{\Gamma(\alpha)^2}{\Gamma(2\alpha)} \int_{\max(0,z-w)}^{\min(z,w)} (z-y)^{\alpha-1}y^{2\alpha-1} \text dy\\ &+ 2\frac{\beta^{3\alpha}}{\gamma(\alpha,w\beta)^3} e^{-z\beta} \frac{\Gamma(\alpha)^2}{\Gamma(2\alpha)} \int_{\max(w,z-w)}^{\min(z,2w)} (z-y)^{\alpha-1}y^{2\alpha-1} \text I_{\frac w y}(\alpha,\alpha) \text dy \\ &- \frac{\beta^{3\alpha}}{\gamma(\alpha,w\beta)^3} e^{-z\beta} \frac{\Gamma(\alpha)^2}{\Gamma(2\alpha)} \int_{\max(w,z-w)}^{\min(z,2w)} (z-y)^{\alpha-1}y^{2\alpha-1} \text dy \\ &= \frac{\beta^{3\alpha}}{\gamma(\alpha,w\beta)^3} e^{-z\beta} \frac{\Gamma(\alpha)^2}{\Gamma(2\alpha)} z^{3\alpha-1} \left[ \mathbb{1}_{0 \leq z < w} \text B_{\frac x z}(2\alpha,\alpha) \middle|_{0}^{z} + \mathbb{1}_{w \leq z < 2w} \text B_{\frac x z}(2\alpha,\alpha) \middle|_{z-w}^{w} \right] \\ &+ 2\frac{\beta^{3\alpha}}{\gamma(\alpha,w\beta)^3} e^{-z\beta} \int_{\max(w,z-w)}^{\min(z,2w)} (z-y)^{\alpha-1}y^{2\alpha-1} \text B_{\frac w y}(\alpha,\alpha) \text dy \\ &- \frac{\beta^{3\alpha}}{\gamma(\alpha,w\beta)^3} e^{-z\beta} \frac{\Gamma(\alpha)^2}{\Gamma(2\alpha)} z^{3\alpha-1} \left[ \mathbb{1}_{w \leq z < 2w} \text B_{\frac x z}(2\alpha,\alpha) \middle|_{w}^{z} + \mathbb{1}_{2w \leq z < 3w} \text B_{\frac x z}(2\alpha,\alpha) \middle|_{z-w}^{2w} \right] \\ &= \mathbb{1}_{0 \leq z < w} \frac{\beta^{3\alpha}}{\gamma(\alpha,w\beta)^3} e^{-z\beta} \frac{\Gamma(\alpha)^3}{\Gamma(3\alpha)} z^{3\alpha-1} \\ &+ \mathbb{1}_{w \leq z < 2w} \frac{\beta^{3\alpha}}{\gamma(\alpha,w\beta)^3} e^{-z\beta} \frac{\Gamma(\alpha)^2}{\Gamma(2\alpha)} z^{3\alpha-1} \left(\text B_{\frac w z}(2\alpha,\alpha) - \text B_{1-\frac w z}(2\alpha,\alpha) \right) \\ &+ 2\frac{\beta^{3\alpha}}{\gamma(\alpha,w\beta)^3} e^{-z\beta} \int_{\max(w,z-w)}^{\min(z,2w)} (z-y)^{\alpha-1}y^{2\alpha-1} \text B_{\frac w y}(\alpha,\alpha) \text dy \\ &- \frac{\beta^{3\alpha}}{\gamma(\alpha,w\beta)^3} e^{-z\beta} \frac{\Gamma(\alpha)^2}{\Gamma(2\alpha)} z^{3\alpha-1} \left[ \mathbb{1}_{w \leq z < 2w} \text B_{\frac x z}(2\alpha,\alpha) \middle|_{w}^{z} + \mathbb{1}_{2w \leq z < 3w} \text B_{\frac x z}(2\alpha,\alpha) \middle|_{z-w}^{2w} \right] \\ \end{aligned}

This is as far as I've gotten.

NB: this is a specific subproblem of the one I posted at Sum of truncated Gammas and degenerate. A good solution to the present problem might help lead to a solution to the previous one, but the two are not exactly the same.

• An analysis similar to that of a sum of uniform distributions suggests you should expect a conditional formula with at least $N$ distinct branches: that is, no neat formula is going to emerge. – whuber Sep 23 '16 at 19:25
• According to the W there is a neat formula, though, for a certain value of neat. I'd love to have something like that. – Pedro Carvalho Sep 23 '16 at 19:35
• You can emulate the techniques in my answer in the thread I referenced to obtain that kind of formula, especially because the formula for the sum of two Gammas of the same scale is so simple. – whuber Sep 23 '16 at 20:24
• For the general case, how large is $N$? It should be reasonably straightforward to find the asymptotic distribution of $N$ via the CLT, which should work nicely other than for small $N$. – wolfies Sep 25 '16 at 12:54
• $N$ has a poisson distribution which I'm marginalising over in my posterior, and the nature of the problem is such that I'm gonna apply this model to cases where $\lambda=5$ and $\lambda=500$ equally, hence why I would like to find a general case. – Pedro Carvalho Sep 25 '16 at 13:21

I'm not sure if the above is correct, or how to calculate the second term in the brackets.

Your solution appears to be correct. And the incomplete Beta does not pose a problem ...

Given: $X_1$ and $X_2$ each have a $\text{Gamma}(a,b)$ distribution truncated above at $w$, with pdf $f(x)$:

Note that your parameter $\beta = \frac{1}{b}$, and that you are using the lower incomplete gamma function, whereas I am using the incomplete gamma. Checking with the development version of mathStatica returns the sum $Y = X_1 + X_2$ has pdf $h(y)$:

where:

• Beta[z,a,c] denotes the incomplete beta function $B_z(a,c)$ and
• Gamma[a,w] is the incomplete gamma function $\Gamma(a,w) = \int _w^{\infty } t^{a-1} e^{-t} d t$

which appears to match your own workings. The inclusion of the incomplete Beta function does not impose any practical problem: it is commonly available in any number of software packages. Here is a plot of the pdf $h(y)$ when $a = 1.2$, $b= 3$, and $w = 4$

Monte Carlo check

Here is a quick Monte Carlo check comparing the 'empirical' pdf of the sum of two truncated Gammas (blue wiggly) to the exact theoretical solution above (dashed red curve), for the same parameter values:

All looks good :)

• That term with the square root of $\pi$ is actually reducible to something prettier, so I'm getting somewhere. – Pedro Carvalho Sep 26 '16 at 17:49
• As I pointed out in a comment to the question, obtaining such a result isn't the challenge: the issue is that the number of distinct lines in your conditional formula is equal to $N$ in general. Now if you could derive a rule for each of those $N$ lines, you would have a closed formula. – whuber Sep 26 '16 at 17:58
• Yeah and I think I might be getting something like that. – Pedro Carvalho Sep 27 '16 at 19:37