1626435387-5d0ec7dd22f55dfcc65f3209f497fb0b (844202), страница 34
Текст из файла (страница 34)
В последнем соотношении использовано свойствогамма-функции (12.11).Следовательно,"#∞XB(i + µ, [ν] + 1) Γ(s + i) ui+µ−1 (1 − u)[ν](µ,ν)fβ(u) =×=B(µ, ν) i! Γ(s)B(i + µ, [ν] + 1)i=0203=∞X(i+µ,[ν]+1)pi fβ(u),(12.19)i=0где величины {pi } представляют собой вероятности и имеют вид=p0 =[ν]!;B(µ, ν)µ(µ + 1) × . . . × (µ + [ν])pi =Γ(i + µ) Γ([ν] + 1) Γ(s + i)=Γ(i + µ + [ν] + 1) B(µ, ν) i! Γ(s)(12.20)[ν]! s (s + 1) × .
. . × (s + i − 1); i = 1, 2, . . . ;B(µ, ν) i! (i + µ)(i + µ + 1) × . . . × (i + µ + [ν])здесь использованы свойства гамма- и бета-функций (12.6), (12.10),(12.11).Из формулы (12.19) получается следующий алгоритм метода дискретной суперпозиции.АЛГОРИТМ 12.2 (см., например, [9, 24]). 1. Моделируем стандартное случайное число α̃0 и, используя наиболее эффективный (экономичный) из методов моделирования целочисленной случайной величины сраспределением (10.5), (12.20) (алгоритмы 10.2, 10.6), получаем значение η0 = k.(µ,ν)2.
Моделируем выборочное значение β0согласно плотности(k+µ,[ν]+1)fβ(u) с использованием соответствующего варианта формулы (12.17):[ν]+1Y 1/(k+µ+i−1)(µ,ν)β0=αi.i=1Недостатком этого метода является то, что выбор номера k плотно(k+µ,[ν]+1)сти fβ(u) весьма трудоемок. Это связано с тем, что таких номеров k бесконечно много (и такие эффективные методы, как алгоритмы 10.5, 10.8, неприменимы) и скорость убывания последовательностисоответствующих им вероятностей pk с ростом k невелика.
Также определенные сложности вызывает то, что нужно вычислять бета-функциюот дробных переменных, причем делать это с высокой точностью [36].Согласно формуле (10.17), трудоемкость стандартного алгоритма 10.2 выбора значения k целочисленной случайной величины η по вероятностям{pi }изсоотношений(12.20)пропорциональнавеличине Eη.204Можно показать, что Eη = +∞ при [ν] = 0.
Поэтому в случаеµ > 1, 0 < ν < 1 следует воспользоваться заменой переменных w = 1−uи формулой (12.12) и строить алгоритм 12.2 для моделирования случайной величины β (ν,µ) (у которой параметры µ и ν меняются местами).В случае 0 < µ < 1, 0 < ν < 1 можно провести дополнительнуюрандомизацию и реализовать соответствующий метод суперпозиции наоснове соотношения(µ,ν)fβ=(µ,ν)(u) = u fβ(µ,ν)(u) + (1 − u) fβ(u) =νµ(µ+1,ν)(µ,ν+1)f(u) +f(u);µ+ν βµ+ν β(12.21)(µ,ν)здесь использован вид плотности fβ(u) (см.
соотношение (12.1)) исвойства (12.10) и (12.11).При моделировании бета-распределения с нецелыми параметрамиможно попытаться воспользоваться формулой (12.9) и использовать соответствующие моделирующие формулы для гамма-распределения.Определенным препятствием здесь является то обстоятельство, чтомоделирование гамма-распределения (12.2) с нецелым параметром ν также является проблематичным.В случае ν ∈/ N представляем параметр в виде суммы целой и дробной части: ν = [ν] + {ν} и, используя утверждение 12.2, представляемслучайную величину γ (λ,ν) в виде суммы двух независимых случайныхвеличинγ (λ,ν) = γ (λ,[ν]) + γ (λ,{ν})(12.22)(при 0 < ν < 1 первое слагаемое в этой сумме отсутствует).Для случайной величины γ (λ,[ν]) из соотношения (12.22) имеется моQ[ν](λ,[ν])делирующая формула вида (12.16): γ0= −(1/λ) × ln i=1 αi .Остается нерешенной проблема моделирования случайной величины γ (λ,{ν}) из соотношения (12.22), имеющей гамма-распределениес параметром, меньшим единицы.Рассмотрим сначала алгоритм Йонка, представляющий собой достаточно редкий пример применения метода интегральной суперпозиции (см.
подраздел 3.1 данного пособия).Заметим, что для любого t > 0 справедлива цепочка равенствZ +∞{ν}−1{ν}−1t=tfγ(1,1−{ν}) (y) dy =0205Z=0+∞t−1 e−y y −{ν} dy=t−{ν} Γ(1 − {ν})Z+∞e−tv0v −{ν}dv;Γ(1 − {ν})здесь использована замена v = y/t.Полагая t = λu, получаем представлениеZ +∞λ{ν} u{ν}−1 e−λuλ e−λuv −{ν}(λ,{ν})fγe−λuv(u) ==dv =Γ({ν})Γ({ν}) 0Γ(1 − {ν})Z=+∞fγ (u|v) fη (v) dv,(12.23)0гдеfη (v) =v −{ν}1×,Γ({ν}) Γ(1 − {ν})v+1fγ (u|v) = λ (v + 1) e−λu(v+1) ,v > 0,(12.24)u > 0.Представление (12.23) дает возможность применить метод интегральной суперпозиции (см.
алгоритм 3.1).АЛГОРИТМ 12.3 (алгоритм Йонка; см., например, [24]). 1. Численно моделируем стандартное случайное число α1 и реализуем выборочное значение η = η0 согласно плотности (12.24).(λ,{ν})2. Моделируем значение γ0согласно экспоненциальной плотности fγ (u|η0 ) = λ(η0 + 1) e−λu(η0 +1) (см. пример 2.1 и соотношение(2.20)):− ln α2(λ,{ν})γ0=.λ (η0 + 1)Для моделирования величины η0 согласно плотности (12.24) в пункте 1 сформулированного алгоритма 12.3 рассмотрим преобразованиеw = g(v) = v/(v + 1).
Обратное преобразование имеет вид v = w/(1 − w).Тогда, согласно утверждению 6.2, случайная величина β̃ = g(η) имеет плотность−{ν}w01−ww1××=fβ̃ (w) =wΓ({ν}) Γ(1 − {ν})1−w+11−w=w−{ν} (1 − w){ν} (1 − w)w(1−{ν})−1 (1 − w){ν}−1=.2Γ(1 − {ν}) Γ({ν}) (1 − w)B(1 − {ν}, {ν})206«Круг замкнулся»: мы получили плотность бета-распределения с параметрами (1 − {ν}) и {ν}, т. е. самый «неприятный» случай, где требуется дополнительная рандомизация вида (12.21) и т. п.Тем не менее, если удается смоделировать выборочное значение β̃0случайной величины β̃, то, учитывая, что η0 = β̃0 /(1 − β̃0 ), получаетсямоделирующая формула(λ,{ν})γ0=(β̃0 − 1) ln α2− ln α2=.λ (η0 + 1)λСуммируя рассуждения данного раздела, констатируем, что оба метода суперпозиции для бета- и гамма-распределений с нецелыми параметрами (алгоритмы 12.2 и 12.3) являются неэффективными (затратными, трудоемкими).
Здесь рекомендуется использовать мажорантные методы исключения, представленные в следующем подразделе 12.4.12.4. Случаи нецелых параметров: методы исключения.Начнем с гамма-распределения. Учитывая соотношения (12.16), (12.22),займемся построением мажорантного метода исключения (алгоритм 11.10) для моделирования случайной величины γ (λ,{ν}) ,имеющей гамма-распределение с параметром {ν}, меньшимединицы.Для функции g(u) = u{ν}−1 e−λ u , пропорциональной плотности (12.2),Г. А.
Михайловым была предложена следующая составная мажоранта(см., например, [23]): {ν}−1uпри 0 < u < 1,g(u) ≤ g (1) (u) =e−λ uпри u ≥ 1.В свою очередь, для моделирования случайной величины ξ (1) с составной плотностью, пропорциональной функции g (1) (u) (см. соотношение (11.18)) в примере 11.3 из подраздела 11.3 данного пособия былпостроен алгоритм 11.5 метода обратной функции распределения.Отсюда получаем следующий алгоритм мажорантного метода исключения.АЛГОРИТМ 12.4 (алгоритм Г. А. Михайлова; см., например,[9, 23, 24]). 1. Моделируем выборочное значение согласно плотностиfξ(1) (u) = C g (1) (u):(1/{ν}(1)α1 (λ + {ν}e−λ )/λпри α1 ≤ λ/(λ + {ν}e−λ ),ξ0 =−λ−1−(1/λ) ln((1 − α1 )(e + λ{ν} )) иначе207(1) (см.
алгоритм 11.5). Моделируем также значение η0 = α2 g (1) ξ0 .(1) (λ,{ν})(1)2. Если выполнено условие η0 < g ξ0 , то полагаем γ0= ξ0 ,иначе повторяем пункт 1 и т. д.Согласно формуле (11.36), трудоемкость алгоритма 12.4 пропорциональна величинеR +∞ (1)g (w) dwλ {ν}−1 + e−λ=s = 0R +∞.λ Γ({ν})g(w) dw0√Например, для {ν} = 1/2 имеем s = (2λ + e−λ )/(λ π). Заметимтакже, что при λ = 1 величина s монотонно растет от s = 1 (при {ν} ↓ 0)до s = 1 + e−1 ≈ 1, 36 (при {ν} = 1). В частности, при {ν} = 1/2 имеемs ≈ 1, 33.Заметим также, что алгоритм Г. А. Михайлова (алгоритм 12.4) заметно эффективнее (экономичнее) алгоритма Йонка (алгоритм 12.3).Это связано в том числе и с тем, что в алгоритме Йонка требуется моделировать бета-распределение (12.1) с нецелыми (а конкретнее,заключенными между нулем и единицей) параметрами µ и ν.В подобных случаях эффективными оказываются различные варианты мажорантного метода исключения (алгоритм 11.10) длямоделирования случайной величины β (µ,ν) , имеющей бета-распределение с нецелыми параметрами µ и ν.Плотность бета-распределения (12.1) пропорциональна функцииg(u) = uµ−1 (1 − u)ν−1 , 0 < u < 1.(12.25)Способы построения мажорант для функции (12.25) весьма разнообразны, причем эти построения и эффективность соответствующих алгоритмов метода исключения существенно зависят от значений параметров µ и ν.Так, для описанного выше «сложного» случая 0 < µ < 1, 0 < ν < 1можно построить мажорантуg(u) ≤ g (1,1) (u) = uµ−1 + (1 − u)ν−1 .Действительно,uµ−1 (1 − u)ν−1 = [u + (1 − u)]uµ−1 (1 − u)ν−1 == uµ (1 − u)ν−1 + uµ−1 (1 − u)ν ≤ (1 − u)ν−1 + uµ−1 .208Здесь использовано то обстоятельство, что при 0 < u < 1 и t > 0выполненоut < 1, (1 − u)t < 1.(12.26)Тогда для случая 0 < µ < 1, 0 < ν < 1 можно реализовать следующий алгоритм метода исключения.АЛГОРИТМ 12.5 (см., например, [9, 24]).
1. Моделируем выборочное(1)значение ξ0 случайной величины ξ (1) , распределенной согласно плотностиµν (1,1)g(u) = p1 f1 (u) + p2 f2 (u);fξ(1) (u) = C g (1,1) (u) =µ+νp1 =νµ, f1 (u) = µuµ−1 ; p2 =, f2 (u) = ν(1 − u)ν−1µ+νµ+ν(здесь используется модифицированный метод суперпозиции – алго1/ν(1)α1 (µ+ν)µ, иначе, то ξ0=ритм 11.3): если α1 < µ+νν1/ν(1)ξ0 = 1 − α1 (µ+ν)−ν.µ(1) 2. Моделируем также величину η0 = α2 g (1,1) ξ0 .(1) (ν,µ)(1)3.
Если η0 < g ξ0 , то β0= ξ0 , иначе повторяем пункты 1 и 2 и т. д.Трудоемкость алгоритма 12.5 пропорциональна величинеR 1 (1,1)g(w) dwµ+ν(1).s = 0R 1=µνB(µ, ν)g(w) dw0Например, для µ = ν ≈ 1/2 имеем s(1) ≈ 4/π ≈ 1, 27.Теперь рассмотрим случаи, когда параметры µ и ν достаточно велики.Для случая µ ≥ ν > 1 с учетом соотношений (12.26) в качествемажоранты функции (12.25) можно взятьg(u) ≤ g (1,2) (u) = uµ−1 (1 − u)[ν]−1 , 0 < u < 1;(12.27)здесь [A] – соответственно целая часть числа A.Плотность fξ(1) (u), пропорциональная функции g (1,2) (u) из соотношения (12.27), является плотностью бета-распределения с целым вторым параметром.
Для моделирования случайной величины ξ (1) следуетиспользовать формулу (12.17).209Отсюда возникает следующий мажорантный метод исключения.(1)АЛГОРИТМ 12.6. 1. Моделируем выборочное значение ξ0 по форQ(1)1/(µ+i−1)[ν]муле ξ0=. Моделируем также величинуi=1 αi(1)(1,2)η0 = α̃0 gξ0 .(1) (1) {µ}(µ,ν)(1)2.
Если η0 < g ξ0(или α̃0 < 1 − ξ0), то β0= ξ0 , иначеповторяем пункт 1 и т. д.Аналогично для случая 1 < µ ≤ ν выбираем мажорантуg(u) ≤ g (1,3) (u) = u[µ]−1 (1 − u)ν−1 , 0 < u < 1.(12.28)Плотность fξ(1) (u), пропорциональная функции g (1,3) (u) из соотношения (12.28), является плотностью бета-распределения с целым первымпараметром. Для моделирования случайной величины ξ (1) следует использовать формулу (12.18).Отсюда возникает следующий мажорантный метод исключения.(1)АЛГОРИТМ 12.7.
1. Моделируем выборочное значение ξ0 по форQ[µ] 1/(ν+i−1)(1)муле ξ0= 1 − i=1 αi. Моделируем также величину(1,3) (1)η0 = α̃0 gξ0 .(1) (1) {µ}(µ,ν)(1)2. Если η0 < g ξ0(или α̃0 < ξ0), то β0= ξ0 , иначеповторяем пункт 1 и т. д.Трудоемкости алгоритмов 12.6 и 12.7 пропорциональны величинамs(2) = B(µ, [ν])/B(µ, ν) и s(3) = B([µ], ν)/B(µ, ν) соответственно.