Hypoexponential distribution is stuck because of $\lambda_i \neq \lambda_j$? The definition of the hypoexponential distribution (HD) requires that:
$$f(x)=\sum_i^d \left(\prod_{j=1,i\neq j}^{d}\frac{\lambda_j}{\lambda_j-\lambda_i}\right)\lambda_i e^{(-\lambda_ix)},\quad x>0
$$
is well-defined, which imposes that all $\lambda_i \neq \lambda_j$. As such, it seems to impose a significant limitation on the underlying exponentially distributed $X_i$'s in that if, say $X_i$ and $X_j$, have $\lambda_i=\lambda_j$, the distribution is no longer HD, neither can be Erlang, which requires that all $\lambda_i$ are the same!
I'm sure I'm missing something!
A second question: I tried to plot the HD using the syntax below for $\lambda = .1:20$ but not sure if it's the right syntax. Maybe you have experience this this R package ir can suggest another one for working with HD (generalized Erlang)?
library(sdprisk)
x <- seq(0.1, 20, 0.05)
y <- dhypoexp(x, rate = r, log = FALSE)
plot(x, y)


 A: The distribution in question is a mixture of Exponential distributions $\mathcal E(\lambda_i)$ with specified weight, i.e. its density is
$$f(x) = \sum_{i=1}^d p_i\,\lambda_i\exp\{-\lambda_i x\}\quad x\ge 0\tag{1}$$
with
$$p_i=\prod_{j=1\\ j\ne i}^d\dfrac{\lambda_j}{\lambda_i-\lambda_j}\quad i=1,\ldots,d$$
There is therefore a single realisation produced by this density (rather than a sequence of $x_i$'s).
Note that, while
$$\sum_{i=1}^d p_i= \sum_{i=1}^d \prod_{j=1\\ j\ne i}^d\dfrac{\lambda_j}{\lambda_i-\lambda_j}=1$$
(by the Cauchy determinant formula), the weights can be negative, which makes (1) a signed mixture.
Now, if $\lambda_1\approx\lambda_2$ (wlog), starting with the case $d=2$ leads to the Gamma $\mathcal G(2,\lambda_1)$ distribution:
$$\lim_{\epsilon\to 0}\lambda_1(\lambda_1+\epsilon)\dfrac{e^{-\lambda_1x}-e^{-[\lambda_1+\epsilon]x}}{\epsilon}=\lambda_1^2xe^{-\lambda_1x}$$
$-$by L'Hospital's rule$-$as expected since this is the distribution of the sum two iid Exponential $\mathcal E(\lambda_1)$ random variates.
In the general case $d>2$, the undefined term in the density
$$\lim_{\epsilon\to 0}\dfrac{e^{-\lambda_1x}}{\epsilon}\underbrace{\prod_{j=3}^d\dfrac{1}{\lambda_j-\lambda_1}}_\rho
-\dfrac{e^{-[\lambda_1+\epsilon]x}}{\epsilon}\prod_{j=3}^d\dfrac{1}{\lambda_j-\lambda_1-\epsilon}$$
is equal to
$$\lim_{\epsilon\to 0}\frac{\rho e^{-\lambda_1x}}{\epsilon\prod_{j=3}^d[\lambda_j-\lambda_1-\epsilon]}
\left\{\prod_{j=3}^d[\lambda_j-\lambda_1-\epsilon]-e^{-\epsilon x}\rho^{-1}\right\}
=\rho e^{-\lambda_1x}\left\{ x - \sum_{j=3}^d (\lambda_j-\lambda_1)^{-1}\right\}$$
by L'Hospital's rule.
A: 
The definition of the hypoexponential distribution (HD) requires that:
$$f(x)=\sum_i^d \left(\prod_{j=1,i\neq j}^{d}\frac{\lambda_j}{\lambda_j-\lambda_i}\right)\lambda_i e^{(-\lambda_ix)},\quad x>0
$$

that expression is only a special case when $\lambda_i \neq \lambda_j$. It is not a requirement for the hypoexponential distribution in general.
The more general hypo-exponential distribution can be expressed as a phase-type distribution
$$f(x) = -\boldsymbol{\alpha}e^{x \boldsymbol{\Theta}} \boldsymbol{\Theta} \mathbf{1} $$
With $$\boldsymbol{\alpha} = (1,0,0,\dots,0,0)$$ and
$$\boldsymbol{\Theta} = \begin{bmatrix} 
-\lambda_1 & \lambda_1 & 0 & \dots & 0 & 0 \\
0 & -\lambda_2 &  \lambda_2 & \dots & 0 & 0 \\
\vdots & \ddots & \ddots & \ddots & \ddots & \vdots\\
0 & 0 & \ddots & -\lambda_{d-2} & \lambda_{d-2} & 0\\
0 & 0 & \dots & 0 & -\lambda_{d-1} & \lambda_{d-1}\\
0 & 0 & \dots & 0 & 0 & -\lambda_{d}\\
\end{bmatrix}$$
That involves the exponentiation of a matrix. And that can be approximated with a sum using a Taylor series.
You can see the matrix as modelling a sort of Markov chain process with non-discrete time steps. A sum of exponential distributed variables is like waiting for several consecutive transitions whose waiting time is each exponential distributed. Those transitions relate to the Markov chain.
Not only when $\lambda_i \neq \lambda_j$ but also when $\lambda_i = \lambda_j$ then the formula can give problems as demonstrated in this question: Why my cdf of the convolution of n exponential distribution is not in the range(0,1)?
