Содержание к диссертации
Введение
Глава 1. Дискретно-стохастические несмещенные оценки в алгоритмах численного интегрирования 13
1.1. Использование функциональных базисов в методах Монте-Карло 13
1.1.1. Моделирование случайных величин 13
1.1.2. Функциональные оценки 16
1.1.3. Выбор функционального базиса 17
1.1.4. Моделируемость аппроксимации Стренга-Фикса 17
1.1.5. Моделируемость аппроксимации Бернштейна 20
1.2. Дискретно-стохастические методы уменьшения дисперсии 21
1.2.1. Дискретно-стохастическая версия метода выборки по важности 21
1.2.2. Дискретно-стохастическая версия метода выделения главной части 25
1.2.3. Сложная многомерная симметризация 27
1.2.4. Дискретно-стохастическая версия метода выборки по группам 30
1.3. Двусторонний геометрический метод Монте-Карло 36
1.3.1. Геометрический метод И. М. Соболя 36
1.3.2. Модификация геометрического метода 36
1.3.3. Дискретно-стохастическая версия двустороннего геометрического метода 37
Глава 2. Дискретно-стохастические состоятельные и асимптотически несмещенные оценки в алгоритмах численного интегрирования .. 41
2.1. Дискретно-стохастическая версия метода взвешенной равномерной выборки.41
2.1.1. Лемма о состоятельных оценках 41
2.1.2. Взвешенная равномерная выборка 41
2.1.3. Использование аппроксимации Стренга-Фикса 42
2.1.4. Зависимость дисперсии от шага сетки 43
2.1.5. Построение доверительных границ и оптимизация оценки вп 43
2.1.6. Результаты тестовых численных экспериментов 44
2.2. Дискретно-стохастическая версия метода Монте-Карло с поправочным множителем 47
2.2.1. Оценка с поправочным множителем 47
2.2.2. Приближение оптимального множителя. Зависимость смещения и дисперсии от шага сетки 47
2.2.3. Построение доверительных границ и оптимизация оценки вп 48
2.2.4. Результаты тестовых численных экспериментов 49
2.3. Рандомизация метода последовательных приближений 53
2.3.1. Итерационный процесс с интегральным оператором 53
2.3.2. Приближения функционала 53
2.3.3. Тестовая задача 55
2.3.4. Использование специального функционала 57
2.3.5. Пример согласованного выбора параметров 59
Глава 3. Стохастическая тестовая система функций 61
3.1. Преобразования спектральных моделей случайных полей 61
3.1.1. Численные спектральные модели гауссовских случайных полей 61
3.1.2. Тестовая спектральная модель 62
3.1.3. Преобразования гауссовских моделей: использование функций многих переменных 63
3.1.4. Использование комбинаций со случайными величинами 66
3.1.5. Функциональная сходимость преобразованных моделей 66
3.1.6. Группировка слагаемых в моделируемой сумме 67
3.2. Тестовая система функций 69
3.2.1. Использование модельных траекторий случайных функций 69
3.2.2. Выполнение требований (0.1а)-(0.1д) 69
3.3. Средние оценки погрешностей простейших квадратурных формул 70
Заключение 76
Литература 77
- Моделирование случайных величин
- Дискретно-стохастическая версия метода выборки по группам
- Взвешенная равномерная выборка
- Численные спектральные модели гауссовских случайных полей
Введение к работе
0.1. Детерминированные и стохастические кубатурные формулы. С развитием вычислительной техники возрастает интерес к численным методам решения прикладных задач. Математическое описание исследуемого процесса сводится, как правило, к рассмотрению неизвестной функции многих переменных, для которой записывается система дифференциальных уравнений (см., например, [1]). Одним из способов приближенного решения этой системы является построение разностной схемы и сведение задачи к решению системы линейных уравнений для значений неизвестной функции в узлах сетки (см., например, [2-4]). Часто также оказывается целесообразным сведение (или постановка) задачи к интегральной форме, когда исследуемая функция представляет собой многократный интеграл, зависящий от параметра, или является решением интегрального уравнения (см., например, [1]). В последнем случае при возникновении интегрального уравнения Фредгольма второго рода соответствующий интегральный оператор предполагается сжимающим, и решение записывается в виде ряда Неймана, представляющего собой сумму параметрических интегралов бесконечно возрастающей кратности (см., например, [5]). Далее используются численные методы приближенного вычисления получаемых многократных интегралов на ЭВМ. При разработке этих методов решается
ЗАДАЧА 0.1. Построить алгоритм вычислегшя интеграла
1 = J д(х) dx = / я(х) dx = J <7(х) dx; (0.1)
здесь X - замыкание тех х Є Rl, для которых <у(х) ф 0.
Для интегралов / малых размерностей I с гладкими (в обычном или обобщенном смыслах) подынтегральными функциями g и относительно простыми областями интегрирования X развита теория квадратурных (для случая I = 1) и кубатурных (для I > 1) формул (см., например, [4, 6]). Кубатурная формула в общем случае имеет вид
/и5п = Х)^(хД (0-2)
j=l
где {xi,... ,хп} - заданные детерминированные (и, как правило, регулярные) узлы сетки в Rl, a {ci,...,c7l} - веса. Вычисление интеграла (0.1) по формуле (0.2) будем называть АЛГОРИТМОМ 0.1. В качестве X в данной работе будем использовать, как правило, простые компактные подмножества в R1 (чаще всего - /-мерный единичный куб Qi = {х = (хъ ..., xi) : 0 < Хі < 1; і = 1,..., I}).
Оптимальный выбор узлов и весов связан с минимизацией погрешности 5п = \I — Sn\ и основан (явно или неявно) на использовании аппроксимаций подынтегральной функции д. Главным достоинством кубатурных формул является возможность получения гарантированной и сравнительно быстрой сходимости 5п к нулю при п —У со для классов гладких подынтегральных функций д.
К недостаткам «классических» (детерминированных) кубатурных формул на классах подынтегральных функций следует отнести:
слабый учет специфики той или иной подынтегральной функции;
необходимость разработки специальных численных алгоритмов поиска оптимальных весов и (или) узлов;
чувствительность к росту размерности и к гладкости начальных данных (подынтегральной функции g и области X);
- трудности в построении показательных тестовых численных примеров и кон
троля точности и затрат при практических вычислениях.
Для существенно многомерных задач 0.1 (т.е. для случая />1и даже / —v со) достаточно эффективным оказывается стандартный метод Монте-Карло (см., например, [7-9]). Этот алгоритм основан на представлении
I = Jq(x)f(x)dx = -EC, C = q(0, g(x) =
Для выполнения (0.3) достаточно потребовать /(х) ф 0 при х Є X. Весовая функция / является вероятностной плотностью в Л', т.е. /(х) > 0 и //(x)dx = 1, а /-мерный случайный вектор распределен согласно плотности /.
АЛГОРИТМ 0.2 [7-9]. 1. Реализуем п значений случайного вектора .- 15... , п.
2. Вычисляем приближенное значение интеграла (0.1):
Если в (0.4) трактовать 1;..., п как набор независимых одинаково распределенных случайных векторов с плотностью распределения /, то случайные величины Сі = jCn = q(n) будут также независимыми одинаково распределенными с математическим ожиданием / (см. соотношение (0.3)) и дисперсией
DC = a2 = Dg(0 = j ?2(х) Дх) dx - /2 = J & dx ~ I2. (0.5)
Если величина (0.5) конечна, то в силу закона больших чисел (см., например, [10]) формула (0.4) верна для достаточно больших п.
Сразу заметим, что для /(х) = 1 (т. е. для равномерного распределения в X = Q{) формула (0.4) является частным случаем формулы (0.2) для с\ = ... — сп = \/п и Xj = j. Таким образом, разница между алгоритмами 0.1 и 0.2 с постоянными весами связана с выбором узлов - детерминированным или стохастическим. Отметим также, что и для «детерминированного» алгоритма 0.1 введение весовой функции типа / и связанного с ней набора узлов позволяет улучшать качество кубатурных формул (0.2).
Алгоритм 0.2 имеет следующие положительные свойства:
возможность уменьшения трудоемкости алгоритма за счет, удачного (согласованного с видом подынтегральной функции д) выбора весовой функции f (алгоритм выборки по важности) или преобразования исходного интеграла (выделение главной части) и весовой функции (лгетоды математического ожидания и расщепления, выборка по группам и т.п.) - см. далее подразд. 0.5;
веса имеют простой вид, а узлы реализуются согласно выбираемому вероятностному распределению с плотностью /;
относительно слабая чувствительность к росту размерности I и к гладкости начальных данных (подынтегральной функции g и области X);
возможность контроля затрат и точности вычислений.
Последнее свойство обусловлено тем, что затраты s алгоритма 0.2 можно подсчитать по простой формуле s = nt, где t - среднее время получения одного выборочного значения случайной величины С (это время, в свою очередь, связано со «сложностью» вычисления функции g и со средними затратами на реализацию одного выборочного
значения вектора ). Далее, в силу центральной предельной теоремы (см., например, [10]) для достаточно больших п имеет место соотношение
p(V-Cn|<7^) =1-е. (0.6)
При фиксированном уровне погрешности число случайных узлов п прямо пропорционально дисперсии а2 случайной величины , и вместо s в качестве величины, отражающей затраты стандартного метода Монте-Карло, можно ввести число
S = a2t, (0.7)
которое называется трудоемкостью алгоритма 0.2 [7-9]. Рассмотрение величины (0.7) вместо s является более удобным при оптимальном выборе весовой функции /, т. к. величина S в явном виде содержит вероятностную характеристику (дисперсию) случайной величины С, которая, в свою очередь, определяется через /. Неизвестное значение (0.5) можно приближенно вычислить по набору выборочных значений l7... ,Л [7-9]:
Формула (0.6) отражает также и главный недостаток алгоритма 0.2 - относительно низкую (порядка 1/у/п) скорость сходимости погрешности к нулю при возрастании числа узлов п. К примеру, для I — In g Є С2(Х) простейшая формула прямоугольников имеет порядок погрешности 1/п2 (см., например, [4]). Это обуславливает использование алгоритма 0.2 только для достаточно больших размерностей /.
0.2. Дискретно-стохастические методы численного интегрирования. Формулы Бахвалова и теория сложности. Эффективными могут оказаться и смешанные, комбинированные процедуры численного интегрирования, сочетающие в себе элементы алгоритмов 0.1 и 0.2. Как правило, эти схемы содержат «детерминированную» составляющую, связанную с регулярной дискретизацией области X, а также «стохастическую» составляющую, связанную с применением метода Монте-Карло. Поэтому, вслед за А.В.Войтишеком [11], мы будем называть такие алгоритмы дискретно-стохастическими.
По-видимому, одна из первых численных схем подобного рода была представлена в работе Н. С. Бахвалова [12] (см. также [4] и подраздел 1.2.4 данной работы). Исходной областью X являлся /-мерный единичный куб Qt, который разбивался на п равных кубов с вершинами в точках равномерной сетки. В каждом j-ом элементарном кубе выбирался узел х^ кубатурной формулы (0.2) случайным образом (согласно равномерному распределению). Представленный алгоритм можно считать как кубатурной формулой (0.2) со случайными узлами (такой термин для подобных конструкций имеется, например, в [8]), так и предельным случаем выборки по группам в методе Монте-Карло [7-9]. Н. С. Бахвалов показал в [12], что его алгоритм является оптимальным (по скорости сходимости к нулю погрешности 5п при п —У со) в пространстве непрерывно дифференцируемых подынтегральных функций СХ{Х). Для этого ему потребовалось получать оценки сверху и снизу для 5п.
Методология работы [12] явилась основой для интенсивного развития так называемой теории сложности (см. [13] и сопутствующие этой монографии работы). Эту теорию можно назвать «философией вычислительных алгоритмов». В ней изучается вопрос о том, каков максимальный порядок стремления к нулю погрешности 5ц для
данного класса вычислительных алгоритмов с количеством вычислительных операций п.
Проблема численного интегрирования (задача 0.1) является наиболее удобной и часто используемой иллюстрацией в теории сложности. В работах по теории сложности рассмотрены вопросы оптимальности кубатурных формул (0.2) для различных пространств В(Х) подынтегральных функций. При построении соответствующих нижних границ для 5п требуется строить конкретную численную схему типа (0.2). В некоторых (достаточно редких) случаях эти алгоритмы имеют вид, вполне пригодный для практических вычислений. Например, для В(Х) = С2(Х) в работе [12] Н. С. Бахвалов построил оптимальный алгоритм, в котором, в дополнение к алгоритму для д Є С1(Х), кроме узла х^ выбирается точка х7, симметричная х3 относительно центра j-ro элементарного куба (см. подраздел 1.2.4 настоящей работы).
С точки зрения теории сложности алгоритм 0.2 решения задачи 0.1 является не слишком интересным, так как, в силу соотношения (0.6), имеет неулучшаемый порядок сходимости 1/у/п. Поэтому в связи с задачей 0.1 в теории сложности распространено обращение к случаю использования так называемых квази-случайных узлов кубатур-ной формулы (0.2) [13]. Здесь теория кубатурных формул смыкается со специальным нетривиальным разделом теории чисел. Наиболее эффективный алгоритм построения квази-случайных узлов представлен в [7] (это ЛПГ-последовательность Соболя). Вопросы использования квази-случайных чисел в численном интегрировании в данной работе не рассматриваются.
В общей теории сложности (в том числе и для задачи 0.1) принимается ряд допущений, которые не всегда выполняются на практике. Например, приравниваются затраты на вычисление значения подынтегральной функции g в точке и одно арифметическое действие. В связи с этим открытым остается вопрос, будут ли теоретически оптимальные схемы иметь на практике минимальную трудоемкость. Забегая вперед, отметим, что в представляемой рабоїе указаны примеры ситуаций, когда оптимальные (с точки зрения теории сложности) кубатурные формулы не являются наилучшими на практике (см. раздел 1.2).
0.3. Цель и структура диссертации. Как уже было сказано выше, возможность эффективной практической реализации оптимальных кубатурных формул является большой редкостью (особенно для многомерных случаев). Поэтому вполне разумными видятся подходы, связанные с построением и исследованием эффективно реализуемых дискретно-стохастических численных алгоритмов, дающих относительно небольшие значения трудоемкости (0.7).
Целью данной работы является разработка и исследование эффективных дискретно-стохастических алгоритмов численного интегрирования, а также тестирование этих алгоритмов на основании построения стохастических систем функций и решения прикладных задач.
Диссертация состоит из введения, трех глав, заключения и списка литературы, содержащего 58 наименований.
0.4. О публикациях по теме диссертации. Основные результаты опубликованы в работах [14-37]. Следует отметить, что работы [14-28] написаны автором совместно с А. В. Войтишеком и Т. Е. Булгаковой. С согласия А. В. Войтишека и Т. Е. Булгаковой совместно полученные в этих работах результаты частично включены в представляемую диссертационную работу. Кроме того, имеется ряд трудов, написанных автором индивидуально и с соавторами-студентами (О. С Герасимовой, Н. В. Лощиной, A. И. Ефремовым, С.В.Бусыгиным); в последнем случае автор диссертации являлся научным руководителем соответствующих исследований [29-37]. Полученные в этих статьях и те-
зисах результаты также включены в диссертационную работу.
0.5. Краткое описание рассматриваемых в диссертации численных схем.
При построении дискретно-стохастических схем численного интегрирования мы будем следовать по пути «улучшения» стандартного алгоритма 0.2, т. е. уменьшения его трудоемкости (0.7). При этом используются как известные методы уменьшения дисперсии DC [7-9], так и способы уменьшения величины t.
В работе рассмотрена, в частности, специальная дискретно-стохастическая версия известного способа уменьшения DC, называемого выборкой по важности [7-9] (см. подраздел 1.2.1). Этот способ основан на утверждении о том, что минимум дисперсии достигается для плотности / вида Дх) — Ci\g(x)\ (здесь С\ - нормирующая константа) [7-9]. Для случая д > 0 (который мы и будем для простоты рассматривать в дальнейшем) при
f(x) = C2g(x) (0.9)
имеем DC = 0. Использование плотности вида (0.9) затруднено или невозможно в силу того, что С*2 = 1/1 [38]. В данной работе изучены алгоритмы метода Монте-Карло, в которых в соотношении (0.9) вместо функции д берется ее аппроксимация на X вида
А/
д(х) ^ LMg{x) = 53w;,-Xi(x), (0.10)
где {хі} _ заданный набор функций (функциональный базис). Как правило, эти функции связаны с дискретной сеткой {у;}. Коэффициенты {wi} выражаются через значения функции д и ее производных в узлах сетки. В большинстве вычислительных схем, рассмотренных в данной работе, коэффициенты {u>i} имеют простой вид Wi — g(yi),
при этом
д(х) и LMg(x) = ]Г д(Уі) Хі(х). (0.11)
Из принципа выборки по важности следует, что если погрешность приближения (0.10) (или (0.11)) мала, то величина DC близка к нулю. Однако в связи с алгоритмом 0.2 возникает следующая
ЗАДАЧА 0.2. Построить алгоритм численного моделирования случайного вектора (для 1 = 1- случайной величины) , имеющего плотность распределения
Дх) = Сз LMg{x), С3 = 1/с, с = I LMg(x) fix; g(x) > 0. (0.12)
Jr!
В случае, когда задача 0.2 эффективно разрешима, будем называть соответствующее приближение (0.10) (или (0.11)) моделируемым. В диссертации проводится исследование известных функциональных базисов с точки зрения сформулированного свойства моделируемости. Забегая вперед, отметим, что наилучшими свойствами моделируемо-сти обладает конечно-элементная аппроксимация Стренга-Фикса (см. [39, 40], а также подраздел 1.1.3). Здесь также уместно отметить, что выбор плотности / вида (0.10) или (0.11) означает сближение «детерминированного» алгоритма 0.1 и «стохастического» алгоритма 0.2, так как оптимальный выбор весов и узлов в формуле (0.2) связан (явно или неявно) с аппроксимацией подынтегральной функции д.
В качестве способа уменьшения дисперсии DC в диссертации рассмотрен также метод выделения главной части (см. [7-9], а также подраздел 1.2.2), состоящий в том, что интеграл (0.1) переписывается в виде
/ = /*0(х) - 0о(х)) dx + /о; /о = У 0о(х) dx, (0.13)
где д0 ~ д, а интеграл /0 берется аналитически. И снова при построении функции до из (0.13) используются приближения (0.10), (0.11).
Еще один способ уменьшения дисперсии - метод противополооїсной переменной или ліетод симметризации подынтегральной функции (в англоязычной литературе для этого приема используется термин «antithetic variates») [7—9j представлен в подразделе 1.2.3. Идею этого приема проще всего объяснить для случая вычисления однократного интеграла I0 = Ja g(x) dx по конечному интервалу а < х < Ь. Если взять f(x) = 1/(6—а) при х Є (а, 6), то согласно формуле (0.3) получаем /0 = Е^0\ где ((0) = {Ь — а)д(а + (6 — а)а) на- стандартное случайное число. Рассмотрим теперь симметризовапную функцию д^Цх) = (д(х) + д(а + Ь- х)) /2 и заметим, что
гь
Io= gw(x)dx^EC{l\ где (^ = (b - a)gw(a+(Ь - а)а). (0.14)
J а
Метод симметризации подынтегральной функции - это алгоритмі 0.2 с оценкой (^ (вместо С^)- Легко показать, что DC(1) < DC(0) [7-9]. Однако для расчета одного значения С(1) из (0.14) надо вычислить два значения функции д(х). Поэтому трудоемкость алгоритма 0.2 с оценкой (^ будет меньше трудоемкости метода противоположной переменной с оценкой () только тогда, когда величина D^1) по крайней мере вдвое меньше, чем D((0). Известно [7-9], что для монотонных функций д(х) это всегда выполнено.
Для уменьшения дисперсии расчетов можно также использовать сложную симметризацию, при которой интервал (а, Ь) разбивается на конечное число частей и для каждой из них используется метод противоположной переменной. Пусть, например, отрезок [а, Ь] разбивается на две равные части. Обозначим с ~ (а + Ь)/2. Тогда
h= I g{x)dx+ g(x)dx = - (д(х) + д(а + с-х)) dx +- (g(x) + g(c + b~x)) dx.
Ja J с w Ja ^ J с
В первом из этих интегралов сделаем замену переменных у = 2х — а, которая преобразует интервал (а, с) в (а,Ь), а во втором - замену у = 2х — Ь, которая преобразует интервал (с, Ь) в (а, Ь). При этом получаем соотношение
h = Ґ 9{2) (У) dy = ЕС(2>, где С(2) = (Ь - a)gW (а + (6 - а)а),
J а
9{2){у)-\
,а + у\ (2а + Ъ-у\ (b + y\^ f2b + a-y
2 * \ 2 * \ 2
9І-7Г-) +9\ ц +9[-тг-]+9'
Процесс сложной симметризации можно продолжить, разделив отрезок [а, Ь] на более чем два равных отрезка и реализовав соответствующее преобразование координат. Предельные свойства этого процесса изучаются в подразделе 1.2.3. Кроме того, в этом подразделе формулируются и исследуются многомерные аналоги процесса сложной симметризации и устанавливается связь этих конструкций с известными дискретно-стохастическими схемами численного интегрирования (в частности, с описанными выше оптимальными кубатурнымл формулами Н. С. Бахвалова).
Наконец, в диссертации относительно подробно изучается и такой способ уменьшения дисперсии DC, как выборка по группам (см. [7-9], а также подраздел 1.2.4). Здесь модификация алгоритма 0.2 состоит в том, что область X разбивается на К непересекающихся подмножеств X = yYi U ... U Хк, и в каждом подмножестве Х{ выбирается щ случайных точек согласно усеченной плотности /j(x) = f(x)/ JXt /(у) dy, при этом
пі +... + пк = п. Здесь уместно заметить, что предельными случаями такого алгоритма являются описанные выше оптимальные кубатурные формулы Н. С. Бахвалова для С1(Х) и С2(Х) из [12] (при /(х) = 1). Заметим, что в работе [11] показаны возможности применения дискретно-стохастических численных схем при реализации других известных способов понижения дисперсии: интегрирования по части области, метода математических ожиданий, метода расщепления и др.
В диссертации представлен ряд методик уменьшения величины t из (0.7). В разделе 1.1 представлены методы моделирования многомерных случайных векторов, плотность распределения которых имеет вид (0.10) или (0.11) (иными словами, способы эффективного решения задачи 0.2). В разделе 1.3 рассмотрен двусторонний геометрический метод, который оказывается весьма эффективным в случае, когда вычисление значений подынтегральной функции д является трудоемким.
В разделах 2.1 и 2.2 рассмотрены возможности применения смещенных состоятельных оценок в численном интегрировании. В частности, в разделе 2.1 рассмотрена следующая модификация стандартного алгоритма 0.2 [7, 41]. Рассмотрим вычисление интеграла по Z-мерному единичному кубу:
1= [ g(x)/(x)dx= [ ... / q(xM,...,x)f(xM,...,xV>)dxw...dx. (0.15)
Jqi Jo Jo
В качестве приближения этого интеграла возьмем величину
п п
I « ] = Х>(аО/(аО/Х)/("0. (-16)
i=i i=i
где вектор а = {&(1\ ..., а^) равномерно распределен в кубе Qi. Такая оценка может оказаться эффективной в случае, когда алгоритм реализации выборочных значений вектора согласно плотности /(х) является трудоемким (т.е. когда величина t из (0.7) достаточно велика). В разделе 2.1 рассмотрена ситуация, когда плотность / представляет собой приближение подынтегральной функции д вида (0.12). В разделе 2.2 рассмотрена оценка с поправочным множителем
e»M = (s«(6))х (;!>«<>) ((Ш)
где h(x) - вспомогательная функция, такая, что Е/г() '= 1. Несложно показать (см., например, [7]), что
Е0 = / - - Г д(х)(1 - А(х))/(х) Ас, (0.18)
п J х
D#> = \ fx (ї(х) - /(2 - h{x))f /(х) dx + O (J^ . (0.19)
Из соотношения (0.18) видно, что если h(x) = 2 — q(x)/I, то главный член выражения для D^2) обращается в нуль и D0^2) = 0(1/п2). При этом смещение в нуль не обращается и, согласно (0.18), равно Е0&2) - / = -~Dq()/(In). Заметим, однако, что оптимальный выбор функции h(x) = 2 — q(x)/I невозможен (величина / неизвестна). В разделе 2.2 предлагается использовать следующее приближение описанного оптимального случая. Для простоты рассмотрен случай вычисления интеграла (0.15) и /(х) = 1. Тогда = а и д(х) = д{к). В качестве функции берется h(x) = 2 - LA/fx LMg(x) dx известен.
В разделе 2.3 сравниваются численные алгоритмы аппроксимации решения интегрального уравнения второго рода на основе рандомизации конечных и бесконечных отрезков ряда Неймана. Исследуются возможности использования отрезков однородных цепей Маркова конечной и случайной длины, а также эффективных дискретно-стохастических методов численного интегрирования. Представлены примеры, подтверждающие целесообразность использования смещенных детерминированно-стохастичес-ких оценок решения.
0.6. Тестирование алгоритмов численного интегрирования. Важнейшим элементом исследования алгоритмов 0.1, 0.2 и их модификаций является тестирование. При этом для параметров задачи 0.1 (подынтегральной функции д и области интегрирования Лг) должны выполняться следующие ТРЕБОВАНИЯ:
(0.1а) чтобы соблюсти «независимость» тестирования, нужно добиваться того, чтобы вид (график) функции g был случаиліьш, заранее непредсказуемым;
(0.16) для рассмотрения случаев «сложных» подынтегральных функций g нуоісно, чтобы имелась возможность варьировать вычислительные затраты на получение одного значения функции;
(0.1 в) нужно, чтобы хотя бы в простейших ситуациях (например, в одномерном случае при I = 1) можно было проверить расчеты аналитически;
(0.1г) в связи с тем, чт,о сходимость многих кубатурных формул обусловлена требованиями вида g Є В(Х), должна присутствовать возможность контроля свойств подынтегральной (функции g (в частности, свойств гладкости);
(0.1д) для изучения теоретической эффективности алгоритмов 0.1 и 0.2 нуоісно, чтобы имелась возможность получать уточненные верхние границы для погрешностей используемых численных методов интегрирования с заданными множествами функций g и областей X.
В связи с требованием (0.1а) возникает идея использования в качестве g модельных траекторий случайных функций'(соответствующие модели описаны, в частности, в работах [8, 9, 42, 43]). При реализации этой идеи (см. [16, 19, 21], а также главу 3 данной работы) пришлось столкнуться с тем, что подходящих (удовлетворяющих требованиям (0.16)-(0.1д), не связанных с конкретными приложениями) моделей случайных процессов и полей не так уж много. Тем не менее, достаточно удачным оказался выбор в качестве тестовых функций траекторий спектральных моделей однородных гауссовских случайных нолей [8, 9, 42, 43] с конечным спектром вида
к
д(х) = Л'Ых), Ых) = >* К1} cos(x, Л,) + 7f sin(x, Afc)], (0.20)
it=i
Здесь 7^, г = 1,2 - независимые в совокупности стандартные нормальные величины; х Є X = (0,1)1. Векторы Ад, выбираются случайно в области Л = {0,А)1. В данной работе рассмотрены два варианта выбора. Первый способ - с разбиением спектра -связан с соотношениями Afc Є Afc и Ад, ~ р{Х)/ри (знак «~» означает «распределен согласно плотности»), где К = т1; Ак — [kih, (ki + l)h) х ... х [kth, (ki + l)h); h = А/т; h = 0,..., т - 1 и рк = jAk р(А) сіЛ = 1/К.
Второй способ выбора Afc - без разбиения спектра - предусматривал численное моделирование Afc в Л согласно плотности равномерного распределения
р(Л) = < -у при А Є [0, А)1; 0 иначе \ . Для обоих способов ак — 1/уК.
Заметим, что тестовые подынтегральные функции (0.20) представляют собой тригонометрические полиномы. Для получения более широкого спектра тестируемых подынтегральных функций в качестве д из (0.1) можно использовать реализации негауссов-ских спектральных моделей. В качестве модели негауссовского поля можно взять сумму вида (0.20) с Afc Є Afc и cHamA^ —> 0. Однако здесь величины \kl и Хк2 при к\ ф к2 должны быть зависимыми.
В связи с этим возможны два пути получения алгоритмов построения спектральных моделей негауссовских случайных функций:
брать наиболее распространенные негауссовские стационарные случайные функции, пытаться находить зависимость между элементами разбиения спектра и строить модели с учетом этой зависимости;
брать такие виды зависимости олеліентов спектра, которые легко моделировать, а затем исследовать, какие случайные функции получаются при моделировании.
Второй путь, который представляется автору диссертации более перспективным, реализован в данной работе (см. раздел 3.1). Негауссовские случайные функции реализуются за счет различных преобразований гауссовских моделей. Исследованы преобразования, основанные на рассмотрении функций на гауссовских траекториях, а также на комбинировании гауссовских спектральных моделей со случайными величинами.
Важным видом сходимости гауссовских и негауссовских моделей вида (0.20) (в том числе, и для приложения, предлагаемого в данной работе, то есть для решения задачи 0.1) является слабая (или функциональная) сходимость в пространстве С(Х). Для гауссовских моделей этот вид сходимости подробно изучен в [44] (см. также [9, 11, 43]). В данной работе рассмотрен вопрос о слабой сходимости в негауссовском случае (см. раздел 3.1).
В разделе 3.2 приведено описание используемой в расчетах тестовой системы подынтегральных функций, а в разделе 3.3 получены средние оценки погрешностей простейших квадратурных формул (а именно, формул прямоугольников, трапеций и Симпсо-на) для рассматриваемых случайных подынтегральных функций. Аналоги этих оценок получены в главах 1 и 2 для различных дискретно-стохастических модификаций алгоритма 0.2.
Модельные (гауссовские и негауссовские) траектории случайных полей можно использовать при построении сложных областей интегрирования X. Например, можно взять X — {х : л-(х) > 0}. Подробного изучения этой возможности в данной работе не представлено.
0.7. Апробация работы. Описанные в диссертации результаты были представлены в докладах на международной конференции по вычислительной математике в Новосибирске (июнь 2002 года) [18]; на Международной конференции по Монте-Карло и квази-Монте-Карло методам (Ульм, Германия, август 2006 года) [23]; на международных семинарах-совещаниях «Кубатурные формулы и их приложения»: в Красноярске (сентябрь 1999 года) [14], в Уфе (июль 2001 года) [17], в Красноярске (август 2003 года) [20, 22], в Улан-Удэ (август 2005 года) [24, 34], в Уфе (июнь 2007 года) [25-27]; на конференциях молодых ученых Института вычислительной математики и математической геофизики СО РАН (март 2001 года [29], март 2002 года [30], апрель 2004 года [32]); на Международных студенческих конференциях «Студент и научно-технический прогресс» (Новосибирск, НГУ; апрель 2003 года [19], апрель 2004 года [31], апрель 2007 года [35-37]) и неоднократно на семинаре отдела статистического моделирования в физике ИВМиМГ СО РАН.
Моделирование случайных величин
Наибольшее внимание в главах 1 и 2 будет уделено использованию приближений вида (1.1) в алгоритмах численного интегрирования (см. далее разделы 1.2, 1.3, 2.1, 2.2). В этом разделе мы отметим ряд других возможностей использования приближений (1.1) в методах Монте-Карло. Прежде всего рассмотрим задачу 0.2 о моделировании случайного вектора (для / = 1 - случайной величины) , имеющего плотность распределения f(x) = CLMg(x), С = 1/с, с= [ LMg(x)dx; д(х) 0; (1.3) см. также формулу (0.12). Приведем примеры ситуаций, в которых возникает необходимость решать задачу 0.2. ПРИМЕР 1.1 (использование гистограммы и полигона частот). Пусть параметр 9 моделируемого физического процесса V случаен и имеется возможность получить выборочные значения в[ ,...,#[ этого параметра с помощью проведения достаточно дорогих экспериментов. Пусть также Wy является математической моделью, описывающей процесс V и включающей параметр 9, причем численная реализация модели Wy требует моделирования большой выборки e[Wv\...,e v\ где К к. (1.4) В этом случае можно разделить интервал [а, 6], a min ,..., ), Ъ = пш( ,... ,6) на полуинтервалы [zt, zl+i); г = 1,..., М — 1; a = zx z2 ... ZM = b и вычислить частоты v\ = тг/к (здесь т1 есть число значений Щ }, попавших в г-й полуинтервал). Пусть уі = (zi + Zj+i)/2 и g(yi) — vl/(zi+x - zt). Тогда аппроксимация вида (1.1) может быть использована для реализации выборки (1.4). Для кусочно-постоянного случая имеем Хг(х) = 1 при .г- Є [zi,zi+i}; Wi — д(уі) и С = 1, и функция (1.1) называется гистограмлюй [45]. Соответствующая кусочно-линейная версия функции (1.1), (1.2) называется полигоном частот [45]. В этих простейших случаях выборочные значения (1.4) могут быть получены с помощью модифицированного метода суперпозиции, который в простейшем случае совпадает с методом обратной функции распределения [9].
Пример 1.2 (приближение «слоокных» плотностей). Рассмотрим ситуацию, в которой функция д является плотностью случайной величины или вектора 0, и не удается построить алгоритм численной реализации выборочных значений этого случайного элемента. В этом случае целесообразно использовать плотность вида (1.3) вместо д и моделировать соответствующую случайную величину , причем приближение (1.1) должно обладать достаточно хорошими аппроксимационными свойствами, т. е. расстояние pB{X)(y,LMg) должно быть мало; здесь рв(х) метрика функционального пространства В(Х). Здесь уместно заметить, что имеется хорошо развитая статистическая теория приближения плотностей [46], в которой в основном изучается случай пространства В(Х) = Li(X), однако метрика этого пространства является недостаточно «сильной» для приложений теории методов Монте-Карло [11]. Описание примера 1.2 закончено.
Рассматриваемая в примере 1.2 ситуация, с одной стороны, аналогична примеру 1.1 (при / = 1), а с другой стороны, применение кусочно-постоянных и кусочно-линейных приближений недостаточно для получения требуемых аппроксимационных свойств приближения (1.1) в представлении плотности (1.3). Поэтому необходимо применение более сложных функциональных базисов Н(Л/). Однако в этом случае, как правило, не работают алгоритмы моделирования случайных векторов, основанные на методе обратной функции распределения [7-9] (даже для I = 1). Поэтому при решении соответствующей задачи 0.2 о моделировании случайной величины согласно плотности (1.3) возникает идея использования методов суперпозиции [7-9], при этом появляются новые требования люделируемости к базису Е л/ (см. далее соотношение (1.5) и подраздел 1.1.3).
Выберем функции Н(Л/) так, что для случайных величин i; распределенных согласно плотностям {/г} из (1.6), имеются эффективные алгоритмы численного статистического моделирования. Тогда возникает следующий АЛГОРИТМ 1.1 [7-9]. Согласно вероятностям {Pi} из (1.6) выбираем номер т : Р(т = г) = Pi. Реализуем согласно плотности /т(х). Количество узлов сетки М может быть достаточно большим, и тогда затраты на поиск номера т могут быть также большими. В этом случае целесообразно использование квантильного метода. АЛГОРИТМ 1.2 [9]. Зададим целое число К и введем на интервале (0,1) равномерную сетку {и)к — к/К; к = 0,1,...,К}. Далее построим массив целых чисел Xk — min{j : Rj = Pi + ..1 + Pj w -i}, к = 1,..., К. Этот массив задает номер j элемента массива {Rz; і = 1,..., п}, с которого следует начинать поиск «вверх» (т. е. вычитать величины Rq, q = j,j+l,... из а до получения первого отрицательного значения) при uik-i а w . Окончательно моделирование дискретной случайной величины выглядит следующим образом: 1. Реализуем значение а равномерно распределенной на интервале (0,1) случайной величины. 2. Вычисляем номер к полуинтервала [wk-i, иік), в который попадает а по формуле к = [К а] + 1; здесь [А] целая часть числа А. 3. Реализуем последовательный поиск «снизу вверх» начиная с Rxh Пример 1.3 {метод исключения). Одним из основных методов реализации выборочных значений случайной величины или вектора 1; распределенного с плотностью /i(x), является мажорантный метод исключения. АЛГОРИТМ 1.3 [7-9]. Рассмотрим неотрицательную функцию g+(x), для которой справедливы соотношения 3"+(х) /i(x) для х Є X, g+(x) = 0 для х 0 X, G+ — I g+(x) dx со, и существует эффективный алгоритм численной реализации выборочных значений случайного вектора /(х), /(х) = /+(х)/ 3+. (1.7) Реализуелі пары ($,у) (здесь вектор - из (1.7), а случайная величина 7 имеет вид 7 = ад+()) до тех пор, пока не будет выполнено неравенство 7 /i( ) f (1-8) и тогда принимаем х = .
Среднее число реализаций (/ + 1)-мерного вектора (, 7) требуемых для получения одного выборочного значения х, равно G+ 1 [7-9]. Это число пропорционально трудоемкости алгоритма 1.3. Поэтому следует выбирать мажоранту д+ так, что реализация вектора достаточно эффективна и число G+ близко к единице, т. е. функция д+ должна быть близка к плотности /і(х). Для построения такой мажоранты целесообразно использовать приближения «сверху» для плотности /і (здесь, особенно в многомерном случае, эффективными оказываются кусочно-постоянные аппроксимации); в этом случае плотность /(х) имеет вид (1.3). ПРИМЕР 1.4 (двусторонний метод исключения). Существует ряд модификаций алгоритма 1.3. В случае, когда вычисление значения функции /і(х) (а значит, и проверка неравенства (1.8)) является трудоемким, возможно применение следующего алгоритма, который носит название двусторонний метод исключения. АЛГОРИТМ 1.4 [9]. Помимо мажоранты д+ выберем миноранту д- такую, что вычисление ее значений не является трудоемким и выполнено неравенство О д (х) /і(х) для х Є X. Снова, как и в алгоритме 1.3, реализуем пары (,7) до тех пор, пока не будет выполнено одно из неравенств 7 #-(0 или (1-8), и затем берем г = . Если функции fi,g+,g близки, то проверка неравенства (1.8) (а значит, и вычисление значений «сложной» функции /i(x)) будет происходить редко. Поэтому в качестве 7_ целесообразно выбирать аппроксимацию «снизу» (например, кусочно-постоянную) функции /i(x). Двусторонняя технология с кусочно-постоянными мажорантой и минорантой может быть применена при реализации геометрического метода Монте-Карло (см. далее раздел 1.3).
Дискретно-стохастическая версия метода выборки по группам
На практике дисперсии (1.55), как правило, неизвестны и выбор {пт} по формуле (1.56) невозможен. Эта формула используется на практике при выборе чисел испытаний {пт} в методе (1 53) Весьма важным является предельный случай расслоенной выборки „ , для которого число делений М области X близко к числу испытаний п. Здесь на классах гладких подынтегральных функций удается получить оптимальные (по порядку t вероятностной погрешности 6п п 1 при п — со) алгоритмы численного интегрирования [4, 9, 52]. Рассмотрим случай М — п, щ — ... — пм = 1; X = Qi и /(х) = 1, х Є Qi (здесь заведомо выполнено условие (1.57)). АЛГОРИТМ 1.8 [52]. Реализуем по одной равномерно распределенной случайной точке т в каждом кубе Хт вида (1.44). Вычисляем интеграл (1.27) согласно фор муле1 йМ) = 1/МЕ 9(и Для функций #(х) из класса C l," {L,..., L; Qi) в [4, 9, 52] для этого случая получено неравенство Р(Є = -ЙМ,І Л;ДЇ7Г) І — Полученный порядок t — 1/2+1/1 является для погрешности 8п n L неулучшаемым. Это дает следующая ЛЕММА 1.4 [52]. Существуїот положительные константы Н(г,Ъ) и Р, удовлетворяющие следующему соотношению. Для любого натурального М и любой расслоенной выборки Q вида (1.53) из алгоритма 1.8 найдется функция д(х) Є Cr(L;Qi), для которой бімНг,Ц = \йМ)-і\ Ш (1-58) То с вероятностью Р; здесь т — (r1;..., r{), L = (Х1;..., L/), \/т = \/г\ + ... + 1/г;. Для рассматриваемого множества функций имеем г = (1,..., 1), L = (L,..., L) и г — 1/7. Следовательно, согласно соотношению (1.58), для класса C l - (L,... ,L;Qi) получаем 5п H/nl 2+1//l.
Заметим, что алгоритмы 1.8, 1.9 допускают несколько названий. Во-первых, это частные случаи метода расслоенной выборки (алгоритм 1.7) для п = М и п = 2М соответственно. Во-вторых, алгоритм 1.9 можно назвать модификацией метода многомерной сложной симметризации (1.51) (для N — 1), в которой розыгрыш случайных точек т и тіЗІт происходит во всех множествах (1.44) (по одному разу). В-третьих, в связи с наличием дискретизации области интегрирования X = Qi, определяемой разбиением на малые кубы Хт вида (1.44), и с выбором одной или двух случайных точек в каждом малом кубе молено отнести алгоритмы 1.8, 1.9 к дискретно-стохастическим алгоритмам численного интегрирования. В-четвертых, алгоритмы 1.8, 1.9 являются частными случаями случайных кубатурных формул (см., например, [9]).
Как уже отмечалось в подразд. 0.2 введения, представленные здесь предельные результаты для g Є C " (L,..., L; Qi) и g Є C2(L\ Qi), полученные в начале шестидесятых годов 20-го века Н.С.Бахваловым, вызвали большой научный резонанс и привели к бурному развитию теории сложности численных алгоритмов (см. монографию [13] и сопутствующие работы). Классической задачей в этой теории является следующая: при фиксированном числе операций п определить максимальный порядок t погрешности 5п п 1 (в обычном или вероятностном смысле) заданного класса вычислительных алгоритмов. Алгоритмы численного интегрирования являются наиболее удачными иллюстрациями конструкций PI методик теории сложности.
В ходе работы над диссертацией (см., например, статьи [21, 22, 29, 33]) были проведены многочисленные расчеты для сравнения эффективности представленных в этом разделе дискретно-стохастических методов уменьшения дисперсии. Тестирование производилось как с помощью набора специально подобранных функций, так и с использованием тестовой системы (3.21) из раздела 3.2. На основе численных экспериментов можно сформулировать следующие выводы. 1. Для размерностей / 5 вычисление интегралов с помощью оптимальных алгоритмов 1.8 и 1.9 от всех тестовых функций более экономично (особенно эффективен здесь алгоритм 1.9 - см. рис. 1.1, 1.2/ 2. При вычислении интегралов от достаточно «гладких» функций из стохастической системы функций (т. е. для небольших А) и размерностей I 7 время работы метода выборки по важности (1.27), (1.28) сравнимо с затратами алгоритмов 1.8 и 1.9. Но для высоких размерностей интеграла и функций малой гладкости этот метод малоэффективен (см. рис. 1.1, 1.3/ 3. При вычислении многократных интегралов I 6 эффективность простейшего метода Монте-Карло (1.27) с /(х) = 1 превосходит эффективность «оптимальных» алгоритмов 1.8 и 1.9 (см. рис. 1.2, 1.3/ 4- Трудоемкость оптимальных 1.8 и 1.9 во всех случаях не превосходит трудоемкости метода прямоугольников (1.60). При небольших размерностях интеграла метод прямоугольников (1.60) эффективнее простейшего метода Монте-Карло. 5. Расчеты подтвердили зависимости от параметров К и А выборочных дисперсий рассматриваемых алгоритмов для случая применения подынтегральных функций g из стохастической тестовой системы (3.21) (т. е. численно подтверждены соотношения (1.35), (1.43), (1.52), (1.61); см. рис. 1.1-1.3/ 1-Ц , Г Ц , , 10 100 с 10QQ 1000 10000 к 100000 Рис. 1.3. Время вычисления интегралов кратности I = 7 и I = 9 методами Монте-Карло (1.27) с /(х) = 1 - МК, прямоугольников (1.60) - гее и с помощью алгоритмов 1.8 - С1 и 1.9 - С2, а также методом выборки по важности (1.27), (1.28) - ImS для случаев А = 0.5 и А = 10 соответственно.
Таким образом, относительно сложные схемы численного интегрирования, имеющие теоретически оптимальные порядки погрешности, могут проиграть по трудоемкости более простым вычислительным конструкциям. Кроме того, подтвержден тезис о том, что для больших размерностей (в нашем случае I 10) следует применять простейший метод Монте-Карло (1.27) с /(х) = 1. Смешанные дискретно-стохастические численные схемы (алгоритм выборки но важности (1.27), (1.28), алгоритмы 1.8 и 1.9) могут быть эффективными для «промежуточных» размерностей (/ — 3 — 9). Для малых размерностей и гладких подынтегральных функций неплохо работают простейшие кубатурные формулы (вплоть до формулы прямоугольников (1.60)). В ходе численных экспериментов установлено также, что методы Бахвалова (алгоритмы 1.8 и 1.9), предложенные им для пространств функций Cl(L;Qi) C2(L\Qi), можно применять и для кусочно непрерывных функций с кусочно непрерывными производными. При этом порядок сходимости методов не изменится, если элементарные области выбирать определенным образом.
Приведем также результаты сравнительных расчетов из нашей работы [27], посвященной исследованию многомерной версии алгоритма сложной симметризации (см. таблицу 1.4). При тестировании здесь сравнивались: алгоритм сложной симметризации (1.51) с N парами симметричных точек в каждом элементе разбиения (алг. 1А), предельный случай метода сложной симметризации для N = 1 (алг. 1Б), модификация выборки по группам (1.53) с N парами симметричных точек в каждом элементе разбиения (алг. 2), алгоритм 1.9 для N = 1 (алг. 3) и метод прямоугольников (1.60) (алг. 4). В таблице 1.4. приводятся средние погрешности рассмотренных методов для функций (3.21) тестовой системы из главы 3 и для различных I, А и К при одинаковом количестве выборочных значений (т.е. при одинаковом времени расчетов).
Взвешенная равномерная выборка
Заметим, что на основании леммы 2.1 можно построить доверительные границы для погрешности 5п —\вп—і\ оценки (2.4). Действительно, лемма 2.1 позволяет утверждать, что случайная величина 9п для больших п имеет распределение, близкое к нормальному с математическим ожиданием I и дисперсией D\. Следовательно, для больших п выполнено приближенное равенство Р($г js/Di) Р(Ы 7)) где 7 - положительная константа, а ш Є iV(0,1) - стандартная нормальная случайная величина.
На основании соотношений (2.9), (2.11) можно провести следующую оптимизацию оценки (2.8) по параметру ц (при фиксированном уровне погрешности 5). Затраты на реализацию этой оценки примерно равны s = Т(/л) + n(ii(//) +t2), где Т(ц) - время на предварительное вычисление величины /, ti{jj) - среднее время вычисления величины Ьмй{оі-і)і a t2 - среднее время для вычисления величины g(cxi) (затраты на сложение этих величин, на умножение на I и на одно деление в (2.8) считаем малыми и не учитываем). Считая, что верхняя граница из соотношения (2.11) воспроизводит реальное поведение погрешности 5п , получаем, что при фиксированном S затраты s равны s{1)=TM + (tM+t2) (2.12) где Ci(fx) - величина из (2.11) при условии, что плотность / выбирается но формуле (2.7). Заметим, что получение точных аналитических зависимостей величин Т, Сі и х от ц затруднено (здесь многое зависит от конкретной реализации оценки (2.8) на ЭВМ). Однако анализ баланса (2.12) приводит к выводу о наличии оптимального значения fj, и увеличении этого значения с уменьшением 5.
Отметим также, что при относительно больших размерностях / и значениях числа делений по каждой координате /л возникают проблемы (по предварительным затратам Т(ц), по использованию оперативной памяти ЭВМ) при построении функции LMQ-Таким образом, с ростом кратности интеграла эффективность применения как самого метода выборки по важности (см. разд. 1.2), так и его модификации (2.8) падает.
Сравнивались: стандартный алгоритм (2.1) для интеграла (2.3) с /(х) = 1 (АЛГОРИТМ 1), алгоритм выборки по важности (2.1), (2.7) (АЛГОРИТМ 2) и его модификация (2.8) (АЛГОРИТМ 3). Для заданного уровня погрешности 6 с помощью серий независимых (по используемым стандартным случайным числам) экспериментов для всех алгоритмов определялось характерное число испытаний п, при котором приближенное значение интеграла гарантированно входило в «5-полосу» около точного значения / = 6.927042. Затем для этого числа испытаний фиксировались временные затраты ЭВМ t (в секундах).
Были также проведены эксперименты для подынтегральных функций стохастической тестовой системы (3.21), представленной далее в разд. 3.2. Для метода выборки по важности на этапе случайного выбора ячейки разбиения куба Qi, в которую попадает очередное выборочное значение, использовался квантильный метод (алгоритм 1.2). В таблицах 2.3 и 2.4, в частности, представлены: доверительные границы е2 = Зу/Щ/п и єз = Зу/ТЭвР для оценок АЛГОРИТМА 2 и АЛГОРИТМА 3 соответственно, средние погрешности до и 53 вычисления интегралов для серий из десяти-одиннадцати независимых (по используемым стандартным случайным числам) экспериментов, времена вычисления интегралов ti,U,h- В квантильном методе число квантилей полагалось равным k\ = fx4 и к2 = /І3 (этим случаям соответствуют времена t\ и t2). Из представленных в таблицах 2.3 и 2.4 результатов можно сделать следующие выводы. 1. При одинаковом уровне допустимой погрешности є и числе испытаний п взвешенная равномерная выборка дает чуть меньшую реальную погрешность ( по сравнению с погрешностью 82 метода выборки по ваоїсиости. 2. Улучшение реализации метода выборки по ваоїсности (в частности, использование квантильного метода с большим числом квантилей) приближает трудоемкость этого метода к трудоемкости метода равномерной выборки. 3. Соотношения трудоемкостей метода выборки по важности и метода взвешенной равномерной выборки качественно не меняются с изменением свойств тестовых подынтегральных функций (т. е. с ростом параметров А и К).
Отметим, что в тестовых расчетах принял участие студент ММФ НГУ А. И. Ефремов. 2.2. Дискретно-стохастическая версия метода Монте-Карло с поправочным множителем Оценка с поправочным множителем. Рассмотрим еще одну модификацию алгоритма (2.1), связанную с применением леммы 2.1. Пусть требуется вычислить интеграл (2.1) с заданной плотностью вероятностей /(х), определенной в X, и 1;..., п - независимые реализации случайного вектора , распределенного согласно плотности /(х). Введем вспомогательную функцию Я(х), х Є X, и рассмотрим случайную величину.
Получение точных аналитических зависимостей величин Т, ( и t\ от fj, затруднено (здесь многое зависит о г конкретной реализации оценки (2.14) на ЭВМ). По аналогии с соотношением (2.12), анализ баланса (2.24) приводит к выводу о наличии оптимального значения /л и увеличении этого значения с уменьшением 8. Этот вывод подтвержден в тестовых экспериментах - см. далее таблицы 2.6-2.10.
Отметим также, что при относительно больших размерностях I и значениях числа делений по каждой координате (л возникают проблемы по предварительному вычислению величины / и по использованию оперативной памяти ЭВМ для многократного вычисления значений функции L g. Таким образом, с ростом кратности интеграла эффективность применения модификации (2.14) падает.
Сравнивались: стандартный алгоритм (2.1) для интеграла (2.3) с /(х) = 1 (АЛГОРИТМ 1) и алгоритм с поправочным множителем (2.14), (2.20). Для заданного уровня погрешности 6 с помощью серий независимых (по используемым стандартным случайным числам) экспериментов для всех алгоритмов определялось характерное число испытаний п, при котором приближенное значение интеграла гарантированно входило в « -полосу» около точного значения /. Затем для этого числа испытаний фиксировались временные затраты ЭВМ t (в секундах). Для демонстрации эффекта уменьшения дисперсии при использовании метода с поправочным множителем в таблицах приведены также величины а: для АЛГОРИТМА 1 это корень из выборочной дисперсии, а для алгоритма (2.14), (2.20) - корень из выборочного приближения величины nD2 (см. формулу (2.16)).
Численные спектральные модели гауссовских случайных полей
Из центральной предельной теоремы для случайных функций [57] следует, что если в формуле (3.3) случайные величины щ независимы (а именно такие величины легче всего реализовать на ЭВМ), то при выполнении условий (3.4) конечномерные распределения модели (3.3) сходятся к гауссовским распределениям. Поэтому модели вида (3.3) используются для моделирования только гауссовских случайных полей о [9, 11]. 3.1.2. Тестовая спектральная модель. Нами в работах [16, 20, 22] исследована и использована численная модель (3.3) однородного гауссовского вещественного случайного поля о (здесь и далее для вещественного случая в (3.1)-(3.4) используются обозначения для поля 0, спектральной плотности / и т.д. - без волны над буквами). Модель (3.5)-(3.7) называется численной моделью однородного гауссовского поля с разбиением спектра [9].
Как отмечено в [8, 9, 56], одномерные распределения модели (3.5) являются стандартными гауссовскимп. Действительно, условные одномерные распределения / -(хо) при фиксированных {А } являются стандартными нормальными. Далее, применяя формулу полной вероятности, получаем, что «глобальное» распределение величины /с(хо;7)- ) в вероятностном пространстве, являющемся прямым произведением соответствующих вероятностных пространств для 7 = {lk \ j = 1,2} и А = {Ajt}, также будет стандартным. Конечномерные распределения модели (3.5), отличные от одномерных, не совпадают с предельными.
Здесь ф - функция (неслучайная), определенная на интервале (—со, +оо), а случайная функция #- является спектральной моделью вида (3.5) - с разбиением или без разбиения спектра. Рассмотрим вопрос о том, как выглядят одномерные распределения случайного поля с . В силу того, что одномерные распределения моделей (3.5) являются стандартными гауссовскими, плотность р= одномерного распределения случайной функции (3.10) совпадает с плотностью Рф(7) распределения случайной величины V(T) гДе 7 Є (0,1) - стандартная нормальная величина. Сразу заметим, что можно добиться того, чтобы случайные функции 3$ из (3.10) имели заданную одномерную функцию распределения .Fxo для х = Хо, если взять фхо{и) = F 0l $(v) (в этом случае моделирование случайной функции согласно формуле (3.10) называется методом обратной функции распределения [9, 43]), где Ф - функция распределения случайной величины.
Рассмотрим функции zx,z2 Є С(Х) и ух = 7 ( ), //2 = 0) В силу равномерной непрерывности -0 функции т/1 и т/г принадлежат С(Х) и, кроме того, существует 6 О такое, что для всех Ші,и 2 Є Я таких, что гиі — w2\ «5, следует, что т/ (гоі)—т/ (гУ2) " ь Пусть теперь рс( ьг2) — suPxexl ifc) г2(х)( 5. Тогда для всех х Є X выполнено zi(x) - z2(x)\ S и, значит, т/і(х) - т/2(х) = ( (х)) - ip(z2(x))\ Sv В силу произвольности х Є X для выбранных нами функций т/і и т/2 выполнено соотношение Рс(УііУ2) і а значит, справедливо неравенство (3.18). Таким образом, для заданного є 0 нашлось 5 0 такое, что для всех i,z2 Є С(Х) таких, что pc(zi,z2) 5, выполнено \F(zi)-F(z2)\ = \F($(zi))-F(ij (z2))\ = \F(yi)-F(y2)\ є, т. е. функционал F является непрерывным в метрике рс- Утверждение 3.1 доказано.
Функция распределения Ф стандартной случайной величины у, несомненно, является равномерно непрерывной на (—оо; +оо), поэтому для сходимости последовательности (3.10) в методе обратной функции распределения достаточно потребовать, чтобы функция F 1 была непрерывной на отрезке [0,1]. Последнему требованию удовлетворяют распределения, сосредоточенные на отрезке.
Из доказанного утверждения 3.1 также следует, что если последовательность (3.5) слабо сходится к моделируемому полю о и функция g(u,v) равномерно непрерывна по и на всей числовой прямой при любом фиксированном v = VQ, то последовательность (3.12) слабо сходится в С(Х) к функции Со(х) = /(о(х),?7)- Заметим, что модели {С/о 1 2} с конечным спектром (3.6) для всех рассмотренных выше примеров функций д, кроме последнего, слабо сходятся к (0.
В связи с применением рассматриваемых здесь моделей для тестирования алгоритмов численного интегрирования (см. разделы 3.2 и 3.3) уместно заметить, что функционал F(K) — fx к-(х) с/х сходится по распределению к F(Q) при К — оо, так как этот функционал непрерывен в метрике рс, ведь для любых т/1,7/2 Є С(Х) выполнено неравенство \F{y{) - F{y2)\ Jx\yi[x) - т/2(х) dx mesXрс(уі,у2) 3.1.6. Группировка слагаемых в моделируемой сумме. С помощью рассмотренных выше приемов получения негауссовских распределений можно строить численные реализации случайных функций, представимых в виде сумм независимых случайных функций.
