Функции от случайных величин

Температура в градусах Цельсия пересчитывается в градусы Фаренгейта, измерение стандартизуется, логарифм цены потенцируется обратно в цену, равномерное случайное число превращается в показательное (экспоненциальное) время ожидания. Каждый раз случайная величина $X$ с известным распределением проходит через функцию $g$, и вопрос в том, каково распределение $Y = g(X)$. Для счётных величин ответ — сумма; для плотностей — формула замены переменных (change-of-variables formula), которая говорит, что вероятность сохраняется, а плотность при этом сжимается или растягивается. Та же формула, прочитанная в обратную сторону, — это способ, которым компьютеры получают случайные значения с любым распределением, и способ, которым нормализующие потоки (normalizing flows) вычисляют правдоподобие.

Обозначения. $Y = g(X)$. Для обратимой $g$ через $x = g^{-1}(y)$ обозначен прообраз $y$, а через $g'$ — производная. $f_X, F_X$ и $f_Y, F_Y$ — плотности и функции распределения $X$ и $Y$.

Дискретный случай#

Если $X$ принимает не более чем счётное число значений, то и $Y$ тоже, а $Y = y$ происходит ровно тогда, когда $X$ попадает в одно из значений, которые $g$ переводит в $y$:

$$P(Y = y) = \sum_{x \,:\, g(x) = y} P(X = x) .$$

Для честной кости $X$ и $Y = (X - 3)^2$ значение $Y = 1$ получается из $X = 2$ и из $X = 4$, поэтому $P(Y = 1) = 2/6$. Больше ничего не нужно: вероятности значений переносятся вместе с ними и складываются там, где значения сливаются.

Метод функции распределения#

Для непрерывной $X$ надёжный общий метод проходит через функцию распределения, потому что $\{Y \le y\}$ — событие, связанное с $X$:

$$F_Y(y) = P(g(X) \le y) ,$$

а плотность $Y$ — производная полученной функции.

Квадрат. Пусть $X \sim U(-1, 1)$ и $Y = X^2$. При $0 \le y \le 1$ условие $Y \le y$ означает $-\sqrt y \le X \le \sqrt y$ — интервал длины $2\sqrt y$ из общей длины 2:

$$F_Y(y) = \sqrt y, \qquad f_Y(y) = \frac{1}{2 \sqrt y} .$$

В половине случаев $X^2$ меньше $1/4$, хотя $1/4$ — лишь четверть области значений; плотность $Y$ в нуле бесконечна: там парабола плоская, и многие значения $X$ теснятся в немногих значениях $Y$.

Линейное отображение. Для $Y = aX + b$ с $a > 0$ имеем $F_Y(y) = P(X \le (y - b)/a) = F_X((y - b)/a)$, поэтому

$$f_Y(y) = \frac{1}{\lvert a \rvert} f_X\Big(\frac{y - b}{a}\Big) ,$$

и то же верно при $a < 0$. Применённая к $X \sim N(\mu, \sigma^2)$ с $a = 1/\sigma$, $b = -\mu/\sigma$, эта формула даёт плотность $N(0, 1)$: это стандартизация $Z = (X - \mu)/\sigma$, обещанная в предыдущей заметке. Верно и обратное: величина $aX + b$ при нормальной $X$ нормальна, так что температура с распределением $N(20, 3^2)$ °C — это температура с распределением $N(68, 5.4^2)$ °F.

Формула замены переменных#

Для строго возрастающей дифференцируемой $g$ имеем $F_Y(y) = F_X(g^{-1}(y))$, и дифференцирование сразу даёт плотность $Y$; для убывающей $g$ производная меняет знак, что и учитывает модуль:

$$f_Y(y) = \frac{f_X(x)}{\lvert g'(x) \rvert}, \qquad x = g^{-1}(y) .$$

Формула говорит, что вероятность сохраняется полоса за полосой. Короткий интервал длины $dx$ около $x$ несёт вероятность $f_X(x)\, dx$, а $g$ переводит его в интервал длины $dy = \lvert g'(x) \rvert\, dx$ около $y$, который должен нести ту же вероятность: $f_Y(y)\, dy = f_X(x)\, dx$. Где $g$ сжимает ($\lvert g' \rvert < 1$), плотность растёт; где растягивает — падает:

Рисунок

Замена переменных: вероятность сохраняется, плотность меняется

Пример — логнормальное распределение (log-normal): у $Y = e^X$ с $X \sim N(0, 1)$ производная $g'(x) = e^x = y$ и

$$f_Y(y) = \frac{\varphi(\ln y)}{y}, \qquad y > 0 .$$

Монотонные отображения переносят квантили, поэтому медиана $Y$ равна $e^0 = 1$, а $P(Y > e) = P(X > 1) = 0.159$. Величины, которые получаются как произведения многих положительных множителей, — доходы, размеры файлов, времена ответа — часто близки к логнормальным: их логарифм — сумма, близкая к нормальной по причинам, которые объясняет центральная предельная теорема.

Если $g$ не монотонна, вклад даёт каждый прообраз: $f_Y(y) = \sum_{x : g(x) = y} f_X(x) / \lvert g'(x) \rvert$. Для $Y = X^2$ из примера выше два прообраза $\pm\sqrt y$ дают по $\frac{1/2}{2\sqrt y}$, и в сумме получается та же $1/(2\sqrt y)$.

Метод обратной функции#

Прочитанный в обратную сторону, метод функции распределения превращает равномерные случайные числа в числа с любым распределением. Если $U \sim U(0, 1)$, а $F$ — непрерывная строго возрастающая функция распределения, то $X = F^{-1}(U)$ имеет функцию распределения $F$, потому что

$$P(F^{-1}(U) \le x) = P(U \le F(x)) = F(x) .$$

Для любой функции распределения, с горизонтальными участками или скачками, то же верно для квантильной функции $Q(u) = \inf\{x : F(x) \ge u\}$ из предыдущей заметки вместо $F^{-1}$: $Q(u) \le x$ ровно тогда, когда $u \le F(x)$, потому что $F$ не убывает и непрерывна справа. Дискретный сэмплер ниже — это и есть $Q$.

Для показательного распределения $F(x) = 1 - e^{-x}$ и $F^{-1}(u) = -\ln(1 - u)$:

Рисунок

Метод обратной функции: от равномерного u к показательному x

Та же идея годится и для сэмплирования из дискретного распределения. Для категориального распределения с вероятностями $p_1, \dots, p_K$ разрежем $[0, 1]$ на идущие подряд отрезки длин $p_1, \dots, p_K$ и вернём номер отрезка, в который попало $U$; это один поиск по накопленным суммам. Языковая модель сэмплирует следующий токен именно так — после того как температура или top-$p$ преобразовали вероятности.

Обратное утверждение — интегральное преобразование вероятности (probability integral transform): если у $X$ непрерывная функция распределения $F$, то $F(X) \sim U(0, 1)$. Так проверяют вероятностные прогнозы, поскольку значения $F(x)$, которые хорошо откалиброванный прогноз даёт в наблюдённых исходах, должны выглядеть равномерно распределёнными; по этой же причине p-значение теста, статистика которого при нулевой гипотезе имеет непрерывное распределение, распределено равномерно, когда эта гипотеза верна. Для дискретной статистики p-значение лишь супер-равномерно, $P(p \le \alpha) \le \alpha$: в двустороннем тесте для 100 бросков монеты из заметки о проверке гипотез честная монета даёт $p \le 0.05$ с вероятностью $0.035$. Тест остаётся корректным, но консервативным.

Нормализующие потоки в одномерном случае#

Нормализующий поток моделирует сложную плотность как простую, пропущенную через цепочку обратимых отображений, $Y = g_K(\cdots g_1(Z))$ с $Z \sim N(0, 1)$. Применив формулу замены переменных по одному разу к каждому отображению, получаем точный логарифм плотности в точке данных:

$$\log f_Y(y) = \log \varphi(z) - \sum_{k=1}^{K} \log \lvert g_k'(z_{k-1}) \rvert, \qquad z_0 = z = g^{-1}(y), \; z_k = g_k(z_{k-1}) .$$

Обучение максимизирует этот логарифм правдоподобия по параметрам отображений. В размерности больше единицы производная превращается в определитель матрицы Якоби (якобиан), и как сделать этот определитель дешёвым в вычислении — задача конструирования, которой посвящена заметка о нормализующих потоках в курсе «Генеративные модели».

Что изменилось#

	До	С формулой замены переменных
Распределение $g(X)$	симулировать и смотреть на гистограмму	точно: $f_X(x) / \lvert g'(x) \rvert$ или $F_Y$ через $F_X$
Стандартизация нормального распределения	принималась на веру	доказана: $(X - \mu)/\sigma \sim N(0, 1)$
Сэмплирование из распределения	для каждого распределения нужен свой приём	$F^{-1}(U)$ для любой функции распределения; накопленные суммы для категориального
Плотность потока	—	точный логарифм правдоподобия, по одному $\log \lvert g' \rvert$ на отображение

Ограничения#

• $F^{-1}$ не всегда выражается в замкнутом виде. Обратная к нормальной функции распределения не выражается через элементарные функции; сэмплеры используют другие преобразования (Бокса–Мюллера, которому нужны две величины) или численное обращение.
• Формуле нужна дифференцируемая, локально обратимая $g$. Функция с горизонтальными участками, например округление или обрезка (clipping), сосредоточивает положительную вероятность в отдельных значениях, и тогда $Y$ не является непрерывной: округление делает её дискретной, обрезка — смесью масс в точках и плотности.
• Для функций от нескольких величин нужно их совместное распределение. Распределение суммы или отношения двух случайных величин зависит от того, как они зависят друг от друга; это тема заметки о совместных распределениях.

Что нужно запомнить#

• Дискретный случай: $P(g(X) = y)$ — сумма $P(X = x)$ по прообразам $y$.
• Непрерывный случай в общем виде: вычислить $F_Y(y) = P(g(X) \le y)$, затем продифференцировать.
• Непрерывный случай, монотонная $g$: $f_Y(y) = f_X(x) / \lvert g'(x) \rvert$ при $x = g^{-1}(y)$; вероятность сохраняется, плотность сжимается или растягивается.
• Линейные отображения сохраняют нормальность; $(X - \mu)/\sigma \sim N(0, 1)$; $e^X$ от нормальной величины имеет логнормальное распределение.
• $F^{-1}(U)$, а в общем случае $Q(U)$, имеет функцию распределения $F$ (метод обратной функции, inverse transform sampling); $F(X)$ распределена равномерно, если $F$ непрерывна (интегральное преобразование вероятности).

Код#

Две полосы с рисунка, квадрат равномерной величины и метод обратной функции, проверенные численно:

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
Z = stats.norm()

# Y = e^X for X ~ N(0, 1): each band keeps its probability over an interval of a different length
for a, b in [(-1.2, -0.7), (0.8, 1.3)]:
    p = Z.cdf(b) - Z.cdf(a)
    print(round(p, 3), round(np.exp(b) - np.exp(a), 3), round(p / (np.exp(b) - np.exp(a)), 3))
# 0.127 0.195 0.649   (squeezed: high density)
# 0.115 1.444 0.08    (stretched: low density)
print(round(Z.sf(1), 4), stats.lognorm(1).median())              # 0.1587 1.0

# Y = X^2 for X ~ U(-1, 1): P(Y <= 1/4) = sqrt(1/4)
y = rng.uniform(-1, 1, 1_000_000) ** 2
print(round((y <= 0.25).mean(), 3))                              # 0.5

# inverse transform: -ln(1 - U) is Exp(1)
u = np.array([0.1, 0.3, 0.5, 0.7, 0.9])
print((-np.log(1 - u)).round(3))                                 # [0.105 0.357 0.693 1.204 2.303]
x = -np.log(1 - rng.random(1_000_000))
print(round(x.mean(), 3), round((x > 2).mean(), 4))              # 1.0 0.1354; exactly 1 and e^-2 = 0.1353

# a categorical distribution sampled through cumulative sums, as for a next token
p = np.array([0.5, 0.3, 0.2])
k = np.searchsorted(np.cumsum(p), rng.random(1_000_000), side='right')
print((np.bincount(k) / 1_000_000).round(3))                     # [0.499 0.301 0.199]

# the probability integral transform: F(X) is uniform
print(np.quantile(Z.cdf(rng.normal(size=1_000_000)), [0.1, 0.5, 0.9]).round(3))   # [0.1   0.499 0.9  ]

Дальше#

Распределение — это целая функция; часто для работы достаточно двух-трёх чисел: где его центр, насколько оно широкое, насколько тяжёлые у него хвосты. Следующая заметка определяет математическое ожидание и дисперсию, показывает, что математическое ожидание $g(X)$ можно вычислить, вовсе не находя распределения $g(X)$, и выводит правило $1/\sqrt n$ для средних из вводной заметки.

Литература#

• J. K. Blitzstein, J. Hwang. Introduction to Probability, 2nd ed. CRC Press, 2019. Глава 5 (универсальность равномерного распределения), глава 8 (преобразования, замена переменных).
• L. Devroye. Non-Uniform Random Variate Generation. Springer, 1986. Глава II (метод обратной функции).
• D. Rezende, S. Mohamed. Variational Inference with Normalizing Flows. ICML, 2015.