lssum: Compute Logarithm of a Sum with Signed Large Summands

Description

Properly compute $\log (x_{1} + \dots + x_{n})$ for given log absolute values lxabs = $l o g (| x_{1} |), . ., l o g (| x_{n} |)$ and corresponding signs signs = $s i g n (x_{1}), . ., s i g n (x_{n})$ . Here, $x_{i}$ is of arbitrary sign.

Notably this works in many cases where the direct sum would have summands that had overflown to +Inf or underflown to -Inf.

This is a (simpler, vector-only) version of copula:::lssum() (CRAN package copula).

Note that the precision is often not the problem for the direct summation, as R's sum() internally uses "long double" precision on most platforms.

Usage

lssum(lxabs, signs, l.off = max(lxabs), strict = TRUE)

Arguments

lxabs

n-vector of values $\log (| x_{1} |), \dots, \log (| x_{n} |)$ .

signs

corresponding signs $s i g n (x_{1}), \dots, s i g n (x_{n})$ .

l.off

the offset to substract and re-add; ideally in the order of max(.).

strict

logical indicating if the function should stop on some negative sums.

Value

$l o g (x_{1} + . . + x_{n}) == l o g (s u m (x)) = l o g (s u m (s i g n (x) * | x |)) == l o g (s u m (s i g n (x) * e x p (l o g (| x |)))) == l o g (e x p (l o g (x 0)) * s u m (s i g n s * e x p (l o g (| x |) - l o g (x 0)))) == l o g (x 0) + l o g (s u m (s i g n s * e x p (l o g (| x |) - l o g (x 0)))) == l . o f f + l o g (s u m (s i g n s * e x p (l x a b s - l . o f f)))$

Examples

Run this code

# NOT RUN {
rSamp <- function(n, lmean, lsd = 1/4, roundN = 16) {
  lax <- sort((1+1e-14*rnorm(n))*round(roundN*rnorm(n, m = lmean, sd = lsd))/roundN)
  sx <- rep_len(c(-1,1), n)
  list(lax=lax, sx=sx, x = sx*exp(lax))
}

set.seed(101)
L1 <- rSamp(1000, lmean = 700) # here, lssum() is not needed (no under-/overflow)
summary(as.data.frame(L1))
ax <- exp(lax <- L1$lax)
hist(lax); rug(lax)
hist( ax); rug( ax)
sx <- L1$sx
table(sx)
(lsSimple <- log(sum(L1$x)))           # 700.0373
(lsS <- lssum(lxabs = lax, signs = sx))# ditto
lsS - lsSimple # even exactly zero (in 64b Fedora 30 Linux which has nice 'long double')
stopifnot(all.equal(700.037327351478, lsS, tol=1e-14), all.equal(lsS, lsSimple))

L2 <- within(L1, { lax <- lax + 10; x <- sx*exp(lax) }) ; summary(L2$x) # some -Inf, +Inf
(lsSimpl2 <- log(sum(L2$ x)))                    # NaN
(lsS2 <- lssum(lxabs = L2$ lax, signs = L2$ sx)) # 710.0373
stopifnot(all.equal(lsS2, lsS + 10, tol = 1e-14))
# }

Run the code above in your browser using DataLab