Numeric precision for log(1-exp(x))

I'm doing some math with really big numbers (I'm using Python, but this question isn't Python specific). For one value, I have a formula that gives me f(t) = Pr(X < t). I want to use this formula to get Pr(X >= t) = 1 - f(t). Because f(t) returns values very close to zero, I've been using a log transform and storing log( f(t) ) instead of f(t). My log( f(t) ) is on the order of -1e5 or so.

For multiplying, this works quite well. log( f(t) * g ) = log( f(t) ) + log(g).

But, it's very difficult to compute log( 1 - f(t) ) using only log( f(t) ); I could, of course, temporarily exponentiate the value I store and compute log( 1 - exp( log( f(t) ) ), but that will return log( 1 - 0.0 ) = 0.0 because log( f(t) ) is so close to zero.

You might ask, "Why do you care? If it's close to zero, then 1 minus it is very close to 1." Well, that's a good point you've made. You're a smart cookie.

The trouble is I want to use this to rank values, so I really care if one is log(0.999) and the other is log(0.9999). You might also ask, "Well, why don't you just rank the log( f(t) ) and then reverse the order to get the rankings for log( 1 - f(t) )." Again, I cannot help but point out how great your questions are. It really has been a pleasure talking to you.

But here's the problem: I don't just want to rank by 1 - f(t); I actually want to rank based on Pr(X >= t) * g(t) = (1 - f(t)) g(t). After taking logs, I get log( 1 - f(t) ) + log( g(t) ); ranking based on f(t) alone won't give the correct answer.

In the past I wrote a little Python function to compute log(a + b) from log(a) and log(b):

def log_add(logA,logB):
    if logA == log(0):
        return logB
    if logA<logB:
        return log_add(logB,logA)
    return log( 1 + math.exp(logB-logA) ) + logA

It helps by first normalizing them so that they're close together and then exponentiating when they're close together.

Unfortunately, I wasn't able to get the same trick to work for my subtraction because there is no normalization factor that will bring log(1) and log( f(t) ) close together because they're so far apart.

Does anyone know how to solve this? It seems like such a classic kind of problem; I'd really expect/hope/pray that there is a clever function that operates on the bit level that can give me log(1-x) from log(x). Also, if you know 开发者_StackOverflowhow it works, I'd really, really love to know.

Cheers! Oliver

If log(f(t)) is indeed -1e5 (or of similar order of magnitude) then 0.0 is the best floating point representation of log(1-f(t)). Indeed, f(t) = exp(-1e5) so, by the Taylor series that dmuir mentioned, log(1-f(t)) = -exp(-1e5) (this is actually not an exact equality, but it is an extremely good approximation). Now, -exp(-1e5) = -3.56e-43430, but there are no floating point numbers between 0 and -4e-324, so the best floating point representation is 0.0.

So, what you want to do is impossible with standard floating point numbers.

Does this matter? You say that want to rank based on Pr(X >= t) * g(t) = (1 - f(t)) g(t), which is equivalent to ranking by log( 1 - f(t) ) + log( g(t) ). We found above that log(1-f(t)) = -3.56e-43430, so this term is only going to make a difference if the different values of log(g(t)) differ by no more than this tiny number and if your computation is accurate enough that it can distinguish by these tiny numbers (if you use standard floating point numbers, then your computation will never be accurate enough). In other words, if log(f(t)) is indeed -1e5 or similar, than you can just rank by g(t).

However, it may be that log(f(t)) is of the order of -1e5, but that it sometimes takes values closer to zero like -10 or -1. In that case, you cannot just ignore it and you have to indeed rank by log(1-f(t)) + log(g(t)). You should write this using the math.log1p function: rank by log1p(-f(t)) + log(g(t)). The reason is that if f(t) is close to zero, then log(1-f(t)) is inaccurate but log1p(-f(t)) is accurate. If f(t) is extremely close to zero, such as when log(f(t)) = -1e5, then log1p(-f(t)) = 0.0 because that is the best it can do using standard floating point numbers.

I am using the expression "standard floating point numbers" for a reason. It is possible to use floating point numbers with more precision, and if you really want to capture tiny numbers like -3.56e-43430 that is what you should do. One possibility is in Python is mpmath (unfortunately, it does not seem to support the log1p function). Be warned that this is much slower than standard floating point numbers, and as I said I don't think you'll need it. However, it's worth a go if you want to understand these issues better.