Читая книгу Нассима Талеба Статистические последствия жирных хвостов, я моделировал в Excel иллюстрации, представленные автором. Особенно мне было интересно генерировать непрерывные случайные величины, заданные различными распределениями. В процессе этой работы я понял, что мне не хватает теоретических знаний в этой области, и нагуглил книгу Геннадия Михайлова и Антона Войтишека Статистическое моделирование. Методы Монте-Карло. Предлагаю вашему вниманию отдельные фрагменты книги. Мои комментарии и примеры в Excel набраны с отступом.
Геннадий Михайлов, Антон Войтишек. Статистическое моделирование. Методы Монте-Карло. – М.: Юрайт, 2024. – 324 с.
Скачать краткое содержание в формате Word или pdf (конспект составляет около 4% от объема книги), скачать архив (содержит Excel-файлы с моделями)
Купить цифровую книгу в ЛитРес, бумажную книгу в Ozon
Введение
Под численным статистическим моделированием обычно понимают реализацию с помощью компьютера вероятностной модели некоторого объекта с целью оценивания изучаемых интегральных характеристик на основе закона больших чисел.
В самом общем виде схема метода Монте-Карло выглядит следующим образом. Пусть требуется вычислить некоторую величину I. Предполагается, что можно построить случайную величину ζ с математическим ожиданием Еζ, равным I, и с конечной дисперсией Dζ, причем выборочные значения ζj случайной величины ζ достаточно просто реализуются на компьютере. Построив большое количество n выборочных значений ζ1, …, ζn, на основе закона больших чисел получаем приближение искомой величины
Примером величин I, допускающими такое представление, является интеграл
Глава 1. МОДЕЛИРОВАНИЕ СЛУЧАЙНЫХ ВЕЛИЧИН
ГЕНЕРАТОРЫ СТАНДАРТНЫХ СЛУЧАЙНЫХ ЧИСЕЛ
Основные свойства стандартного случайного числа
Основным инструментом для моделирования случайных величин на ЭВМ является подходящий генератор стандартных случайных чисел, дающий выборочные значения αi случайной величины α, равномерно распределенной в интервале (0,1). Распределение α является абсолютно непрерывным с плотностью
Функция распределения равна
Математическое ожидание и дисперсия
Утверждение 1.2. Пусть имеется интервал (а, b) ⊆ (0, 1). Тогда случайная величина α равномерно распределена в (а, b) при условии попадания в этот интервал и Р(α ∈ (а, b)) = b – а.
Два типа генераторов стандартных случайных чисел
Существует немало специалистов-вычислителей, которые предпочитают физические датчики случайных чисел, поэтому работы по конструированию таких устройств продолжаются. Большинство же расчетов по методу Монте-Карло проведено и проводится с помощью генераторов псевдослучайных чисел, представляющих некоторые вычислительные программы. Аргументами в пользу применения псевдослучайных чисел являются следующие: возможность воспроизводить расчеты, быстрота получения чисел, отсутствие внешних устройств и необходимости многократной проверки качества получаемых чисел, малая загруженность памяти ЭВМ.
Большинство известных алгоритмов реализации псевдослучайных чисел имеет вид
где начальное число α0 задано. Областью значений функции ψ должен являться интервал (0,1).
Одно из соображений о том, каким образом следует выбирать функцию ψ из (1.8), состоит в следующем. Пары точек
должны, с одной стороны, располагаться на кривой y = ψ(x), а с другой – быть равномерно распределены в квадрате Q2 = {(x,y): 0 < х <1, 0 < y < 1}. Поэтому график функции ψ(x) должен достаточно плотно заполнять квадрат Q2. Примером такой функции ψ(x) может служить
для большого множителя М; здесь выражение в фигурных скобках обозначает дробную часть произведения Мх.
Алгоритм (1.8) с функцией (1.10) называется мультипликативным методом вычетов и является одним из наиболее часто употребляемых алгоритмов при моделировании псевдослучайных чисел.
Свойства мультипликативного метода вычетов
Замечание 1.1. Как правило, генераторы псевдослучайных чисел, представленные в современных версиях языков программирования (FORTRAN, PASCAL, СИ++ и др.), достаточно хорошо протестированы и дают статистически удовлетворительные результаты вычислений по методу Монте-Карло. В дальнейшем мы будем полагать, что используемый в расчетах генератор стандартных случайных чисел дает «настоящие» (теоретические) выборочные значения α.
Я спросил у Chat GPT, как обстоят дела с генератором случайных чисел в Excel. В Excel используется алгоритм Mersenne Twister. Для него существует огромный период ~219937. В Excel вы не столкнетесь с ситуацией, когда случайные числа начнут повторяться.
Тестирование и модификация генераторов случайных и псевдослучайных чисел
Окончательный вывод о качестве того или иного генератора случайных или псевдослучайных чисел следует из результатов тестирования этого генератора. Существует множество тестов. Например, исследуют случайность первой цифры стандартных чисел α, проверяют частоту появления различных цифр в числах, реализуемых генератором (тест «проверка частот»); частоту различных двузначных чисел среди пар цифр, реализуемых подряд генератором (тест «проверка пар») и т.д.
В Excel я взял 500 выборок по 1М значений, нашел среднее по каждой выборке, и построил распределение средних (см. приложенный файл «01. Равномерное распределение.xlsx»):
Рис. 1. Распределение средних значений 500 выборок по 1М значений в каждой выборке; чтобы файлы меньше весили перед их сохранением я установил n = 10k
МОДЕЛИРОВАНИЕ ДИСКРЕТНОГО РАСПРЕДЕЛЕНИЯ (стандартный алгоритм)
Наличие надежного генератора стандартных случайных чисел α позволяет получать выборочные значения случайных величин (одномерных и многомерных) с произвольными законами распределения. Для ряда важнейших (с прикладной точки зрения) распределений стандартные алгоритмы либо не реализуемы, либо не эффективны. В этих случаях предпринимаются (часто довольно успешные) попытки построить специальные алгоритмы моделирования, учитывающие свойства заданного распределения.
Рассмотрим вопрос о моделировании дискретной случайной величины ξ с конечным числом значений х1, …, хN и распределением вероятностей
В силу соотношений (1.18) реализация того или иного значения хm случайной величины ξ означает розыгрыш события с вероятностью рm. Кроме того, из соотношений (1.18) следует, что интервал (0,1) можно разбить на полуинтервалы Δm длины рm:
для m = 1 полагаем Rm-1 = R0 = 0. Используя соответствующий генератор, реализуем выборочное значение α стандартного случайного числа. В силу утверждения 1.2 имеем Р(α ∈ Δm) = рm. Таким образом, если α ∈ Δm, то для данного испытания полагаем ξ = хm.
Для примера возьмем дискретное распределение, заданное таблицей вероятностей:
Рис. 2. Дискретное распределение, заданное таблицей вероятностей
Разыграем стандартное случайное число α. В зависимости от значения α вернем значение ξ (левая часть рис. 3). Далее разыграем n случайных чисел (столбец Е), и подсчитаем частоту появления ξi в выборке из n значений (столбец I):
Рис. 3. Розыгрыш значений ξ на основе α
Пример 1.1. Геометрическое распределение с параметром р:
Соответствующая случайная величина η определяет количество испытаний Бернулли γ до получения первого успеха. Отметим, что во многих изданиях в качестве случайной величины, имеющей геометрическое распределение, берется η’ = η – 1. Эта величина равна числу испытаний Бернулли, предшествующих первому успеху. Для нужд численного моделирования удобнее использовать версию (1.30) геометрического распределения.
Для генерирования случайной величины η в Excel удобнее воспользоваться не стандартным алгоритмом (1.18а), а специальным алгоритмом, предложенным авторами позже:
где u – случайная величина равномерно распределенная в интервале (0,1), ⌈∙⌉– функция округления вверх.
Рис. 4. Розыгрыш случайной величины x, распределенной геометрически с вероятностью успеха в одном испытании р = 0,5
Пример 1.2. Биномиальное распределение с параметрами р и N:
Соответствующая случайная величина
определяет количество успехов в N испытаниях Бернулли.
У биномиального распределения нет простой формулы (типа 1.30а) для генерации случайной величины методом инверсии. Поэтому используется стандартный алгоритм моделирования дискретного распределения (1.18а). Например, для 8 испытаний с вероятностью р = 0,5 получаем распределение:
Рис. 5. Розыгрыш случайной величины x с биноминальным распределением, число испытаний N = 8 вероятность успеха в одном испытании р = 0,5
Пример 1.3. Распределение Пуассона с параметром λ > 0:
Согласно теореме Пуассона, это распределение является предельной формой биномиального распределения при
У распределения Пуассона тоже нет простой формулы для генерации случайной величины методом инверсии. Поэтому используется стандартный алгоритм моделирования дискретного распределения. Например, для λ = 1,5 получаем распределение:
Рис. 6. Розыгрыш случайной величины x с распределением Пуассона для λ = 1,5
СТАНДАРТНЫЙ АЛГОРИТМ МОДЕЛИРОВАНИЯ НЕПРЕРЫВНОЙ СЛУЧАЙНОЙ ВЕЛИЧИНЫ
Метод обратной функции распределения
Рассмотрим алгоритмы моделирования случайной величины ξ, которая принимает значения в интервале (a, b), где ‑∞ < а < b < +∞, и ее функция распределения F(x) = Р(ξ < х) непрерывна и строго возрастает при х ∈ (a, b). При этом
Для случая а = -∞ и b = +∞ соотношения (1.41) приобретают вид
В отличие от дискретного случая, отдельное значение x0 ∈ (a, b) имеет нулевую вероятность. Функция распределения позволяет вычислять вероятности того, что ξ принадлежит некоторому интервалу:
Сформулируем стандартный алгоритм (метод обратной функции распределения).
Алгоритм 1.14. Для численной реализации (моделирования) выборочного значения ξ ∈ (a, b) используем формулу
где α – случайное число, равномерно распределенное на интервале (0, 1).
Алгоритм 1.14, на первый взгляд, закрывает вопрос о моделировании случайных величин ξ ∈ (a, b). Однако остается одна важная «техническая» проблема, связанная с использованием формул вида (1.43) в реальных вычислительных программах.
Задача 1.1. Представить зависимость φ(х) = F–1(x) в виде простой композиции элементарных функций так, чтобы вычисление значения φ(х) могло быть эффективно реализовано на ПК.
В случае, когда задача 1.1 разрешима, будем называть распределение случайной величины ξ и соответствующие формулу (1.43) и алгоритм 1.14 элементарными (с точки зрения возможности численного моделирования). Сразу заметим, что практически для всех распределений, для которых задача 1.1 неразрешима, удается построить альтернативные алгоритмы численной реализации (моделирования) выборочных значений. Однако для случайных величин, имеющих элементарные распределения, алгоритм 1.14 является, как правило, наиболее эффективным (экономичным).
Решение задачи 1.1 будем рассматривать для случая, когда распределение величины ξ ∈ (a, b) является абсолютно непрерывным, что означает существование неотрицательной функции f(u) такой, что для любого интервала (c,d) ⊆ (а, b) выполнено
Это аналог соотношения (1.42). Функция f(u) называется плотностью распределения. Свойствами плотности являются:
Из соотношений (1.42), (1.44) следует, что
и почти всюду имеет место равенство f(u) = dF(u)/du. Из соотношения (1.47) следует, что абсолютно непрерывное распределение можно задавать не функцией распределения F(x), а плотностью f(u).
Итак, пусть имеется непрерывная случайная величина ξ ∈ (a, b), распределенная согласно плотности f(u). С учетом того, что величины ξ и F–1(α) принадлежат интервалу (a, b), а функция F(x) является возрастающей на этом интервале, перепишем (1.43) в эквивалентной форме F(ξ) = α. В свою очередь, в силу соотношений (1.46), (1.47), последнее равенство можно переписать в виде
Распределение случайной величины ξ является элементарным, если решение уравнения (1.48) представимо в виде ξ = ψ(α), где ψ(α) — простая композиция элементарных функций.
Уравнение (1.48) может быть неразрешимым по двум причинам. Во-первых, интеграл в левой части равенства (1.48) не берется (т. е. соответствующая первообразная не выражается в элементарных функциях); примером может служить стандартное нормальное распределение с плотностью
Во-вторых, даже если интеграл из (1.48) берется, получаемое уравнение может быть неразрешимым (в элементарных функциях) относительно ξ.
Продемонстрируем простейшие примеры элементарных распределений и особо отметим важность этих распределений в приложениях метода Монте-Карло и теории вероятностей.
Пример 1.4. Рассмотрим экспоненциальное распределение с плотностью
Сфера применения этого распределения весьма широка. На основе этого распределения формируются пуассоновские потоки, используемые в теории массового обслуживания, в простейших моделях теории переноса излучения, при моделировании случайных полей и т. д.
Решив уравнение вида (1.48)
получим формулу
Заметим, что величина α’ = 1 – α равномерно распределена в (0,1). Действительно, в силу того, что α ∈ (0, 1), имеем Fα’(x) = 0 при х ∈ (–∞, 0) и Fα’(x) = 1 при х ∈ (1, +∞). Наконец, для х ∈ (0, 1) выполнено
т. е. случайная величина α’ имеет функцию распределения
Обращаясь к датчику типа RANDOM, можно считать, что реализуется выборочное значение α’, и тогда моделирующая формула для экспоненциального распределения приобретает вид
Последнее замечание о замене (1 – α) на α’ является весьма важным с прикладной точки зрения, так как во многих практических расчетах количество обращений n к формуле (1.53) очень велико (n ≫ 1), и небольшая экономия ε, связанная с ликвидацией одного вычитания, может дать ощутимый выигрыш в эффективности на величину nε. При практическом применении тех или иных моделирующих соотношений в трудоемких прецизионных расчетах следует весьма тщательно выверять эти формулы на предмет их эффективности.
Рис. 7. Моделирование случайной экспоненциально распределенной величины
Пример 1.5. Рассмотрим степенное распределение с плотностью
Решив уравнение (1.48) для плотности (1.54), получим
Пример 1.6. Похожая на (1.55) формула получается при моделировании случайной величины ξ с распределением Парето
Распределение (1.56) встречается в задачах экономической статистики. Решением уравнения
является
Учитывая, что случайная величина α’ = 1 – α равномерно распределена (см. (1.52 и 1.52a)), целесообразно использовать формулу
Рис. 8. Моделирование случайной величины распределенной по Парето; с = 1,13 соответствует классическому правилу 80/20
МЕТОД СУПЕРПОЗИЦИИ
В этом и следующем разделах изучаются два метода реализации выборочных значений одномерных случайных величин, для которых не удается построить эффективный стандартный алгоритм 1.14. Это методы суперпозиции и исключения, в которых предусматривается моделирование дополнительных вероятностных распределений.
Наибольший практический интерес представляет метод дискретной суперпозиции, когда одномерная целочисленная случайная величина η задается с распределением Р(η = i) = pi, i = 1, 2, … , М; М < ∞. В этом случае
Для реализации этого алгоритма должны существовать эффективные методы получения выборочных значений случайной величины η для любого i. Таким образом, соотношение (1.79) представляет собой взвешенную сумму (смесь) эффективно моделируемых плотностей {fi(u)}. Метод суперпозиции для плотности (1.79) формулируется следующим образом.
Алгоритм 1.21. Реализовав выборочное значение α стандартного случайного числа, согласно вероятностям {pi}, выбираем номер m. Реализуем выборочное значение случайно величины ξ ∈ (a,b) по формуле ξ = ψm(β) , полученной с помощью решения уравнения
относительно переменной ξ.
Пример 1.12. Пусть требуется построить алгоритм моделирования случайной величины имеющей плотность распределения:
Эта функция не является плотностью элементарного распределения, так как соотношение
равносильно уравнению
которое неразрешимо относительно ξ.
Справедливо представление
Моделирующие формулы метода суперпозиции выглядят следующим образом:
Для алгоритма (1.21) величина β равна (6/5)α при m = 1 и β = 6α – 5 при m = 2 и, следовательно,
Для моделирования в Excel использовал формулу с let
Рис. 9. Моделирование случайной величины ξ, отвечающей распределению (1.82a)
МЕТОДЫ ИСКЛЮЧЕНИЯ
В теории статистического моделирования широко используются методы исключения (иногда применяются термины методы отбора и методы отказов, которые также соответствуют, хотя и в меньшей степени, английскому термину rejection technique). Суть этих методов состоит в следующем. Пусть случайный вектор (случайная точка) ξ распределен в некотором множестве X и дано подмножество A ⊆ X.
Алгоритм 1.23. Проводится некоторое статистическое испытание Т и считается, что Т состоялось, если численная реализация вектора ξ ∈ А, и Т не состоялось, если ξ ∉ А.
Мажорантный метод исключения
В подавляющем числе случаев алгоритм 1.23 применяется в следующей ситуации. Пусть требуется моделировать случайный вектор (случайную величину) распределенный в области U ∈ Rl согласно плотности f(u), которая пропорциональна заданной неотрицательной функции g(u), т. е.
Предполагается, что ни один из рассмотренных ранее методов не дает эффективного алгоритма моделирования вектора ξ. Возможность построить моделирующий алгоритм для случайного вектора с плотностью (1.87) дает следующее утверждение.
Утверждение 1.14. Пусть точка (ξ, η) равномерно распределена в области
т. е. в «подграфике» функции g(u); при этом ξ ∈ U и η ∈ (0, g(ξ)). Тогда случайный вектор ξ распределен согласно плотности (1.87).
Таким образом, если получить выборочное значение точки (ξ, η), равномерно распределенной в G, то значение ξ будет распределено согласно плотности (1.87).
Возникает вопрос: каким образом можно реализовать точку, равномерно распределенную в «подграфике» заданной функции? Ответ на этот вопрос дает следующее утверждение, которое является по сути обратным к утверждению 1.14.
Утверждение 1.15. Пусть случайный вектор распределен согласно плотности
а условное распределение при фиксированном значении ξ1 случайной величины η является равномерным на интервале (0, g(ξ)). Тогда случайная точка (ξ1, η) равномерно распределена в «подграфике»
функции g1(u).
Если в утверждении 1.15 выбрать ξ1 = ξ, то в совокупности с утверждением 1.14 получаем логический круг: нам нужно разыграть случайную точку (ξ, η), равномерно распределенную в «подграфике» G (и тогда координата ξ имеет требуемое распределение с плотностью (1.87)), но для этого нужно реализовать вектор ξ по плотности (1.87). Имеется еще утверждение 1.1, из которого следует, что если погрузить «подграфик» G в область G1 в системе координат (u, v) (т. е. G ⊆ G1) и реализовать выборочное значение случайного вектора (ξ1, η), равномерно распределенного в G1, то при условии (ξ1, η) ∈ G пара (ξ1, η) равномерно распределена в G. Тогда, согласно утверждению 1.14, вектор имеет требуемое распределение с плотностью (1.87).
Конструирование области G1 связано с расширением «подграфика» G1 в направлении оси v. Другими словами, рассматривается мажоранта g1(u) функции g(u), такая, что g(u) ≤ g1(u) при u ∈ U. Первое требование к мажоранте g1(u) таково, что для плотности (1.89) имеется эффективный алгоритм (формула) вида ξ1 = ψ1(α̅1) для реализации выборочного значения случайного вектора ξ1 согласно одному из вариантов алгоритма 1.16 (здесь α̅1 – соответствующий набор стандартных случайных чисел). Это дает мажорантный метод исключения.
Алгоритм 1.24. 1. Реализуем выборочное значение ξ1 согласно плотности (1.89): ξ1 = ψ1(α̅1), а также значение η = α2g1(ξ1)
- Если
то в качестве искомого выборочного значения принимаем ξ = ξ1. В случае, когда неравенство (1.91) не выполнено, повторяем п. 1 данного алгоритма и т. д.
Согласно утверждению 1.15, точка (ξ1, η), реализуемая в п. 1 алгоритма 1.24, равномерно распределена в области G1. Если выполнено условие (1.91), то пара (ξ1, η) принадлежит области G и, согласно утверждению 1.1, равномерно распределена в этой области, и тогда, согласно утверждению 1.14, величину можно принять в качестве искомого выборочного значения ξ.
Мажоранту g1(u) функции g(u) следует подбирать так, чтобы объемы G̅1 и G̅ были близки; это выполнено при g1(u) ≈ g(u).
Мажорантный метод исключения используется намного чаще других методов отбора, поэтому в литературе алгоритм 1.24 часто называется просто методом исключения.
Пример 1.13. Рассмотрим случайную величину ξ, имеющую плотность распределения f(u), пропорциональную функции
Интегрируя по частям, получаем:
Т.е.
В Excel я построил графики плотности вероятности f(u) (1.94a) и интегральной функции F(u):
Для нахождения аналитического выражения для F(u) я воспользовался ChatGPT:
Рис. 10. ChatGPT интегрирует F(u)
Рис. 11. Плотность вероятности f(u) и интегральная вероятность F(u)
Функция f(u) не является плотностью элементарного распределения, т. е. уравнение метода обратной функции распределения…
…неразрешимо относительно ξ. В силу того, что |sin u| ≤ 1, в качестве мажоранты функции (1.93) можно взять функцию g1(u) = Зе–u/2. Легко вычислить интеграл…
Следовательно, f1(u) = е–u, u > 0. Это частный случай экспоненциальной плотности (1.51) для λ = 1. Отсюда получаем следующий алгоритм метода исключения:
- Согласно формуле (1.53), получаем следующее выборочное значение: ξ1 = –lnα1. Реализуем также величину η = α2g1(ξ1) = 3α2exp(–ξ1)/2.
- Проверяем неравенство η < g(ξ1) или Зα2 < 2 + sin ξ1. Если оно выполнено, то полагаем, что выборочное значение случайной величины ξ равно ξ1, иначе повторяем п. 1 и т.д.
В Excel формула фильтрации случайного массива работала нестабильно, поэтому пришлось сначала сгенерировать случайные массивы, а затем сохранить их, как значения:
Рис. 12. Разыгрывание случайной величины ξ1, имеющую плотность распределения f(u) (1.94а)
Это лишь малая часть численных методов моделирования, изложенных в книге.