On 2/17/2011 3:49 PM, William Dunlap wrote:
-----Original Message-----
From: r-help-boun...@r-project.org
[mailto:r-help-boun...@r-project.org] On Behalf Of Petr Savicky
Sent: Wednesday, February 16, 2011 10:46 PM
To: r-help@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@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@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@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.
