# Probability determined from sum of exponentials numerically stable log 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)$$

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

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]



• Thanks @JimB. I'm afraid I'm only using Mathematica to work things out to transfer into R and Stan though! Commented Aug 12 at 20:08
• @kjetilbhalvorsen I appreciate any edits someone does to any of my answers but in this case I just can't tell what was modified. It looks like all of the original coded was deleted but then put all back in. I must be missing something obvious.
– JimB
Commented Aug 12 at 21:32
• I added code fence, look after mathematica. That effects layout and code coloring Commented Aug 12 at 21:34
• More on what @kjetilbhalvorsen did: Language-specific code blocks versus default code blocks. Commented Aug 13 at 1:36
• @JimB, if you are not sure what was edited, you can always check the "Side-by-side markdown" which would show you the minuscule change that would otherwise be not noticed. Commented Aug 13 at 5:34

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.