Sums of exponentials joint probability

Question

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

Answer 1

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

Answer 2

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. $$