1626435387-5d0ec7dd22f55dfcc65f3209f497fb0b (844202), страница 31
Текст из файла (страница 31)
формулу (2.23)).(1)Согласно утверждению 11.5, точка ξ 0 , η0 равномерно распреде-186(1)лена в области G(1) , т. е. ξ 0 , η0 ∈ U G(1) .2. Если(1) η0 < g ξ 0 ,(11.35)(1)то пара ξ 0 , η0 принадлежит области G и, согласно утвержде(1)нию 2.3, равномерно распределена в этой области: ξ 0 , η0 ∈ U (G),(1)и тогда, согласно утверждению 11.4, величину ξ 0 можно принятьв качестве искомого выборочного значения, имеющего нужную плот(1)ность распределения (11.31): ξ 0 = ξ 0 .В случае, когда неравенство (11.35) не выполнено, повторяемпункт 1 данного алгоритма и т. д.Рис. 11.3.
Схема мажорантного метода исключенияСогласно формуле (11.30) и утверждению 2.3, трудоемкость (затраты) s̃ алгоритма 11.10 пропорциональна величинеs=Ḡ(1)1. (1) =ḠP ξ ,η ∈ G(11.36)Таким образом, мажоранту g (1) (u) функции g(u) следует подбирать187так, чтобы объемы (площади) Ḡ(1) и Ḡ были близки, т. е.
близкимидолжны быть сами функции: g (1) (u) & g(u).ПРИМЕР 11.6 (Г; 1 балл; см., например, [9, 13]). Пусть требуется построить алгоритм моделирования выборочного значения ξ0 случайнойвеличины ξ, имеющей плотность распределения fξ (u), пропорциональную функцииsin ue−u , u > 0.(11.37)g(u) = 1 +2R +∞Заметим, что плотность fξ (u) = g(u)/Ḡ; Ḡ = 0 g(u) du не является элементарной. Действительно, если записать соответствующееRξуравнение типа (2.16) 0 0 fξ (u) du = α0 и проинтегрировать по частям,то получатся уравненияξ0(cos u + sin u) e−u 1−u−e − = α0 или4Ḡ04e−ξ0+ (sin ξ0 + cos ξ0 )e−ξ0+ 4Ḡα0 − 5 = 0,которые неразрешимы относительно ξ0 .Из последних соотношений, заменяя ξ0 на +∞, а α0 – на единицу,получаем Ḡ = 5/4.В силу того, что | sin u| ≤ 1, в качестве мажоранты функции (11.37)можно взять функцию g (1) (u) = 3 e−u /2, u > 0.R +∞Легко вычислить интеграл Ḡ(1) = 0 g (1) (u) du = 3/2.Следовательно, пропорциональная мажоранте g (1) (u) плотность случайной величины ξ (1) (см.
соотношения (11.33)) имеет видfξ(1) (u) = e−u , u > 0. Это частный случай экспоненциальной плотности(2.18) с параметром λ = 1 и моделирующей формулой(1)ξ0 = − ln α1(11.38)(см. соотношение (2.20)). Отсюда получаем следующий алгоритм метода исключения.АЛГОРИТМ 11.11. 1. Моделируем выборочное значение (11.38), атакже величину(1) (1) (1) η0 = α2 g (1) ξ0 = 3α2 exp − ξ0 /2 ∈ U 0, g (1) ξ0(1)(см. формулу (2.23)). Точка ξ0 , η0 равномерно распределена в «под(1)графике» функции g (1) (u): ξ0 , η0 ∈ U G(1) .188(1) 2. Проверяем неравенство η0 < g ξ0или(1)3α2 < 2 + sin ξ0 .(11.39)(1)Если неравенство (11.39) выполнено, то точка ξ0 , η0 принадлежит «подграфику» функции g(u) из (11.37) и является равномерно(1)распределенной в этом множестве: ξ0 , η0 ∈ U (G).
Тогда в качестве(1)выборочного значения ξ0 случайной величины ξ берем ξ0 = ξ0 .Если неравенство (11.39) не выполнено, то повторяем пункт 1 ит.д.Трудоемкость s этого алгоритма (т. е. среднее число попыток розыг(1)рыша пар ξ0 , η0 до выполнения неравенства (11.39)) равна величинеs = Ḡ(1) /Ḡ = 3/2 : 5/4 = 1, 2 (см. соотношение (11.30)). Эта величинаблизка к единице, и поэтому алгоритм 11.11 можно считать достаточноэффективным (экономичным).Описание примера 11.6 закончено.11.6. Замечания о построении мажорант. Двустороннийметод исключения. Моделирование усеченных распределений.Наиболее простой вариант мажорантного метода исключения получается в случае, когда функция g(u) определена на ограниченном множествеX в Rd и существует (известна) константа H, такая, что g(u) ≤ H приu ∈ X.
В качестве мажоранты здесь можно выбрать g (1) (u) ≡ H. При(1)этом на первом этапе алгоритма 11.10 выборочное значение ξ 0 вектораξ (1) моделируется согласно равномерному распределению в множествеX.В одномерном случае для X = (a, b) ⊂ R; −∞ < a < b < +∞алгоритм 11.10 выглядит здесь следующим образом.АЛГОРИТМ 11.12 (см., например, [9]).
1. Моделируем(1)ξ0 = a + α1 (b − a) и η0 = α2 H (оба раза использована формула вида(2.23)).(1) 2. Если η0 ≤ g ξ0 , то в качестве выборочного значения ξ0 слу(1)чайной величины ξ ∈ (a, b) принимаем ξ0 = ξ0 , иначе повторяемпункт 1 и т.д.В свое время Дж. Нейман первым предложил мажорантный методисключения именно в такой форме, и часто этот простой частный случай называют методом Неймана (мы также будем придерживатьсяэтой терминологии).
В литературе по методам Монте-Карло иногда методами Неймана называют все методы исключения.189Важное обобщение алгоритма 11.12 связано с использованием кусочно-постоянных мажорант (см. подразделы 11.3, 11.4), когда интервал(a, b) разбит на полуинтервалы ∆i = (ui−1 , ui ], i = 1, . . . , M точкамиa = u0 < u1 < . . . < uM −1 < uM = b, и для каждого ∆i известнаконстанта Hi такая, что g(u) ≤ Hi при u ∈ ∆i .Если g (1) (u) ≡ Hi при u ∈ ∆i , то на первом этапе алгоритма 11.10(0)требуется моделировать выборочное значение ξ1 случайной величиныξ (1) согласно кусочно-постоянной плотности (здесь следует использовать алгоритм 11.6).Следующая модификация алгоритма 11.10 эффективна в достаточно распространенном случае, когда требуется моделировать выборочноезначение случайного вектора ξ, плотность которого пропорциональнафункции g(u), вычисление значений которой весьма трудоемко.В этом случае помимо мажоранты g (1) (u) строим миноранту g (2) (u),такую, что0 ≤ g (2) (u) ≤ g(u) ≤ g (1) (u); u ∈ X.(11.40)АЛГОРИТМ 11.13 (см., например, [9, 24]).
1. Моделируем выбороч(1)ное значение ξ 0 = ψ (1) ᾱ1 согласно плотности (11.33), а также зна(1) чение η0 = α2 g (1) ξ 0 .2. Вместо неравенства (11.35) проверяем сначала соотношение(1) η0 < g (2) ξ 0 .Если оно выполнено, то пара (ξ 1 , η) принадлежит «подграфику»функции g (2) (u), а значит, и области G. Тогда можно положить, чтовыборочное значение ξ 0 случайного вектора ξ, распределенного согласно(1)плотности (11.31), равно ξ 0 = ξ 0 .(1)В случае же η0 ≥ g (2) ξ 0 проверяем неравенство (11.35).
Если оно(1)выполнено, то ξ 0 = ξ 0 , иначе повторяется пункт 1 данного алгоритма и т. д.В связи с соотношением (11.40) алгоритм 11.13 называют двусторонним методом исключения.В случае, когда все три функции из неравенства (11.40) близки:g (2) (u) . g(u) . g (1) (u), а миноранта g (2) (u) и мажоранта g (1) (u) легко вычислимы, проверка (11.40), связанная с трудоемким вычислением(1) значения g ξ 0 , будет происходить относительно редко, и двусторонний метод может дать существенный выигрыш по сравнению с «односторонним» алгоритмом 11.10.190В качестве функций g (2) (u) и g (1) (u) часто используются кусочнопостоянные приближения снизу и сверху для функции g(u).Теперь напомним алгоритм 5.1 из подраздела 5.3, с помощью которого можно моделировать случайные векторы (случайные величины) сусеченными распределениями.Рассмотрим случайный вектор ξ (1) , распределенный в областиY ∈ Rd согласно плотности fξ (1) (u).ОПРЕДЕЛЕНИЕ 11.1 (см., например, [9]).
Случайный вектор ξ имеет усеченное распределение вектора ξ (1) , если он распределен вподобласти X ⊂ Y и его плотность распределения fξ (u) пропорциональна в X плотности fξ (1) (u):fξ (u) = H fξ (1) (u) = Rfξ (1) (u), u ∈ X ⊂ Y.f(w) dwX ξ (1)(11.41)В случае, когда имеется эффективный алгоритм моделирования вы(1)борочного значения ξ 0 случайного вектора ξ (1) , можно использоватьследующий алгоритм исключения для моделирования выборочного значения вектора ξ, имеющего усеченное распределение (11.41).АЛГОРИТМ 11.14 (см., например, [9, 24], а также алгоритм 5.1 изподраздела 5.3 данного пособия). 1.
Моделируем выборочное значение(1)ξ 0 случайного вектора ξ (1) в области Y согласно плотности fξ (1) (u).(1)(1)2. Если ξ 0 ∈ X, то принимаем ξ 0 = ξ 0 , иначе повторяетсяпункт 1 данного алгоритма и т. д.Несложно понять, что алгоритм 11.14 является частным случаем алгоритма 11.10 в области Y , в котором для функции (11.41) рассмотренамажоранта H fξ (1) (u) для всех u ∈ Y . В этом случае при u ∈ X имеемfξ (u) = H fξ (1) (u), а при u ∈ Y \ X выполнено fξ (u) = 0 < H fξ (1) (u).Моделирование значения случайной величины η0 не требуется (см.(1)пункт 1 алгоритма 11.10), так как при ξ 0 ∈ X неравенство (11.35)(1)заведомо выполнено, а при ξ 0 ∈ Y \ X – заведомо не выполнено.Трудоемкость алгоритма 11.14 пропорциональна величине s = H(см.
соотношения (11.30), (11.41)).Заметим, что во многих случаях для моделирования выборочногозначения ξ 0 случайного вектора ξ с усеченным распределением вида(11.41) удается построить более эффективную, чем алгоритм 11.14, чис-191ленную процедуру, не связанную с включением области X в множествоY (см., в частности, пример 5.2).Более того, для одномерного случая удается доказать следующееутверждение.УТВЕРЖДЕНИЕ 11.6. Если плотность fξ(1) (u) случайной величины ξ (1) ∈ (A, B) (здесь −∞ ≤ A < B ≤ +∞) является элементарной,то для любых a и b, таких, что A ≤ a < b ≤ B, усеченная плотностьраспределенияfξ (u) = H fξ(1) (u) = R bafξ(1) (u),u ∈ (a, b)(11.42)fξ(1) (w) dwслучайной величины ξ ∈ (a, b) также является элементарной.ДОКАЗАТЕЛЬСТВО.
Элементарность плотности fξ(1) (u) означает,что при решении основного уравнения (2.16) метода обратной функцииR ξ(1)распределения вида A0 fξ(1) (u) du = α0 удается перейти к уравнениюξ(1)с первообразной F̃ (u)A0 = α0 и разрешить его относительно верхнего(1)предела: ξ0 = ψ (1) (α0 ); здесь функции F̃ (u) и ψ (1) (u) представляютсобой композиции элементарных функций.Теперь рассмотрим уравнение (2.16) для плотности (11.42):Z ξ0Z ξ0Z afξ (u) du = α0 , илиH fξ(1) (u) du = α0 −H fξ(1) (u) du, илиaAZAξ0fξ(1) (u) du = α0 [F̃ (b) − F̃ (a)] − F̃ (a) + F̃ (A), илиAξ0 = ψ (1) α0 [F̃ (b) − F̃ (a)] − F̃ (a) + F̃ (A)Rb(здесь использовано соотношение 1/H = a fξ(1) (u) du = F̃ (b) − F̃ (a)),т.
е. для выборочного значения ξ0 случайной величины ξ с плотностью(11.42) удалось вывести моделирующую формулу, представляющую собой композицию элементарных функций от стандартного случайногочисла α0 ∈ U (0, 1). Утверждение 11.6 доказано.ПРИМЕР 11.7 (А; 0,5 балла; [9, 24]). Пусть требуется построить алгоритм моделирования выборочного значения ξ0 случайной величиныξ, имеющей плотность распределенияfξ (u) =λ e−λu, 0 < u < A, λ > 0.1 − e−λA192(11.43).С учетом примера 2.1 и определения 11.1, распределение (11.43)можно назвать усеченным экспоненциальным распределением.Это распределение широко используется в приложениях (например, онопозволяет реализовать так называемое блуждание без вылета при моделировании переноса частиц – см.
раздел 6 данного пособия).Алгоритм 11.14 здесь выглядит следующим образом.(1)АЛГОРИТМ 11.15. 1. Моделируем выборочное значение ξ0 случай(1)ной величины ξ (1) согласно табличной формуле (2.20): ξ0 = − lnλα0 .(1)2. Если ξ0 ≤ A, то в качестве выборочного значения ξ0 случайной(1)величины ξ принимаем ξ0 = ξ0 , иначе повторяем пункт 1 и т. д.Трудоемкость этого алгоритма пропорциональна величине s = 1−e1−λ A .При малых A это значение может быть достаточно большим.С другой стороны, непосредственно решая уравнение (2.16) методаобратной функции распределения, по аналогии с примером 2.1 получаем следующую моделирующую формулу:ln 1 − α0 1 − e−λ Aξ0 = −.λ−λ AПроверка 2.1 для α0 = 0 дает ξ0 = − ln 1 −0×1−e/λ = 0, а−λ A−λ Aдля α0 = 1 имеем ξ0 = − ln 1 − 1 × 1 − e/λ = − ln e/λ = A.Последнее соотношение для ξ0 ненамного сложнее формулы (2.20),и не требует процедуры исключения, как в алгоритме 11.15.