Probability determined from sum of exponentials numerically stable log probability

Question

In this question, I asked:

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

and I received an answer from Ben:

$$ \sum_{i=1}^{n-1} m_{n,i} \cdot [ \exp (- \lambda_i t) - \exp (- \lambda_n t) ], $$

where the terms $m_{n,j}$ are defined as:

$$ m_{n,i} = \frac{\lambda_i}{\lambda_n} \prod_{j=1 \\ j \neq i}^{n} \frac{\lambda_j}{\lambda_j-\lambda_i}. $$

I have found that when $n$ is moderately big (even above say 10), that the probabilities calculated by the above rule can become negative computationally. I.e. since the $m_{n,i}$ can be both positive and negative, the precision of their calculation can really matter to ensure the probabilities come out positive.

I would like to know if there is a computationally stable way of calculating the log-probability?

Note: I have realised that a part of $m_{n,i}$ is what is known as a Lagrange polynomial:

$$ l_j(0) := \prod_{j=1 \\ j \neq i}^{n} \frac{\lambda_j}{\lambda_j-\lambda_i}, $$

which does have a stable form defined here. I have tried this though and I can still get negative probabilities for moderately sized n. This is partly why I am hoping I can find a more stable expression for the log-probability; it is also more convenient to work with the log-probability because I am doing MCMC that uses this.

Note: it is not straightforward to use some of the clever solutions for stabilising a log of a sum devised here because the $m_{n,i}$ can be negative.

Answer 1

I know you've asked for a "log stable" approach but because in your previous question you mentioned that you have access to Mathematica, there are straightforward steps to take to obtain more precise results.

Below are the results with no adjustments, then rationalizing the values of $\lambda$, and then increasing the working precision of the Plot function.

(* Define probability function *)
p[t_, \[Lambda]_, n_] := 
  Sum[(\[Lambda][[i]]/\[Lambda][[n]])*
     Product[If[i != j, \[Lambda][[j]]/(\[Lambda][[j]] - \[Lambda][[i]]),  1], {j, 1, n}]*
    (Exp[-\[Lambda][[i]]  t] - Exp[-\[Lambda][[n]] t]), {i, 1, n - 1}];

(* Take a random sample of \[Lambda] values *)
n = 30;
SeedRandom[12345];
\[Lambda] = RandomVariate[ChiSquareDistribution[5], n]
(* {1.61194, 4.18716, 3.48782, 6.66293, 3.07288, 8.31177, 4.93457, 
    6.6087, 20.4747, 2.47034, 1.54648, 6.67345, 3.94997, 3.54765, 
    2.65146, 4.54139, 5.62278, 1.57728, 7.04541, 9.54187, 6.49925, 
    2.14693, 6.00469, 1.20595, 6.44199, 8.14599, 4.01718, 3.32529, 
    3.51255, 7.0308} *)

(* Plot results but with just machine precision numbers *) 

Plot[p[t, \[Lambda], n], {t, 0, 20}, PlotRange -> {All, {-0.03, 0.06}},
 Frame -> True, FrameLabel -> {"t", "Probability"}]

(* Rationalize the \[Lambda] values and then plot *)
\[Lambda] = Rationalize[\[Lambda], 0];
Plot[p[t, \[Lambda], n], {t, 0, 20}, PlotRange -> {All, {-0.03, 0.06}},
 Frame -> True, FrameLabel -> {"t", "Probability"}]

(* Increase the working precision in the Plot function *)
Plot[p[t, \[Lambda], n], {t, 0, 20}, PlotRange -> {All, {-0.03, 0.06}},
 Frame -> True, FrameLabel -> {"t", "Probability"}, 
 WorkingPrecision -> 30]

Answer 2

Something similar to my answer above can be done in R using arbitrary precision arithmetic. Here I use the bignum library but there are others that can work just as well or better. I also use the same lambda values as in my other answer.

# Load bignum library
  library(bignum)

# Define function to produce the probability
  p <- function(t, lambda, n) {
    t0 <- as_bigfloat(t)
    total <- as_bigfloat(0)
    for (i in 1:(n-1)) {
        product <- lambda/(lambda - lambda[i])
        product[i] <- 1  # This is the element that is inf in the above vector
        total <- total + (lambda[i]/lambda[n]) * prod(product) *
          (exp(-lambda[i]*t0) - exp(-lambda[n]*t0))
    }
    total
  }

# set lambda values
  n <- 30
  lambda <- as_bigfloat(c(1.61194, 4.18716, 3.48782, 6.66293, 3.07288, 8.31177, 
    4.93457, 6.6087, 20.4747, 2.47034, 1.54648, 6.67345, 3.94997, 3.54765, 
    2.65146, 4.54139, 5.62278, 1.57728, 7.04541, 9.54187, 6.49925, 2.14693, 
    6.00469, 1.20595, 6.44199, 8.14599, 4.01718, 3.32529, 3.51255, 7.0308))
  
# Plot results for values of t from 0 to 20
  t <- c(0:200)/10
  prob <- NULL
  for (i in 1:length(t)) {
      prob[i] <- p(t[i], lambda, n)
  }
  plot(t, prob, type="l", xlab="t", ylab="Probability", las=1,
    ylim=c(-0.03, 0.06), lwd=3)
  lines(c(0, 20), c(0, 0))

I note that for this set of lambda values, values of $t \leq 0.04$ are slightly negative on the order of -1e-37. The library bignum appears to only use up to 50 decimal digits of precision. I suspect you'll need some use of arbitrary precision arithmetic and that tricks using logs just won't cut it. (But I've been wrong many times before.)

In practice for this particular probability function, if you see a smooth curve, you've got enough precision. If you see sections of the cure with wildly fluctuating values (even if none are negative), you need more precision.