Исследование и уменьшение дисперсии весовых оценок в методе Монте-Карло тема диссертации и автореферата по ВАК РФ 01.01.07, кандидат физико-математических наук Медведев, Илья Николаевич

• Медведев, Илья Николаевич
• кандидат физико-математических наук
• 2005

Оглавление диссертации кандидат физико-математических наук Медведев, Илья Николаевич

Введение

1. Оптимизация весовых методов Монте-Карло по вспомогательным переменным

1.1. Модификация фазового пространства и весовой оценки

1.2. Оптимизация моделирования по части переменных

2. Ценностное моделирование

2.1. Эффективность ценностного моделирования начального распределения

2.2. Моделирование начального распределения процесса для некоторых задач теории переноса.

2.3. Частичное ценностное моделирование элементов траектории

2.4. Метод исследования дисперсии весовой оценки.

3. Исследование дисперсии весовых оценок в методе "экспоненциального преобразования"

3.1. Асимптотическая оптимизация распределения длины свободного пробега.

3.2. Исследование дисперсии весовых оценок.

3.2.1. Одномерный вариант.

3.2.2. Сферический вариант.

3.2.3. Численные расчеты для сферического варианта.

Введение диссертации (часть автореферата) на тему «Исследование и уменьшение дисперсии весовых оценок в методе Монте-Карло»

Разработка алгоритмов численного статистического моделирования в настоящее время имеет особое значение в связи с возможностью их идеального распараллеливания путем распределения численных статистических испытаний по отдельным процессорам. Практически статистическое моделирование наиболее часто используется для решения задач физики и техники, в основе которых лежат вероятностные модели, связанные с некоторыми цепями Маркова [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<p* -h) + К;9л (0.6) где К* - оператор с ядром к2(х,х')/р(х,х').

Введем спектральный радиус р(К) по формуле р{К)= р{К*) = Ы || К*11 Ц1/", п

Известно [2], что если р(Кр) < 1 и /2/7г G N\(X), то D£ < +оо. Также известно [17],[2], что, если h{x) > 0 и к(х>,х)<р*(х) f(x)<p*(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. Проведен детальный анализ конечности дисперсий весовых оценок для сферического варианта экспоненциального преобразования.

Список литературы диссертационного исследования кандидат физико-математических наук Медведев, Илья Николаевич, 2005 год

1. Михайлов Г.А. Построение весовых методов Монте-Карло на основе увеличения размерности фазового пространства // Доклады РАН - 2003. - Т. 389, Л*2 4 -С.461-464.

2. Ермаков С.М., Михайлов Г.А. Статистическое моделирование. М, Наука, 1982.

3. Фролов А.С., Ченцов Н.Н. Решение трех типичных задач теории переноса методом Монте-Карло.- В кн.: Метод Монте-Карло в проблеме переноса излучения. М., Атомиздат, 1967, С. 25-52.

4. Михайлов Г.А. Оптимизация весовых методов Монте-Карло. М.: Наука,1987 Engl. transl.: Springer-Verlag, 1992].

5. Марчук Г.И., Михайлов Г.А., Назаралиев М.А. и др. Метод Монте-Карло в атмосферной оптике. Новосибирск: Наука, 1976.

6. Дэвисон Б. Теория переноса нейтронов. М.: Атомиздат, 1960.

7. Михайлов Г.А., Медведев И.Н. Оптимизация весовых методов Монте-Карло по вспомогательным переменным // Сиб. матем. жури. 2004 - Т. 45, №2 - С. 399409.

8. Михайлов Г.А., Медведев И.Н. Эффективность "ценностного" моделирования в методе Монте-Карло. // Доклады РАН. 2005 - Т.401, Л* 1. - С. 16-20.

9. Михайлов Г.А., Медведев И.Н. "Метод исследования дисперсии весовой оценки численного статистического моделирования". // Доклады РАН. в печати.

10. Михайлов Г.А. Весовые методы Монте-Карло. Новосибирск: Издательство СО РАН, 2000. - 248 С.

11. Medvedev I.N., Mikhailov G.N. Optimization of weighted methods for part of variables // Russ. J. Numer. Anal. Math. Modelling. 2003. - Vol.18, 5, - P. 397-412

12. Medvedev I.N., Mikhailov G.N. Study of "value" modelling efficiency in the Monte-Carlo Method // Russ. J. Numer. Anal. Math. Modelling. 2005. - Vol.20, № 2, - P. 185-207

13. Medvedev I.N. The variance of weight estimate in the case of partial "value" modelling // Proceedings of the 5th St.Petersburg Workshop on Simulations. N11 Chemistry St.Petersburg University Publishers, 2005, P.465-470.

14. Михайлов Г.А. Весовые Алгоритмы статистического моделирования. -Новосибирск: Изд. ИВМиМГ СО РАН, 2003.

15. Ермаков С.М., Некруткин В.В., Сипин А. С. Случайные процессы для решения класических уравнений математической физики. М, Наука, 1984.

16. Соболь И.М. Численные методы Монте-Карло. М, Наука, 1973.

17. Kalos М.Н. Monte-Carlo integration of the adjoint gamma-ray transport equation // Nuclear Sci. and Eng., 1968, P. 284-290

18. Медведев И.Н. Сравнение эффективности "ценностного" моделирования начального распределения в методе Монте-Карло. // Труды конференции молодых ученых. Новосибирск, 2004. - С.111-118

19. Марчук Г.И. Методы расчета ядерных реакторов. М.: Атомиздат, 1961.