1626435387-5d0ec7dd22f55dfcc65f3209f497fb0b (844202), страница 27
Текст из файла (страница 27)
10.3. Схема алгоритма УолкераПроведем на высоте h = 1/M линию, параллельную основанию диаграммы. Часть столбцов диаграммы имеет высоту, большую h, а часть– меньшую h.Далее проводим перераспределения частей столбцов диаграммы длятех столбцов, для высот которых выполнены неравенства11или p̃i < h =;(10.20)p̃i > h =MM163столбцы с высотой p̃i = h остаются неизменными. Перераспределениепроисходит таким образом, чтобы все столбцы имели высоту h и в каждом из них были не более чем две доли исходных столбцов (рис.
10.3).Процедура включает порядка M шагов. Если на данном шаге имеются столбцы, высоты которых удовлетворяют соотношениям (10.20),то выбирается столбец (например, s-й), для которого p̃s < h (такойPMстолбец должен найтись, иначе i=1 p̃i > 1). Если же столбцов, удовлетворяющих неравенствам (10.20), нет, то процедура перераспределениязаканчивается.Кроме s-го столбца, для которого p̃s < h, выбираем еще один столбец(например, l-й), высота которого p̃l > h (такой столбец тоже долженPMнайтись, иначе i=1 p̃i < 1).В силу того, что p̃l + p̃s > p̃l > h, мы можем дополнить s-й столбецдо высоты h частью l-го столбца (и делаем это). При этом p̃l – высотаl-го столбца – может остаться больше, чем h, или стать меньше, чем h(т. е.
для него останется выполненным одно из соотношений (10.20), ион участвует в дальнейших шагах процедуры перераспределения), илиравно h (в этом случае столбец не участвует в процедуре перераспределения).В результате на описываемом шаге формируется ячейка Es = (Fs , ls ),включающая число Fs = M p̃s и номер ls того столбца, из которого произошло дополнение s-го столбца до высоты h.После описанного шага s-й столбец выбывает из процедуры перераспределения вероятностей, так как для него уже не выполнены строгиенеравенства (10.20).После не более чем M шагов все ячейки Ei , i = 1, . . . , M , будутсформированы (для столбцов, имевших исходную высоту h, вторые компоненты этих ячеек пусты), и подготовительная работа заканчивается(рис. 10.3).
Суммарный объем требуемой оперативной памяти равен 2M(отсюда и ограничение M . N0 /2).АЛГОРИТМ 10.8 (см., например, [9, 24]). 1. Моделируя первое стандартное случайное число α1 ∈ U (0, 1) и используя формулуm = [α1 M + 1] (это аналог соотношения (10.18)) выбираем номер ячейки Em = Fm ; lm .2. Моделируем второе стандартное случайное число α2 ∈ U (0, 1).Если α2 < Fm , то ξ0 = xm , иначе ξ0 = xlm .Обоснование алгоритма 10.8 можно провести на основании утверждения 2.3.
То, что происходит в алгоритме 10.8, можно трактовать164(1)(2)как розыгрыш равномерной точки β0 = M α1 , β0 = hα2 в прямоугольнике ∆ единичной площади с длинами сторон M и h, получающемся на диаграмме после процедуры перераспределения (рис. 10.3).Для каждого i (i = 1, . . . , M ) суммарная площадь частей исходногоi-го столбца, являющихся подмножествами ∆, равна pi , и поэтому(вследствие утверждения 2.3) в алгоритме 10.8 воспроизводится распределение (10.1): P{ξ = xi } = pi .ЗАМЕЧАНИЕ 10.2 [9].
Алгоритм 10.8 допускает следующую модификацию: на втором шаге алгоритма вместо α2 ∈ U (0, 1) можно взятьвеличинуα1 − (m − 1)/M= M α1 − m + 1.(10.21)β0 =1/MОбоснование такой замены следует из того, что случайная величина β0равномернораспределенана полуинтервале [0, 1) для любогоα1 ∈ (m − 1)/M, m/M , что, в свою очередь, следует из утверждения 2.4 (см. также обоснование модифицированного метода суперпозиции из подраздела 11.2 – см. замечание 11.3).Целесообразность применения формулы (10.21) связана с тем, что реализация нового выборочного значения α2 с помощью обращения к генератору типа RAN DOM является относительно трудоемкой операцией(см.
замечание 9.2).Что касается сравнения квантильного метода (алгоритм 10.6) и метода Уолкера (алгоритм 10.8), то здесь многочисленные тестовые и прикладные вычисления показали, что, без учета подготовительной работы(формирования массива {Xj }Kj=1 вида (10.19) в одном случае и массиMва {Ej = (Fj , lj )}j=1 – в другом), алгоритм Уолкера является болееэкономичным. Тем не менее, к несомненным достоинствам квантильного алгоритма следует отнести относительную простоту формированиямассива (10.19) и определенную «универсальность», связанную с возможностью применения метода для больших M (вплоть до M = ∞).10.8.
Анализ использования специальных алгоритмов (напримере моделирования геометрического распределения).В подразделе 10.6 показано, что специальный алгоритм 10.5 моделирования равномерного дискретного распределения превосходит по эффективности (экономичности) стандартный алгоритм 10.1.Специальные алгоритмы моделирования выборочных значений имеются и для целого ряда других дискретных распределений (10.1). Однако далеко не всегда эти алгоритмы являются заведомо более эффек-165тивными по сравнению со стандартным алгоритмом.Для иллюстрации этого тезиса в качестве примера рассмотрим целочисленную случайную величину η из (10.5), имеющую геометрическое распределение (это частный случай биномиального распределения; см., например, [14]):pi = P{η = i} = pq i−1 , 0 < p < 1; q = 1 − p; i = 1, 2, 3, .
. .(10.22)Соответствующая случайная величина η определяет количество испытаний Бернулли (с вероятностью «успеха» p) до получения первого«успеха».Соотношение (10.22) представляет собой пример распределения дискретной случайной величины с бесконечным числом значений (см. подраздел 10.4). Стандартным методом моделирования случайной величины η является алгоритм 10.4.
При его реализации для пересчета вероятностей целесообразно использовать мультипликативную рекуррентнуюформулу (10.14) при r(i + 1) = pi+1 /pi ≡ q.Согласно формуле (10.17), затраты алгоритма 10.4 для распределения (10.22) равны величине t1 = a + (b + c)/p.Это, в свою очередь, следует из того, чтоEη =∞Xi=1i pi =∞Xi p q i−1 =i=11.p(10.23)Действительно, разлагая функцию 1/(1 − z)2 в ряд Тейлора в точкеz0 = 0 и полагая z = q, получаем∞∞XXp11i p q i−1 =iq i−1 , и тогда Eη === .22(1 − q)(1−q)pi=1i=1В попытках уменьшить величину t1 можно сформулировать и обосновать следующий специальный алгоритм моделирования геометрического распределения.АЛГОРИТМ 10.9 (см., например, [9, 24]). Моделируем стандартноеслучайное число α0 ∈ U (0, 1) и полагаемln α0η0 =+ 1;(10.24)ln(1 − p)здесь, как и ранее, [A] обозначает целую часть числа A.166Докажем, что случайная величина η̃ = [ln α/ ln(1 − p)] + 1 имеет распределение (10.22), т.
е. pk = P{η̃ = k} = p q k−1 для k = 1, 2, . . ..Используя определение целой части числа и соотношение (9.6), получаемln αP{η̃ = k} = P k − 1 ≤<k =ln(1 − p)= P k ln(1 − p) < ln α ≤ (k − 1) ln(1 − p) == P ln(1 − p)k < ln α ≤ ln(1 − p)k−1 = P (1 − p)k < α ≤ (1 − p)k−1 == (1 − p)k−1 − (1 − p)k = (1 − p)k−1 [1 − (1 − p)] = p q k−1 ,что и требовалось доказать.Сравнивая алгоритмы 10.4 и 10.9 для геометрического распределения (10.22), можно отметить, что несмотря на определенную «компактность» формулы (10.24) неверно говорить, что алгоритм 10.9 всегда экономичнее алгоритма 10.4.
Для вероятности «успеха» p, близкойк единице, затраты t1 из (10.17) относительно невелики, а в формуле(10.24) для любого p применяется трудоемкая операция логарифмирования. Поэтому для выбранного компьютера и данного языка программирования можно экспериментально найти число p0 , для которого приp ≥ p0 более экономичным является алгоритм 10.4, а при p < p0 – алгоритм 10.9.
Проведенные нами расчеты дали p0 ≈ 0, 75.Описанная выше связь случайной величины η с комбинацией незави(s)симых бернуллиевскихвеличин {γ} с вероятностью «успе (s) случайных(s)ха» p (здесь P γ = 1 = p; P γ = 0 = q для всех s = 1, 2, ...)позволяет моделировать выборочное значение η0 с помощью реализа(s)ции соответствующей комбинации выборочных значений {γ0 }; такойалгоритм называется методом браковки.Так, для геометрического распределения(s)η0 = min s : γ0 = 1 ,(10.25)т. е. метод браковки состоит в последовательной проверке неравенстваαs < p (см.
алгоритм 10.3) до тех пор, пока оно не окажется верным.В силу замечания 9.2, методы браковки неэффективны из-за необходимости реализации большого количества стандартных случайных чисел {αs }, поэтому метод браковки (10.25) следует «забраковать».16711. Альтернативные алгоритмы широкогоприменения: метод дискретной суперпозициии мажорантный метод исключения11.1. Метод дискретной суперпозиции. В подразделе 3.1 былпредставлен алгоритм метода интегральной суперпозиции (рандомизации) (алгоритм 3.1) для моделирования выборочного значенияξ 0 случайного вектора (случайной величины) ξ, имеющей плотностьраспределения видаZfξ (u) =fη (v)fξ (u|v) dv.(11.1)VВ алгоритме 3.1 сначала моделировалось выборочное значение η 0вспомогательного вектора (величины) η, имеющего плотность распределения fη (u), согласно эффективному(экономичному) алгоритму (формуле) η 0 = ψη ᾱ1 , а затем получалось требуемое значениеξ 0 = ψξ ᾱ2 ; η 0 , распределенное согласно плотности fξ uη 0 (здесьтакже предполагалось наличие соответствующего эффективного алгоритма или формулы для каждого фиксированного η 0 ); при этомᾱ1 , ᾱ2 ∈ U (0, 1) – нужные наборы стандартных случайных чисел.В замечании 3.2 было сказано о том, что наиболее частым являетсяприменение алгоритма 3.1 в случае, когда вспомогательный случайныйвектор η представляет собой одномерную целочисленную случайнуювеличину η с распределением (10.5).В этом случае, формально подставляя в формулу (11.1) аналог соPMотношения (10.2) вида fη (v) = i=1 pi δ(v − i), получаемfξ (u) =MXpi fi (u); fi (u) = fξ (u|η = i);M ≤ ∞.(11.2)i=1Согласно требованиям, при которых применим алгоритм метода суперпозиции, должны существовать эффективные алгоритмы получениявыборочных значений целочисленной случайной величины η (в рассматриваемом частном случае следует применять наиболее эффективный изалгоритмов, описанных в разделе 10, – как минимум, стандартный алгоритм 10.2) и вектора ξ согласно плотности fi (u) для любого i.Таким образом, соотношение (11.2) представляет собой взвешеннуюсумму (смесь) эффективно моделируемых плотностей {fi (u)}.168Метод суперпозиции для плотности (11.2) формулируется следующим образом.АЛГОРИТМ 11.1 (см., например, [9, 13]).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
