Probability density function of a summation of continuous random variables Let $Z_{i} = \tau + X_{i}$, where $X_{i}$ is a exponential random variable ($X_i \sim \varepsilon(\lambda)$),  $0<\tau, \lambda < \infty$
Assume $X_i$ are independent random variables.
Suppose $T_{n} = Z_{1} + Z_{2} + \cdots + Z_{n}$, find the PDF of $T_{n}$, $h(t_n)$.
My approach:


*

*$T_{n} = \sum_{i=1}^{n}\left ( \tau + X_{i}   \right )$

*$T_{n} = n\tau+\sum_{i=1}^{n}\left (X_{i}  \right )$

*Let $Y$ be $\sum_{i=1}^{n}\left (X_{i}  \right )$, I find the PDF of $Y$, $g(y)$

*$T_{n} = n\tau+ Y$

*$Y = T_{n} - n\tau$

*Using the Jacobian of transformation, I find PDF of $T_{n}$: $h(t_n) = \left | \frac{\mathrm{d} }{\mathrm{d} t_{n}}(t_{n} - n\tau) \right |\times g(t_{n} - n\tau)$
I am stuck on step 3, as I have no idea how to find the PDF of $\sum_{i=1}^{n}\left (X_{i}  \right )$. Any ideas?
By the way, if there's any wrong with my approach, please say so. 
Thanks a lot.
 A: The main thing you need to know here is that the probability distribution of a sum of independent exponentially distributed random variables with the same scale parameter (or, equivalently, with the same rate parameter) has a gamma distribution, whose "shape parameter" is the number of random variables being added.  That can be proved by at least (estimating conservatively) two different methods.
(I don't know whether you intended $\lambda$ to be the scale parameter and $1/\lambda$ the rate parameter or vice-versa.  It's usually a good idea to specify that when writing about such things, but in this case the paragraph above stands either way.)
PS: Here is a proof.  The gamma distribution is the distribution whose density is
$$
\text{constant}\cdot x^{n-1} e^{-x/\lambda} \text{ on the interval }0\le x<\infty
$$
where $\lambda>0$ is the scale parameter and $n>0$ is sometimes called the shape parameter because it determines the shape of the graph of the density function.
Notice that an exponential distribution is a gamma distribution with $n=1$.
Suppose two independent random variables $X$ and $Y$ have gamma distributions with the same scale parameter $\lambda>0$ and shape parameters $m$ and $n$.  Then their joint density is
$$
\text{constant}\cdot x^{n-1} e^{-x\lambda} y^{m-1} e^{-y/\lambda}.
$$
in the quadrant $0\le x$, $0\le y$.  Then
$$
\Pr(X+Y\le w) = \iint\limits_{\{(x,y)\,\mid\,x+y\,\le\,w\}} \text{constant}\cdot x^{n-1} e^{-x\lambda} y^{m-1} e^{-y/\lambda} \, dx\,dy.\tag1
$$
We will use the substitution
\begin{align}
u & = x+y, \\
v & = y/(x+y), \\[10pt]
\text{hence }
u\,du\,dv & = dx\,dy, \\
x & = u(1-v), \\
y & = uv.
\end{align}
Then $0\le u<\infty$ and $0\le v\le 1$.  The integral $(1)$ becomes
$$
\begin{align}
& \phantom{={}}\text{constant}\cdot\int_0^w \int_0^1 u^{n-1}(1-v)^{n-1} u^{m-1} v^{m-1} e^{-u/\lambda}   (u\,dv \,du) \\
& = \text{constant}\cdot\int_0^1 (1-v)^{n-1}v^{m-1} \, dv \cdot \int_0^w u^{n+m-1} e^{-u/\lambda} \,du \\
& = \text{constant} \cdot \int_0^w u^{n+m-1} e^{-u/\lambda} \,du.
\end{align}
$$
The integral involving $v$ is a "constant" since it does not depend on $w$.  What we get is therefore a gamma distribution with the same scale parameter $\lambda$ and shape paramter $n+m$.
Thus we have: The distribution of the sum of two independent gamma-distributed random variables with the same scale parameter is another gamma-distributed random variable with the same scale paremeter, and the shape parameter is just the sum of the two shape parameters.
Now apply this repeatedly adding one more independent gamma-distributed random variable at each step.
A: \begin{align*}
    Y &= X_1 +  X_2.\\
    {\rm cdf}_Y(y) &=\Pr\{X_1+X_2\leq y\},\\
        &\overset{a)}{=}\int_{-\infty}^\infty\int_{-\infty}^{y-x_1}{\rm pdf}_{X_1,X_2}(x_1,x_2)\,{\rm d} x_2 {\rm d} x_1,\\
    {\rm pdf}_Y(y)&=\frac{{\rm d}}{{\rm d} y}{\rm cdf}_Y(y),\\
    &\overset{b)}{=}\int_{-\infty}^\infty \left[ {\rm pdf}_{X_1,X_2}(x_1,y-x_1)\frac{{\rm d} (y-x_1)}{{\rm d} y} 
    -{\rm pdf}_{X_1,X_2}(x_1,-\infty)\frac{{\rm d} (-\infty)}{{\rm d}y}+0\right]\,{\rm d} x_1,\\
    &=\int_{-\infty}^\infty {\rm pdf}_{X_1,X_2}(x_1,y-x_1)\,{\rm d} x_1,\\
    &\overset{c)}{=}\int_{-\infty}^{\infty}{\rm pdf}_{X_1}(x_1){\rm pdf}_{X_2}(y-x_1)\,{\rm d} x_1.
\end{align*}
a) -- Integrate over sufficient 2D area, where $x_1 + x_2 \leq y$.
b) -- for derivative under integral see
https://en.wikipedia.org/wiki/Leibniz_integral_rule
c) -- assume RV $X_1$ and $X_2$ independent.
so probability density function of sum of random variables is convolution of their pdf-s.
Define $Y_n = \sum_{i=1}^n X_i$.
We will compute ${\rm pdf}_{Y_n}$ for some $n$ and will observe the pattern, formula follows.
Note: exponential RV can have nonnegative values only, so their sum also has nonnegative values only.
For a nonnegative $x$
$${\rm pdf}_{X}(x):={\rm pdf}_{X_i}(x)=\lambda {\rm e}^{-\lambda x}.$$
\begin{align}
{\rm pdf}_{Y_2}(y)&=\int_{-\infty}^{\infty}{\rm pdf}_X(x){\rm pdf}_X(y-x)\,{\rm d}x.
\end{align}
The integration limits srunk:
\begin{align}
x\geq 0 &\quad \text{and} \quad y-x\geq 0,\\
0\leq &x \leq y.
\end{align}
\begin{align}
{\rm pdf}_{Y_2}(y)&=\int_0^y\lambda {\rm e}^{-\lambda x}\lambda {\rm e}^{-\lambda(y-x)}\,{\rm d}x,\\
  &=\lambda^2{\rm e}^{-\lambda y}\int_0^y 1\,{\rm d}x = \lambda^2 {\rm e}^{-\lambda y}y.\\
{\rm pdf}_{Y_3}(y)&=\int_{-\infty}^{\infty}{\rm pdf}_{Y_2}(x){\rm pdf}_X(y-x)\,{\rm d}x,\\
  &=\int_0^y\lambda^2{\rm e}^{-\lambda x}x\lambda {\rm e}^{-\lambda (y-x)}\,{\rm d}x,\\
  &=\lambda^3{\rm e}^{-\lambda y}\int_0^y x\,{\rm d}x = \lambda^3{\rm e}^{-\lambda y}\frac{y^2}{2}.\\
{\rm pdf}_{Y_4}(y)&=\int_{-\infty}^{\infty}{\rm pdf}_{Y_3}(x){\rm pdf}_X(y-x)\,{\rm d}x,\\
  &=\int_0^y\lambda^3{\rm e}^{-\lambda x}\frac{x^2}{2}\lambda{\rm e}^{-\lambda(y-x)}\,{\rm d}x,\\
  &=\lambda^4{\rm e}^{-\lambda y}\int_0^y\frac{x^2}{2}\,{\rm d} x=\lambda^4{\rm e}^{-\lambda y}\frac{y^3}{2\!\cdot\! 3}.
\end{align}
Observing the pattern, we have:
$$ {\rm pdf}_{Y_n}(x) = \lambda^n \frac{x^{n-1}}{(n-1)!}{\rm e}^{-\lambda x}.$$
Appendix -- relation to the Poisson distribution
Using
$$\int x^n{\rm e}^{-x}\,{\rm d}x = -n!{\rm e}^{-x}\sum_{k=0}^n\frac{x^k}{k!}.$$
we get 
$${\rm cdf}_{Y_n}(x) = 1-\sum_{k=0}^{n-1}\frac{(\lambda x)^k}{k!}{\rm e}^{-\lambda x}.$$
We can define elementary event $A_k$ -- that exactly $k$ events occur in time interval of length $y$, i.e. $\lambda y$ is an average number of events in the unit interval. ${\rm cdf}_{Y_n}$ can be than written as:
\begin{align}
{\rm cdf}_{Y_n}(y)&={\rm Pr}[Y_n \leq y] = 1-{\rm Pr}[\lor_{i=0}^{n-1} A_i],\\
  &={\rm Pr}[\neg \lor_{i=0}^{n-1} A_i] = {\rm Pr}[\land_{i=0}^{n-1}(\neg A_i)].
\end{align}
I.e. the probability that time between two subsequent events will be less than $y$ is equal to the probability that at least $n$ events occur in time-interval of length $y$.
I don't understand the interpretation. Maybe, somebody will find a better one.
