# Sums of exponentials joint probability

If we have that: $$\tau_i \overset{\text{independent}}{\sim} \exp(\lambda_i)$$, for $$i=1,2,3,...,n$$, where $$\lambda_i\neq \lambda_j, \forall i\neq j$$ then I would like to find a general form for the probability:

$$\text{Pr}(\sum_{i=1}^{n-1} \tau_i \leq t, \sum_{i=1}^{n} \tau_i > t)$$

Note: I have managed to find some specific results in Mathematica but the expressions seem to become unwieldy as $$n$$ grows. For example, if $$n=2$$, then the above probability becomes:

$$\frac{\left(e^{-t \lambda_1} - e^{-t \lambda_2}\right) \lambda_1}{\lambda_2 - \lambda_1}.$$

For $$n=3$$, it becomes:

$$\frac{e^{-t (\lambda_1 + \lambda_2 + \lambda_3)} \lambda_1 \lambda_2 \left(e^{t (\lambda_1 + \lambda_2)} (\lambda_1 - \lambda_2) + e^{t (\lambda_2 + \lambda_3)} (\lambda_2 - \lambda_3) + e^{t (\lambda_1 + \lambda_3)} (-\lambda_1 + \lambda_3)\right)}{(\lambda_1 - \lambda_2) (\lambda_1 - \lambda_3) (\lambda_2 - \lambda_3)}$$

Note: I realise that:

$$\sum_{i=1}^{n} \tau_i \sim \text{Hypoexponential}(\lambda_1,\lambda_2,...,\lambda_n),$$

but I am not quite sure how to use this to derive a general result (or even a recursive result or good approximation).

As you point out in your note, the sum of exponential random variables with different rate parameters follows a hypoexponential distribution. This distribution has density function:

$$\text{HypoExp}(r|\lambda_1,...,\lambda_{k}) = \sum_{i=1}^{k} \bigg( \prod_{j=1 \\ j \neq i}^{k} \frac{\lambda_j}{\lambda_j-\lambda_i} \bigg) \lambda_i \exp(- \lambda_i r).$$

Now, applying the law of total probability gives the resulting probability (which I will call the HEProb, short for hypoexponential probability):

\begin{align} \text{HEProb}(t|\boldsymbol{\lambda}) &\equiv \mathbb{P} \bigg( \sum_{i=1}^{n-1} \tau_i \leqslant t, \sum_{i=1}^{n} \tau_i > t \bigg) \\[6pt] &= \int \limits_0^t \mathbb{P} \bigg( \sum_{i=1}^{n} \tau_i > t \bigg| \sum_{i=1}^{n-1} \tau_i = r \bigg) p \bigg( \sum_{i=1}^{n-1} \tau_i = r \bigg) \ dr \\[6pt] &= \int \limits_0^t \mathbb{P}(\tau_n > t-r) \cdot \text{HypoExp}(r|\lambda_1,...,\lambda_{n-1}) \ dr \\[6pt] &= \int \limits_0^t \exp(- \lambda_n (t-r)) \cdot \bigg( \sum_{i=1}^{n-1} \bigg( \prod_{j=1 \\ j \neq i}^{n-1} \frac{\lambda_j}{\lambda_j-\lambda_i} \bigg) \lambda_i \exp(- \lambda_i r) \bigg) \ dr \\[6pt] &= \sum_{i=1}^{n-1} \bigg( \prod_{j=1 \\ j \neq i}^{n-1} \frac{\lambda_j}{\lambda_j-\lambda_i} \bigg) \lambda_i \exp (- \lambda_n t) \int \limits_0^t \exp ((\lambda_n - \lambda_i) r) \ dr \\[6pt] &= \sum_{i=1}^{n-1} \bigg( \prod_{j=1 \\ j \neq i}^{n-1} \frac{\lambda_j}{\lambda_j-\lambda_i} \bigg) \lambda_i \exp (- \lambda_n t) \cdot \frac{1}{\lambda_n - \lambda_i} [ \exp ((\lambda_n - \lambda_i) t) - 1] \\[6pt] &= \frac{1}{\lambda_n} \sum_{i=1}^{n-1} \bigg( \prod_{j=1 \\ j \neq i}^{n} \frac{\lambda_j}{\lambda_j-\lambda_i} \bigg) \lambda_i [ \exp (- \lambda_i t) - \exp (- \lambda_n t) ] \\[6pt] &= \sum_{i=1}^{n-1} m_{n,i} \cdot [ \exp (- \lambda_i t) - \exp (- \lambda_n t) ], \\[6pt] \end{align}

where we use the multipliers $$m_1,...,m_{n-1}$$ defined by:

$$m_{n,i} \equiv m_{n,i}(\boldsymbol{\lambda}) \equiv \frac{\lambda_i}{\lambda_n} \prod_{j=1 \\ j \neq i}^{n} \frac{\lambda_j}{\lambda_j-\lambda_i} = \bigg( \prod_{j=1}^{n-1} \lambda_j \bigg) \Bigg/ \bigg( \prod_{j=1 \\ j \neq i}^{n} (\lambda_j-\lambda_i) \bigg). \\[6pt]$$

This gives you an explicit closed form for the probability of interest. It is notable that you can precompute the multipliers from the rate parameters without use of the time variable, which aids in computation of the probability function. The forms of these multipliers for $$n=2,3,4,5$$ are given by:

\begin{align} m_{2,1} &= \frac{\lambda_1}{\lambda_2-\lambda_1}, \\[6pt] m_{3,1} &= \frac{\lambda_1 \lambda_2}{(\lambda_2-\lambda_1)(\lambda_3-\lambda_1)}, \\[6pt] m_{3,2} &= \frac{\lambda_1 \lambda_2}{(\lambda_1-\lambda_2)(\lambda_3-\lambda_2)}, \\[6pt] m_{4,1} &= \frac{\lambda_1 \lambda_2 \lambda_3}{(\lambda_2-\lambda_1)(\lambda_3-\lambda_1)(\lambda_4-\lambda_1)}, \\[6pt] m_{4,2} &= \frac{\lambda_1 \lambda_2 \lambda_3}{(\lambda_1-\lambda_2)(\lambda_3-\lambda_2)(\lambda_4-\lambda_2)}, \\[6pt] m_{4,3} &= \frac{\lambda_1 \lambda_2 \lambda_3}{(\lambda_1-\lambda_3)(\lambda_2-\lambda_3)(\lambda_4-\lambda_3)}, \\[6pt] m_{5,1} &= \frac{\lambda_1 \lambda_2 \lambda_3 \lambda_4}{(\lambda_2-\lambda_1)(\lambda_3-\lambda_1)(\lambda_4-\lambda_1)(\lambda_5-\lambda_1)}, \\[6pt] m_{5,2} &= \frac{\lambda_1 \lambda_2 \lambda_3 \lambda_4}{(\lambda_1-\lambda_2)(\lambda_3-\lambda_2)(\lambda_4-\lambda_2)(\lambda_5-\lambda_2)}, \\[6pt] m_{5,3} &= \frac{\lambda_1 \lambda_2 \lambda_3 \lambda_4}{(\lambda_1-\lambda_3)(\lambda_2-\lambda_3)(\lambda_4-\lambda_3)(\lambda_5-\lambda_3)}, \\[6pt] m_{5,4} &= \frac{\lambda_1 \lambda_2 \lambda_3 \lambda_4}{(\lambda_12-\lambda_4)(\lambda_2-\lambda_4)(\lambda_34-\lambda_4)(\lambda_5-\lambda_4)}. \\[6pt] \end{align}

Computational implementation: This probability can be programmed as a vectorised function in R using the code below. This HEProb function allows the user to input the vector of rates and a vector of values for the time variable t. The function returns the corresponding HEProb values (either as probabilities or log-probabilities depending on the log.p option) for those time variables using the stipulated parameters.

#Set function
HEProb <- function(t, rates, log.p = FALSE) {

#Check input t
if (!is.vector(t))          stop('Error: Input t should be a vector')
if (!is.numeric(t))         stop('Error: Input t should be a numeric vector')

#Check input rates
if (!is.vector(rates))      stop('Error: Input rates should be a vector')
if (!is.numeric(rates))     stop('Error: Input rates should be a numeric vector')
if (min(rates) <= 0)        stop('Error: Input rates should be positive')
n <- length(unique(rates))
if (n < 1)                  stop('Error: Input rates should have at least one value')
if (length(rates) != n)     stop('Error: Input rates should be distinct')

#Check input log.p
if (!is.vector(log.p))      stop('Error: Input log.p should be a vector')
if (!is.logical(log.p))     stop('Error: Input log.p should be a logical value')
if (length(log.p) != 1)     stop('Error: Input log.p should be a single logical value')

#Deal with trivial special case where n = 1
if (n == 1) {
if (log.prob) {
return(-rates[1]*pmax(t,0)) } else {
return(exp(-rates[1]*pmax(t,0))) } }

#Deal with remaining cases where n > 1
#Compute multipliers
MULTMAT <- matrix(1, nrow = n-1, ncol = n)
MULTVEC <- rep(0, n-1)
for (i in 1:(n-1)) {
for (j in 1:n) {
if  (i != j) {
MULTMAT[i, j] <- rates[j]/(rates[j]-rates[i]) } }
MULTVEC[i] <- (rates[i]/rates[n])*prod(MULTMAT[i, ]) }

#Compute probabilities
T <- length(t)
PROBS <- rep(0, T)
TERMS <- matrix(1, nrow = n-1, ncol = T)
for (i in 1:(n-1)) {
for (r in 1:T) {
TERMS[i, r] <- exp(-rates[i]*t[r]) - exp(-rates[n]*t[r]) } }
for (r in 1:T) {
PROBS[r] <- sum(MULTVEC*TERMS[, r]) }

#Return output
if (log.p) { log(PROBS) } else { PROBS } }


Here is an example of this function using the parameters $$\boldsymbol{\lambda} = (3, 2, 5, 4)$$ showing the time variable ranging over $$t=0,...,4$$ in increments of $$\Delta t = 0.1$$. As can be seen, the function is quasi-concave with respect to the time variable.

RATES <- c(3, 2, 5, 4)
PROBS <- HEProb(t = 0:40/10, rates = RATES)
NAMES <- c('0', rep('', 9), '1', rep('', 9), '2', rep('', 9), '3', rep('', 9), '4')
names(PROBS) <- NAMES
barplot(PROBS, col = 'blue', ylim = c(0, 0.2),
xlab = 'Time', ylab = 'Probability')


• thank you so much for this comprehensive answer. Commented Aug 12 at 8:00

This problem gets easier if you think about your process in terms of continuous-time Markov chains. You are looking for the probability of being in state $$n$$ after time $$t$$ for the CTMC with rate (or generator) matrix $$Q = \begin{bmatrix} -\lambda_1 & \lambda_1\\ & -\lambda_2 & \lambda_2\\ & & \ddots & \ddots\\ & & & -\lambda_n & \lambda_n\\ & & & & 0 \end{bmatrix}$$ starting from initial distribution $$\pi_0 = e_1$$. So the probability can be written as $$p = e_{n}^\top \exp(Qt)e_1$$.

Note also that you can omit the last row and column, since the matrix is upper triangular, and reduce to $$p = e_n^T\exp\left(t \begin{bmatrix} -\lambda_1 & \lambda_1\\ & -\lambda_2 & \lambda_2\\ & & \ddots & \ddots\\ & & & -\lambda_n\\ \end{bmatrix}\right)e_1.$$

• Thanks @Federico. I actually was originally using this method but I wondered if the approach based on the hypoexponential distribution might speed things up (especially when n is big). It does speed things up when n is small. But when n gets even moderately large, there's numerical underflow in the multiplier terms which causes issues. So, ultimately, it looks like the CTMC approach will be better for my problem. (Sorry, I should have said I was using a CTMC approach previously.) Commented Aug 12 at 12:26
• For anyone interested, I tried stabilising the weights as per here but this did not help either: en.wikipedia.org/wiki/Lagrange_polynomial Commented Aug 12 at 13:12
• What do you mean with "stabilizing"? Anyhow, those weights require lots of finite differences, which are infamous for being unstable, so I am not surprised. It seems a tricky problem to solve numerically, especially for large $n$ and nearby-but-not-equal $\lambda_i$. Commented Aug 12 at 17:56
• Numerically, algorithms for the matrix exponential might be the best choice. There are specialized algorithms for exponentials of matrices with the sign structure of that $Q$, for instance doi.org/10.1137/120894294 . They give explicit error bounds from above and below. Commented Aug 12 at 17:58