Численно устойчивая оценка log интеграла функции с чрезвычайно малыми значениями

Question

Если у меня есть случайное число Z это определяется как сумма двух других случайных чисел, X а также Y, то распределение вероятностей Z является сверткой распределений вероятностей для X а также Y, Свертка в основном является интегралом от произведения функций распределения. Часто не существует аналитического решения интеграла в свертке, поэтому его необходимо вычислять с помощью базового квадратурного алгоритма. В псевдокоде:

prob_z(z) = integrate(lambda t: prob_x(t) * prob_y(z-t), -inf, inf)

Для конкретного примера, сумма Z нормально распределенной переменной X и лог нормально распределенная переменная Y можно рассчитать с помощью следующего кода Python/Scipy:

from scipy.integrate import quad
from scipy.stats import norm, lognorm
from scipy import log

prob_x = lambda x: norm.pdf(x, 0, 1)  # N(mu=0, sigma=1)
prob_y = lambda y: lognorm.pdf(y, 0.1, scale=10)  # LogN(mu=log(10), sigma=0.1)
def prob_z(z):
    return quad(lambda t: prob_x(t)*prob_y(z-t), -inf, inf)

Теперь я хочу вычислить логарифмическую вероятность. Наивным решением является простое выполнение:

def log_prob_z(z):
    return log(prob_z(z))

Тем не менее, это численно нестабильно. Приблизительно после 39 стандартных отклонений распределение вероятностей численно равно 0,0, поэтому, даже если логарифмическая вероятность имеет какое-то конечное значение, ее нельзя вычислить, просто взяв логарифм вероятности. сравнить norm.pdf(39, 1, 0) от 0,0 до norm.logpdf(39, 1, 0) что около -761. Очевидно, Сципи не вычисляет logpdf как log(pdf)- он находит какой-то другой путь - потому что иначе он вернется -infнеполноценный ответ. Точно так же я хочу найти другой способ решения моей проблемы.

(Вы можете задаться вопросом, почему меня волнует логарифмическая вероятность значений, далеких от среднего. Ответ - подгонка параметров. Алгоритмы подбора могут сближаться, когда логарифмическая вероятность представляет собой какое-то чрезвычайно отрицательное число, но с этим ничего нельзя поделать -inf или же nan.)

Вопрос: кто-нибудь знает, как я могу переставить log(quad(...)) так что я не вычисляю quad(...) и тем самым избежать создания 0.0 в журнале?

Answer 1

Проблема в том, что значения интегрируемой функции слишком малы, чтобы представлять их с двойной точностью, что хорошо только до 1e-308 или около того.

спасение

Когда двойной точности недостаточно для численных вычислений, вызывается mpmath, библиотека для операций с плавающей запятой произвольной точности. Имеет свой quad рутина, но вам нужно реализовать свои функции pdf, чтобы они работали на уровне mpmath (иначе интегрировать некуда). Есть много встроенных функций, включая обычный pdf, поэтому я собираюсь использовать это для иллюстрации.

Здесь я сворачиваю два обычных PDF-файла на расстоянии 70 друг от друга с помощью SciPy:

z = 70
p = quad(lambda t: norm.pdf(t, 0, 1)*norm.pdf(z-t, 0, 1), -np.inf, np.inf)[0]

К сожалению, р точно 0,0.

И здесь я делаю то же самое с mpmath, после import mpmath as mp:

z = 70
p = mp.quad(lambda t: mp.npdf(t, 0, 1)*mp.npdf(z-t, 0, 1), [-mp.inf, mp.inf])

Теперь p - это объект mpmath, который печатается как 2.95304756048889e-543, что намного превышает масштаб двойной точности. И его логарифм, mp.log(p), -1249.22086778731.

SciPy-альтернатива: логарифмическое смещение

Если по какой-то причине вы не можете использовать mpmath, вы можете, по крайней мере, попытаться "нормализовать" функцию, переместив ее значения в диапазон двойной точности. Вот пример:

z = 70
offset = 2*norm.logpdf(z/2, 0, 1)
logp = offset + np.log(quad(lambda t: np.exp(norm.logpdf(t, 0, 1) + norm.logpdf(z-t, 0, 1) - offset), -np.inf, np.inf)[0])

Здесь logp печатает -1264.66566393, что не так хорошо, как результат mpmath (поэтому мы потеряли часть функции), но это разумно. То, что я сделал, было:

• рассчитать логарифм максимального значения логарифма нашей функции (это переменная смещения)
• вычесть это смещение из логарифма PDF; это часть norm.logpdf(t, 0, 1) + norm.logpdf(z-t, 0, 1) - offset
• возвести в степень результат, так как мы не можем просто поместить логарифм в интеграл. Алгебраически это будет то же самое, что произведение pdfs на exp(-offset). Но численно это число с меньшей вероятностью переполнится; действительно, при t = z/2 это exp(0)=1.
• нормально интегрироваться; возьмите логарифм, добавьте смещение к логарифму. Алгебраически, результатом является просто логарифм интеграла, который мы хотели взять.