CC BY
29
Вычислительные технологии
Том 16, № 4, 2011
Модификации геометрических программных модулей, связанные с построением моделируемых вероятностных плотностей*
М.В. Бессмельцев, A.B. Войтишек Институт вычислительной математики и математической геофизики,
Новосибирск, Россия e-mail: [email protected]
Представлен программный модуль А1Тпск,$ ОеотИапс1от, позволяющий получать выборочные значения случайных векторов, распределенных согласно заданной плотности в геометрически сложных областях. Приведены полезные для пользователей модуля методики по конструированию вероятностных плотностей, допускающих эффективную численную реализацию выборочных значений. Рассмотрены численные схемы, связанные с применением полиномиальных и кусочно-полиномиальных приближений вероятностных плотностей (сопряженные с использованием метода дискретной суперпозиции при реализации выборочных значений). В качестве примера применения предлагаемых технологий рассмотрен итерационный дискретно-стохастический алгоритм построения адаптивных сеток.
Ключевые слова: численное моделирование случайных величин и векторов, технологии конструирования моделируемых плотностей, моделируемые приближения плотностей, итерационный дискретно-стохастический алгоритм построения адаптивных сеток.
Введение
В настоящее время широкое применение находят алгоритмы и пакеты программ, в которых реализуются множества точек, распределенных в сложных областях согласно заданной плотности распределения. К таким задачам относятся визуализация объектов на ЭВМ, построение адаптивных сеток и др. В данной работе предложен программный модуль А1Тпскв СеотНапйот (см. [1] и раздел 9), в котором использован ряд принципиальных положений по оптимизации соответствующих алгоритмов численного моделирования. В частности, существенным является то обстоятельство, что при применении кусочно-полиномиальных (чаще всего кусочно-постоянных) аппроксимаций функции плотности возможна реализация эффективных модификаций соответствующего метода суперпозиции (квантильного метода, метода Уолкера и др. — см., например, [2], а также раздел 6 настоящей работы).
Во многих ситуациях для более "быстрого" получения выборочных значений случайно распределенных точек вместо приближения плотности целесообразно выбирать моделируемые вероятностные распределения, допускающие построение эффективных
*Работа выполнена при финансовой поддержке РФФИ (гранты № 10-01-00040 и 09-01-00035).
алгоритмов численного моделирования (эта возможность также предусмотрена в программном модуле А1Тпскз СеотКапёот). В разделах 1-4 работы представлены соответствующие технологии "искусственного" конструирования вероятностных распределений, допускающих эффективное моделирование методами: обратной функции распределения, моделирования двумерного случайного вектора с зависимыми компонентами, интегральной и дискретной суперпозиции, исключения (соответствующие алгоритмы описаны, например, в [2]), Это является основой создания "банка" моделируемых вероятностных распределений с целью использования его при построении эффективных алгоритмов численного статистического моделирования.
Далее (в разделах 5-7) сделан обзор подходов, реализованных, в том числе, в программном модуле А1Тпскз СеотКапйот и связанных с приближением (как правило, полиномиальным или кусочно-полиномиальным) вероятностных плотностей, сопряженных с применением соответствующего метода дискретной суперпозиции при реализации выборочных значений. Подробнее рассмотрены практически важные случаи кусочно-постоянного приближения плотности и моделирования равномерного распределения в сложных областях и на сложных поверхностях,
В качестве примера применения представленных технологий в разделе 8 рассмотрен итерационный дискретно-стохастический алгоритм построения адаптивных сеток [3], Наконец, в заключение сформулированы основные результаты работы,
1. Технология вложенных замен
Стандартный алгоритм численной реализации выборочного значения £0 непрерывной случайной величины распределенной па интервале (а, Ь) согласно непрерывной, монотонно возрастающей на (а, Ь) функции распределения ^(х) (метод обратной функции распределения — см., например, [2]), основан на соотношении
где а0 — "стандартное случайное число" (т.е. выборочное значение случайной величины а, равномерно распределенной на интервале (0,1)), На ЭВМ выборочные значения а случайной вели чины а реализуются с помощью соответствующих генераторов (подпрограмм) типа RAND или RANDOM.
При использовании алгоритма (1) возникает "программистская" проблема выражения обратной функции F-1(x) в элементарных функциях, С учетом того что эта проблема рассматривается для практически важных случаев абсолютно непрерывных распределений, описываемых кусочно-непрерывными плотностями распределения f (u) случайных величин £ (см., например, [2]), встает вопрос о разрешимости уравнения
относительно верхнего предела интеграла в элементарных функциях. Здесь могут возникнуть трудности, связанные со взятием интеграла в левой части, а также с аналитическим разрешением уравнения (2) (когда первообразная интеграла из левой части (2) выражается в виде композиции элементарных функций), В случае существования
£о = F-1(a о)
(1)
(2)
относительно простого для программирования выражения для решения уравнения (2) вида (1)
Со = ф(ао) (3)
плотность Д (и) и сама формула (3) называются элементарными [2], В связи с тем, что элементарные плотности нужны как в других общих алгоритмах численной реализации случайных величин и векторов (имеются в виду методы суперпозиции, исключения и специальные методы), так и в многочисленных приложениях метода Монте-Карло (см., например, [2], а также разделы 2-4), возникает проблема расширения спектра таких плотностей. Практически неограниченные возможности конструирования элементарных плотностей дает следующая технология "вложенных замен".
Технология 1 [2]. Пусть Дп(у) — плотность случайной величины п, имеющей элементарное распределение в интервале (с, в), т. е. из соотношения типа, (2)
По
! Д(у) ву = ао
с
для, соответствующего выборочного значения п0 случайной величины, п можно получить формулу типа, (3): по = Фп(а0), где фп(т) — простая композиция элементарных функций. Рассмотрим взаимно-однозначное преобразование, задаваемое монотонно возрастающей дифференцируемой функцией у(х), переводящей интервал (а,Ь) в интервал (с, в); в частноети у (а) = с, у(Ь) = в. Полагаем также, что функцию у(х) и обратную к ней функцию у-1 (у) можно представить в виде простой композиции элементарных функций. Пусть случайная, величина, С имеет плотность распределения,
f (и) = Д(у(и)) у'(и), и е (а, Ь). (4)
При сделанных предположениях можно утверждать, что f (и) является, плотностью элементарного распределения, т. е. уравнение (2) разрешимо относительно Со в элементарных функциях и справедлива, формула Со = у-1(фп(а0)).
Действительно, записывая уравнение (2) для плотности (4), имеем
У Д(у(и)) у'(и) ви = ао, или У Д(у) ву = ао,
а <р(а)
или у(Со) = Фп (ао), ИЛИ Со = У-1(Фп (ао)). (5)
Название технология вложенных замен для технологии 1 связано с тем, что полученную плотность (4) можно взять в качестве исходной плотности fп (у) и осуществить еще одно взаимно-однозначное преобразование типа у (и). С помощью таких вложенных замен можно получать неограниченное количество новых плотностей элементарных распределений.
Пример 1. В методах численного статистического моделирования широкое применение находит формула
1п ао , ,
•По =--д-, (6)
соответствующая экспоненциальному распределению с плотностью
Д(у) = Ае-Хь, у> 0, А > 0. (7)
На основе формулы (6) формируются и реализуются пуассоновские потоки, используемые в теории массового обслуживания, в простейших моделях теории переноса излучения, при моделировании случайных полей и т.д. (см., например, [2]),
Рассмотрим также случайную величину имеющую плотность распределения
f (и) = ехри х ехр(— ехри), —то < и < +то. (8)
Это плотность экстремального (точнее, минимального) распределения (см., например, [4]), описывающая одно из трех возможных асимптотических распределений линейных комбинаций вида ап шт{п(1),... ,п(п)} + Ьп при ап = 0, п ^ то; здесь ап,Ьп — числовые последовательности, а {п(г)} — независимые одинаково распределенные случайные величины. Применение распределения (8) связано с множественными сравнениями в сложных процедурах принятия решений (таких, в частности, как ранжирование средних).
Функция (8) может быть получена из плотности (7) согласно технологии 1 с преобразованием <^(я) = ехр я, переводящим интервал (а,Ь) = (—то, +то) в интервал (с, д) = (0, +то). Согласно соотношению (5) моделирующая формула для распределения (8) имеет вид £0 = 1п(— 1па0).
Приведенный пример показывает, что применение технологии 1 позволяет получать плотности распределения и соответствующие моделирующие формулы для различных разделов теории вероятностей и связанных с ними приложений.
2. Технология взвешенного параметра
Во многих приложениях численного статистического моделирования (чаще всего при реализации траекторий цепи Маркова, а также при использовании метода двойной рандомизации, см., например, [2]) требуется конструировать "моделируемые" плотности распределения f (и, V) двумерных случайных векторов (£,п) с зависимыми компонентами, Справедливы два представления (см., например, [2])
= ДЫ = J ¡(и,ь)(1ь, Ь(ь\и) = ^^у (9)
= /^ЩиИ, ^(у) = J f(u,v)du, = (10)
Представлению (9) соответствует следующий алгоритм численного моделирования вектора (£, п): сначала реализуется выборочное значение £0 согласно плотности Д(и), затем моделируется выборочное значение п0 согласно плотности f(^о^)//^(£0). Аналогично для представления (10) — сначала реализуется выборочное значение п0 согласно плотности (V), далее моделируется выборочное значение £0 согласно плотности f (и^^/^(п0)- Сформулированные алгоритмы могут быть далеко не равнозначными с точки зрения их эффективной реализации на ЭВМ,
Пример 2. Пусть требуется построить эффективный алгоритм моделирования двумерного случайного вектора (£, п) с плотностью распределения
/(и, г;) = ^ ье~иь, и> 0, 0 < V < 2.
Рассмотрим представление (10)
/»= I ¿и =1-, 0 < г; < 2; ¡¿и\у) = = уе~™, и > 0.
о п
Плотность /пявляется плотностью равномерного распределения на интервале (0, 2)-, соответствующая моделирующая формула щ = 2а1. Функция /(и1п0) является плотностью экспоненциального распределения с параметром Л = щ (см. соотношение (7)), и, следовательно, С = — (1п а2)/п0 (см. формулу (6)). Теперь рассмотрим представление (9). Интегрируя по частям, имеем
^ 2 1 Л 1 — (2и + !)е-2и
/?(«) = / 2 ^ =- 2и2 -> и >
0
Полученная функция, очевидно, не является элементарной плотностью распределения, и поэтому для этого примера представление (9) является заведомо худшим (с точки зрения реализации на ЭВМ) по сравнению с представлением (10).
Примеры плотностей, для которых по крайней мере одно из разложений (9) или (10) соответствует эффективному алгоритму численной реализации выборочных значений (£0,П0); дает следующая технология "взвешенного параметра".
Технология 2. Рассмотрим элементарную плотность распределения, /(и; Л), и Е (а, Ь), зависящую от па,ра,м,етра, Л, допустимые значения которого принадлежат интервалу (С, Б). Элементарность распределения, означает существование простой (элементарной) формулы £0 = (а\; Л) для, получения выборочного значения случайной величины, Рассмотрим также еще одну элементарную плотность /п (у) случайной величины, п, принимающей значения в интервале (с, в) С (С, Б); при этом, имеется соответствующая элементарная моделирующая формула п0 = Фп(а2)- Теперь поставим, задачу построения эффективного алгоритма реализации выборочных значений (С0,щ) двумерного случайного вектора (£,ц), принимающего значения в прямоугольнике С = {(и, у) : а < и < Ь; с < у < в} и имеющего плотность распределения,
/(и, у) = /п(у) х /(и; у), (и, у) Е С. (11)
Это результат формального умножения плотностей /п(у) и / (и; у) (здесь происхо-
уЛ ности (11) имеем, /(и1у) = /%(и; у). Для, этого представления получаем эффективный алгоритм
П0 = Фп ^^ & = ф (а2; П0). (12)
Для, представления (9) плотности (11) эффективных формул типа, (12) построить, как правило, не удается.
В частности, пример 2 описывает ситуацию, в которой применена технология 2. В качестве исходной плотности с параметром использована функция (7) (здесь (С, Б) = (0, и на подмножестве (с, в) = (0, 2) С (С; Б) выбрана плотность равномерного
распределения.
3. Технология формирования смеси
Частным случаем применения метода двойной рандомизации (или метода интегральной суперпозиции), в котором при моделировании случайной величины £ вводятся вепомо-
п
f (и) = / /п^ШиИ ¿V, (13)
причем соответствующее разложение (10) дает эффективный алгоритм реализации пары (£0, п0) (см., например, [2]), является случай, когда вспомогательная величина п берется дискретной и целочисленной с распределением Р(п = г) = рг, г = 1, 2,... Здесь плотность (13) имеет вид
f (и) = X] Л(и) = Л (и|п = г) (14)
г
а моделирующий алгоритм (метод дискретной суперпозиции — см., например, [2]) состоит в выборе номера п0 = т согласно стандартному методу реализации дискретной случайной величины или какой-либо его модификации (см. [2] и раздел 6 данной статьи) с последующим моделированием значения £0 согласно плотности ^(и). Достаточно большой спектр примеров эффективного применения метода дискретной суперпозиции можно получить для небольшого количества М номеров г, г = 1,..., М (в частности, для М = 2). Сформулируем соответствующую технологию "формирования смеси". Технология 3. Возьмем две плотности, элементарных распределений Д(и) м ^(и),
(а, Ь)
ними коэффициентами
f (и) = Р1 Л(и)+ Р2 f2(u), и € (а, Ь), Р1 > 0, Р2 > 0, Р1 + Р2 = 1 (15)
не является, плотностью элементарного распределения. Такие плотности, Д(и) и f2(u) можно получить, в частности,, с помощью разнородных замен в технологии, 1. Длд выборочных значений £0г), реализуемых согласно плотностям ^(и), выписываются моделирующие формулы £0г) = Ф(а0), г = 1, 2. Длд плотности (15) можно построить экономичный алгоритм дискретной, суперпозиции: если а1 < р1; то значение п0 вспо-
п
н,и,е £0 случайной вели,чины, £ реализуется по формуле £0 = ф1(а2^ иначе £0 = ф2(а2).
Пример 3 [2]. Пусть требуется построить алгоритм моделирования случайной величины £, имеющей плотность распределения
¡■(и) = ^(1+и2), -1<и<1. (16)
Соотношение (16) представляет так называемый закон, Релея молекулярного рассеяния, фотонов в атмосфере, используемый в теории переноса излучения. Функция (16) не
являстся плотпос™ Элемептарп„г„ распределен, так как уравнение („) А, = а°
-1
сводится к соотношению ^ + 3£0 — 8а0 — 4 = 0, не позволяющему получить элементарную формулу моделирования случайной величины Плотность (16) предетавима в виде смеси (15):
^ N 3 1 1 3 2 ^=4Х2 4Х2И' -1<и<1,
т.е. pl = 3/4, fl(u) = 1/2, p2 = 1/4; f2(u) = 3u2/2. Функция fl(u) является плотностью равномерного распределения на интервале (—1,а функция f2(u) — элементарной (степенной). Алгоритм дискретной суперпозиции здесь выглядит следующим образом: если ai < 3/4, то = 2a2 — 1, иначе = \/2а2 — 1,
Обобщение технологии 3 может быть связано с увеличением числа слагаемых M в сумме (15) (вплоть до рассмотрения функциональных рядов), а также с переходом к моделированию многомерных случайных величин (случайных векторов)
4. Технология "порчи" моделируемой плотности
К дальнейшему расширению спектра вероятностных распределений, допускающих эффективное численное моделирование выборочных значений, ведет применение специальных вариантов мажорантного метода исключения, суть которого состоит в следующем [2], Пусть требуется численно получать выборочные значения ^ случайного вектора (случайной величины) распределенного в облаети U £ Rd согласно плотности f (u), которая пропорциональна заданной неотрицательной функции g(u), т.е.
/(u) = ^, G = jg(u)du. (17)
и
Предполагается, что ни один из известных стандартных и специальных методов не дает эффективного алгоритма реализации значений j Рассматривается мажоранта g(1)(u) функции g(u) такая, что g(u) < g(1)(u) при u £ U. Первое требование к мажоранте g(1)(u) таково, что для плотности
= G^ = Jg^(u)du, (18)
и
имеется эффективный алгоритм (формула) вида = ф(1) (аа1) для реализации выборочного значения случайного вектора (здееь а1 — соответствующий набор стандартных случайных чисел).
Алгоритм мажорантного метода исключения состоит в том, что реализуется выборочное значение согласно плотности (18), а также значение no = a2g(1)(^01)). Несложно показать (см., например, [2]), что пара (£0^,^) равномерно распределена в "подграфике" G(1) = {u £ U, 0 <v< g(1)(u)} функции g(1)(u), Если
no <д(й1]), (19)
то реализованная точка попадает в "подграфик" G = {u £ U, 0 < v < g(u)} функции g(u). Так как в этом случае пара (£0^, П0) равномерно распределена в области G,
то в качестве искомого выборочного значения £0 вектор а £ принимаем £0 = ^с^- Если неравенство (19) не выполнено, снова разыгрываем пару (£(1),П0); проверяем неравенство (19) и т.д. Несложно показать (см., например, [2]), что реализованная таким образом точка £0 распределена согласно плотности (17),
Среднее время реализации выборочного значения £0 пропорционально математическому ожиданию количества реализуемых пар (£01), По); которое равно в = (51/(5 (см., например, [2]), Близость величины в к единице характеризует эффективность метода исключения.
Примеры эффективного применения метода исключения можно построить с помощью следующей технологии "порчи" моделируемой плотности.
Технология 4. Конструируем сначала плотность /(1)(и) (и € и С Д^) вектора £(1), для которого существует эффективный алгоритм (формула) численной реализации £01) = (этот алгоритм используется затем в первом пункте алгоритма метода исключения). Для, построения, функции /(1)(и) можно использовать весь арсенал конструирования моделируемых плотностей, (в частности,, технологии, 1-Ъ). Далее преобразуем плотность /(1)(и) таким образом,, чтобы она трансформировалась в функцию д(и), пропорциональную "немоделируемой" плотности, /(и) (по сути мы "портим" моделируемую плотность /(1)(и)у). Одним из простейших преобразований является умножение плотности / (1)(и) на, мало меняющуюся ф ункцию У (и):
д(и) = / (1)(и) У (и), и € и; где 0 < А < У (и) < В (20)
и (В — А) — близкая, к нулю положительная величина. В качестве мажоранты тогда, можно взять д(1)(и) = В/(1)(и). Плотность, пропорциональная этой, функции, очевидно, равна, /(1)(и). Интегрируя неотрицательные функции д(1)(и) м д(и) по области и с учетом соотношения А/(1)(и) = Ад(1)(и)/В < д(и), получаем А(5(1)/В < (5, и тогда, в < В/А, т. е. при А « В величина в невелика (близка к единице) и соответствующий алгоритм метода исключения может считаться эффективным, (экономичным).
Для удобства выкладок в равенстве (20) вместо плотности /(1)(и) можно рассматривать пропорциональную ей функцию д(1)(и) (опуская, к примеру, нормирующую константу) ,
Пример 4. Пусть требуется построить алгоритм моделирования случайной величины имеющей плотность распределения /(п), пропорциональную функции
, . / агсвтп\ , #(«)=( 2Н----1 и, 0<и<1.
/(п)
что д(п) = У(п) х д(1)(п), где д(1)(п) = п3 и У(п) = 2 + (агсвт п)/(5п), причем в силу монотонности функции агсвтп на интервале (0,1) выполнено неравенство 2 < У(п) < 2.1. Тогда д(п) < д(1)(п) = 2.1 п3, Плотность, пропорциональная мажоранте д(1)(п), равна /^(и) = 4и3, 0 < и < 1; соответствующая моделирующая формула ¿¡д1"1 = а,о [2]. Алгоритм метода исключения содержит следующие пункты,
1. Реализуются выборочное значение ^ по формуле ^ = а также величина
По = д(1)(£01)) = 2.1 а2 (С01))3- Точка (С01),П0) равномерно распределена в "подграфике" мажоранты д(1)(п).
2. Проверяется неравенство п0 < д^о^) или
10.5 7ГО!2 < 107Г + агсвт (21)
Если это неравенство выполнено, то точка (Со^, По) принадлежит "подграфику" функции д(п) и является, равномерно распределенной в данном множестве. Тогда, в качестве выборочного значения, С0 случайной вели чины С берется Со = Со^ • Если же неравенство (21) не выполнено, то повторяется, п. 1, и т. д.
Трудоемкость в (среднее число попыток розыгрыша пар (Со1),П0) Д° выполнения неравенства (21)) оценивается сверху величиной в < 2.1/2 = 1.05,
5. Использование приближения плотности и метода суперпозиции
Описанные в предыдущих разделах технологии связаны с применением "теоретически точных" численных методов моделирования выборочных значений случайных величин и векторов. Практическая реализация этих методов на ЭВМ может быть сопряжена с погрешностями (как правило, незначительными), связанными с применением датчиков псевдослучайных чисел и с машинными ошибками [2].
Как указывалось выше, в прикладных задачах часто приходится иметь дело с ситуацией, когда моделируемое распределение случайного вектора £ в области О с Яа задается "извне" (здесь прежде всего следует упомянуть различные вариации метода выборки по важности — см., например, [2]), При этом можно, в том числе, выбрать распределение, близкое к требуемому, из "банка" плотностей, созданного согласно технологиям 1-4,
Однако в ряде приложений (в том числе в представленном в разделе 8 диекретно-ето-хаетичееком алгоритме построения адаптивных сеток) более универсальным видится следующий подход. Вместо вектора £ будем моделировать вектор £, плотность распределения которого f(x) пропорциональна некоторому (как правило, полиномиальному или кусочно-полиномиальному) приближению функции f (х) вида
м
/(х) = СЬмf (х) = СУш({)хг(х), С = 1/с, с = Ьмf (у) йу. (22)
г=1
С
Аппроксимация Ьм f (х) из формулы (22) строится па сетке X(м) = {х1,..., хм}, Базисные полиномиальные функции 2(м) = {х1(х),... ,хм(х)} и коэффициенты Ш(м) = {и^),...,им (£)} определенным образом связаны с узлами сетки X(м \ В частности, коэффициенты Ш(м) являются, как правило, комбинациями значений { = ^(х1),..., f (хм)); чаще всего
иг({) = f (хг), г = 1,... ,М. (23)
С позиций теории приближения функций (см., например, [5]) функция Ьм f (х) должна обладать хорошими аппроксимационными свойствами, т.е. величина рв^,Ьмf) должна быть мала; здесь рв(д1,д2) = ||д1 —д2 Нв _ расстояние в соответствующем нормированном функциональном пространстве В (О). Достаточно хорошо развита Ь1-теория приближения плотностей [6], в которой рассматриваются также вопросы приближения
вида (/, Л) линейных функционалов (/, Л) = J / (х)Л(х) ¿х = ЕЛ(£), Однако в целом
о
ряде приложений численного статистического моделирования требуется рассматривать более сильные (по сравнению с р^) метрики пространств СР(С) и ЖР(С) (см., например, [2, 7]).
В нашем случае приближение Ьм / (х) должно обладать свойством моделируемоети [2, 7], что означает наличие эффективного алгоритма численного моделирования вектора £ согласно плотности (22). Поскольку функция / (х) представляет собой линейную комбинацию заданных базисных функций, то это наводит на мысль о применении метода дискретной суперпозиции (см., например, [2], а также раздел 3). Перепишем плотность (22) в виде
м ( )
/~(х) = ^ВДх); /г(х) = ^, Уг =
i=1 г
Выберем функциональный базис 2(м) таким образом, чтобы неотрицательные коэффициенты Ш(м) обеспечивши малость величины рв (/, Ьм/) и для случайных векторов £(г), распределенных согласно плотпостям fi(x) из (24), имелись эффективные алгоритмы численного статистического моделирования. В частности, должны быть выполнены соотношения
Хг(х) > 0 для х е Я4, (25)
) > 0, г =1,...,М. (26)
Тогда по аналогии с алгоритмом из раздела 3 можно выбрать номер т согласно вероятностям {Рг} го (24) и реализовать выборочное значение согласно плотности /т(х).
В работе [7] проведен подробный анализ известных функциональных базисов с точки зрения одновременного наличия у них свойств аппроксимации и моделируемоети. Оказалось, что далеко не все "классические" аппрокеимационные базисы являются моделируемыми. Так, функции базиса Лагранжа (см., например, [5])
м
Хг(х) = Д (х - х^-)/(хг - х^-), х е Я,
з=1,з=г
являются знакопеременными (т.е. требование (25) не выполнено). Аналогичные недостатки имеют тригонометрические базисы [8]. Отметим также, что в ряде прикладных задач кроме аппрокеимационных свойств и свойств моделируемоети требуется свойство устойчивости приближения (22) (см., например, [2, 5]).
В работе [7] показано, что наилучшими (с точки зрения сочетания свойств аппроксимации, моделируемоети и устойчивости) являются простейшие конечно-элементные приближения Стренга —Фикса [9, 10]: кусочно-постоянная и .мулы н. пшенная аппроксимации. Здесь порядок аппроксимации по шагу Л сетки X(м) относительно невелик (для функций из С2 (С) он второй), однако он обеспечивает коэффициенты вида (4), для которых условие (6) заведомо выполнено. Попытка увеличить порядок аппроксимации до четвертого предпринята в работе [11], в которой в качестве базисных функций из 2(м) рассмотрены кубические сплайны. Здесь обнаружились существенные трудности при нахождении коэффициентов Ш(м\ обеспечивающих оптимальный порядок аппроксимации и выполнение условия (26).
/хг(у) ¿у, Рг = С^)Уг. (24)
Так же, как в работе [7], можно отметить, что простейший вариант аппроксимации Стренга —Фикса — кусочно-постоянная аппроксимация — дает наиболее универсальный и относительно просто реализуемый вариант описанного выше алгоритма метода дискретной суперпозиции.
При реализации предлагаемой методологии приближение области С происходит следующим образом, В пространстве КЛ рассматривается прямоугольная система координат и в ней — равномерная сетка ^ с шагом К, состоящая из узлов = О^Н,... здесь ^ — целые числа. Выбираются элементарные кубы ("воксели" [3])
= (О^Н (О1 + 1)Н) х ... х К, № + 1)Н) , г = 1,... ,М, (27)
имеющие непустое пересечение с областью С. В каждом кубе (27) выбирается узел х. сетки X(м\ Если куб полностью принадлежит области С, то х _ центр этого куба:
X = (О? + 1/2)Н,..., ^ + 1/2)Н). (28)
Сложнее ситуация при неполной принадлежности множества области С. Здесь возможны варианты. Если шаг К сетки Z мал, то можно либо продолжить каким-либо образом функцию / (х) на все множество либо, наоборот, исключить соответствующий куб из множества (27); в обоих случаях область С заменяется на приближение С = и^А, а сетка X(м" формируется из узлов вида (28), Для таких подходов про-
С
свойство гладкости этой границы). Если подобное нарушение недопустимо, то следует определить способ выбора узла х внутри куба Далее полагаем
Ьм/(х) = /(х.) при X е 3., (29)
что определяет кусочно-постоянную аппроксимацию (интерполяцию) функции /(х). При этом функции х.(х) из представления (22) тождественно равны единице при х е и нулю иначе, В свою очередь, все плотности /.(х) из формулы (24) тождественно равны 1/Нл при х е и нулю иначе. Это означает, что на втором шаге метода дискретной суперпозиции выборочное значение £0 вектор а £ реализуется в соответствующем кубе согласно равномерной плотности распределения
£о = т + а1Н,... + ааЬ) . (30)
Если выбранный на первом шаге метода суперпозиции куб Цт является граничным, то при моделировании вектора £ можно использовать следующий простой факт (см., например, [2]),
Утверждение 1. Если ¿-мерная точка а равномерно распределена в области А1 С * ^«ш Л =1 Л, » оиа Рааи0М,Рио „.^еиа . пртаа0М-
ной подобласти, А С А1 объем,а А при условии, попадания в эту подобласть; при этом,
Р(а е А) = А/Ах.
Для граничного Цт реализуется выборочное значение вектора £0 по формуле (30), и в случае выполнения условия £0 е С оно принимается, иначе алгоритм метода суперпозиции реализуется заново. При программной реализации алгоритма метода суперпозиции можно не различать граничные и "внутренние" кубы 3т, а проводить проверку условия £0 е С для каждого выборочного значения, реализованного по формуле (30) на втором шаге алгоритма.
6. Применение квантильного метода
Для того чтобы приближение (22), (29) обладало хорошими аппрокеимационными свойствами, следует выбирать Л достаточно малым. При этом параметр М может быть весьма большим. Это ведет к определенным проблемам, связанным с хранением значений (29) в оперативной памяти ЭВМ. Кроме того, возникают трудности реализации номера т
Напомним (см., например, [2]), что стандартным алгоритмом численного моделирования выборочного значения т целочисленной случайной величины п с распределением Р(п = г) = Рг (в нашем случае Рг = С/(хг)Л4) является следующий.
Алгоритм 1. Реализуем значение В := а0 и полагаем г := 1. Производим переприсваивание
В := В - Рг. (31)
Если новое значение В не положительно, то полагаем, т = г, в противном случае производим переприсваивания г := г + 1 и (31) и вновь производим проверку В на положительность и т. д.
м
Средняя трудоемкость в алгоритма 1 равна 5 = а + ( ^ гРг I х Ь; здесь а — затраты
на моделирование стандартного случайного числа а0, Ь — затраты на одно сравнение ВМ
Трудоемкость в можно существенно уменьшить, если применить так называемый квантильный метод моделирования дискретных случайных величин (см., например, [2]), который состоит в следующем.
Зададим целое число К и разобьем интервал (0,1) на К равных частей [(? - 1)/К, ?/К), ? = 1,..., К, Далее построим массив целых чисел {X, }К=1 такой, что
X,- = ш1п{к : Я = Р1 + Р2 + ... + Р* > (? - 1)/К},
называемый массивом нижних квантилей. Этот массив задает номер к элемента массива {Яг; г = 1, 2,..., М}, с которого следует начинать поиск "вверх" (т.е. как и в алгоритме 1, вычитать величины , д = к, к + 1,..., из а0 с помощью переприсваивания (31) до получения первого отрицательного значения) при (? - 1)/К < а < ?/К.
Окончательно моделирование дискретной случайной величины выглядит следующим образом.
Алгоритм 2 (квантильный метод).
1. Реализуем выборочное значение а0 равномерно распределенной, в интервале (0,1) случайной величины, а.
2. Вычисляем, номер п полуинтервала, [(п - 1)/К, п/К), в который попадает а0, по формуле моделирования равномерного дискретного распределения (см., например, [1]^):
п =[Ка0] + 1, (32)
здесь [Т] — целая, часть числа Т.
3. Реализуем последовательный поиск "снизу вверх" начиная с ЯХп.
Алгоритм 2 требует использования дополнительной (по сравнению с алгоритмом 1) оперативной памяти ЭВМ для хранения массивов номеров {X,} и сумм }■ Тестовые вычисления показали, что следует выбирать число К квантилей [(? - 1)/К, ?/К) так,
чтобы выполнялось соотношение M/K « 2 (при этом трудоемкость s1 алгоритма 2 с ростом M практически не меняется).
Отметим, что можно построить эффективную реализацию алгоритма 2 и в случае, когда число M относительно мало и D0 > r =1 + 1/ mini=1,...,M Pj (здееь D0 — максимально допустимый размер массива в выбранном языке программирования). Тогда можно взять число квантилей как целую часть K = [r]; при этом в каждом квантиле [(j -1)/K,j/K) будет не более одного значения Rk и в пункте 3 алгоритма 2 потребуется не более одного вычитания.
Заметим также, что в целом ряде ситуаций целесообразно использовать вместо кван-тильного метода алгоритм Уолкера или даже бинарный поиск, однако алгоритм 2 является более универсальным и простым для реализации на ЭВМ [2].
7. Случай равномерного распределения
Представленные выше алгоритмы допускают эффективные модификации в случае, когда f (x) = 1/G при x б G и f (x) = 0 при x </ G (G — объем компактной области G), т. e, случайный вектор £ равномерно распределен в G. Здесь кусочно-постоянное приближение LM f (x), описанное в разделе 5, дает одинаковые вероятности Pj = hd/G1, i = 1,..., M (в этом случае G1 — область G, дополненная теми частями граничных кубов Qm, которые не входят в G), Поэтому при выборе номера m в алгоритме метода суперпозиции вместо алгоритма 1 (или его модификации — квантильного алгоритма 2) можно применить простую формулу (типа (32)) моделирования дискретного равномерного распределения (см., например, [2]): m = [Ma0] + 1.
Как уже отмечалось, при малом шаге h вспомогательной сетки Z могут возникнуть проблемы с сохранением информации об узлах (28) и значениях (29) в оперативной памяти ЭВМ, Для равномерного распределения можно обходить эту трудность, пользуясь, в частности, тем обстоятельством, что все значения (29) одинаковы. Кроме того, можно использовать утверждение 1, например, из элементарных кубов Qj из (27) сформировать прямоугольный параллелепипед
максимального объема П, полностью принадлежащий области С (будем называть такой параллелепипед вписанным в область С), и с вероятностью П/С разыгрывать факт принадлежности выборочного значения £0 этому параллелепипеду. Если £0 е П, то реализация его выборочного значения происходит по формулам типа (30):
иначе в области С \ П применяется алгоритм метода суперпозиции.
Можно (а чаще всего нужно) строить "вложенную" процедуру вписывания параллелепипедов максимального объема, для которой очередной параллелепипед строится в множествах типа С \ П, В пределе получаются конструкции, аналогичные описанным в разделе 5, с той разницей, что только подмножества вида (27) и (33) получаются разной формы и разного объема.
Заметим также, что параллелепипеды (33) малого объема можно использовать вместо кубов (27) в качестве элементов разбиения области (вокселей) с заменой формулы
П = (a(1), b(1)) х ... х (a(d),
(33)
£о = (a(1) + (b(1) - a(1))a1,..., a(d) + (b(d) - a(d))ad)
(34)
(30) на соотношение (34) в соответствующем методе суперпозиции; это предусмотрено в программном модуле AITricks GeomRandom (см, раздел 9 и сайт [1]),
Кроме того можно пытаться вписывать в область G (и использовать в качестве вокселей) другие стандартные фигуры, для которых имеются формулы реализации равномерно распределенной точки, В частности, в случае использования триангуляции двумерной области (или поверхности) требуется применение эффективного алгоритма моделирования случайного вектора £ (двумерного или трехмерного), равномерно распределенного в треугольнике с вершинами r2, r3 (см., например, [12]), Имеем £ = ri + (r2 — ri)n(1) + (r3 — г2)п(2), где вектор (П(1),П(2)) равномерно распределен в "стандартном" треугольнике {(u,v) :0 < u < 1,0 <v< u}, Для реализации выборочных значений (По^По^) наиболее эффективным является следующий алгоритм:
(1) (2) (1) ^ (2) (1) 1 (1) (2) 1 " (2) / а „
По := аь по := «2/ если По < По , m0 По := 1 — ПО , По := 1 — По (здесь := — знак переприсваивания). Этот алгоритм является более экономичным, чем формулы
По1"1 := л/аГ, По2"* := рекомендованные в [12],
Известны также алгоритмы моделирования случайных точек, равномерно распределенных в других "примитивных" фигурах (круг, сфера, усеченный конус, тор и др.). Так, формулы
¿¡о1"1 = Ry/cn sin 27го!2, ¿2) = iV^cos27ra2 (35)
дают реализацию вектора £ = (£(1),£(2)); равномерно распределенного в круге радиуса R, а формулы
Со1) = sin ezcos фс, ^2) = rzsin ezsin фс, ¿3) = rzcos ez, (36)
где cos вz = 1 — 2а1, r^ = R(a2)1/3, ф^ = 2па3, дают реализацию вектора £ =
(£(1),£(2),£(3)); равномерно распределенного в шаре радиуса R (см., например, [2]), Следует, однако, отметить, что для фигур Н, отличных от параллелепипеда (33), гораздо труднее строить приближения дополнений G \ Ни реализовывать на них алгоритм метода суперпозиции,
£
G
G
G С П) и реализуется случайная точка (34), Если
£о е G, (37)
то соответствующее выборочное значение принимается, иначе снова применяется формула (34) вплоть до того момента, когда будет выполнено соотношение (37), Согласно утверждению 1, реализованное таким образом выборочное значение соответствует рав-
G
Метод исключения будет тем эффективнее, чем раньше (в среднем) будет выполнено соотношение (37), Известно (см., например, [2]), что среднее число обращений к формуле (34) до выполнения (37) равно отношению объемов s = П/G, Таким образом,
П
добиваться того, чтобы величина s была близка к единице:
s = П/G > 1.
(38)
В большом числе случаев добиться выполнения соотношения (38) не удается, и поэтому описанные "прямые" методы, связанные с "вокеелизацией" и с применением алгоритма метода суперпозиции и его модификаций, оказываются более экономичными (здесь, однако, метод исключения может пригодиться при моделировании точек вблизи границы области — см, раздел 5), Аналогичное замечание можно сделать и для методов исключения, связанных с "погружением" области G в другие стандартные области (например, в окружность или круг с применением формул (35), (36)),
8. Итерационный дискретно-стохастический метод построения адаптивных сеток
Описанные в предыдущих разделах технологии конструирования вероятностных плотностей, допускающих точную или приближенную численную реализацию выборочных значений, находят широкое применение в прикладных задачах,
В качестве примера такого применения рассмотрим проблему построения адаптивных сеток на сложных многомерных компактных областях. Использование адаптивных сеток в задачах численного моделирования позволяет повысить точность приближенного решения задачи без существенного увеличения числа узлов (см., например, [3, 13-17]), "Сгущения" узлов (т. е, их распределение согласно заданной плотности f (x)) часто требуются в подобластях, где либо само решение, либо его градиент принимают большие значения. Кроме того, важные приложения адаптивных сеток (такие, как обработка изображений, визуализация данных и т, п.) связаны с учетом сложной геометрии объемных областей и их поверхностей.
Среди всех видов методов построения сеток можно выделить важный класс, в котором адаптивные сетки получаются в результате отображения "вычислительной" области Q С RQ на "физическую" область X, Как правило, k = d (для "объемных" сеток) или k = d — 1 (для "поверхностных" сеток), К указанному классу относятся метод эк-вираепределения [16], метод Томпсона [15], эллиптические методы [17] и др. Все они, а также алгебраические методы и метод конформных отображений [13, 14] требуют решения сложных систем нелинейных дифференциальных уравнений с частными производными (при этом приходится накладывать довольно жесткие условия гладкости на граничные и начальные условия, а также на способы задания области). Еще один недостаток перечисленных методов связан с трудностями автоматизации и распараллеливания соответствующих численных схем.
Преодолеть указанные недостатки позволяет подход, предложенный в работе [3] и подразумевающий использование стохастических нейросетевых моделей самоорганизации, таких как самоорганизующиеся карты Т. Кохонена (Self-Organizing Maps — SOM), растущий нейронный газ (Growing Neural Gas — GNG), растущие клеточные структуры (Growing Cell Structures — GCS) [18],
Опишем соответствующую постановку задачи, В "физической" области X (или на ее поверхности Gx) требуется построить сетку X(M) = |xb..., xM}, распределение узлов которой соответствует заданной плотности f (x), x G Rdx. Структура этой сетки (порядок и структура расположения узлов) задается "картой" Q(M) = (qi,..., qM} и системой так называемых нейронов E(M) = (e1,... ,eM} (где ei = (q^ x^)), определяющих соответствие между сетками X(M) и Q(M\
Приближение нейронов осуществляется с помощью процедуры самообучения, представляющей собой итерационный процесс, основанный на последовательном форми-
ровании обучающего множества S(T) = {£0(1),..., £0(T)} в виде выборки из вероятностного распределения случайного вектора имеющего плотноеть f (x); здееь T — количество итераций, £0(t) G X (или £0(t) G GX; t = 1,...,T), Кроме того, с помощью специальных коэффициентов обучения 94j (t, q^) G [0,1) на каждом шаге итерации устанавливаются латеральные связи между нейронами ei и е.,-. В результате получается последовательность сеток X(M^(t) = (xi(t),...,xM(t)}; при этом требуется выполнение соотношения
X(M) « X(M) = lim X(M)(t) « X(M)(T). (39)
t—y^o
Знаки приближенных равенств в (39) означают не только малость расстояний р(хг(сх>),хг), р(х*(го), Xi(T)) и p(xi(T), Xi), где
р(х, у) = г(1) - у(1))2 + . . . + (а**) - у х=(Ж(1),...,^)), у = (?/(1),---,Л
но и воспроизведение (с хорошей точностью) требуемых свойств сетки X(M) (например, свойств прямоугольноети, отсутствия граничного эффекта и др.) при реализации следующего итерационного процесса.
Алгоритм 3 [3, 18]. 1. Устанавливаются, начальные положения, узлов сетки X(M )(0) = {xi(0),..., xm (0)}.
2. На каждой итерации с номером t = 1,..., T выполняются следующие действия:
а) выбирается очередной, элемент £0(t) выборки S(T);
б) вычисляются расстояния p(£0(t), x^t — 1)) от точки £0(t) до всех узлов x^t — 1) и выбирается, ближайший к £0(t) узел, xm(t — 1) в соответствии, с условием
т = arg rnn р ^С^ xi(t — 1));
i=1,...,M
такой узел, xm(t — 1) называют победителем;
ej проводится корректировка положений всех узлов в соответствии, с формулой
xi(t) = xi(t — 1) + 0qm (t, qi) (^(t) — xi(t — 1)) (40)
для, всех i = 1,..., M.
На каждой итерации алгоритма 3 узлы сетки сдвигаются к случайной точке £0(t). Поэтому в местах высокой концентрации элементов выборки постепенно скапливается все больше узлов, за счет чего достигаются сгущения сетки. В работе [3] показано, что при T — оо алгоритм 3 ведет к выполнению аналога принципа эквираепределения, т. е. к получению требуемой плотности сетки, определяемой функцией f (x). Весьма нетривиальной (с точки зрения обоснованного применения алгоритма 3) является проблема выбора коэффициента обучения 04m (t, qi) из формулы (40).
9. Программный модуль AITricks GeomRandom
При реализации алгоритма 3 определенную трудность представляет собой пункт 2а, в котором требуется численно моделировать многомерную (чаще всего двумерную или трехмерную) случайную точку (случайный вектор) £0, распределенную согласно заданной плотности f (x). Здесь наиболее перспективным видится применение технологий приближения плотностей, описанных в разделах 5-7. Соответствующие рекомендации использованы при создании программного модуля AITncks GeomRandom [1].
GeorriRandom представляет собой кросс-платформенную библиотеку С -классик и шаблонов для моделирования случайных точек внутри и на границе одно-, двух- и трехмерных фигур. Библиотека поддерживает реализацию случайных точек как на фигурах, полученных с помощью булевых операций над примитивными фигурами (куб, параллелепипед, круг, сфера, усеченный конус, тор и др.), так и на сложных объектах, заданных поверхностной триангуляцией, У пользователя имеется возможность загрузить файл с CAD .моделью в одном из стандартных форматов, тем самым обеспечивая совместимость с такими распространенными (MD-пакетами, как Autodesk AutoCAD/Inventor и др. Реализована экономичная версия вокселизации (разбиения области на простые подобласти) с целью применения соответствующего метода суперпозиции, В случае большого количества подобластей разбиения предусмотрено использование квантиль-ного метода (см, раздел 6),
Особо отметим, что описанный программный модуль сначала разрабатывался для преимущественного применения при реализации алгоритма 3 и его модификаций. Следует, однако, подчеркнуть, что этот модуль может быть эффективно использован при решении других задач, связанных с моделированием случайных точек. Здесь в первую очередь следует упомянуть задачи, решаемые методами численного статистического моделирования (см., например, [2]), Полученные с помощью библиотеки GeomRandom точки можно либо использовать внутри соответствующей программы, либо экспортировать в один из удобных форматов, в том числе в стандартный формат PTS.
Заключение
В работе сформулированы технологии конструирования вероятностных плотностей распределения, допускающих эффективную численную реализацию выборочных значений: технология вложенных замен (для использования метода обратной функции распределения) , технология взвешенного параметра (для метода моделирования двумерного случайного вектора с зависимыми компонентами), технология формирования смеси (для метода дискретной суперпозиции), технология "порчи" моделируемого распределения (для мажорантного метода исключения).
Показано, что в ряде практически важных случаев весьма перспективным является применение приближений (в частности, полиномиальных и кусочно-полиномиальных, а конкретнее — кусочно-постоянных) вероятностных плотностей, В качестве примера рассмотрен итерационный дискретно-стохастический метод построения адаптивных сеток, Представлен программный модуль AITricks GeomRandom, позволяющий получать выборочные значения случайных векторов, распределенных согласно заданной плотности в геометрически сложных областях.
Список литературы
fl] HTTP://aitricks.com/products/geomrandom/
[2] Михайлов Г.А., Войтишек A.B. Численное статистическое моделирование. Методы Монте-Карло. М.: Изд. центр "Академия", 2006.
[3] Нечаева О.И. Нейросетевые модели, алгоритмы и комплекс программ для построения адаптивных сеток. Дисс. ... канд. физ.-мат. наук. Новосибирск, НГУ, 2007.
[4] Дейвид Г. Порядковые статистики. М.: Наука, 1979.
[5] Бахвалов Н.С. Численные методы. М.: Наука, 1975.
[6] Деврой Л., Дерфи Л. Непараметрическое оценивание плотности (Li). М.: Мир, 1988.
[71 voytishek A.v., Kablukova E.G. Using the approximation functional bases in Monte Carlo methods // Rus. J. of Numerical Analysis and Math. Modelling. 2003. Vol. 18, No. 6. P. 521-542.
[8] Березин И.С., Жидков Н.П. Методы вычислений. М.: Физматгиз, 1962.
[9] Стренг Г., Фикс Дж. Теория метода конечных элементов. М.: Мир, 1977.
[10] Марчук Г.И., Агошков В.И. Введение в проекционно-сеточные методы. М.: Наука, 1981.
[11] Милосердое В.В. Дискретно-стохастические численные алгоритмы со сплайн-восполнениями. Дисс. ... канд. физ.-мат. наук. Новосибирск, НГУ, 2006.
[12] Махоткин О.А. Моделирование точек, равномерно распределенных в многоугольниках // Сиб. журн. вычисл. математики. 2002. Т. 5, № 4. С. 331-350.
[13] Годунов С.К., Прокопов Г.П. О расчетах конформных отображений и построении разностных сеток // Журн. вычисл. математики и матем. физики. 1967. Т. 7. С. 1031-1059.
[14] Gordon W.J., Thiel L.C. Transfinite mappins and their applications to grid generations // Numerical Grid Generation. Appl. Math, and Comput. 1982. Vol. 2/3. P. 171-192.
[15] Thompson J.F., Warsi Z.U.A., Mastin C.W. Numerical Grid Generation: Foundations and Applications. Amsterdam (Netherlands): North-Holland, 1985.
[16] Хакимзянов Г.С., Шокин Ю.И., Барахнин В.Б., Шокина Н.Ю. Численное моделирование течений жидкости с поверхностными волнами. Новосибирск: Изд-во СО РАН, 2001.
[17] Лисейкин В.Д., Лебедев А.е., Китаева И.А. Универсальный аналитический метод построения разностных сеток. Новосибирск: НГУ, 2004.
[18] KOHONEN Т. Self-Organizing Maps. Springer-Verlag, 2001.
Поступила в редакцию 13 сентября 2010 г., с доработки — 20 апреля 2011 г.
