[R] logsumexp function in R

[William Dunlap]

> -----Original Message-----
> From: r-help-bounces at r-project.org 
> [mailto:r-help-bounces at r-project.org] On Behalf Of Petr Savicky
> Sent: Wednesday, February 16, 2011 10:46 PM
> To: r-help at r-project.org
> Subject: Re: [R] logsumexp function in R
> 
> On Wed, Feb 16, 2011 at 04:14:38PM -0500, Brian Tsai wrote:
> > Hi,
> > 
> > I was wondering if anyone has implemented a numerically 
> stable function to
> > compute log(sum(exp(x))) ?
> 
> Try the following
> 
>    log(sum(exp(x - max(x)))) + max(x)

Sometimes the log1p(x) function, which gives log(1+x)
accurately for small abs(x), helps a little more.
Compare the following 3 functions, which I think give
the same thing for 'ordinary' values of x:

f0 <- function(x) log(sum(exp(x)))
f1 <- function(x) log(sum(exp(x - max(x)))) + max(x)
f2 <- function(x) { x <- sort(x) # mainly so max(x)==x[length(x)]
                    n <- length(x)
                    log1p(sum(exp(x[-n] - x[n]))) + x[n]
                  }
But for the following x f2 gives a more accurate result
than f1, which in turn often gives a more accurate result
than f0:
  > x <- c(0, -50)
  > exp(x)
  [1] 1.000000000000000e+00 1.928749847963918e-22
  > f0(x)
  [1] 0
  > f1(x)
  [1] 0
  > f2(x) # log(1+epsilon) =~ epsilon
  [1] 1.928749847963918e-22
I don't think f2 should ever be less accurate.

expm1(x) is the inverse of log1p(x): it gives exp(x)-1
accurately for small abs(x).

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com  

> 
> If this is not sufficient, consider using CRAN package Rmpfr with
> extended arithmetic.
> 
> Petr Savicky.
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>