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: sta...@nus.edu.sg Singapore 117546 http://www.stat.nus.edu.sg/~statba
Index: src/nmath/pgamma.c =================================================================== --- src/nmath/pgamma.c (revision 48409) +++ src/nmath/pgamma.c (working copy) @@ -231,7 +231,11 @@ */ double logspace_sub (double logx, double logy) { - return logx + log1p (-exp (logy - logx)); + register double diff = logy - logx ; + + return diff < 0.1 + ? logx + log (-expm1 (diff)) + : logx + log1p (-exp (diff)); }
______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel