A large chunk of the function below could be replaced with a call to the Reduce function. I don't know if it would be faster, slower, or depend on the situation, but it might make it a little more readable.
-- Gregory (Greg) L. Snow Ph.D. Statistical Data Center Intermountain Healthcare greg.s...@imail.org 801.408.8111 > -----Original Message----- > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r- > project.org] On Behalf Of Berwin A Turlach > Sent: Saturday, May 09, 2009 10:18 AM > To: spencerg > Cc: r-help@r-project.org; joaks1 > Subject: Re: [R] Beyond double-precision? > > G'day all, > > On Sat, 09 May 2009 08:01:40 -0700 > spencerg <spencer.gra...@prodsyse.com> wrote: > > > The harmonic mean is exp(mean(logs)). Therefore, log(harmonic > > mean) = mean(logs). > > > > Does this make sense? > > I think you are talking here about the geometric mean and not the > harmonic mean. :) > > The harmonic mean is a bit more complicated. If x_i are positive > values, then the harmonic mean is > > H= n / (1/x_1 + 1/x_2 + ... + 1/x_n) > > so > > log(H) = log(n) - log( 1/x_1 + 1/x_2 + ... + 1/x_n) > > now log(1/x_i) = -log(x_i) so if log(x_i) is available, the logarithm > of the individual terms are easily calculated. But we need to > calculate the logarithm of a sum from the logarithms of the individual > terms. > > At the C level R's API has the function logspace_add for such tasks, so > it would be easy to do this at the C level. But one could also > implement the equivalent of the C routine using R commands. The way to > calculate log(x+y) from lx=log(x) and ly=log(y) according to > logspace_add is: > > max(lx,ly) + log1p(exp(-abs(lx-ly))) > > So the following function may be helpful: > > logadd <- function(x){ > logspace_add <- function(lx, ly) > max(lx, ly) + log1p(exp(-abs(lx-ly))) > > len_x <- length(x) > if(len_x > 1){ > res <- logspace_add(x[1], x[2]) > if( len_x > 2 ){ > for(i in 3:len_x) > res <- logspace_add(res, x[i]) > } > }else{ > res <- x > } > res > } > > R> set.seed(1) > R> x <- runif(50) > R> lx <- log(x) > R> log(1/mean(1/x)) ## logarithm of harmonic mean > [1] -1.600885 > R> log(length(x)) - logadd(-lx) > [1] -1.600885 > > Cheers, > > Berwin > > =========================== Full address ============================= > Berwin A Turlach Tel.: +65 6515 4416 (secr) > Dept of Statistics and Applied Probability +65 6515 6650 (self) > Faculty of Science FAX : +65 6872 3919 > National University of Singapore > 6 Science Drive 2, Blk S16, Level 7 e-mail: sta...@nus.edu.sg > Singapore 117546 http://www.stat.nus.edu.sg/~statba > > ______________________________________________ > 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.
