1626435387-5d0ec7dd22f55dfcc65f3209f497fb0b (844202), страница 9
Текст из файла (страница 9)
также алгоритм 1.2) и ζ (K) = ζ.(1) (2) В этом случае выборочные значения ξ i = ξ i , ξ iполучаются со(1)гласно многомерной версии алгоритма 2.1: значения ξ i моделируют(2)ся согласно плотности fξ (1) (u), а значения ξ i – согласно плотностям (1) fξ (2) vξ i .Равенство I = Eζ (K) (см. соотношение (3.21)) может быть обосновано следующим образом (см. также[9]).Рассмотрим вектор ~ξ(2)=(2)ξ1 , . . .
, ξ(2)(2) , где ξ k – независимые,ξодинаково распределенные (как ξ (2) ) случайные векторы (величины).Используя формулу полного математического ожидания и свойство аддитивности математического ожидания (см., например, [14]),получаемEζ (K) = Eξ (1) E~ (2) ζ (K) ξ (1) =ξ(1)PK ξ(1) (2) (1)Eqξ,ξ(2)ξk=1k ξ~== Eξ (1) K ξ (1)KhiEξ (2) q ξ (1) , ξ (2) ξ (1) × K ξ (1)== Eξ (1) K ξ (1)hi= Eξ (1) Eξ (2) q ξ (1) , ξ (2) ξ (1) = Eq(ξ) = I.Теперь рассмотрим вопрос о выборе оптимального числа K ξ (1) ,минимизирующего трудоемкость S (K) = t(K) × Dζ (K) (см.
также [9]).51Согласно утверждению 3.1 получаемDζ (K) = Dξ (1) E~ (2) ζ (K) ξ (1) + Eξ (1) D~ (2) ζ (K) ξ (1) .ξξИспользуя свойства математического ожидания и дисперсии (см.,(2)например, [14]), а также независимость компонент вектора ~ξ , получаемDξ (2) q(ξ)ξ (1).Dζ (K) = Dξ (1) Eξ (2) q(ξ)ξ (1) + Eξ (1) K ξ (1)Пусть t1 обозначает среднее время моделирования одного выбороч(1)ного значения ξ i случайного вектора (случайной величины) ξ (1) , a(1)t2 ξ– среднее время получения на ЭВМ одного выборочного значе(2)ния ξ i,k случайного вектора (случайной величины) ξ (2) при фиксиро(1)ванном ξ (1) = ξ i . Тогда среднее время получения выборочного значе(K)ния ζiслучайной величины ζ (K) равноhit(K) = t1 + Eξ (1) K ξ (1) t2 ξ (1) .Для простоты дальнейших выкладок найдем оптимальное K ξ (1)для весьма распространенного на практике простейшего варианта метода расщепления, для которого K ξ (1) ≡ K = const.
В этом случаеt(K) = t1 + Kt2 , Dζ (K) = A1 +A2,Kгде t2 = Eξ (1) t2 ξ (1) иA2 = Eξ (1) Dξ (2) q(ξ)ξ (1) .A1 = Dξ (1) Eξ (2) q(ξ)ξ (1) ,Найдем минимум функцииA2S(r) = (t1 + t2 r) × A1 +rдля r > 0.Вычислим производнуюA1 t2S (r) = 2r052A2 t 1r −A1 t 22,(3.22)учтем положительность величин A1 , A2 , t1 , t2 и переменной r и получимточку минимумаrA2 t 1rmin =.(3.23)A1 t 2Таким образом, в качестве оптимального параметра K для простейшего варианта метода расщепления целесообразно использоватьближайшее к rmin целое число:rA2 t 1Kopt ≈.(3.24)A1 t 2Другая возможность выбора оптимального K – считать это значениецелочисленной дискретной случайной величиной κ с двумя значениямиK (0) и K (0) + 1, такими, что K (0) < rmin < K (0) + 1, и распределениемP{κ = K (0) } = K (0) + 1 − rmin и P{κ = K (0) + 1} = rmin − K (0) .
Такаявеличина κ имеет наименьшую дисперсию среди всевозможных целочисленных случайных величин κ̃, имеющих математическое ожиданиеEκ̃ = rmin .Соответствующий алгоритм моделирования выборочных значенийи обоснование указанного свойства случайной величины κ приведеныдалее в подразделе 10.3.Некоторая проблема, связанная с использованием величин (3.23) и(3.24), состоит в том, что величины t1 , t2 , A1 и A2 заранее не известны.Однако они могут быть оценены с помощью предварительных расчетов следующим образом.Реализуя алгоритм 3.3 для двух параметров K1 , K2 и относительнонебольшого числа испытаний n = n̂, получаем приближения t̃(K1 ) , t̃(K2 )времен t(K1 ) , t(K2 ) (по формулам вида (1.19)) и приближения D̃ζ (K1 ) ,D̃ζ (K2 ) дисперсий Dζ (K1 ) , Dζ (K2 ) (по формулам вида (1.20) или (1.22)).Тогда с учетом соотношений (3.19) можно рассмотреть следующуюсистему приближенных линейных уравнений относительно неизвестныхвеличин t1 , t2 , A1 и A2 : (K1 )t̃≈ t1 + K 1 t2D̃ζ (K1 ) ≈ A1 + A2 /K1t̃(K2 ) ≈ t1 + K2 t2D̃ζ (K2 ) ≈ A1 + A2 /K2Решая эту систему, получаем приближения требуемых значенийt1 ≈t̃(K1 ) − t̃(K2 )K2 t̃(K1 ) − K1 t̃(K2 ); t2 ≈;K2 − K1K2 − K1531K2 D̃ζ (K2 ) − K1 D̃ζ (K1 ) ;K2 − K1K1 K2 (K1 )A2 ≈D̃ζ− D̃ζ (K2 ) .K2 − K1На практике метод расщепления применяется также в случае, когдаξ (1) и ξ (2) (а значит, и q(ξ)) имеют дискретное распределение, при этомсоответствующие интегралы превращаются в суммы.
Приведем характерный пример из теории переноса частиц (излучения).ПРИМЕР 3.3 [9]. Пусть требуется оценить вероятность P прохождения частицы через слой вещества {x, y, z : 0 ≤ z ≤ H}. Моделируемна ЭВМ n траекторий блуждания частиц в слое, полагая, что каждаячастица движется в веществе прямолинейными «пробегами» случайнойдлины; в конце каждого пробега с некоторой вероятностью она можетпоглотиться или рассеяться по случайному закону (см.
далее раздел 6данного пособия). Источник частиц расположен на плоскости z = 0.Рассмотрим бернуллиевскую случайную величину ζ, которая равнаединице, если частица вылетает из слоя через плоскость z = H (приэтом P{ζ = 1} = P ), и нулю, если частица поглощается или вылетаетиз слоя через плоскость z = 0. Стандартный алгоритм приближениявеличины P = Eζ описывается формулой (1.1).Модификация стандартного алгоритма (простейший вариант методарасщепления) состоит в том, что фиксируется точка первого пересечения частицей плоскости z = z0 ; 0 < z0 < H (которая называется плоскостью расщепления); из этой точки «испускается» K «новых» независимых частиц, для которых результат прохождения слоя учитываетсяс «весом» 1/K.Рассмотрим случайную величину ξ (это аналог ξ (1) в двумерном случае), которая равна единице, если частица пересекла плоскость расщепления (вероятность этого события обозначим p1 ), и нулю иначе.Введем также случайную величину η (это аналог ξ (2) ), которая равна единице, если «новая» частица, выпущенная из точки расщепления,достигает плоскости z = H (вероятность этого события обозначим p2 ),и нулю иначе.
Заметим, что P = p1 × p2 .Для введенных таким образом случайных величин ξ и η выполняются соотношенияA1 ≈P{ξ = 1} = p1 , P{ξ = 0} = 1 − p1 , P{η = 1|ξ = 0} = 0,P{η = 0|ξ = 0} = 1, P({η = 1|ξ = 1} = p2 , P{η = 0|ξ = 1} = 1 − p1 .(3.25)54Рассмотрим также функцию0 при v = 0,q(u, v) =1 при v = 1.Несложно понять, что P = p1 × p2 = Eq(ξ), где ξ = (ξ, η), и длявычисления этого математического ожидания можно применить методрасщепления (3.21).Исследуем вопрос о выборе оптимальных параметров z0 и K, используя построенную выше теорию оптимизации простейшего вариантаметода расщепления.
Применяя формулы (3.25), вычислимA1 = Dξ Eη q(ξ)|ξ = p22 p1 (1 − p1 ) = P × p2 (1 − p1 ),A2 = Eξ Dη q(ξ)|ξ = p2 (1 − p2 )p1 = P × (1 − p2 ).Предположим теперь, что t2 = t1 p1 , то есть среднее время моделирования одной траектории «новой» частицы равно среднему времениблуждания до расщепления.В этом случае величина трудоемкости S (K) равна1 − p2(K).S= t1 (1 + Kp1 ) × P × p2 (1 − p1 ) +KОптимальное значение равноsrA2 t 11 − p2;=Kopt ≈A1 t 2P × (1 − p1 )(здесь, как и ранее, знак «≈» означает, что берется ближайшее к этомузначению натуральное число).Подставляя полученное значение в S (K) , получаемp2p(K)S0 = t1 × P ×p2 (1 − p1 ) + p1 (1 − p2 ) .(K)Рассмотрев S0как функцию от p1 > 0 (с учетом того, чтоp2 = P/p1 ), несложно√получить, что минимум этой функции достигается при p1 = p2 = P .
Такой выбор p1 , p2 и соотношение t2 = t1 p1определяют оптимальный параметр z0 .Таким образом, минимум трудоемкости равен√√ (K)Sopt = 4t1 P P 1 − P .55Максимальный выигрыш от введения расщепления равенS (1)(K)Sopt√ 2√ t1 1 + P P (1 − P )1+ P√√ =√=;4t1 P P 1 − P4 Pздесь S (1) – трудоемкость алгоритма блуждания частиц без расщепления.Если вероятность P мала, то√K ≈ 1/ P иS (1)(K)Sopt1≈ √ ,4 Pт. е. выигрыш от расщепления может быть весьма ощутимым. Описаниепримера 3.3 закончено.Для ряда приложений (прежде всего для задачи переноса излучения – см., в частности, только что рассмотренный пример 3.3) иногдаиспользуется метод многократного расщепления, когда переменное xразделяется на более чем две части (рис. 3.2). Однако для этого методаэкономия времени вычисления не всегда окупает значительное усложнение расчетных программ и трудности в определении параметров.Рис.
3.2. Пример использования многократного расщепления в теориипереноса частиц (излучения)564. Принцип выборки по важности4.1. Теорема о минимальной дисперсии. Рассмотрим снова весовой метод Монте-Карло для приближения интеграла"#ZZg(x)1 g(ξ 1 )g(ξ n )fξ (x) dx = Eζ ≈+ ... +;I=g(x) dx =n fξ (ξ 1 )fξ (ξ n )XX fξ (x)(4.1)здесь X ⊆ Rd и fξ > 0 при g(x) 6= 0 (см. также алгоритм 1.2 и формулу(3.12)).Напомним, что при оптимизации алгоритма (4.1) следует минимизировать величину трудоемкости S = t × Dζ (см. соотношение (1.18)), гдеt – среднее время получения выборочного значения ζi = g(ξ i ) fξ (ξ i )оценивателя (монте-карловской оценки) ζ = q(ξ) = g(ξ)/fξ (ξ) на ЭВМ,а Dζ – это дисперсия:ZZg 2 (x) dx− I 2 (4.2)Dζ = Eζ 2 − (Eζ)2 =q 2 (x)fξ (x) dx − I 2 =XX fξ (x)(см.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
