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.
