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

Экономичность алгоритмов метода Монте-Карло

Алгоритмы численного (компьютерного) статистического (вероятностного) моделирования (или алгоритмы метода Монте-Карло) находят весьма широкое применение. В этих численных схемах ключевым элементом является многократное применение алгоритмов (формул) вида

$\xi _{0}=\psi _{\xi }(\alpha _{1},...,\alpha _{k})$

для компьютерного моделирования выборочных значений $\xi _{i}$ одномерных непрерывных случайных величин $\xi \in (a,b);\;-\infty \leq a<b\leq +\infty$, имеющих заданные (выбранные) плотности распределения

$f_{\xi }(u);\;\;u\in (a,b)$

(последняя запись означает, что $f_{\xi }(u)>0$ при $u\in (a,b)$ и $f_{\xi }(u)=0$ при $u\notin (a,b)$; здесь $\alpha _{i}\in U(0,1)$ – стандартные случайные числа, т. е. выборочные значения случайной величины $\alpha \in U(0,1)$, равномерно распределенной в интервале $(0,1)$ это выборочное значение реализуется на компьютере с помощью соответствующего генератора – специальной подпрограммы, именуемой в языках программирования как RAND или RANDOM).

Многократность применения алгоритмов (формул) вида для компьютерной реализации выборочных значений $\xi _{i}$ обусловлена низкой скоростью сходимости (порядка $1/{\sqrt {n}}$ , где $n$ – число обращений к этим формулам) алгоритмов метода Монте-Карло. Поэтому формулы (алгоритмы) вида $\xi _{0}=\psi _{\xi }(\alpha _{1},...,\alpha _{k})$ следует выбирать максимально экономичными.

Метод обратной функции распределения и его обоснование

Стандартным алгоритмом моделирования непрерывной случайной величины $\xi$ является метод обратной функции распределения, основанный на применении формулы

$\xi _{0}=F_{\xi }^{-1}(\alpha _{0});\;\;\alpha _{0}\in U(0,1);$

здесь $F_{\xi }(x)={\bf {P}}\{\xi <x\}=\int _{-\infty }^{x}f_{\xi }(u)\,du$ – функция распределения случайной величины $\xi$, монотонно возрастающая (а значит, имеющая обратную функцию) на интервале распределения $(a,b)$.

Для обоснования формулы метода обратной функции распределения следует показать, что случайные величины $\xi$ и ${\tilde {\xi }}=F_{\xi }^{-1}(\alpha )$ , где $\alpha \in U(0,1)$ , имеют одинаковую функцию распределения. Так как $F_{\xi }^{-1}(\alpha )\in (a,b)$ , то при $x\leq a$ получаем $F_{\tilde {\xi }}(x)={\bf {P}}\{{\tilde {\xi }}<x\}\equiv 0=F_{\xi }(x)$, а при $x\geq b$ имеем $F_{\tilde {\xi }}(x)\equiv 1=F_{\xi }(x)$. Из-за монотонности функции $F_{\xi }(x)$ на интервале $(a,b)$ при $a<x<b$ верны равенства: $F_{\tilde {\xi }}(x)={\bf {P}}\{F_{\xi }^{-1}(\alpha )<x\}={\bf {P}}\{\alpha <F_{\xi }(x)\}=F_{\xi }(x)$; в последнем случае использовано свойство стандартной случайной величины ${\bf {P}}\{\alpha \in (c,d)\subseteq (0,1)\}=d-c$ для $c=0$ и $d=F_{\xi }(x)$ . Таким образом, доказано равенство $F_{\tilde {\xi }}(x)=F_{\xi }(x)$ для всех $x\in (-\infty ,+\infty )$, а значит, $\xi ={\tilde {\xi }}=F^{-1}(\alpha );\;\alpha \in U(0,1)$.

Элементарные плотности распределения

При заданной плотности распределения $f_{\xi }(u)$ для получения вычислимой (т.е. представляющей собой композицию элементарных функций от стандартного случайного числа $\alpha _{0}$) формулы метода обратной функции распределения нужно разрешить уравнение

$\int _{a}^{\xi _{0}}f_{\xi }(u)\,du=\alpha _{0};\;\;\alpha _{0}\in U(0,1)$

относительно верхнего предела интегрирования $\xi _{0}$ в элементарных функциях. Если такое решение удается получить, то плотность $f_{\xi }(u)$ называется элементарной.

ПРИМЕР 1. Рассмотрим случайную величину $\xi ^{(\lambda )}\in (0,+\infty )$, распределенную согласно плотности экспоненциального распределения

$f_{\xi ^{(\lambda )}}(u)=\lambda e^{-\lambda u};\;u>0,\lambda >0.$

Эта плотность является элементарной. Действительно, решая уравнение $\int _{0}^{\xi _{0}^{(\lambda )}}\lambda e^{-\lambda u}\,du=\alpha _{0}^{\prime };\;\alpha _{0}^{\prime }\in U(0,1)$, последовательно получаем

$-e^{-\lambda u}|_{0}^{\xi _{0}^{(\lambda )}}=\alpha _{0}^{\prime };\;1-e^{-\lambda \xi _{0}^{(\lambda )}}=\alpha _{0}^{\prime }\;\;$и, наконец, $\;\;\xi _{0}^{(\lambda )}=-{\frac {\ln(1-\alpha _{0}^{\prime })}{\lambda }}.$

Учитывая, что $1-\alpha _{0}^{\prime }\in U(0,1)$ последнюю формулу (из приведенных выше соображений экономичности вычислений) целесообразно использовать в виде $\xi _{0}^{(\lambda )}=-{\frac {\ln \alpha _{0}}{\lambda }};\;\alpha _{0}\in U(0,1)$.

ПРИМЕР 2. Рассмотрим случайную величину $\xi ^{(0,1)}\in (-\infty ,+\infty )$, распределенную согласно плотности стандартного гауссовского (нормального) распределения

$f_{\xi ^{(0,1)}}(u)={\frac {e^{-u^{2}/2}}{\sqrt {2\pi }}};\;\;-\infty <u<+\infty .$

Эта плотность не является элементарной: в левой части уравнения $\int _{-\infty }^{\xi _{0}^{(0,1)}}{\frac {e^{-u^{2}/2}\,du}{\sqrt {2\pi }}}=\alpha _{0};\;\alpha _{0}\in U(0,1)$ интеграл не берется аналитически (первообразная не может быть выражена в элементарных функциях).

ПРИМЕР 3. Рассмотрим случайную величину $\xi ^{(pol)}\in (0,1)$, распределенную согласно плотности, представляющей собой полином

$f_{\xi ^{(pol)}}(u)=\sum _{m=0}^{M}c_{i}u^{i};\;\;0<u<1.$

В частных случаях плотность этого распределения является элементарной; например, при $c_{p}=(p+1)$ и $c_{i}=0$ для $i\neq p$ получается плотность степенного распределения $f_{\xi ^{(p)}}(u)=(p+1)u^{p};\;0<u<1$, которая является элементарной (несложно вывести соответствующую моделирующую формулу метода обратной функции распределения $\xi _{0}^{(p)}=\alpha _{0}^{1/(p+1)};\;\alpha _{0}\in U(0,1)$). Однако в общем случае при решении уравнения $\int _{0}^{\xi _{0}^{(pol)}}\sum _{m=0}^{M}c_{i}u^{i}\,du=\alpha _{0};\;\alpha _{0}\in U(0,1)$ после применения формулы Ньютона–Лейбница получается соотношение $\sum _{m=0}^{M}{\frac {c_{i}\left[\xi _{0}^{(pol)}\right]^{i+1}}{i+1}}=\alpha _{0}$ , из которого далеко не всегда удается получить аналитическое выражение для $\xi _{0}^{(pol)}$.

Конструирование элементарных плотностей: технология вложенных замен

Существует простой аналитический способ – технология вложенных замен – для получения неограниченного количества элементарных плотностей, которую можно описать следующим образом.

Рассмотрим случайную величину $\eta \in (c,d)$ с элементарной плотностью распределения $f_{\eta }(v);\;c<v<d$ ; это означает, что из уравнения $\int _{c}^{\eta _{0}}f_{\eta }(v)\,dv=\alpha _{0};\;\alpha _{0}\in U(0,1)$ можно получить вычислимое на компьютере выражение $\eta _{0}=\psi _{\eta }(\alpha _{0})$ (в виде композиции элементарных функций).

Возьмем взаимнооднозначное преобразование интервала $(a,b)$ в интервал $(c,d)$, задаваемое монотонной (возрастающей или убывающей) дифференцируемой функцией $\varphi (w)$; причем и сама функция $\varphi (w)$, и ее производная $\varphi ^{\prime }(w)$, и обратная к ней функция $\varphi ^{-1}(s)$ могут быть представлены в виде композиций элементарных функций.

Рассмотрим новую случайную величину $\xi \in (a,b)$, распределенную согласно плотности

$f_{\xi }(u)=f_{\eta }[\varphi (u)]\times |\varphi ^{\prime }(u)|;\;\;u\in (a,b).$

При сделанных предположениях можно утверждать, что плотность $f_{\xi }(u)$ является элементарной; соответствующая формула метода обратной функции распределения имеет вид

$\xi _{0}=\varphi ^{-1}[\psi _{\eta }(\alpha _{0})];\;\;\alpha _{0}\in U(0,1).$

Действительно, для возрастающей функции $\varphi (w)\uparrow$ (для случая убывающей функции $varphi(w)\downarrow$ выкладки аналогичны) с помощью правила замены переменных $v=\varphi (u)$ под знаком интеграла получаем, что уравнение вида $\int _{a}^{\xi _{0}}f_{\xi }(u)\,du=\alpha _{0};\;\alpha _{0}\in U(0,1)$ последовательно преобразуется следующим образом:

$\int _{a}^{\xi _{0}}f_{\eta }[\varphi (u)]\varphi ^{\prime }(u)\,du=\alpha _{0};\;\;\int _{\varphi (a)}^{\varphi (\xi _{0})}f_{\eta }(v)\,dv=\alpha _{0};\;\;\varphi (\xi _{0})=\psi _{\eta }(\alpha _{0}),$

что дает выписанную выше формулу для $\xi _{0}$.

ПРИМЕР 4. Рассмотрим случайную величину $\xi ^{(min)}\in (-\infty ,+\infty )$, распределенную согласно плотности экстремального (точнее, минимального) значения, известной из теории порядковых статистик:

$f_{\xi ^{(min)}}(u)=\exp(-e^{u})\times e^{u};\;\;-\infty <u<+\infty .$

Эта элементарная плотность получается из рассмотренной в примере 1 плотности экспоненциального распределения для $\lambda =1$ с помощью замены $\varphi (u)=e^{u}$; поэтому соответствующая моделирующая формула имеет вид $\xi _{0}^{(min)}=\ln(-\ln \alpha _{0});\;\alpha _{0}\in U(0,1)$.

С помощью технологии вложенных (последовательных) замен можно получать неограниченное количество новых плотностей элементарных распределений, так как в описанном выше процессе в качестве исходной плотности $f_{\eta }(v)$ можно взять преобразованную плотность $f_{\xi }(v)$ и т.д.

Трудоемкость формул метода обратной функции распределения. Компьютерная система NMPUD

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

Например, выведенная в примере 4 формула $\xi _{0}^{(min)}=\ln(-\ln \alpha _{0})$ для распределения минимального значения является в 1.7 раза более трудоемкой, чем исходная формула $\xi _{0}^{(\lambda =1)}=-\ln \alpha _{0};\;\alpha _{0}\in U(0,1)$ из примера 1 для плотности $\lambda =1$, из которой получалась и плотность, и моделирующая формула для распределения минимального значения с помощью преобразования $\varphi (u)=e^{u}$.

Последнее соотношение между трудоемкостями моделирующих формул получено с использованием доступной компьютерной системы NMPUD (Numerical Modelling of Probabilistic Univariate Distributions) – см. https://nmpud.netlify.app, предназначенной для исследования формул метода обратной функции распределения.

Для изучения того или иного примера плотности $f_{\xi }(u);\;u\in (a,b)$ и моделирующей формулы $\xi _{0}=\psi _{\xi }(\alpha _{0})=F_{\xi }^{-1}(\alpha _{0});\;\alpha _{0}\in U(0,1)$, пользователю системы NMPUD демонстрируется основная страница – см. рисунок.

Внизу страницы (в левой части) демонстрируется (если плотность взята из банка системы) или заносится (с помощью естественного языка ввода математических выражений) формула плотности $f_{\xi }(u)$ и соответствующий интервал распределения $(a,b)$ . При этом в центре страницы система строит график этой функции (кривая белого цвета).

Моделирующая формула (алгоритм) $\xi _{0}=\psi _{\xi }(\alpha _{0})$ демонстрируется (или заносится) внизу в правой части основной страницы системы. Система с помощью многократного обращения к формуле строит гистограмму, которая отображается в центре страницы желтыми столбиками; этот процесс можно наблюдать поэтапно. Правильность моделирующей формулы $\xi _{0}=\psi _{\xi }(\alpha _{0})$ определяется визуальной близостью построенной гистограммы и графика плотности $f_{\xi }(u)$.

Также вычисляется и показывается с правой стороны экрана (крупно, зеленым цветом) среднее время одного обращения к формуле $\xi _{0}=\psi _{\xi }(\alpha _{0})$ в наносекундах.

Отметим, что компьютерная система NMPUD является удобным инструментом для того, чтобы определить, какой вклад вносят основные математические операции и функции в трудоемкость формул метода обратной функции распределения. В ходе экспериментов выяснилось, что компьютерные затраты на сложение, умножение, деление и взятие квадратного корня невелики и примерно одинаковы; "умеренно большими" можно назвать затраты на вычисление тригонометрических функций (и обратных к ним) и логарифмов; "экстремально большими" являются затраты на вычисление степенной и показательной функций.

Модификации метода обратной функции распределения. Альтернативные экономичные методы с уравниванием вероятностей

В ряде научных источников упоминается метод обратной функции распределения для дискретных случайных величин с распределением

${\bf {P}}\{\xi =x_{i}\}=p_{i};\;p_{i}>0;\;\sum _{i=1}^{M}p_{i}=1;\;i=1,...,M;\;M\leq \infty ;$

при этом имеется в виду стандартный алгоритм моделирования $\xi$, в котором принимается значение $\xi _{0}=x_{m}$ , где номер $m$ получается вычитанием из реализованного на компьютере стандартного случайного числа $\alpha _{0}\in U(0,1)$ сумм вероятностей $\sum _{j=1}^{m^{\prime }}p_{j}$ до первого получения такого номера $m^{\prime }=m$, для которого $\alpha _{0}-\sum _{j=1}^{m}p_{j}<0$.

Отметим, что:

– название «метод обратной функции распределения» оправдано только для случая $x_{1}<x_{2}<...<x_{M-1}<x_{M}$ (последнее условие существенно ограничивает сферу применения стандартного метода моделирования дискретной случайной величины);

– вместо обратной функции в записи вида $\xi _{0}=F_{\xi }^{-1}(\alpha _{0});\;\alpha _{0}\in U(0,1)$ для дискретной случайной величины нужно использовать запись $\xi _{0}=G(\alpha _{0})=\inf\{y:\alpha _{0}<F_{\xi }(y)\};\;\alpha _{0}\in U(0,1)$;

– запись $\xi _{0}=G(\alpha _{0})$ лишь затушевывает простую идеологию приведенного выше стандартного алгоритма.

Таким образом, на практике следует различать стандартные алгоритмы для моделирования дискретных и непрерывных случайных величин. Несложно построить аналоги (или комбинации) этих алгоритмов и для смешанных "дискретно-непрерывных" случайных величин и/или для непрерывных случайных величин, распределенных на нескольких отстоящих друг от друга интервалах, но эти случаи являются крайне редкими, "экзотическими".

Отметим также нецелесообразность предлагаемого в некоторых статьях и книгах, описывающих различные применения методов Монте-Карло, табулирования функции $F_{\xi }^{-1}(w)$ (это фактически соответствует рассмотрению вместо непрерывной случайной величины $\xi$ "близкой" к ней дискретной случайной величины ${\tilde {\xi }}$ с большим числом значений) и использования соответствующих алгоритмов для полученного дискретного распределения: помимо того, что здесь происходит искажение моделируемого распределения, моделирующие численные схемы получаются трудоемкими.

В заключение упомянем – в качестве альтернативных трудоемким формулам метода обратной функции распределения – разработанные недавно алгоритмы метода исключения с уравниванием вероятностей – двусторонний алгоритм с кусочно-постоянными мажорантой и минорантой и новый, модифицированный зиккурат-метод, которые имеют относительно малые, одинаковые для большого класса распределений (в том числе, для распределений с элементарными плотностями) времена получения выборочных значений; соответствующие описания и коды приведены в библиотеке компьютерной системы PrEMA (Probability Equalization Modelling Algorithms) – см. https://prema.andronix1.ru.

Отметим также, что система PrEMA включает в себя диалоговую систему, позволяющую наглядно (с помощью соответствующих диаграмм) сравнивать затраты на использование формул метода обратной функции распределения и алгоритмов с уравниванием вероятностей. Например, расчеты, проведенные с помощью системы PrEMA, показали, что компьютерные затраты на использование широко применимой формулы метода обратной функции распределения $\xi _{0}^{(p)}=\alpha _{0}^{1/(p+1)};\;\alpha _{0}\in U(0,1)$ для моделирования степенного распределения $f_{\xi ^{(p)}}(u)=(p+1)u^{p};\;0<u<1$ (см. выше пример 3) в 2.5 раза больше затрат на применение двустороннего алгоритма с кусочно-постоянными мажорантой и минорантой и примерно в 3 раза больше затрат на применение зиккурат-метода для того же распределения.

Ме́тод обра́тного преобразова́ния.

См. также

• Преобразование Бокса — Мюллера
• Алгоритм Зиккурат
• Метод Монте-Карло

Литература

1. Metropolis N., Ulam S. The Monte Carlo method // Journal of American Statistical Association. 1949. V. 44. No 249. P. 335–341.

2. Бусленко Н. П., Голенко Д. И., Соболь И. М., Срагович В. Г., Шрейдер Ю. А. Метод статистических испытаний (метод Монте-Карло). М.: Физматгиз, 1962.

3. Hammmersley J. M., Handscomb D. C. Monte Carlo Methods. New York: Jonh Wiley and Sons, 1964.

4. Spanier J., Gelbard E. Monte Carlo Principles and Newtron Transport Problems. Addison–Wesley, Reading, 1969.

5. Соболь И. М. Численные методы Монте-Карло. М.: Наука, 1973.

6. Михайлов Г. А. Некоторые вопросы теории методов Монте-Карло. Новосибирск: Наука, 1974.

7. Ермаков С. М. Метод Монте-Карло и смежные вопросы. М.: Наука, 1974.

8. Марчук Г. И., Михайлов Г. А., Назаралиев М. А., Дарбинян Р. А., Каргин Б. А., Елепов Б. С. Метод Монте-Карло в атмосферной оптике. Новосибирск: Наука, 1976.

9. Ермаков С. М., Михайлов Г. А. Статистическое моделирование. М.: Наука, 1982.

10. Kalos M. H., Whitlock P. A. Monte Carlo Methods. New York: Jonh Wiley & Sons, 1986.

11. Михайлов Г. А. Оптимизация весовых методов Монте-Карло. М.: Наука, 1988.

12. Сабельфельд К. К. Методы Монте-Карло в краевых задачах. М.: Наука, 1989.

13. Михайлов Г. А., Войтишек А. В. Численное статистическое моделирование. Методы Монте-Карло. М.: Академия, 2006.

14. Войтишек А. В. Лекции по численным методам Монте-Карло. Новосибирск: НГУ, 2018.

15. Михайлов Г. А., Войтишек А. В. Статистическое моделирование. Методы Монте-Карло. М.: Юрайт, 2025.