[R] Beyond double-precision?
I need to perform some calculations with some extremely small numbers (i.e. likelihood values on the order of 1.0E-16,000). Even when using the double() function, R is rounding these values to zero. Is there any way to get R to deal with such small numbers? For example, I would like to be able to calculate e^-1 (i.e. exp(-1)) without the result being rounded to zero. I know I can do it in Mathematica, but I would prefer to use R if I can. Any help would be appreciated! Many Thanks in Advance! -- View this message in context: http://www.nabble.com/Beyond-double-precision--tp23452471p23452471.html Sent from the R help mailing list archive at Nabble.com. __ 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.
Re: [R] Beyond double-precision?
Yes, all of the numbers are positive. I actually have a Bayesian posterior sample of log likelihoods [i.e. thousands of ln(likelihood) scores]. I want to calculate the harmonic mean of these likelihoods, which means I need to convert them back into likelihoods [i.e. e^ln(likelihood)], calculate the harmonic mean, and then take the log of the mean. I have done this before in Mathematica, but I have a simulation pipeline written almost entirely in R, so it would be nice if I could do these calculations in R. If R cannot handle such small values, then perhaps there's a way to calculate the harmonic mean from the log likelihood scores without converting back to likelihoods? I am a biologist, not a mathematician, so any recommendations are welcome! Thanks! -Jamie spencerg wrote: > > Are all your numbers positive? If yes, have you considered using > logarithms? > > I would guess it is quite rare for people to compute likelihoods. > Instead I think most people use log(likelihoods). Most of the > probability functions in R have an option of returning the logarithms. > > Hope this helps. > Spencer > > joaks1 wrote: >> I need to perform some calculations with some extremely small numbers >> (i.e. >> likelihood values on the order of 1.0E-16,000). Even when using the >> double() function, R is rounding these values to zero. Is there any way >> to >> get R to deal with such small numbers? >> >> For example, I would like to be able to calculate e^-1 (i.e. >> exp(-1)) without the result being rounded to zero. >> >> I know I can do it in Mathematica, but I would prefer to use R if I can. >> Any help would be appreciated! >> >> Many Thanks in Advance! >> > > __ > 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. > > -- View this message in context: http://www.nabble.com/Beyond-double-precision--tp23452471p23457955.html Sent from the R help mailing list archive at Nabble.com. __ 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.
Re: [R] Beyond double-precision?
Thanks Berwin, Spencer, and Gabor!!!
Berwin A Turlach wrote:
>
> G'day all,
>
> On Sat, 09 May 2009 08:01:40 -0700
> spencerg 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 TurlachTel.: +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 117546http://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.
>
>
--
View this message in context:
http://www.nabble.com/Beyond-double-precision--tp23452471p23466589.html
Sent from the R help mailing list archive at Nabble.com.
__
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.
