[R] Beyond double-precision?

[Greg Snow]

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.snow at imail.org
801.408.8111


> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
> project.org] On Behalf Of Berwin A Turlach
> Sent: Saturday, May 09, 2009 10:18 AM
> To: spencerg
> Cc: r-help at r-project.org; joaks1
> Subject: Re: [R] Beyond double-precision?
> 
> G'day all,
> 
> On Sat, 09 May 2009 08:01:40 -0700
> spencerg <spencer.graves at 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: statba at nus.edu.sg
> Singapore 117546                    http://www.stat.nus.edu.sg/~statba
> 
> ______________________________________________
> 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.