1626435387-5d0ec7dd22f55dfcc65f3209f497fb0b (844202), страница 11
Текст из файла (страница 11)
Здесь функции fj (x) являются плотностямираспределения на гиперповерхностях Sj , а числа {pj } являются вероятPAностями (т. е. pj > 0 и j=1 pj = 1).64Учитывая то, что для x ∈ Sj выполнено равенство fξ (x) = pj fj (x),перепишем интеграл (4.13) в виде (4.1):ZAXg(x)δ(Y(x))jjI = G(x) pj fj (x) dx =pf(x)jjj=1Z=AXG(x)gj (x)δ (Yj (x))pj fj (x)j=1 f (x) dx = Eζ.ξЗдесьζ=AXG(ξ)gj (ξ)δ (Yj (ξ))pj fj (ξ)j=1,(4.14)и случайный вектор ξ распределен согласно плотности fξ (x).Таким образом, можно построить алгоритм (4.1) (или алгоритм 1.2),причем для моделирования выборочных значений ξ i следует использовать метод дискретной суперпозиции (см., например, [9, 13], а такжераздел 11 данного пособия).АЛГОРИТМ 4.2.
Моделируем номер mi согласно вероятностям {pj };при этом выборочное значение ξ i численно реализуется на гиперповерхности Smi согласно плотности fmi (x). Соответствующий вклад встатистическую оценку (4.1) равенζi =G(ξ i )gmi (ξ i ).pmi fmi (ξ i )По аналогии с утверждением 4.1 несложно показать, что минимальная дисперсия Dζ оценивателя (4.14) достигается, когда плотность fξ (x)берется в видеfξ (x) = HAX|G(x)| gj (x) δ (Yj (x)) ,j=1где H – нормирующая константа.655. Методы уменьшения дисперсии весовогооценивателя интеграла5.1.
Метод выделения главной части. Следующие рассуждения дают один из наиболее эффективных методов уменьшения дисперсии Dζ из соотношения (4.2); здесь ζ – весовой оцениватель (монтекарловская оценка) интеграла I из соотношения (4.1); см., например,[9, 13].Предположим, что существует функцияRg (0) (x), близкая к g(x) дляx ∈ X ⊆ Rd и такая, что интеграл I (0) = X g (0) (x) dx может бытьвычислен аналитически.Тогда для повышения эффективности алгоритма метода Монте-Карло(4.1) (или алгоритма 1.2) можно использовать стандартную технологиювыделения главной части (см., например, [9, 13]), основанную на соотношенииZ(0)I=I +g(x) − g (0) (x) dx.XПрименим стандартный весовой алгоритм (4.1) (или алгоритм 1.2)для второго слагаемого последней суммы. Строим соответствующийоцениватель и реализуем приближение интеграла:ng ξ (0) − g (0) ξ (0)1 X (0)(0)(0)(0),ζ , где ζ = I +I = Eζ ≈n i=1 ifξ (0) ξ (0)а случайный вектор ξ (0) распределен согласно плотности fξ (0) (x).Последняя плотность fξ (0) (x) может весьма значительно отличатьсяот плотности fξ (x) из исходного алгоритма (4.1).Например, в случае, когда интегрирование ведется по многомерномуединичному кубу X = ∆(d) = (0, 1)d и g(x) ≈ g (0) (x) (это означает, чторазность g(x) − g (0) (x) равномерно мала, т.
е. близка к малой константев кубе ∆(d) ), вместо плотности (4.7) для получения малой дисперсииоценивателя ζ (0) можно взять плотность случайной величины ξ (0) = α,равномерно распределенной в кубе ∆(d) , выборочные значения которойлегко моделируются на ЭВМ.Так или иначе, при g(x) ≈ g (0) (x) дисперсия случайной величины(0)ζможет быть достаточно малой.ПРИМЕР 5.1 [5, 9]. Рассмотрим следующую тестовую задачу. Исследуем алгоритмы приближенного вычисления интеграла по единичному66двадцатимерному кубу ∆(20) с известным значением:Z 1Z 1(1)(20)I=...ex ...x− 1 dx(1) . .
. dx(20) ≈ 9, 54 · 10−7 .00Возьмем в качестве fξ (x) плотность равномерного распределения вдвадцатимерном единичном кубе ∆(20) : fξ (x) ≡ 1, x ∈ ∆(20) . Тогда,согласно формуле (4.1), имеемI = Eζ,ζ = eα(1)...α(20)− 1,где α(i) ∈ U (0, 1); i = 1, ..., 20 – стандартные случайные числа.Несложно подсчитать, что Dζ ≈ 3, 0 · 10−10 .Учитывая вид разложения экспоненты в ряд Тейлора в окрестностинуляu2eu = 1 + u ++ o u2 ,(5.1)2целесообразно выделить главную частьg (0) x(1) , . . . , x(20) = x(1) × . . . × x(20) ,при этом(1)I (0) = 2−20 и ζ (1) = 2−20 + eα...α(20)− 1 − α(1) × .
. . × α(20) ;здесь fξ (0) (x) ≡ fξ (x).Заметим, что, с учетом соотношения (5.1),Dζ (0) ≈14Z1Z..01x(1) × .. × x(20)40dx(1) ..dx(20) =5−20≈ 2, 5 · 10−15 ,4то есть дисперсия уменьшается на пять порядков. Описание примера 5.1 закончено.При построении интегрируемой функции g (0) (x) можно использовать классические (чаще всего – кусочно-полиномиальные) приближения подынтегральной функции g(x) (в этом случае метод выделенияглавной части можно причислить к дискретно-стохастическим алгоритмам численного интегрирования [9, 17]). Здесь, в отличие от подходов, связанных с выборкой по важности (см.
подраздел 4.2), не нужныдополнительные требования «моделируемости» (т. е. функция g (0) (x) не67должна быть пропорциональной плотности распределения случайноговектора ξ (0) , для которого имеются эффективные алгоритмы численного моделирования).5.2. Интегрирование по части области. Рассмотрим следующийаналог метода выделения главной части.Пусть, как и в соотношении (4.1), интеграл представлен в видеZZg(x)I=.g(x) dx =q(x) fξ (x) dx; q(x) =fξ (x)XXПредположим, что можно аналитически вычислить следующие интегралы по некоторому подмножеству X2 области интегрированияX ⊆ Rd :ZZq(x)fξ (x) dx = I2 иfξ (x) dx = i2 ,X2X2при этом 0 < i2 < 1; здесь мы считаем, что fξ (x) = 0 при x ∈/ X.В этом случае целесообразно представить исходный интеграл I ввидеZZI=Xq(x)fξ (x) dx = I2 +Z= I2 + i1 ×X1X1q(x)fξ (x) dx =q(x)fξ (1) (x) dx = Eζ (1) ,гдеX1 = X\X2 ,i1 = 1 − i2 ,ζ (1) = I2 + i1 q ξ (1)и ξ (1) – случайный вектор, распределенный в подобласти X1 согласноплотности усеченного распределения fξ (1) (x) = fξ (x)/i1 ; x ∈ X1 .Соответствующий алгоритм метода Монте-КарлоI ≈ I2 +ni1 X(1) q ξin i=1(5.2)называется методом интегрирования по части области (см.,например, [9, 13]).Следующее утверждение показывает целесообразность использования метода (5.2) вместо (4.1) для уменьшения дисперсии оценивателя(и, как правило, трудоемкости) весового метода Монте-Карло.68УТВЕРЖДЕНИЕ 5.1 (см., например, [5, 9]).
Справедливо следующеенеравенство Dζ (1) ≤ i1 × Dζ.ДОКАЗАТЕЛЬСТВО. Используя соотношения вида (4.2), а такжесвойства интеграла и дисперсии (см., например, [14]), имеемZi1 × Dζ = i1q 2 (x)fξ (x) dx − I 2 =XZ= i1q (x)fξ (x) dx + i1X1Dζ(1)=i21 DqξZ2(1)=i21ZZX12 Zq (x)fξ (1) (x) dx− i12X1= i1q 2 (x)fξ (x) dx − i1 I 2 ,X2X1q 2 (x)fξ (x) dx −q(x)fξ (1) (x) dx=2ZX1q(x)fξ (x) dx.Рассмотрим разностьi1 × Dζ − Dζ(1)Z2= i1X22Zq (x)fξ (x) dx − i1 I +X12q(x)fξ (x) dxи покажем, что она неотрицательна. Заметим, что2ZX1где ∆ =q(x)fξ (x) dxRX2q(x) −I2i2= (I − I2 )2 иZX2q 2 (x)fξ (x) dx = ∆ +I22,i22i1 × Dζ − Dζ (1)fξ (x) dx ≥ 0. ПоэтомуI22= i1 ∆ +− i1 I 2 + (I − I2 )2 = i1 ∆+i2I2+ 2 + i2 I 2 − 2II2 = i1 ∆ +i2√I√2 − i2 Ii22≥ 0.Утверждение 5.1 доказано.Если подобласть X2 близка к области интегрирования X, то методинтегрирования по части области превращается в аналог метода выделения главной части (здесь g (0) (x) = g(x) для x ∈ X2 и g0 (x) ≡ 0иначе).69Тем не менее выгода от применения метода интегрирования по частиобласти имеется и в случае, когда множество X2 значительно меньше,чем X, но в этом случае уменьшение дисперсии Dζ может быть не такимзначительным.Если имеется алгоритм моделирования (реализации на ЭВМ) выборочного значения ξ 0 случайного вектора ξ ∈ X согласно плотностиfξ (x), то для случайного вектора ξ (1) ∈ X1 ⊂ X можно использоватьследующий алгоритм метода исключения для моделирования усеченного распределения.АЛГОРИТМ 5.1 (см., например, [9], а также алгоритм 11.14 из подраздела 11.6 данного пособия).
1. Моделируем ξ 0 согласно плотностиfξ (x).(1)2. Если ξ 0 ∈ X1 , то ξ 0 = ξ 0 , иначе снова повторяем пункт 1 ит. д.При этом метод интегрирования по части области практически недает выигрыша, так как уменьшение дисперсии вида i1 ×Dζ сочетается снеобходимостью моделирования в среднем P{ξ ∈ X}/P{ξ ∈ X1 } = 1/i1выборочных значений вектора ξ.Более предпочтительным является моделирование выборочного зна(1)чения ξ 0 случайного вектора ξ (1) непосредственно в области X1 согласно плотности fξ (1) (x) (см. подраздел 11.6 данного пособия).Приведем пример применения метода интегрироваания по части области.ПРИМЕР 5.2 [5, 9, 17]. Пусть требуется вычислить объем Ḡ трехмерной фигуры G, ограниченной поверхностью (в сферических координатах (r, φ, θ))r = a + h t(φ, θ), где − 1 ≤ t(φ, θ) ≤ 1 и a − h > 0.Обозначим через X2 и X вписанный в G и описанный около G шары,радиусы и объемы которых равны соответственно a − h, a + h и X̄2 , X̄(рис.
5.1). Заметим, что искомая величина равна интегралуZḠ =χ(G) x(1) , x(2) , x(3) dx(1) dx(2) dx(3) ,Xгде χ(G) x(1) , x(2) , x(3) – индикатор множества G.Рассмотрим сначалавесовой алгоритм 1.2 со случайным векторомξ = ξ (1) , ξ (2) , ξ (3) , имеющим плотность равномерного распределения70Рис. 5.1.
Иллюстрация к примеру 5.2в области X:31=, x = x(1) , x(2) , x(3) ∈ X.fξ (x) ≡34π(a + h)X̄Этот алгоритм дает приближениеḠ = Eζ ≈nX̄ X (G) (1) (2) (3) χξi , ξi , ξi ; ζ = X̄χ(G) (ξ).n i=1(5.3)Моделирующие формулы для выборочных значений компонент случайного вектора ξ имеют вид (см., например, [9], а также подраздел 6.4данного пособия):(1)ξi= ρi sin ϕi sin θ̃i ,(2)ξi= ρi cos ϕi sin θ̃i ,1/3(3)ξi= ρi cos θ̃i ;(5.4)ϕi = 2πα2,i , cos θ̃i = 1 − 2α3,i ; αj,i ∈ U (0, 1). (5.5)(1) (2) (3) для выЗаметим, что при получении значений χ(G) ξi , ξi , ξiчисления приближения (5.3) переход к декартовым координатам (5.4)можно не осуществлять, а лишь проверять условиеρi ≤ a + h t ϕi , θ̃i ,(5.6)ρi = (a + h) α1,i ,71(1) (2) (3) при выполнении которого χ(G) ξi , ξi , ξi=(G) (1) (2) (3)χξi , ξi , ξi= 0.Выделим теперь объем X̄2 = 34 π(a − h)3 шара X2 :ZḠ = X̄2 +χ(G) (x) dx,1,аиначеX1n22где X1 = X\X2 =x(1) , x(2) , x(3) : (a − h)2 < x(1) + x(2) +o2+ x(3) < (a + h)2 .Рассмотрим соответствующую модификацию метода интегрирования по части области:nX̄1 X (G) (1,1) (1,2) (1,3) (1)χξi , ξi , ξi; ζ = X̄2 +X̄1 χ(G) ξ (1) ,n i=1(5.7)где случайная точка ξ (1) = ξ (1,1) , ξ (1,2) , ξ (1,3) распределена равномернов X1 согласно плотностиḠ = Eζ (1) ≈ X̄2 +fξ (1) (x) =13, x ∈ X1 .=34π [(a + h) − (a − h)3 )]X̄1По аналогии с реализацией равномерно распределенной случайнойточки в шаре с помощью перехода к сферическим координатам можнополучитьформулыдлямоделированияслучайныхточек(1,1) (1,2) (1,3) ξi , ξi , ξi, равномерно распределенных в шаровом слое X1 .
Этиформулы совпадают с соотношениями (5.4) и (5.5) с той лишь разницей,что по-другому вычисляется компонента ρi :q(1)(1)(1) 3ρi = α1,i (a + h)3 + 1 − α1,i (a − h)3 .(5.8)(1,1) (1,2) (1,3) При подсчете χ(G) ξi , ξi , ξiиз (5.7) можно не использовать(1)преобразование (5.4), а лишь проверять условие (5.6) для ρi = ρi .Из формул (5.3) – (5.8) следует, что вычислительные затраты на реализацию приближений (5.3) и (5.7) практически одинаковы, в то времякакDζ = D X̄χ(G) (ξ) > Dζ (1) = D X̄2 + X̄1 χ(G) (ξ (1) ) ,так какD X̄χ(G) (ξ) = X̄ Ḡ − Ḡ2 = Ḡ(X̄ − Ḡ),72D X̄2 + X̄1 χ(G) (ξ (1) ) = X̄1 (Ḡ − X̄2 ) − (Ḡ − X̄2 )2 = (Ḡ − X̄2 )(X̄ − Ḡ);здесь использовано равенство X̄1 + X̄2 = X̄.Таким образом,D X̄2 + X̄1 χ(G) (ξ (1) )Ḡ − X̄2X̄ − X̄2=≤=(G)D(X̄χ (ξ))ḠX̄24π336h(a + h)23 (a + h) − (a − h)≤.=4π3(a − h)33 (a − h)Последняя величина имеет асимптотику 6ha при h → 0.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
