Accurately computes the logarithm of the sum of exponentials, that is, log(sum(exp(lx))). If lx = log(x), then this is equivalently to calculating log(sum(x)).
This function, which avoid numerical underflow, is often used when computing the logarithm of the sum of small numbers (|x| << 1) such as probabilities.
This is function is more accurate than
when the values of x = exp(lx) are |x| << 1.
The implementation of this function is based on the observation that
log(a + b) = [ la = log(a), lb = log(b) ] = log( exp(la) + exp(lb) ) = la + log ( 1 + exp(lb - la) )
Assuming la > lb, then |lb - la| < |lb|, and it is
less likely that the computation of 1 + exp(lb - la) will
not underflow/overflow numerically. Because of this, the overall
result from this function should be more accurate.
Analoguously to this, the implementation of this function finds the
maximum value of
lx and subtracts it from the remaining values
This method is optimized for correctness, that avoiding underflowing. It is implemented in native code that is optimized for speed and memory.
 R Core Team, Writing R Extensions, v3.0.0, April 2013.
 Laurent El Ghaoui, Hyper-Textbook: Optimization Models and Applications, University of California at Berkeley, August 2012. (Chapter 'Log-Sum-Exp (LSE) Function and Properties', https://inst.eecs.berkeley.edu/~ee127a/book/login/def_lse_fcn.html)
 R-help thread logsumexp function in R, 2011-02-17. https://stat.ethz.ch/pipermail/r-help/2011-February/269205.html
To compute this function on rows or columns of a matrix,
For adding two double values in native code, R provides
the C function
For properties of the log-sum-exponential function, see .
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
## EXAMPLE #1 lx <- c(1000.01, 1000.02) y0 <- log(sum(exp(lx))) print(y0) ## Inf y1 <- logSumExp(lx) print(y1) ## 1000.708 ## EXAMPLE #2 lx <- c(-1000.01, -1000.02) y0 <- log(sum(exp(lx))) print(y0) ## -Inf y1 <- logSumExp(lx) print(y1) ## -999.3218 ## EXAMPLE #3 ## R-help thread 'Beyond double-precision?' on May 9, 2009. set.seed(1) x <- runif(50) ## The logarithm of the harmonic mean y0 <- log(1/mean(1/x)) print(y0) ## -1.600885 lx <- log(x) y1 <- log(length(x)) - logSumExp(-lx) print(y1) ##  -1.600885 # Sanity check stopifnot(all.equal(y1, y0))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.