Числовая точность для log(1-exp(x))

Я занимаюсь математикой с очень большими числами (я использую Python, но этот вопрос не специфичен для Python). Для одного значения у меня есть формула, которая дает мне f(t) = Pr(X < t). Я хочу использовать эту формулу, чтобы получить Pr(X >= t) = 1 - f(t). Поскольку f(t) возвращает значения, очень близкие к нулю, я использовал логарифмическое преобразование и сохранял log( f(t) ) вместо f(t). Мой log( f(t) ) порядка -1e5 или около того.

Для умножения это работает очень хорошо. log( f(t) * g ) = log( f(t) ) + log(g).

Но очень сложно вычислить log( 1 - f(t) ), используя только log( f(t) ); Я мог бы, конечно, временно возвести в степень сохраняемое значение и вычислить log( 1 - exp( log( f(t) ) ), но это вернет log( 1 - 0.0 ) = 0.0, потому что log( f(t) ) так близко к нулю.

Вы можете спросить: «Какое вам дело? Если оно близко к нулю, то 1 минус очень близко к 1». Что ж, это хорошее замечание, которое вы сделали. Ты умница.

Проблема в том, что я хочу использовать это для ранжирования значений, поэтому мне действительно важно, будет ли одно из них log(0.999), а другое — log(0.9999). Вы также можете спросить: «Почему бы вам просто не ранжировать log( f(t) ), а затем изменить порядок, чтобы получить ранжирование для log( 1 - f(t) )». Опять же, я не могу не отметить, насколько хороши ваши вопросы. Было действительно приятно поговорить с вами.

Но вот проблема: я не просто хочу ранжироваться на 1 - f(t); На самом деле я хочу ранжироваться на основе Pr(X >= t) * g(t) = (1 - f(t)) g(t). После сбора логов получаю log( 1 - f(t) ) + log( g(t) ); ранжирование, основанное только на f(t), не даст правильного ответа.

В прошлом я написал небольшую функцию Python для вычисления log(a + b) из log(a) и 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

Это помогает сначала нормализовать их, чтобы они были близко друг к другу, а затем возводить в степень, когда они близко друг к другу.

К сожалению, я не смог использовать тот же трюк для моего вычитания, потому что нет коэффициента нормализации, который сближал бы log(1) и log( f(t) ), потому что они так далеко друг от друга.

Кто-нибудь знает, как это решить? Это похоже на классическую проблему; Я действительно ожидаю/надеюсь/молюсь, что есть умная функция, работающая на уровне битов, которая может дать мне log(1-x) из log(x). Кроме того, если вы знаете, как это работает, мне бы очень хотелось узнать.

Ваше здоровье! Оливер


statistics probability numerical-methods
person user    schedule 26.07.2011    source источник
comment
Ряд Тейлора для x->log(1-x) равен -(x + x^2/2 + .. x^n/n + ..). С имеющимся у вас диапазоном x, чтобы отличить -x от log(1-x), вам нужно будет использовать ~ 1,4e5 бит, поэтому, возможно, вы могли бы аппроксимировать log (1-x) на -x.   -  person dmuir    schedule 27.07.2011

Ответы (1)


arrow_upward
2
arrow_downward

Если log(f(t)) действительно равно -1e5 (или аналогичного порядка величины), то 0,0 является лучшим представлением log(1-f(t)) с плавающей запятой. В самом деле, f(t) = exp(-1e5) итак, по ряду Тейлора, о котором упоминал дмьюир, log(1-f(t)) = -exp(-1e5) (на самом деле это не точное равенство, но очень хорошее приближение). Теперь -exp(-1e5) = -3.56e-43430, но между 0 и -4e-324 нет чисел с плавающей запятой, поэтому наилучшее представление с плавающей запятой — 0,0.

Итак, то, что вы хотите сделать, невозможно со стандартными числами с плавающей запятой.

Имеет ли это значение? Вы говорите, что хотите ранжироваться на основе Pr(X >= t) * g(t) = (1 - f(t)) g(t), что эквивалентно ранжированию на log( 1 - f(t) ) + log( g(t) ). Мы обнаружили выше, что log(1-f(t)) = -3.56e-43430, поэтому этот термин будет иметь значение только в том случае, если разные значения log(g(t)) отличаются не более чем на это крошечное число, и если ваши вычисления достаточно точны, чтобы различать эти крошечные числа (если вы используете стандартные чисел с плавающей запятой, то ваши вычисления никогда не будут достаточно точными). Другими словами, если log(f(t)) действительно равно -1e5 или подобному, то вы можете просто ранжировать по g(t).

Однако может случиться так, что log(f(t)) имеет порядок -1e5, но иногда принимает значения ближе к нулю, например -10 или -1. В этом случае вы не можете просто игнорировать это, и вам действительно нужно ранжировать на log(1-f(t)) + log(g(t)). Вы должны написать это, используя функцию math.log1p: rank by log1p(-f(t)) + log(g(t)). Причина в том, что если f(t) близко к нулю, то log(1-f(t)) неточное, а log1p(-f(t)) точное. Если f(t) очень близко к нулю, например, когда log(f(t)) = -1e5, то log1p(-f(t)) = 0.0, потому что это лучшее, что можно сделать, используя стандартные числа с плавающей запятой.

Я использую выражение «стандартные числа с плавающей запятой» не просто так. Можно использовать числа с плавающей запятой с большей точностью, и если вы действительно хотите захватить крошечные числа, такие как -3.56e-43430, это то, что вам нужно сделать. Одна из возможностей в Python — это mpmath (к сожалению, похоже, что он не поддерживает функцию log1p) . Имейте в виду, что это намного медленнее, чем стандартные числа с плавающей запятой, и, как я уже сказал, я не думаю, что вам это понадобится. Тем не менее, стоит попробовать, если вы хотите лучше разобраться в этих вопросах.

person Jitse Niesen    schedule 28.07.2011
comment
+1 Спасибо, отличные идеи. Я мог бы использовать высокую точность только для этого шага. Или, возможно, просто используйте расширение Тейлора для log (1-exp (x)). - person user; 31.07.2011