Разработка алгоритмов численного статистического моделирования в настоящее время имеет особое значение в связи с возможностью их идеального распараллеливания путем распределения численных статистических испытаний по отдельным процессорам. Практически статистическое моделирование наиболее часто используется для решения задач физики и техники, в основе которых лежат вероятностные модели, связанные с некоторыми цепями Маркова [2],[15],[16]. В принципе, такие задачи можно решать непосредственно численно моделируя траектории этих цепей. На основе исходного «феноменологического» описания проблемы можно даже строить весовые модификации алгоритма, домножая вспомогательный «вес» после каждого элементарного перехода на отношение соответствующих (может быть обобщенных) плотностей исходного и моделируемого распределений. С другой стороны, моделирование траекторий частиц можно интерпретировать как алгоритм оценки функционалов от соответсвующего интегрального уравнения второго рода, ядро которого совпадает с плотностью перехода базовых цепей Маркова [2],[14]. Однако зачастую прямое моделирование не позволяет оценивать изучаемые величины с требуемой точностью. В этом случае можно использовать весовые методы Монте-Карло, состоящие в том, что на ЭВМ моделируется подходящая цепь Маркова, а требуемые функционалы оцениваются с помощью веса, который после очередного перехода в цепи домножается на отношение ядра интегрального уравнения к переходной плотности. Несомненный плюс рассмотрения соответствующих интегральных уравнений — возможность детального изучения эффективности различных весовых модификаций. В частности, на основе таких уравнений можно разрабатывать «ценностные» алгоритмы с малыми вероятностными погрешностями, что особенно важно для оценок малых вероятностей [2]. В последнее время было выяснено [1], что для построения весовых модификаций может быть особенно эффективным увеличение размерности фазового пространства путем включения в число фазовых координат моделируемых вспомогательных случайных величин. Соответствующая факторизация ядра базового интегрального уравнения и вспомогательной плотности перехода иногда дает требуемую весовую модификацию.
Изложим более подробно введенные понятия. Математическая модель ряда прикладных задач строится на основе рассмотрения некоторого скачкообразного обрывающегося с вероятностью единица однородного марковского процесса (см. например, [2]). При этом траектория процесса вполне определяется ее состояниями в моменты скачков, т. е. фактически можно рассматривать обрывающуюся однородную цепь Маркова с заданной переходной функцией P (x', S), где х G X — т-мерное евклидово пространство, S С X — измеримое по Лебегу множество. Для построения весовых алгоритмов моделирования целесообразно рассматривать соответствующую условной мере P (x', S) обобщенную субстохастическую плотность перехода к (х', х). Обобщенная плотность распределения к (х', х) определяется для любого х' 6 X равенством.
J h{x)P{x', dx) = J k{x', x) h{x)dx V/г € C0(X),.
X X где Cq (X) — множество непрерывных ограниченных функций. Отметим, что здесь и далее рассматриваются и ненормированные (т.е. не обязательно вероятностные) распределения.
Настоящая работа ориентирована на приложения в теории переноса частиц, в которой кроме измеримых плотностей, соответствующих абсолютно непрерывным распределениям, используются также «дельта-функции», означающие интегрирование по некоторым многообразиям меньшей сравнительно с т размерности (см., например, [2] раздел 2.1.1). Использовать обобщенные плотности (вместо интегрирования по соответствующим мерам) в теории статистического моделирования предложил Н. Н. Ченцов [3] в связи с тем, что такой подход упрощает построение и реализацию модификаций моделирования. Это важно для настоящей работы, так как в ней рассматриваются дополнительные фазовые переменные, причем базовые переменные — координаты и скорости — связаны с дополнительными функционально, т. е. их условные распределения определяются соответствующими «дельта-функциями». Предполагается, что.
J к (х', х) dx = q (x') <1−5, S> 0.
Величина q (x') имеет смысл вероятности необрыва траектории в заданной точке х'. Вследствие последнего неравенства цепь обрывается с вероятностью единица и среднее число состояний конечно. Пусть x0, xi, ., xN.
— однородная обрывающаяся цепь Маркова, определяемая плотностью f (x) распределения начального состояния xq и субстохастической обобщенной плотностью перехода к (хх). Здесь N — номер состояния, в котором реализуется обрыв траектории (иначе, момент остановки). Ясно, что обобщенная плотность распределения состояний, непосредственно следующих за начальным, выражается равенством J f (x')k (x', x) dx'= [Kf](x).
Следовательно, обобщенная плотность распределения среднего числа фазовых состояний цепи оо п=0 где фп (х) — плотность распределения состояний номера п, представляет собой ряд Неймана для следующего интегрального уравнения 2-го рода: ф) = J к{х', х)<�р (х') dx' + f{x), (0.1) или (р — К<�р + /. Это уравнение будем рассматривать в пространстве Ni (X) обобщенных плотностей мер ограниченной вариации, так как ряд.
I Knf, h I в силу условия q (x') < 1 — 5 сходится для любой h G Cq (X). n=o /.
Для заданной функции h G Со (Х) рассмотрим также сопряженное уравнение р* = К*(р* + h, (0.2) где [К*(р*](х) = f k (x, x')ip*(x') dx'. Предполагается, что.
К* е [С0(Х) — Со (Х)}. Справедливо представление оо оо п=0 п=0 где 5(-) — плотность локализованного в соответствующей точке источника.
Методы Монте-Карло обычно используются для оценки линейных функционалов вида оо tp (x)h (x) dx = ^(iT1/, h), he CQ (X). x.
Если реализуется прямое моделирование исходной цепи Маркова, то для оценки величины Ih используется соотношение n.
4 = Е ?>(*") 0.
Однако можно моделировать и другую, вспомогательную, цепь Маркова с плотностью перехода р (х', х), взаимно регулярной с к (х', х), и начальной плотностью п (х), взаимно регулярной с f (x). Кроме того, предполагается, что р (х', х) ф 0 и тх{х) ф 0 на носителях функций к (х', х) и f (x) соответственно. При выполнении этих «общих условий несмещенности» отношения к (х', х)/р (х', х) и J (x)/, k{x) имеют смысл и определяются отношениями измеримых сомножителей рассматриваемых функций. Это позволяет ввести вспомогательные веса, но формулам:
QoW = -г, Qn = Qn-1−7-г- (0.3).
Полагая n.
С =2Qnh (xn), п=0 имеем Ih — Е£ (см., например, [2]). Величина? представляет собой стандартную весовую оценку «по столкновениям». Справедливы следующие соотношения:
С = ix = h (x) + У) Qnh (xn). (0.4).
Заметим, что в теории переноса, если в точке столкновения величина h (xn) рассчитывается после выбора нового направления пробега частицы, то оценку (0.4) принято называть оценкой «по рассеяниям» [4].
Известно [2], что для оценки функционала можно построить большое количество различных оценок, среднее значение которых равно /д (например оценка по поглощениям). Общим недостатком таких оценок является то, что они включают в себя, как правило величины, которые сильно меняются в реальных задачах, либо трудно вычисляются. Практические расчеты показывают, что для задач теории переноса значительно более эффективно разрабатывать модификации основной оценки по столкновениям, чем применять другие виды оценок (за исключением задач о вычислении вероятности поглощения какого-либо типа) [2].
Функцию ц>*(-) в теории весовых методов Монте-Карло принято называть «функцией ценности» в связи с ее вероятностным представлением: n ip*(x) = = h{x) + E^Qn/i (xn), xQ = x, (0.5) n=l причем для оценки (0.4) величина Е^ определяется рядом Неймана для следующего уравнения [4] g = h (2.
Введем спектральный радиус р (К) по формуле р{К)= р{К*) = Ы || К*11 Ц1/", п.
Известно [2], что если р (Кр) < 1 и /2/7г G N (X), то D? < +оо. Также известно [17],[2], что, если h{x) > 0 и к (х>, х)<�р*(х) f (x).
Р[Х'Х) [К*<�р*](?') ' 1 ^ (/>*) ' 1 j то D? = 0. Отметим, что реализация оценки с нулевой дисперсией невозможна по следующим причинам: 1) изначально функция ip* неизвестна- 2) вследствие того, что вероятность обрыва тождественно равна нулю, необходимо моделировать бесконечные цепочки. Поэтому на практике используют приближенное моделирование по ценности, основанное на замене в характеристиках цепи (0.7) функции ip* на функцию р ^ </?* (замечательным свойством получаемого алгоритма является его независимость от постоянного множителя функции р). Кроме того, начиная с т — го состояния вводят поглощение с вероятностью, обеспечивающей обрыв траектории с вероятностью 1. Например, если вместо <�р*(х) здесь использовать функцию: т (х) = Ср*(х)[ 1 + б (я)], Ия)| < е> то при малом б величина D? мала [4]. Соответствующие алгоритмы объединяются под названием «моделирование по ценности» [2],[4].
Как правило, при увеличении размерности фазового пространства, переход х —> х' осуществляется в результате выбора совокупности значений вспомогательных случайных величин, причем, если соответствующие функции ценности определяются формулами типа (0.5), то ценностное моделирование всех элементарных вспомогательных переходов дает оценку с нулевой дисперсией [1]. Однако на практике осуществляется весовая модификация моделирования только части вспомогательных переменных, например в теории переноса наиболее часто используется модификация моделирования длины пробега на основе так называемого «экспонециального преобразования» [2]. При этом оказалось, что для некоторых задач частичная ценностная модификация вспомогательных переменных может увеличивать дисперсию по сравнению с прямым моделированием [12], а в некоторых случаях возникает вопрос о конечности дисперсии и, как следствие, величины трудоемкости s = D?, где t (?) — среднее время расчетов на ЭВМ для получения одного выборочного значения ?
Далее следует краткое содержание диссертационной работы по главам и параграфам.
Заключение
.
Сформулируем основные результаты диссертационной работы.
1. Решена задача минимизации дисперсий оценок линейных функционалов путем вариации плотности моделируемого распределения одной (первой по порядку) вспомогательной, может быть векторной, случайной величины. На примере практически важной задачи теории переноса представлена асимптотическая оптимизация распределения длины свободного пробега.
2. Получены условия, при выполнении которых ценностное моделирование начального распределения уменьшает дисперсию по сравнению с прямым. Используя полученные условия, проведен анализ эффективности ценностного моделирования начального распределения для некоторых задач теории переноса.
3. На примере модельной задачи теории переноса с «дельта» рассеянием и задачи оценки решения системы линейных алгебраических уравнений показано, что частичное ценностное моделирование может увеличивать дисперсию весовой оценки по сравнению с прямым моделированием.
4. Доказано, что дисперсия весовой оценки в случае «частично-ценностного» статистического моделирования конечна. На этой основе предложен метод решения вопроса о конечности дисперсии весовой оценки с использованием мажорантного сопряженного уравнения.
5. Рассмотрено применение предложенного метода исследования конечности дисперсии для анализа классического «метода экспоненциального преобразования» моделирования длины свободного пробега частицы.
6. Проведен детальный анализ конечности дисперсий весовых оценок для сферического варианта экспоненциального преобразования.
