Does a hypo-exponential distribution become deterministic as the number of phases approach $\infty$?

Question

A hypo-exponential distribution is a more general case of the Erlang distribution. The former has different transition rates $\lambda_i,$ $i=1,2,\dots,N$ while the latter uses the same transition rates $\lambda_i = \lambda,$ $i=1,2,\dots,N$. Furthermore, it is easy to show that the Erlang distribution becomes deterministic as the number of phases increases to $\infty$. This requires making use of its squared coefficient of variation $CV(X_N) = \operatorname{Var}(X_N)/(\mathbb{E}[X_N])^2$ where $\operatorname{Var}(X_N) = N/\lambda^2$ and $\mathbb{E}[X_N] = N/\lambda$. Hence, $CV(X_N) = 1/N$ such that $\lim_{N\to \infty } 1/N = 0$. Thus the distribution becomes deterministic.

I conjecture that the same happens to the hypo-exponential distribution. If this conjecture is correct, I would like to show this similarly as was done with the Erlang distribution. I start with the expected value $$ \mathbb{E}[X_N] = \sum_{n=1}^N\frac{1}{ \lambda_n} $$ and variance $$ \operatorname{Var}(X_N) = \sum_{n=1}^N \frac{1}{\lambda_n^2} $$ where it is understood that $\lambda_n > 0, $ $n=1,2,\dots,N$. The coefficient of variation thus follows as $$ CV(X_N) = \frac{\sum_{n=1}^N\frac{1}{ \lambda_n^2}}{\left( \sum_{n=1}^N\frac{1}{ \lambda_n} \right)^2} $$ where the denominator can be expanded $$ CV(X_N) = \frac{\sum_{n=1}^N\frac{1}{ \lambda_n^2}}{\sum_{n=1}^N\frac{1}{ \lambda_n^2} + \sum_{n=1}^N \frac{1}{\lambda_n}\sum_{\substack{j=1\\j\neq n}}^N \frac{1}{\lambda_j}}. $$ This makes it clear that $CV(X_N) < 1$. The next step then results in

$$ CV(X_N) = \frac{1}{1 + \frac{\sum_{n=1}^N \frac{1}{\lambda_n}\sum_{\substack{j=1\\j\neq n}}^N \frac{1}{\lambda_j}}{\sum_{n=1}^N \frac{1}{\lambda_n^2}} } $$ or more succinctly $$ CV(X_N) = \frac{1}{1+ \Psi(X_N)} $$ where $$ \Psi(X_N) = \frac{\sum_{n=1}^N \sum_{\substack{j=1\\j\neq n}}^N \frac{1}{\lambda_n\lambda_j}}{\sum_{n=1}^N \frac{1}{\lambda_n^2}}. $$ Now if we were to show that $\Psi(X_N)$ is a monotonically increasing function of $N$ then we are sorted. Unfortunately, I have not been able to do so yet. I will probably resort to plotting $\Psi(X_N)$ as to determine this. However, that is not rigorous at all.

Can you please help show me if it is a monotonically increasing function?

Answer 1

I am posting an answer to my own question. I have chosen not to place it as an edit in the question as I would like this post to reflect that there is at least one answer. Furthermore, the question is still open to alternative solutions.

TL;DR: $\Psi(X_N)$ is not monotonically increasing in $N$ but seems to be increasing in $N$ or so simulations would suggest. As such, one cannot conclude that the hypo-exponential distribution becomes deterministic as one adds infinitely more phases. However, in practice we do observe that $CV(X_N)$ decreases close to zero.

Long answer: we are going to have a look at the following: $$ \Delta_N := \Psi(X_{N+1}) - \Psi(X_N). $$ The goal is to determine/characterise the sign of $\Delta_N$. Note that in the above, $\Psi(X_N)$ makes use of $\vec{\lambda}_{N} = [\lambda_1,\dots,\lambda_N]^\top$ while $\Psi(X_N)$ shares the same $N$ values but with the addition of $\lambda_{N+1}$, i.e. $\vec{\lambda}_{N+1} = [\lambda_1,\dots,\lambda_N,\lambda_{N+1}]^\top$. Furthermore, we need to define the combination function $C(\vec{\lambda},n)$ that returns all $n$ length combinations $\vec{c}$ that can be made from the elements of $\vec{\lambda}$. These combinations will be returned as $\vec{c}$. For example $C([\lambda_1,\lambda_2,\lambda_3],2) = \{ [\lambda_1,\lambda_2],[\lambda_1,\lambda_3],[\lambda_2,\lambda_3] \}$

So $\Psi(X_N)$ changes due to the addition of $\lambda_{N+1}$. Hence, it makes sense to solve for $\lambda_{N+1}$ in $\Delta_N = 0$. Call this solution $\lambda_{N+1}^*$. It follows as: $$ \lambda_{N+1}^* = \frac{\left( \sum_{\vec{c} \in C\left(\vec{\lambda_N},N-2\right)} \prod_{\lambda_i \in \vec{c}} \lambda_i\right) \left( \prod_{i=1}^N \lambda_i^2 \right) }{\left( \sum_{\vec{c} \in C\left(\vec{\lambda_N},N-1\right)} \prod_{\lambda_i \in \vec{c}} \lambda_i\right)\left( \sum_{\vec{c} \in C\left(\vec{\lambda_N},N-1\right)} \prod_{\lambda_i \in \vec{c}} \lambda_i^2\right)}, \quad N \geq 3. $$ If $N=2$ then replace the first product of the numerator with 1. If $N=1$ then the solution is trivial/obvious as $\lambda_{2}^* = 0$. This solution was obtained through the help of Sympy.

We now have the following three cases: $$ \Delta_N \in \begin{cases} \mathbb{R}_{>0}, & \lambda_{N+1} > \lambda_{N+1}^*\\ \{0\}, & \lambda_{N+1} = \lambda_{N+1}^*\\ \mathbb{R}_{<0}, & \lambda_{N+1} < \lambda_{N+1}^*\end{cases}. $$ So depending on the value of $\lambda_{N+1}$ relative to $\lambda_{N+1}^*$, it becomes clear that $\Psi(X_N)$ can increase or decrease in $N$. However, we need it to be monotonically increasing such that $\lim_{N\to \infty }\Psi(X_N) = \infty$ as this would prove that $\lim_{N\to \infty }CV(X_N) = 0$.

There might of course exist special cases. The first one that jumps to mind is setting all the $\lambda_i$ to a single value $\lambda$. This of course gives the Erlang distribution for which we know $\lim_{N\to \infty }CV(X_N) = 0$. I would like to just verify our solution $\lambda_{N+1}^*$ through showing that for the Erlang distribution we have $$ \lambda_{N+1}^* = \frac{\lambda}{N}. $$ Hence, $\lambda_{N+1} = \lambda \geq \lambda_{N+1}^*$, $N=1,2,3,\dots$ such that $\Delta_N \geq 0$, $N=1,2,3,\dots$

Two factors that I have considered are:

1. The ordering of $\vec{\lambda}$ as increasing, decreasing or unordered (randomly ordered).
2. The range/dispersion of $\vec{\lambda}$. As noted with the Erlang distribution, if the range is zero then we have guaranteed convergence to a deterministic distribution.

In regard to the first point, increasing order seems to produce increasing behaviour for $\Psi(X_N)$ even if the range is large. This is shown in the figure below. I guess this is because it increases the likelihood that $\lambda_{N+1} > \lambda_{N+1}^*$.

With decreasing order, it would seem that the function inevitably starts to decrease. A wider range seems to just accelerates this. Perhaps a small enough range can ensure that it is always increasing.

Here is what the $CV(X_N)$ looks like for the wide range, decreasing order case. Have a look at its tail. It starts to increase thus ruling out convergence to a deterministic distribution. With that said, the $CV(X_N)$ value is tiny such that one could consider it near deterministic at this point. Now we get to the random ordering. I present here only the plots for a wide range as we have already seen that small ranges lead to higher changes of $\Psi(X_N)$ increasing.

As we've seen, this is not always the whole story. Here is what the resultant $CV(X_N)$ plot looks like.

Again, it does decrease to a point where we can consider it deterministic. We must thus conclude that the hypo-exponential distribution is not guaranteed to converge to a deterministic distribution as more phases are added. However, it typically does become more deterministic; especially if increasing $\lambda_i$ values are added.

Python code Here is some code if you would like to replicate some of the plots.

import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm

def psi(lambdas):
    num = 0
    den = 0
    for lambda0 in lambdas:
        for lambda1 in lambdas:
            if lambda0 == lambda1:
                continue
            num += 1/(lambda0*lambda1)
        den += 1/lambda0**2
    return num/den
    
def cv(lambdas):
    lambdas = np.array(lambdas)
    var = (1/lambdas**2).sum()
    mean = (1/lambdas).sum()
    return var/mean**2
    
    
def experiment(N,func,
               lambda_range=(0,100),
               seed=123):
    np.random.seed(seed)
    lambdas = np.random.uniform(0.9,100,N) # no order
    lambdas = np.sort(lambdas)  # increasing
    lambdas = lambdas[::-1]  # decreasing
    results = []
    for n in tqdm(range(2,N+1)):
        results.append(func(lambdas[:n]))
    fig,ax = plt.subplots(1,1)
    ax.plot(range(2,N+1),results)
    ax.set_xlabel('$N$')
    return fig,ax
    
def psi_experiment(N,lambda_range=(0,100),seed=123):
    fig,ax = experiment(N,psi,lambda_range,seed)
    ax.set_ylabel('$\Psi(X_N)$')
    return fig,ax
    
def cv_experiment(N,lambda_range=(0,100),seed=123):
    fig,ax = experiment(N,cv,lambda_range,seed)
    ax.set_ylabel('$CV(X_N)$')
    return fig,ax
    
    
    
if __name__ == '__main__':
    N = 300
    lambda_range = (0,10)
    seed = 123
    test_psi = True
    if test_psi:
        fig,ax = psi_experiment(N,lambda_range,seed)
    else:
        fig,ax = cv_experiment(N,lambda_range,seed)