G'day all,
I am working on problems where I have to calculate the logarithm of a
sum or difference from the logarithms of the individual terms; so the
functions logspace_add and logspace_sub which are part of R's API come
in handy.
However, I noticed that logspace_sub can have problems if both
arguments are (very) small or the difference between the arguments are
vary small. The logic behind logspace_sup, when called with arguments
lx and ly, is:
log( exp(lx) - exp(ly) ) = log( exp(lx) * ( 1 - exp(ly - lx) ) )
= lx + log( 1 - exp(ly - lx) )
= lx + log1p( - exp(ly - lx) )
and it is the last expression that is evaluated by logspace_sub.
However the use of log1p for additional precision is appropriate if
exp(ly-lx) is small, i.e. when lx >> ly. If |ly-lx| << 1, then
exp(ly-lx) is close to one; if this term becomes numerically one, then
log1p will return -Inf and "large handfuls of accuracy" are thrown
away. In such circumstances, it would be better to use the following
evaluation scheme:
log( exp(lx) - exp(ly) ) = lx + log( - ( exp(ly - lx) - 1 ) )
= lx + log( - expm1(ly-lx) )
The following code, using the equivalent commands and evaluation schemes
at the R level, illustrates my points:
R> lx <- 2e-17
R> ly <- 1e-17
R> lx + log1p(-exp(ly-lx)) ; lx + log(-expm1(ly-lx))
[1] -Inf
[1] -39.1439465808988
R> lx <- 2e-17
R> ly <- 1e-20
R> lx + log1p(-exp(ly-lx)) ; lx + log(-expm1(ly-lx))
[1] -Inf
[1] -38.4512995253805
R> lx <- 1e-16
R> ly <- 1e-20
R> lx + log1p(-exp(ly-lx)) ; lx + log(-expm1(ly-lx))
[1] -36.7368005696771
[1] -36.8414614929051
In all these cases the output from the second evaluation scheme compares
favourably with the output of yacas or bc.
The current version of R-devel, with this patch applied, passes "make
check FORCE=FORCE" on my machine. The cut-off point for switching
between the two evaluation schemes was chosen arbitrarily.
I hope this patch will proof to be useful and will be applied. :)
Cheers,
Berwin
PS; If use of the equivalent commands and evaluation schemes at the R
level are not convincing enough, then the following code can be used to
verify that logspace_sub has the same behaviour. But one would need two
versions of R for to do so, one without the patch and one with it.
R> install.packages("inline") ## if necessary
R> library(inline)
R> sig <- signature(lx="double", ly="double",
res="double")
R> code <- "*res = logspace_sub(*lx, *ly);"
R> incl <- "#include <Rmath.h>"
R> fn <-
cfunction(sig,code,convention=".C",language="C",includes=incl)
R> lx <- 1e-16
R> ly <- 1e-20
R> fn(lx, ly, res=0)$res
[1] -36.7368005696771
=========================== Full address ============================Berwin A
Turlach Tel.: +65 6516 4416 (secr)
Dept of Statistics and Applied Probability +65 6516 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
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: R-patch
URL:
<https://stat.ethz.ch/pipermail/r-devel/attachments/20090427/bba56f7d/attachment.pl>
