Sum of truncated Gammas

Question

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.

Answer 1

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 :)