ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2007, том 47, < 1, с. 162-173
УДК 519.676
СТАТИСТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ОДНОГО ТИПА ПАР СЛУЧАЙНЫХ ВЕЛИЧИН С ИСПОЛЬЗОВАНИЕМ "ФИКТИВНЫХ СКАЧКОВ"
© 2007 г. А. И. Хисамутдинов
(630090 Новосибирск, пр-т акад. Коптюга, 3, Ин-т геофиз. СО РАН)
e-mail: khisam@emf.ru
Поступила в редакцию 26.05.2005 г. Переработанный вариант 12.11.2005 г.
Рассматривается статистическое моделирование посредством техники отбора пары случайных величин (T, <), T е [0, < е ж", d > 1. Задано совместное распределение пары в форме, которая объединяет родственные задачи из разных областей; оно имеет вид
P{T е dt, < е du} = f(t, u)exp
-J J f( tU)m (dU) dt'
v 0
dtm( du),
ж"
где/и т - некоторые функция и мера соответственно. Первая из величин Т - хорошо известное случайное время ожидания. Конструируется метод розыгрыша пары (Т, ^1) вводом реализации вспомогательной марковской последовательности пробных пар, и рассматриваются его применения в переносе частиц и кинетике разреженных газов. Библ. 18.
Ключевые слова: статистическое моделирование, пара случайных величин, совместное распределение, техника отбора, марковская последовательность пробных пар, трудоемкость алгоритма, розыгрыш соударений частиц, метод Монте-Карло.
1. ВВЕДЕНИЕ
Идеи "отбора" являются основой многих эффективных методов для моделирования различных случайных величин (см., например, [1]-[3]). В настоящей работе эти идеи и представления используются для конструирования методов розыгрыша пары случайных величин (Т, Т е [0,
°и, е й > 1. Совместное распределение пары имеет вид
/ t
(i)
P{T е dt, < е du} = w (t, u)p(t, u)exp
J w (t') dt'
m( du )dt, (1.1)
где
™(*) = | t, и)р(t, и)т(йи), (1.2)
(
^(1)(', ■) - некоторая неотрицательная измеримая функция на [0, <8> (, т(йи) есть с-конечная мера на с-алгебре борелевых множеств из ( , а р(*, и) по переменной и, и е ( , является плотностью относительно меры т(йи) V* е [0, и есть измеримая функция по * \1и е ( . Важными в данной работе частными случаями т являются лебегова мера (т(йт) = йи) и считающая мера; последняя является атомической, и каждый атом наделен массой 1. Первая из величин Т - случайное время ожидания с показательным распределением
t
P{T > t} = exp
-J w (t') dt'
0 < t < (1.3)
а условное распределение d-мерной случайной величины Ч имеет вид
P{Ч е du\T = t} = w(l\t, и)p(t, и)m(du)/w(t), и e ?Ad. (1.4)
Например, пара (T, Ч) характерна для марковских скачкообразных процессов, в том числе в их применениях в переносе частиц и в кинетике разреженных газов; причем распределение (1.4) для Ч входит в вероятностную "смесь", задающую распределение значений скачка в момент T = t. Следует заметить, что о многомерных случайных величинах нередко говорят, как о случайных векторах или о случайных точках; данное замечание в случае d > 2 касается как Ч, так и других случайных величин, которые вводятся ниже. Эффективность методов розыгрыша распределения (1.1) зависит от свойств функций w(1) и p. В физических применениях функции w(1) и w можно охарактеризовать как парциальную и полную частоты взаимодействий.
Методы отбора и их разновидности применяются для моделирования распределений вида (1.3) и (1.4). В физике переноса частиц известен метод (дополнительного) 5-рассеяния, или "выравнивания полной частоты (полного сечения)", в применении к (1.3); как правило, его применяют, чтобы обойти трудности, связанные с вычислением интеграла J^w (к) dK. При этом сама
функция w считается легко вычислимой. Рассмотрение и формальное обоснование метода отбора для (1.3), в котором функция w мажорируется постоянной и разновидностью которого является названный метод выравнивания w, было дано в [4]. Описания в связи с переносом частиц метода выравнивания w и вообще приема "введение дополнительного 5-рассеяния" (и истоков последнего) имеются в [5]-[7] и др. Пример моделирования (1.4) в переносе частиц в случае, когда w(1) - функция вида
w(1)(t, и) = const] v0- и\, d = 3, и е v0 е (1.5)
а p является плотностью нормального (максвелловского) распределения, дается работой [8]; v0 -некоторый заданный вектор. Здесь посредством того самого метода отбора (в узком смысле) или исключения, которое было определено в [1], обходится трудность вычисления w, т.е. интеграла (1.2). В настоящей работе рассматривается общий абстрактный случай, когда функция w(-) не является "легко вычислимой" и поэтому последовательное моделирование T согласно (1.3), а затем Ч согласно условному распределению (1.4) не является эффективным; мы рассматриваем случай, когда более эффективными, т.е. менее трудоемкими, оказываются алгоритмы методов, связанных с непосредственным моделированием совместного распределения (1.1) с помощью техники отбора посредством марковской последовательности "пробных" пар.
В переносе частиц и кинетике разреженных газов способ дополнительного 5-рассеяния как метод выравнивания уже парциальной частоты w(1) применяется для непосредственного розыгрыша совместного распределения (1.1). Данный метод имеет ясную физическую интерпретацию, и именно на основе физических соображений он рассматривался в ряде работ, включая [9], [10]. Отметим характерные черты способов дополнительного 5-рассеяния для розыгрыша как (1.3), так и (1.1). В этих способах соответствующая частота взаимодействий является ограниченной и в качестве верхней границы для этой частоты взаимодействий назначается постоянная; разность между верхней границей и частотой трактуется как частота "фиктивных рассеяний", которые не изменяют скоростей, энергий и положений соударяющихся частиц, т.е. не изменяют соответствующей фазовой плотности частиц. Будем трактовать оба названных метода выравнивания как разновидности методов отбора (см. п. 2.3).
В настоящей работе, рассматривая указанный выше абстрактный случай, мы вводим некоторую общую конструкцию - вспомогательный марковский скачкообразный процесс и на его основе рассматриваем розыгрыш пары (T, Ч). Конструируется метод отбора, использующий некоторые функции wM и W, посредством которых задается скачкообразный процесс (и последовательность пробных пар); w(1)(-) < wM) (•), а w(-) < W(-). Заметим, что функция wM) может быть и
неограниченной. В разд. 3 и 4 выделены два подслучая вида функции w(1). В обоих подслучаях рассматриваются конкретные приложения
t
I p(t, и)m(du)< + ^, w(t)< const < + lim I w(к)dK = .
J t ^
md 0
Символ Р {•} обозначает везде вероятность соответствующего события. Далее множество интегрирования не указывается, если оно совпадает со всем или полным множеством определения подынтегрального выражения. Термины "легко вычислимая" функция, трудность вычисления и менее трудоемкий используются в данной работе в контексте вычислительной математики и математического моделирования и не будут более уточняться и формализовываться. В тексте их истолкование не вызывает неоднозначности.
2. ОБЩАЯ КОНСТРУКЦИЯ
2.1. Пусть мМ (г, и) и Ж(г) - измеримые функции на [0, ® ^ и [0, соответственно и такие, что
г, и) < мМ'(г, и), мм(г) = |г, и)р(г, и)т(ёи) < Ж(г) < Сш< +
где СЖ - некоторая постоянная. Введем теперь (вспомогательный) марковский скачкообразный
процесс {(аг, Уг) }0° с состояниями в {/, рН } ® 9Лё, где / и рН - символы, соответственно, фиктивного и истинного ("физического") скачков и где аг означает тип скачка, который совершается последним в интервале [0, г) (или непосредственно предшествует моменту г + 0). Процесс считается непрерывным справа, и полагается, что в начальный момент г = 0 задано состояние после
"нулевого" скачка г0 = 0 и аг = /. Пусть V и ^ - соответственно, случайные номер и момент скачка такие, что V > 1, > 0 и впервые а^ = рН, т.е. аг = ... = а^ = f, аtv-0 = /, а аtv = рН. В
процессе {(аг, Уг) }0° нас будут интересовать именно этот случайный момент ^ и значение У^, т.е. пара (гу, Ytv). Будем использовать далее обозначение У, = Уг , j > 1. Нам достаточно задавать процесс {(аг, Уг) }0° как марковскую последовательность (цепь) случайных пар
(го,Уо ),...,( ^У*), V > 1, (2.1)
тем самым фактически трактуя состояния со значением рН как поглощающие состояния. Переходная плотность вероятностей цепи (2.1), а именно плотность вероятностей события (г', аг,
Уг —► г, аг, Уг) при условии, что аг■ = /, конструируется в форме произведения плотностейр1 ир2:
р1(гг|а, = /)р2(а„ У^г).
Следствием этого определения является такое свойство: если аг = /, то моделирование Уг можно не производить, поскольку распределение следующей после (аг, Уг) пары не зависит от Уг. Учитывая данное свойство, случайную величину У0 в (2.1) можно, в частности, не доопределять.
Введем (вспомогательную) последовательность случайных величин { Т) 1 такую, что
О = 0_1 + т), j >1
Теперь в случае использования этих величин последовательность (2.1) и, главное, пара (гу, Уу) в ней определяются следующими положениями.
1. Все Tj, 1 <j, являются независимыми и одинаково распределенными случайными величинами,
/ Аг
V/ > 1 P{Tj > At|tj-1 = t' > 0, atj-1 = f} = exp
2. Vj > 1
-J W (t' + к) dK
v 0
At > 0.
P{v> j\tj = t } = P{at/ = at/-o = f\tj = t} = c (t) + c2 (t),
J J J J
ci(t) = 1 - wM(t)/W(t), c2(t) = WM(t)/W(t) - w(t)/W(t).
3. Все Ту, ] > 1, являются независимыми и одинаково распределенными случайными величинами,
WM( 1) и)
р{Т/е^ = 1} = 1т^мжр(и)т( ).
4а. У] > 1 P{v = ] = 1, Т = и} = Р{ а = рк = 1, Т, = и} = w(1)(t, и)/м>Ц) (1, и).
46. Р^ >]\] = 1, Т = и} = 1 - w(1)(t, u)/w{M) (1, и). Нетрудно видеть, что
с2( 1) = | P{v>/|1/ = 1,Т/ = и } Р{Т / е = 1} = |( 1- w(1)( и и)) р^ и) т^и)
и что
Г(0
Р{Т, е = 1} = ^ = 1- с (1).
Шаг 4а. С вероятностью w(1)(t, и)/wM1) (1, и) цепь обрывается, = рк и V = ]; пара (1У, Ту) прини-
Алгоритм моделирования перехода
(1/ = 1 ',а, = /,Т/ 1/,а1/,Т/), ] > 1, (2.2)
предлагаемый здесь, имеет следующий вид.
Шаг 1. = + Ту, Т' разыгрывается согласно ее плотности вероятностей + А()ехр[1' + + к)а?к]; полагается = 1 = 1' + А1.
Шаг 2. С вероятностью с1(1) полагается а1 = / (и V > ]), т.е. последовательность продолжается; с п. 1) начинается моделирование следующего (]
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.