Диссертация (1145371), страница 26
Текст из файла (страница 26)
.i=1При полуцелом ν = m + 1/2 — по формуле!mYγ(m + 1/2) = − lnαi + (Ξ[1])2 /2,m = 0, 1, . . . .i=1В общем случае, при 0 < δ < 1 используется формула γ(m + δ) = γ(m) + γ(δ), вправой части которой случайные величины должны быть независимыми. Случайную величину γ(δ) обычно моделируют методом отбора, используя мажорантуe−x , x > 1;g(x) = xδ−1 , 0 < x 6 1;0,x 6 0.Эффективность моделирования равна Γ(δ)/(1/δ + 1/e).Для моделирования случайной величины можно использовать формулуЙонка:γ(δ) = (ζν − 1)α,где величина ζν имеет бета-распределение с параметрами 1 − ν и ν.Заметим, что случайная величина γ(ν)/θ имеет при x > 0 плотностьраспределения (двухпараметрическое гамма-распределение)f (x) = θν xν−1 e−θx /Γ(ν).207Пусть p > 0,m > 0 — параметры.Z1B(p, m) =xp−1 (1 − x)m−1 dx =Γ(p)Γ(m)Γ(p + m − 1)0— бета-функция.
Бета-распределение имеет на отрезке (0, 1) плотностьxp−1 (1 − x)m−1fp,m (x) =.B(p, m)Моделирование произвольного бета-распределения выполняется в несколькоэтапов и подробно описано в [11]. Мы рассмотрим лишь те случаи, которыетребуются в данной работе. Отметим, что если β(p, m) имеет распределение спараметрами p, m, то 1 − β(p, m) имеет распределение с параметрами m, p.Случай 1.
m = 1. Моделирование выполняется по явной формуле:β(p, 0) = α1/p .Случай 2. p = 1. Моделирование выполняется по явной формуле:β(0, m) = 1 − α1/m .Случай 3. p, m — целые. Моделирование выполняется с помощью порядковых статистик при выборке объема p + m − 1 из равномерного распределенияна отрезке (0, 1). А именно, β(p, m) = αp — p-й член вариационного ряда.Случай 4. m — целое.β(p, m) =mY1/(p+k−1)αkk=1Случай 5. Целая часть [m] = 1.В этом случае распределение является смесьюfp,m (x) =∞Xpk fp+k,2 ,k=0где a = 2 − m иpk =1a(a + 1) . . .
(a + k − 1).B(p, m) k!(k + p)(k + p + 1)208Пусть дискретная случайная величина ν имеет распределение pk , k = 0, 1, . . . ,тогда1/(p+ν) 1/(1+p+ν)α2.β(p, m) = α1Случай 6. m < 1, p < 1. Этот случай сводится к предыдущему, так какfp,m =mp(x)fp+1,m (x) +fp,m+1 (x).p+mp+mИменно этот случай наиболее интересен. Он позволяет промоделировать гаммараспределение с нецелым показателем. В свою очередь можно моделироватьбета- распределение через гамма-распределение по формуле:β(p, m) =γ(p)γ([p]) + γ({p})=.γ(p) + γ(m) γ([p]) + γ({p}) + γ([m]) + γ({m})Входящие в данную формулу случайные величины должны быть независимыв совокупности.A.3Краевые задачи для параболического уравненияМоделирование всех случайных величин, необходимых для реализациирассмотренных алгоритмов, кроме геометрического распределения, уже обсуждалось.A.3.1Моделирование геометрического распределенияГеометрическое распределение на множестве натуральных чисел задает-ся рядом вероятностей pk = q(1 − q)k−1 , k = 1, 2, .
. .. Используя формулу длясуммы первых m членов геометрической прогрессии, легко получить следующую формулу для моделирования этого распределения:ln(α)N=+1 .ln(1 − q)209A.4Первая краевая задача для эллиптическогоуравненияA.4.1Моделирование величины ρ в блуждании по сферамДля моделирования случайной величины ρ, имеющей на отрезке [0, 1]плотность распределения pn (r) = 2n r − rn−1 /(n − 2), можно предствить ее ввиде смеси плотностей порядковых статистикpn (r) =n−2Xk=12nqk (r),(n − 2)(k + 1)(k + 2)где qk (r) = (k + 1)(k + 2)rk (1 − r) — плотность k + 1-й порядковой статистикипри выборке объема k + 2 из равномерного распределения на отрезке [0, 1].Суммы коэффициентов смеси находятся явно: S0 = 0, S1 = n/3(n − 2), .
. . , Sl == nl/(n − 2)(l + 2). Поэтому, по α, распределенному равномерно на отрезке[0, 1], выбираем номер ν порядковой статистики ν = [2(1 − α(1 − 2/n))−1 ] −1 и моделируем порядковую статистику ν + 1-ю порядковую статистику привыборке объема ν + 2. Для этого достаточно следить за двумя последнимичленами вариационного ряда.
Приведем фрагмент программы, выполняющейнеобходимые действия.α := random;β := random;if α > β then begin γ := α; α := β; β := γ end;for k := 1 to n − 2 dobegin γ := random;if γ > β then begin α := β; β := γ endelse if γ > α then α := γ;end;ρ := α;210A.4.2Моделирование блуждания по эллипсоидамИспользуя обозначения § 3.2.1, определим в области T (x) следующиефункции:p(x, y) = Ny L(y, x),p1 (x, y) = Ny L(y, x), при c(x) ≡ 0,p2 (x, y) = µ(x)(n − 2)p(σ)σ1−nn XnXaij (y)i=1 j=1∂σ ∂σ,∂yi ∂yj(A.4.1)µp3 (x, y) = 2 µ(x)(n − 2)p(σ)σ 1−n ,νp4 (x, y) = µ(x)(n − 2)p(σ)σ 1−n .Эти функции неотрицательные и связаны неравенствами:p(x, y) 6 p1 (x, y) 6 2p2 (x, y) 6 p3 (x, y).(A.4.2)Первое неравенство следует из равенства p(x, y) − p1 (x, y) = c(x)L(y, x)и неравенств c(x) 6 0, L(y, x) > 0.Второе неравенство является следствием выбора функции p(ρ), а именно,p(ρ) следует выбирать исходя из условияp2 (x, y) > |p1 (x, y) − p2 (x, y)| + kckµ(x)σ 2−nZRp(ρ)dρ,σкоторое обеспечивает требуемое неравенство.
Для этого в доказательстве теоремы 3.2.3 постоянные κ1 , κ2 , κ3 cледует выбирать при c(x) ≡ 0, а затем увеличитьκ3 на величину kck = maxx∈D |c(x)|.Третье неравенство является следствием оценкиn XnXi=1 j=1aij (y)∂σ ∂σ6 µ| grady σ(y, x)|2 =∂yi ∂yj|A−1 (x)(y − x)|2µ= µ −16(A (x)(y − x), A(x)A−1 (x)(y − x))ν211и второго неравенства.При c(x) ≡ 0, функция u(x) ≡ 1 является решением уравненияM u(x) = 0, поэтомуZ1 = u(x) =Zp1 (x, y)u(y)dy =T (x)p1 (x, y)dy.T (x)Следовательно, p1 (x, y) — плотность распределения по переменной y. В силулеммы 2.1.2,RZZpp4 (x, y)dy = σn det(A(x)) Ep4 (x, x + rΩ)rn−1 dr =0T (x)1=q(R)ZRp(r)dr = 1,0поэтому p4 (x, y) — плотность распределения по переменной y.Для получения следующей точки Y марковской цепи по текущей X выполняем следующие действия:– моделируем Y с плотностью p1 (X, y);– моделируем случайную величину α, распределенную равномерно на [0, 1];– если αp1 (X, Y ) < p(X, Y ), то переходим в точку Y ;– если αp1 (X, Y ) > p(X, Y ), то переходим в поглощающее состояние.Неравенства (A.4.2) позволяют использовать для моделирования точкиY метод отбора.
В каждом конкретном случае выбирается своя мажоранта.Опишем процедуру моделирования Y с помощью мажоранты p3 (x, y). Пустьфункция next(x) возвращает реализацию случайной величины Z с плотностьюраспределения p4 (x, z), тогда величина Y , имеющая плотность p1 (X, y) реализуется в циклеrepeat Y := next(X); α := random;until αp3 (X, Y ) < p1 (X, Y );212Эффективность данной процедуры равна ν/(2µ) и может оказаться весьма низкой.
Примеры более эффективных процедур моделирования Y имеютсяв [12].Отметим, что в силу леммы 2.1.2, вектор Z = X + rpA(X)Ω1 име-ет плотность распределения p4 (X, z), если r распределено на [0, R] с плотностью p(ρ)/q(R), а Ω1 — изотропный вектор. Аналогичная формула Y =p= X + r A(X)Ω1 верна и для моделирования случайного вектора Y с плотностью L(y, x)/h(x), но теперь величина r распределено на [0, R] с плотностьюrn−1pe(r) =h(x)(n − 2)q(R)ZRr2−n − ρ2−n p(ρ)dρ.rЛегко проверить, что такую плотность имеет произведение двух независимыхслучайных величин % и ϑ, распределенных на [0, R] и [0, 1], соответственно, сRRплотностями ρ2 p(ρ)/ ρ2 p(ρ)dρ и 2n(s − sn−1 )/(n − 2).0Построение несмещенной оценки для интегралаZF (x) =f (y)L(y, x)dyT (x)можно упростить, представив его в видеZRf (x + rF (x) = Epn−1A(x)Ω1 )r(n − 2)q(R)ZRr2−n − ρ2−n p(ρ)dρdr,r0поменяв порядок интегрирораванияZRF (x) = Ep(ρ)(n − 2)q(R)0Zρpr2−n − ρ2−n f (x + r A(x)Ω1 )rn−1 drdρ,0и выполнив линейную замену переменной r = sρZRF (x) = E0ρ2 p(ρ)2nq(R)Z10p2ns − sn−1 f (x + sρ A(x)Ω1 )dsdρ.n−2213Теперь видно, что несмещенной оценкой для F (x) является случайная величинаp%2 f (x + ϑ% A(x)Ω1 )/(2n), где случайные величины ϑ и % – независимые ираспределены с плотностями 2n(s − sn−1 )/(n − 2) и p(ρ)/q(R).Алгоритм моделирования величины ϑ представлен в §A.4.1.Выбор константы κ в блуждании по эллипсоидамДля реализации алгоритма блуждания по эллипсоидам необходимоиметь оценку снизу для константы κ в неравенствеZRp(σ) > κp(ρ)dρ.σВыполним необходимые для этого вычисления.
Запишем Ny L в виде1Ny L(y, x) =q(R)ZRσnn∂H ∂σp(σ) X X−aik (y)p(ρ)dρ · Ny H −q(R) i=1∂yi ∂ykk=1ZR− µ(x)C(y)ρ2−n p(ρ)dρ =σ1q(R)ZRp(ρ)dρ · (Ny H − C(y)H) + µ(x)(n − 2)p(σ)σ1−nnn XXi=1 r=1σZR+ µ(x)C(y)aik (y)∂σ ∂σ+∂yi ∂yk(σ 2−n − ρ2−n )p(ρ)dρ. (A.4.3)σОбозначим z = grady σ = A−1 (y − x)/σ. Отметим, чтоσ 2 = (A−1 (y − x), y − x) = (A−1 (y − x), A(x)A−1 (y − x)) 6 ν1 |A−1 (y − x)|2 ,√√поэтому |z| ν1 > 1. Аналогично, |z| ν 5 1.Запишем неравенствоn XnXi=1 r=1aik (y)∂σ ∂σ= (A(y)z, z) > ν|z|2 .∂yi ∂yk214Выполним вычисленияNy H − C(y)H =n XnXi=1 k=1nX∂H∂ 2HBk (y)−.aik (y)∂yi ∂yk∂ykk=1Пусть βk = maxx∈D |Bk (y)|, тогда верны неравенства nnX∂H X ∂H Bk (y)βk 6 |β||gradH|,6∂yk ∂yk k=1k=1√|gradH| = q(R(x))µ(x)(n − 2)σ 1−n |z| 6 q(R(x))µ(x) ν1 (n − 2)σ 1−n |z|2 ,Pв которых |β| = ( nk=1 βk2 )1/2 .Заметим, чтоn XnXaik (x)i=1 k=1поэтомуnn XXi=1 k=1∂ 2H= 0,∂yi ∂yknnXX∂ 2H∂ 2Haik (y)(aik (y) − aik (x))=.∂yi ∂yk∂y∂yiki=1k=1По теореме Лагранжаaik (y) − aik (x) =nX∂aik (y ∗ )l=1∂yl(yl − xl ).Таким образом, справедливы равенстваnn XXnnnX X X ∂aik (y ∗ )∂ 2H∂ 2Haik (y)=(yl − xl )=∂yi ∂yk∂y∂y∂yliki=1i=1 k=1n Xn XnX=i=1 k=1 l=1k=1 l=1∂aik (y ∗ )Aik (x)zi zk(yl − xl )q(R(x))µ(x)(n − 2)σ 1−n (−+n·)=∂ylσσ= q(R(x))µ(x)(n − 2)σ 1−n ×× −(grady (Sp(A(y ∗ )A−1 (x))), A(x)z) + n(A(x)z, grady (z > A(y ∗ )z) ,где при дифференцировании по переменной y зависимость z от y не учитывается.215Пусть матрица A0l (y ∗ ) составлена из частных производных матрицы A(y)по переменной yl и kA0l (y ∗ )k 6 k l .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
