[Rd] Patch proposal for logspace_sub

[Berwin A Turlach]

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>