[R] logsumexp function in R

[Spencer Graves]

On 2/17/2011 3:49 PM, William Dunlap wrote:
>> -----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]
>                    }


       That's great.  The following seems to be equivalent but takes 
less time by replacing an expensive sort with "which.max":


f2.0 <- function(x){
   xmax <- which.max(x)
   log1p(sum(exp(x[-xmax]-x[xmax])))+x[xmax]
}


  set.seed(1)
 > system.time(f2(rnorm(1e7)))
    user  system elapsed
    8.33    0.25    8.69
 > system.time(f2.0(rnorm(1e7)))
    user  system elapsed
    3.41    0.56    4.09
 > system.time(f2(rnorm(1e7)))
    user  system elapsed
    7.99    0.27    8.35
 > system.time(f2.0(rnorm(1e7)))
    user  system elapsed
    3.55    0.33    4.07


       Spencer

> 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.
>>
> ______________________________________________
> 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.
>