2
$\begingroup$

To subtract a small probability from another, this answer has constraint on log probabilities l1 > l2: Subtracting very small probabilities - How to compute? but I need a function that works for general l1, l2 (l1=log(p1), l2=log(p2)).

If p1-p2<0, that's fine and let's record that.

For instance, I'll need exp(logdiff(-Inf,log(0.01)))=-0.01 which mirrors 0-0.01=-0.01.

Another case is exp(logdiff(log(0.01),log(0.02)))=-0.01 which mirror 0.01-0.02=-0.01.

The first equality in the question is true:

> exp(log(0.01))-exp(log(0.02))
[1] -0.01
> exp(log(0.01))*(1-exp(-(log(0.01)-log(0.02))))
[1] -0.01

But the second equality isn't holding:

log(0.01)+log(1-exp(-(log(0.01)-log(0.02))))
[1] NaN
Warning message:
In log(1 - exp(-(log(0.01) - log(0.02)))) : NaNs produced

So I'm not sure if the deductions in the above answer are still valid.

Improve this question
$\endgroup$
1
  • $\begingroup$ Based on your answer and comments, perhaps what you really want here is an implementation of the function $\log(\exp(\max(\ell_1,\ell_2)) - \exp(\min(\ell_1,\ell_2)))$ (which does not involve logarithms of negative numbers)? $\endgroup$
    – Ben
    Commented May 19 at 1:41

1 Answer 1

Reset to default
5
$\begingroup$

This requires fixing the logarithmic functions in R to handle negative numbers

The problem here is that the base R functions for logarithms can't handle logarithms of negative numeric values. This is a particularly frustrating deficiency because these functions are actually programmed to give the complex logarithm for a complex number input, but if you put in a negative numeric input the functions do not recognise that this should give a complex output. (The core team should have fixed this long ago but they haven't.) This means that every function that is built on top of the base logarithmic functionality in R has this problem that it can't handle logarithms of negative numbers, and will throw an error when this occurs.

There has been a small step towards solving this in the utilities package, which reprograms the functions log, log2 and log10 so that they work for negative inputs (giving complex numbers as the outputs). In order to get the function you want you would also need to reprogram the other relevant logarithmic functions accordingly. The best way to do this would be to reprogram log1p or log1mexp from scratch in a way that allows large arbitrary numeric or complex inputs. However, here is a "quick and dirty" method that should work for this case most of the time:

logdiff <- function(l1, l2) { 
  
  if (l1 == -Inf) { 
    utilities::log(-exp(l2)) } else {
    l1 + utilities::log(1-exp(-(l1-l2))) } }
Improve this answer
$\endgroup$
2
  • 2
    $\begingroup$ The utilities package is no longer on CRAN. I tend to think (though that may be misguided) that packages removed from CRAN are no longer supported, so shouldn't be recommended. $\endgroup$
    – dipetkov
    Commented May 18 at 13:25
  • 6
    $\begingroup$ This is tricky. What's your guess? Mine is that >99% of the time someone feeds a logarithm function negative numbers they just don't understand what they're doing and/or they haven't noticed negatives as contaminants in their data, whereas <1% of the time they really want a complex result. A complete solution would handle zeros too. Does your code do that? I am not a routine R user. In software I know much better log of zero or log of negative arguments just returns missing, but there is a complex logarithm function you can use it if you want it. $\endgroup$
    – Nick Cox
    Commented May 18 at 14:59

Not the answer you're looking for? Browse other questions tagged or ask your own question.