Numerically stable log-odds of intersection of independent events Suppose we have independent events $E_1,\dots,E_n$, and define $E = E_1 \land \cdots \land E_n$ to be the conjunction (logical-and) of these events.  Let $p_i=\Pr[E_i]$, $p=\Pr[E]$.
If we are given the probabilities $p_1,\dots,p_n$, it is easy to compute $p$ via the formula $p=p_1 \cdots p_n$.
But what if we are working in log-odds space, instead of probability space?  Define $\sigma(\ell)=e^\ell/(e^\ell+1)$ to be the function that maps from log-odds to probability, and $\sigma^{-1}(p)=\log(p/(1-p))$ the function that maps in the opposite direction.  Let $\ell_i=\sigma^{-1}(p_i)$ and $\ell=\sigma^{-1}(p)$.
If we are given the log-odds $\ell_1,\dots,\ell_n$ for the events $E_1,\dots,E_n$, is there a numerically-stable formula to compute the log-odds $\ell$ for $E$?
Of course I know that we can convert to probability space, multiply, and convert back: $\ell=\sigma^{-1}(\sigma(\ell_1) \cdots \sigma(\ell_n))$.  However, I am not sure whether this will be numerically stable.

For dealing with probabilities on a log-scale, it is often very useful to use the hyperbolic functions:
$$\begin{align}
\text{log1pexp}(x) 
&= \log(1 + \exp(x)), \\[12pt]
\text{log1mexp}(x) 
&= \log(1 - \exp(-x)). \\[6pt]
\end{align}$$
(These functions are programmed into most statistical software; in R they are available in the VGAM package as log1pexp and log1mexp.)  Using these functions you can write the joint probability and its negation as:
$$\begin{align}
p 
&= p_1 \times \cdots \times p_n \\[10pt]
&= \frac{e^{\ell_1}}{1+e^{\ell_1}} \times \cdots \times \frac{e^{\ell_n}}{1+e^{\ell_n}} \\[6pt]
&= \frac{\exp(\sum \ell_i)}{\prod (1+\exp(\ell_i))} \\[6pt]
&= \frac{\exp(\sum \ell_i)}{\exp(\sum \log(1+\exp(\ell_i)))} \\[6pt]
&= \frac{\exp(\sum \ell_i)}{\exp(\sum \text{log1pexp}(\ell_i))}, \\[14pt]
1-p 
&= \frac{\exp(\sum \text{log1pexp}(\ell_i) - \exp(\sum \ell_i)}{\exp(\sum \text{log1pexp}(\ell_i))}. \\[6pt]
\end{align}$$
You can then write the log-odds as:
$$\begin{align}
\ell 
&= \log \Bigg( \frac{p}{1-p} \Bigg) \\[6pt]
&= \log \Bigg( \frac{\exp(\sum \ell_i)}{\exp(\sum \text{log1pexp}(\ell_i)) - \exp(\sum \ell_i)} \Bigg) \\[8pt]
&= \log \Bigg( \frac{1}{\exp(\sum \text{log1pexp}(\ell_i))} \cdot \frac{\exp(\sum \ell_i)}{1 - \exp(- \sum [\text{log1pexp}(\ell_i) - \ell_i])} \Bigg) \\[12pt]
&= \sum \ell_i - \sum \text{log1pexp}(\ell_i) - \text{log1mexp} \bigg( \sum [\text{log1pexp}(\ell_i) - \ell_i] \bigg). \\[6pt]
\end{align}$$
You can write this more succinctly as:
$$\ell = -\sum r_i - \text{log1mexp} \Big(\sum r_i \Big)
\quad \quad \quad \quad \quad 
r_i = \text{log1pexp}(\ell_i) - \ell_i.$$

Programming the joint log-odds function: Here is a simple R function to compute the joint log-odds from a vector of the individual log-odds of the events.
logodds.joint <- function(logodds) {
  r <- VGAM::log1pexp(logodds) - logodds
  - sum(r) - VGAM::log1mexp(sum(r)) }

This gives us a numerically stable version of the log-odds calculation, where everything is computed in log-space.  Below I test this function against direct calculation for a small example where direct calculation is feasible.  As you can see, our function works properly and gives the correct result.
#Create a vector of probabilities pp and joint probability p
set.seed(1)
n  <- 10
pp <- runif(n)
p  <- prod(pp)

#Compute the individual and joint log-odds
ll <- log(pp/(1-pp))
l  <- log(p/(1-p))

#Check results from direct calculation
l
[1] -8.394771

#Check results against our function
logodds.joint(ll)
[1] -8.394771

