On Sat, May 9, 2009 at 12:17 PM, Berwin A Turlach
<ber...@maths.uwa.edu.au> wrote:
> log(H) = log(n) - log( 1/x_1 + 1/x_2 + ... + 1/x_n)
>...But we need to calculate the logarithm of a sum from the logarithms of the
>individual terms.
> ...The way to calculate log(x+y) from lx=log(x) and ly=log(y) ...
> max(lx,ly) + log1p(exp(-abs(lx-ly)))
Agreed completely so far. But instead of calculating the logsum
pairwise, you can do it all in one go, which is both more efficient
and more accurate.
Here are some timing and accuracy measurements of the one-shot logsum
compared to the loop and the Reduce versions. (Full code at the bottom
of this email.) The vector version is much faster and much more
accurate in general. There must be cases where the log1p method
increases accuracy, but I couldn't find them.
-s
Large examples to test accuracy and speed
Test case: runif(1e+06)
function. time error
1 logsum 0.22 9.31e-16
2 logsum_s 0.15 9.31e-16
3 logsum_r 9.75 3.10e-13
Test case: rexp(1e+06)
function. time error
1 logsum 0.21 -1.40e-15
2 logsum_s 0.15 -1.40e-15
3 logsum_r 10.13 -1.38e-14
Test case: abs(rnorm(1e+06))
function. time error
1 logsum 0.24 -4.38e-16
2 logsum_s 0.14 -4.38e-16
3 logsum_r 10.01 -8.74e-14
Test case: rep(1, 1e+05)
function. time error
1 logsum 0.01 1.46e-16
2 logsum_s 0.01 1.46e-16
3 logsum_r 0.96 6.24e-14
Test case: rep(10^-(1:10), each = 10000)
function. time error
1 logsum 0.02 6.14e-16
2 logsum_s 0.01 6.14e-16
3 logsum_r 0.95 -6.96e-12
More accurate even for small cases
Test case: 1:100
function. time error
1 logsum 0 -3.60e-16
2 logsum_s 0 -3.60e-16
3 logsum_r 0 3.24e-15
Test case: abs(rnorm(100))
function. time error
1 logsum 0 -3.48e-16
2 logsum_s 0 -3.48e-16
3 logsum_r 0 -2.09e-15
##########################
# Fast, accurate sum in log space
#
logsum <- function(l) {
maxi <- which.max(l)
maxl <- l[maxi]
maxl + log1p(sum(exp(l[-maxi]-maxl))) }
##########################
##########################
# Simpler, perhaps less accurate sum in log space
#
logsum_s <- function(l) {
maxl <- max(l)
maxl + log(sum(exp(l-maxl))) }
##########################
# Pairwise reduction
logsum_r <- function(x) Reduce( function(lx, ly) max(lx, ly) +
log1p(exp(-abs(lx-ly))), x )
function_names <- c("logsum","logsum_s","logsum_r")
logsum_test <- function(l) {
cat("\nTest case:",deparse(substitute(l)),"\n")
realsum <- sum(l)
logl <- log(l)
results <- times <- list()
lapply( function_names, function(f) times[[f]] <<- system.time(
results[[f]] <<- getFunction(f)(logl))[1])
data.frame(`function`=function_names,
time=as.numeric(times),
error=(exp(as.numeric(results))-realsum)/realsum )
}
set.seed(1)
cat("\n\nLarge examples to test accuracy and speed\n\n")
logsum_test(runif(1000000))
logsum_test(rexp(1000000))
logsum_test(abs(rnorm(1000000)))
logsum_test(rep(1,100000))
logsum_test(rep(10^-(1:10),each=10000))
cat("\n\nMore accurate even for small cases\n\n")
logsum_test(1:100)
logsum_test(abs(rnorm(100)))
______________________________________________
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.
