Re: [R] Beyond double-precision?
joaks1 wrote: I need to perform some calculations with some extremely small numbers The R package Brobdingnag does this. It uses logarithmic representation, as recommended by others in this thread, but wraps it all up for you in a custom numeric class so you can use them as any other numeric type. http://cran.r-project.org/web/packages/Brobdingnag/vignettes/brobpaper.pdf __ 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?
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 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.
__
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?
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. timeerror
1logsum 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
1logsum 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
1logsum 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. timeerror
1logsum 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 = 1)
function. time error
1logsum 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
1logsum0 -3.60e-16
2 logsum_s0 -3.60e-16
3 logsum_r0 3.24e-15
Test case: abs(rnorm(100))
function. time error
1logsum0 -3.48e-16
2 logsum_s0 -3.48e-16
3 logsum_r0 -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(100))
logsum_test(rexp(100))
logsum_test(abs(rnorm(100)))
logsum_test(rep(1,10))
logsum_test(rep(10^-(1:10),each=1))
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.
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 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 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.
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?
Dear Jamie: The harmonic mean is exp(mean(logs)). Therefore, log(harmonic mean) = mean(logs). Does this make sense? Best Wishes, Spencer joaks1 wrote: 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. __ 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?
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 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.
Re: [R] Beyond double-precision?
Dear Berwin: Thanks for the elegant correction. Spencer
Berwin A Turlach wrote:
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 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.
Re: [R] Beyond double-precision?
The following packages support high precision precision arithmetic (and the last two also support exact arithmetic): bc - interface to bc calculator http://r-bc.googlecode.com gmp - interface to gmp (gnu multiple precision) http://cran.r-project.org/web/packages/gmp rSymPy - interface to sympy computer algebra system http://rsympy.googlecode.com Ryacas - interface to yacas computer algebra system http://ryacas.googlecode.com On Fri, May 8, 2009 at 4:54 PM, joaks1 joa...@gmail.com 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! -- 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. __ 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] 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?
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.
