Стандартный метод моделирования непрерывной случайной величины

Стандартный метод моделирования непрерывной случайной величины также называют Методом обратной функции распределения.

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

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

$\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)$ называется элементарной.

ПРИМЕР: Рассмотрим случайную величину $\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)$.

Для получения неограниченного количества элементарных плотностей используется технология (последовательных) вложенных замен.

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

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

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

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

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

${\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 }}$ с большим числом значений) и использования соответствующих алгоритмов для полученного дискретного распределения: помимо того, что здесь происходит искажение моделируемого распределения, моделирующие численные схемы получаются трудоемкими.

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

См. также

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

Литература

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.*

*Некоторых книг нет в открытом доступе, и их прийдётся покупать или скачивать со сторонних ресурсов. Для нахождения литературы воспользуйтесь названием книги.