CC BY
509
2015 ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА Сер. 1. Том 2(60). Вып. 1
МАТЕМАТИКА
УДК 519.644.2
АЛГОРИТМ ЧИСЛЕННОГО ИНТЕГРИРОВАНИЯ МЕТОДОМ КВАЗИ МОНТЕ-КАРЛО С АПОСТЕРИОРНОЙ ОЦЕНКОЙ ПОГРЕШНОСТИ*
А. А. Антонов
Санкт-Петербургский государственный университет,
Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7-9
В работе рассматривается возможность применения теоретико-вероятностных методов к детерминированной процедуре оценки погрешности метода квази Монте-Карло вычисления многократных интегралов. Существующие методы оценки упомянутой погрешности являются неконструктивными. Ранее полученный результат, выведенный при помощи теории случайных кубатурных формул, дополнен новыми теоретическими результатами. Алгоритм Qint, представленный в пошаговом виде, удобном для непосредственного применения, обсуждается в контексте параметризации. Продемонстрированы важные свойства Qint: монотонность оценки по количеству разбиений и зависимость точности от правила этих разбиений.
Предлагаются новые численные результаты, полученные с использованием последовательностей Соболя и подынтегральных функций из программной библиотеки TESTPACK. Демонстрируется преимущество по сравнению с оценкой, полученной традиционным методом Монте-Карло. Библиогр. 14 назв. Ил. 4.
Ключевые слова: метод Монте-Карло, квази Монте-Карло, доверительный интервал, случайные кубатурные формулы, функции Хаара, последовательности Соболя.
1. Введение. Как известно, методы Монте-Карло основаны на результатах теории вероятностей и математической статистики. Методы же квази Монте-Карло, имеющие формальное сходство с методами Монте-Карло, основаны на детерминистских (теоретико-числовых, в основном) результатах. Вопрос о сходстве и различии этих методов представляет несомненный интерес как с теоретической, так и с прикладной точек зрения. Имеется ряд работ, где теоретико-вероятностные методы применяются для анализа детерминированных последовательностей чисел. Достаточно заметить, что так называемые датчики псевдослучайных чисел являются детерминированными алгоритмами и законность их использования подтверждается обычно системой статистических тестов (см. Кнут, [8]).
Рассмотрение квазислучайных последовательностей с вероятностной точки зрения показывает, что члены этих последовательностей статистически зависимы, а при
* Работа выполнена под руководством д-ра физ.-мат. наук проф. С. М. Ермакова и при финансовой поддержке РФФИ (грант №14-01-00271).
вычислении с их помощью интеграла для оценки погрешности не может быть использована центральная предельная теорема так, как это делается в методе Монте-Карло. Оценка погрешности метода квази Монте-Карло теоретически решается с помощью неравенства Коксмы—Хлавки [9]. Практически же его использование связано с рядом трудностей, что отмечалось многими авторами (см., например, работу Нидер-рейтера [11], где указаны более эффективные подходы к анализу погрешности интегрирования). Таким образом, если при вычислении интегралов поведение средних
N
77 гДе -^з —квазислучайные векторы в единичном в-мерном гиперкубе II
3=1
N
во многом сходно с поведением средних ^ гДе Щ ~ равномерно распреде-
3 = 1
лённые случайные векторы в П8, то средние значения сумм квадратов /2(Х3) и / 2(23) пригодны для оценки погрешности во втором случае и не пригодны в первом.
Задачи настоящей работы состоят в том, чтобы выяснить, что может служить формальным аналогом дисперсии при использовании метода квази Монте-Карло, и указать эффективные средства оценки погрешности интегрирования с использованием такого аналога.
Первая из задач использует связи квазислучайных последовательностей со случайными квадратурными формулами и расслоенной выборкой. Пусть гиперкуб Ц
разбит на Ь непересекающихся подобластей Ц1,..., Ць. Расслоенная выборка пред-
ь
полагает, что в каждой области выбирается к3 случайных точек, так что ^ к3 = N.
3=1
N
Известно выражение для дисперсии суммы ^ /(^3) в этом случае (см. [5] в кон-
3=1
тексте генерации выборок из распределений, [4] — в контексте численного интегрирования). Если объемы всех Ц равны, а к\ = к2 = • • • = кь, то такая выборка обладает улучшенными свойствами равномерности. Подход, основанный на теории случайных кубатурных формул, связывает расслоенную выборку с улучшенными свойствами равномерности с формулами, точными для кусочно-постоянных функций с постоянными значениями в областях ... ь и с ортогональными разложениями по кусочно-постоянным функциям. В работе получен новый результат, обобщающий частный случай Ь = 2, исследованный ранее в [4].
Квазислучайные последовательности, предложенные И. М. Соболем, очевидным образом связаны с расслоением и ортогональными разложениями по кусочно-постоянным функциям (функции Хаара, см. [3]). Показано, что с квазислучайными последовательностями может быть связана дисперсия случайной квадратурной формулы. Приведённые в работе численные примеры показывают, что основанные на этой дисперсии доверительные интервалы могут оценивать остаток интегрирования.
Однако приведённые выше результаты не являются строго обоснованным методом оценки остатка. Строгое обоснование может быть получено на пути рандомизации квазислучайных последовательностей, который обсуждается далее.
Рассмотрим задачу численного интегрирования вещественнозначной функции / по единичному гиперкубу размерности в: ^ /(и)йи. При этом мы предполагаем, что размерность в достаточно высока, а функция / устроена достаточно сложно1. Как
1 В нашей работе мы сознательно ограничиваемся контекстом Монте-Карло и квази Монте-Карло. Сравнение с традиционными методами вычислительной математики остаётся за рамками
уже отмечено ранее, в подобных случаях на практике часто используется рандомизированный метод квази Монте-Карло (RQMC), который позволяет строить доверительные интервалы на основе оценки дисперсии. Такой подход ведет к вопросу выбора достаточного количества рандомизаций (см. [10]).
Так, пусть мы используем k рандомизаций {Pn,j}j=i одного и того же множества квази Монте-Карло Pn мощности п. При этом дисперсия может иметь порядок вплоть до (при существенных ограничениях на функцию /, см. [13]), в
то время как традиционный Монте-Карло имеет дисперсию порядка Поэто-
му для получения более ощутимого выигрыша необходимо увеличивать п, что при фиксированных вычислительных затратах ведет к уменьшению k. Это, в свою очередь, приводит к недостоверной оценке дисперсии, поэтому для практических задач рекомендуется выбирать k ^ 10.
Пусть мы имеем фиксированное k, тогда конечное RQMC-множество узлов интегрирования представляет собой объединение всех рандомизаций: P = uk=iPn,j. Если исходный набор Pn обладает некими хорошими квази Монте-Карло свойствами, а выбранный способ рандомизации их не нарушает, то результирующий набор такими свойствами, вообще говоря, уже не обладает. Например, если Pn есть (t, m, в)-сеть для некоторых параметров t и m, а в качестве рандомизации выбран digital shift (см., напр., [6]), сохраняющий структуру сети для каждого Pnj-, то в общем случае нельзя утверждать, что P является (t*, m*, в)-сетью для каких-то t* и m*. Вопрос о том, насколько сильно нарушается структура конечного RQMC-множества, исследован мало; очевидно лишь, что чем больше количество рандомизаций k, тем потенциально сильнее это нарушение.
Предлагаемый в статье подход позволяет, в частности, построить альтернативную оценку дисперсии рандомизированного метода квази Монте-Карло. Используя аппарат случайных квадратурных формул из [4], мы получаем формулу Монте-Карло, точную для ортонормированной системы Хаара в Us. Полученная формула, по сути, эквивалентна одному частному случаю расслоения. Полученный результат, однако, не является тривиальным, поскольку он
1) показывает связь аппарата случайных квадратур с традиционными методами понижения дисперсии;
2) имеет естественную интерпретацию в терминах интегралов по подмножествам Us определённого вида;
3) имеет дисперсию, не превосходящую дисперсии традиционного метода Монте-Карло с тем же количеством точек;
4) устанавливает свойства точности для хааровских систем, что дает возможность рассмотреть предложенный подход в рамках квазислучайных множеств, в частности, (t, m, в)-сетей.
В заключительной части статьи мы приводим численные примеры, свидетельствующие о перспективности предложенного гибридного метода.
2. Теоретические основы Qint. Рассмотрим произвольное вероятностное пространство (X, где X — непустое множество, B — а-алгебра для подмножеств X, снабжённая вероятностной мерой Выберем некоторое число P G No и возьмём N = 2P. Разбив X на N дизъюнктных частей равной меры Xi,X2,..., Xn, покрывающих собой всё X, построим ортонормированную в L2(d^) систему из N функций
данной статьи; среди подобных работ можно выделить [14], откуда, в частности, и взяты функции для вычислительных примеров, приведённых в заключительной части работы.
Хаара, порождённую Х4,Х2,..., Хм. Здесь и далее при проведении подобных построений неявно предполагается, что такое разбиение возможно для данного N. Свойства такого рода систем подробно описаны в [3], а явный вид функций и процесс построения указан в [1].
Чтобы построить оценку интеграла I = / f (х)^(^х) для функции / £ Ь2(^,),
Jx
воспользуемся аппаратом случайных квадратурных формул, изложенным в [4], в особенности теоремой Ермакова—Грановского. Эта теорема позволяет построить ^точечную формулу, обладающую двумя важными свойствами: свойством несмещённости для искомого интеграла I и свойством точности для рассмотренной системы Хаара. Такая формула (см. [1], теорема 1) имеет вид
1 м
= ^ £/(*•), (1)
¿=1
при этом узлы формулы являются случайными величинами, совместная плотность распределения которых задаётся как
( \ _ ) Ж' если (ж!,ж2,... ,ждг) € £а1(«1,г2, • • • ,«дг),
^>(х1 , х2, • • • , хМ) Л
I 0, иначе.
Здесь использовано определение латинского множества, задаваемого перестановкой (¿1, ¿2, .. ., ¿м):
(х1, Х2,..., х*) £ £а{.(«1, ¿2, • • •, ¿м) V? £ {1, 2,..., N} х, € . Далее, в той же работе ([1], теорема 2) показано, что дисперсия формулы равна
DSN[f] = 1[f{x)^dx) - 1 ( [f(xMdx) ] - 1 - aj)\ (2)
i<j
где ai = / f (x)^(dx), i = 1, 2,..., N. Эта дисперсия ни при каких условиях не меньше
JXi
нуля и не больше дисперсии обычного Монте-Карло: 0 ^ DSN[f] ^ Dmc•
Хотя полученная формула обладает рядом важных характеристик, её использование на практике затруднено из-за сложности моделирования распределения с плотностью у>(ж1, Х2,..., xn), особенно в случае произвольного пространства и произвольной схемы его разбиения. Поэтому мы будем использовать полученную формулу в варианте квази Монте-Карло с применением последовательностей Соболя.
В самом деле, пусть здесь и далее X есть s-мерный куб [0,1]s, снабжённый мерой Лебега. Известно, что s-мерная последовательность Соболя для некоторого известного t(s) является (t(s), s ^последовательностью2 по основанию 2. Таким образом,
2 Напомним, последовательность {хп} называется (Ь, з)-последовательностью по основанию Ь, если для произвольных натуральных к ^ 0, т ^ Ь отрезок последовательности {хк-ьт,..., х(к+1)-Ьт — 1} содержит в точности Ь точек внутри произвольного элементарного подмножества вида 01=1 11 (а.,, (¿1 € N0, а,- < Ь^э для всех ]), имеющего объём Ь4—т. Это 1Ь з Ь 3 J
замечательное свойство впервые установлено Соболем [2] и впоследствии обобщено Нидеррейтером [12] на случай произвольного основания.
если строить Xi,X2,..., Xn на основе элементарных подмножеств, мы гарантируем, что для достаточно большого P любой отрезок последовательности Соболя длины N = 2P, имеющий вид {хт n +i, • • •, Х(т +i) N}, принадлежит объединению латинских множеств, то есть в точности тех областей, которые не обращают в ноль плотность ¥>(xi, Х2, . . . ,XN).
Резюмируем выше приведённое рассуждение: мы будем использовать формулу Sn [f ] в сочетании с последовательностью Соболя, а оценивая дисперсию в соответствии с (2), мы получим оценку погрешности интегрирования, аналогичную доверительному интервалу метода Монте-Карло. Таким образом, на каждом шаге алгоритма Qint мы будем получать одновременно оценку интеграла I и оценку погрешности е = |/ - Sn[f]|.
Опишем шаги алгоритма Qint и исследуем вопрос выбора параметров. Пусть дана s-мерная функция f G L2(d^), определённая и принимающая конечные значения в единичном гиперкубе [0,1]s (рассмотрение функций с особенностями остаётся за рамками данной статьи; отметим лишь, что требование конечности значений можно существенно ослабить).
1. Выбираем параметр подразбиений (partition parameter) P G N0.
2. Выбираем параметр повторений (repetition parameter) R G N.
3. Строим разбиение [0,1]s на N = 2P подмножеств равного объёма на основе элементарных подмножеств.
4. Для r = 1, 2,..., R повторяем следующий вычислительный цикл:
(a) берём N следующих3 точек s-мерной последовательности Соболя {x(r-i).N+i,...,xr.N}; вычисляем значения подынтегральной функции в них и дописываем в массив уже имеющихся значений длины rN;
rN
(b) обновляем оценку интеграла: I = ^ f{xi)>
i=i
(c) обновляем оценку стандартного отклонения S в соответствии с (2), где коэффициенты a.j оцениваются по г точкам каждый: a.j = ^ f(xi)',
ieXj
(d) для количества вычислений подынтегральной функции, равного rN, местонахождение истинного значения оценивается при помощи квазидоверительного интервала (I — 3S; I + 3S).
Важно понимать, что чем больше параметр P, тем больше N, и тем реже будет обновляться оценка, поскольку она определена только для {N, 2N, 3N,..., kN,...}. Для промежуточных значений {kN + 1, kN + 2, kN + 3,..., (k + 1)N — 1} мы будем использовать последнюю доступную оценку, полученную на итерации k. Таким образом, с ростом P оценка будет носить всё более ярко выраженный ступенчатый характер.
Отметим отдельно, что оценка интеграла I при помощи выборочного среднего I есть стандартный подход квази Монте-Карло, поэтому требование сходимости оценки к истинному значению не требует никакого формального доказательства: оно выполнено автоматически. При этом основной результат алгоритма Qint состоит в том,
3 Вообще говоря, можно брать N любых подряд идущих не использованных ранее точек; мы пренебрегаем этой возможностью для простоты.
что мы имеем альтернативную оценку погрешности I — I . Отсутствие такого инструмента всегда являлось важным недостатком использования квази Монте-Карло в прикладных задачах; мы надеемся, что распространение и дальнейшее усовершенствование Qint позволит расширить область применения квази Монте-Карло.
3. Параметризация РШ;. 3.1. Двойственность параметров подразбиений и повторений. Зачастую в прикладных вычислениях явный вид функции / неизвестен, а вычисление её значения может занимать значительное время. Поэтому для сравнения различных алгоритмов численного интегрирования часто используется ограничение на количество вычислений функции. Мы будем использовать этот подход для сравнения вариантов Qint с разными наборами параметров.
Пусть мы ограничены М вычислениями подынтегральной функции, причём М = 2т, где т — достаточно большое натуральное число. Такое ограничение несущественно и вводится только для удобства демонстрации связи параметров Р и Д. Если представить М в виде М = Д • 2Р, то допустимые значения для этих параметров исчерпываются парами вида {Р = 0, Д = 2т}, {Р = 1, Д = 2т-1}, {Р = 2, Д = 2т-2},..., {Р = г, Д = 2т-®},..., {Р = т, Д = 1}. Мы уже упоминали, что с ростом Р оценка становится более ступенчатой, но, что куда более важно, дисперсия Бм [/] монотонно убывает по Р.
Теорема 1. Пусть для параметров {Р, Д}, Д ^ 2, есть 'разбиение Х1;Х2,..., Х2р. Если для параметров {Р + 1,Д/2} взять такое разбиение ф1;ф2,...,ф2р+1, что X® = ф2®-1 и ф2® для любого г, то ОБМ{Рд}[/] ^ ЮБМ,{Р+1,Д/2}[/].
Доказательство. Пусть для исходного разбиения определены коэффициенты а® = ^ /(ж)^(^ж), г = 1,2,..., 2Р, участвующие в (2). Для нового разбиения пусть вj = /(х)^(^ж), ] = 1,2,..., 2Р+1. Достаточно доказать только
^га<за («»- — а3с )2 гв<оц (вгв — )2. Этот факт докажем почленно: рассмотрим сначала га = 1, = 2 и соответствующие им гр = 1, 2, = 3,4. Поскольку новое разбиение получено на основе исходного, верны соотношения «1 = в1 + в2, «2 = вз + в4. Если неравенство выполнено, то аналогичное неравенство будет верно и для произвольных га,)а, которыми исчерпывается исходная сумма, и соответствующих им непересекающихся наборов г в, . Итак, для фиксированных индексов имеем следующую цепочку неравенств:
(«1 — «2)2 <
< (в1 — в2)2 + (в1 — вз)2 + (в1 — в4)2 + (в2 — вз)2 + (в2 — в4)2 + (вз — в4)2;
«1 — 2«1«2 + «2 ^
< 3в2 + 3в2 + 3вз2 + 3в4 — 2в1в2 — 2в1вз — 2в1в4 — 2в2вз — 2в2в4 — 2взв4;
2в1в2 + 2взв4 < 2в2 + 2в2 + 2вз2 + 2в4 — 2в1в2 — 2взв4; 0 < 2(в1 — в2)2 +2(вз — в4)2.
Итак, наращивая параметр подразбиений и уменьшая параметр повторений, мы вправе ожидать более точной оценки погрешности. Это довольно естественный результат, потому что вырожденный случай P = 0 соответствует обычному Монте-Карло, в котором вместо псевдослучайных чисел используется последовательность Соболя. С другой стороны, чем больше P, тем меньше Д, и тем потенциально более неустойчивы оценки ai (для них используются только те точки, которые лежат в Xi, которых в точности Д). Практика показывает, что максимально большое P (случай {P = m, Д = 1}) приводит к недостоверной оценке дисперсии, так что мы будем рассматривать только такие P, для которых Д ^ 2.
Подведём итог рассуждениям о двойственности параметров Qint. Для постоянного количества вычислений подынтегральной функции мы можем варьировать соотношение между параметром подразбиений P и параметром повторений Д. В силу монотонности по P чем больше P, тем более точной и ступенчатой будет оценка погрешности. С другой стороны, чем больше Д, тем более оценка погрешности будет стабильной и плавной. Эти эффекты мы увидим в вычислительных примерах, приведённых в заключительной части статьи.
3.2. Параметризация формы разбиений. Алгоритм Qint не описывается однозначно парой параметров (P, Д), поскольку необходимо дополнительно указать правила формирования разбиений Xi,X2,..., для каждого r = 1, 2,... Д. Вообще говоря, для разных r могут быть совершенно разные наборы, но мы будем рассматривать только такие универсальные способы построения, для которых характерна естественная итеративность по P. В частности, в [1] описан один из таких способов, который устроен чрезвычайно просто: разбиение всегда проводится по первой координате гиперкуба.
В этой статье мы используем более нетривиальный способ, который мы назовём правилом кубической формы (cubic shape rule). Мы не будем его подробно описывать; скажем лишь, что гиперкуб разбивается равномерно по всем имеющимся размерностям. Шаги правила кубической формы в трёхмерном случае представлены на рис. 1.
4. Численные примеры. В предыдущей работе [1] рассмотрена функция довольно простого вида. В качестве более продвинутых демонстрационных примеров мы возьмём три функции: одну функцию специального вида, близкую к первой ха-аровской и две из набора TESTPACK (набор функций для тестирования алгоритмов многомерного интегрирования, составленный Аланом Генцем, [7]) —рис. 2-4 соответственно.
4.1. Кусочно-линейная функция. Рассмотрим функцию, задаваемую прямым произведением покоординатных функций, для произвольной размерности i задаваемых как
где с^ = 27^10' Такая функция обладает интересным свойством: в низких размерностях она очень близка к первой функции Хаара, в высоких становится практически линейной /(х) = х. Интегрируя такую функцию при помощи Qint, мы рассчитываем на очень точную оценку дисперсии. Кроме того, мы сможем оценить преимущество, получаемое при помощи правила кубической формы.
Размерность пространства в = 8, истинное значение интеграла I = 1.
x е (0.5 - ci, 0.5 + ci], x > 0.5 + ci,
(d) Шаг 4
(е) Шаг 5 Шаг 6
Рис. 1. Первые шесть шагов правила кубической формы, я = 3.
Модуль ошибки Оценка Монте-Карло
10000 20000 30000 40000 50000 Количество вычислений функции
Рис. 2. Кусочно-линейная функция в размерности 8.
Стоит отметить тот факт, что качественный скачок в точности оценки дисперсии происходит тогда, когда параметр подразбиений становится равен размерности пространства. Это происходит ровно в тот момент, когда правило кубической формы задействует разбиения по всем имеющимся размерностям. Такого скачка не наблюдалось бы с разбиениями только по первой координате.
4-2. Осциллирующая функция из TESTPACK. В качестве второго примера возьмём функцию под кодовым именем «Oscillatory» из набора TESTPACK [7]:
fosc(x) = CCS £ Xi \i=i
Такого рода функции часто встречаются в практических задачах и дают представление о преимуществах использования Qint в тех случаях, когда функция не подобрана специальным образом, как это было в первом примере.
Размерность пространства в = 20, истинное значение интеграла I « -0.362.
Модуль ошибки Оценка Монте-Карло
Количество вычислений функции
Рис. 3. Осциллирующая функция в размерности 20.
а Ш"
В
Н К
3 I
о
« <L> Ю Ö
Is
О
О 10000 20000 30000 40000 50000 Количество вычислений функции
Рис. 4. Функция с пиком в размерности 8.
На этом примере хорошо видны все особенности оценок Qint: растущий вместе с параметром подразбиений выигрыш по сравнению с доверительным интервалом Монте-Карло; всё более выраженная ступенчатость оценок; монотонность по параметру подразбиений.
4-3. Функция с явно выраженным пиком из TESTPACK. Заключительным примером будет функция под кодовым именем «Corner peak» из того же набора TESTPACK:
-s-1
Она интересна тем, что у неё есть область быстрого роста, и значения в этой области значительно выше, чем среднее значение по всему гиперкубу5. Это означает,
5 Так, для рассматриваемого далее случая « = 8 минимум функции составляет 1x10 9, значение
в центре куба примерно 2.2 X 10-7, а пик находится в точке f (0,..., 0) = 1.
что дисперсия, посчитанная методом Монте-Карло, будет большой, а доверительный интервал будет нестабилен. Продемонстрируем, является ли это препятствием для Qint.
Размерность пространства s = 8, истинное значение интеграла I « 2.8 х 10-6. Безусловно, скачки в оценке дисперсии портят не только доверительный интервал Монте-Карло, но и оценки Qint. Однако увеличение параметра подразбиений и связанное с этим усреднение помогает справиться с этой проблемой и получить стабильные оценки, пожертвовав при этом начальным отрезком. Так, для P = 14 мы не сможем построить оценку раньше, чем при N, близком 16000; тем не менее, для остальных N оценка стабильна и значительно более точна по сравнению с обычным Монте-Карло.
Литература
1. Антонов А. А., Ермаков С.М. Эмпирическая оценка погрешности интегрирования методом квази Монте-Карло // Вестн. С.-Петерб. ун-та. Сер. 1. Т. 1(59). 2014. C.3-11.
2. Соболь И. М. О распределении точек в кубе и приближенном вычислении интегралов // Ж. вычисл. матем. и матем. физ. Т. 7(4). 1967. C. 784-802.
3. Соболь И. М. Многомерные квадратурные формулы и функции Хаара. М.: Наука, 1969.
4. Ермаков С. М. Метод Монте-Карло и смежные вопросы. М.: Наука, 1975.
5. Cochran W. G. Sampling Techniques, 3rd Edition. John Wiley, 1977.
6. Dick J., Pillichshammer F. Digital Nets and Sequences. New York: Cambridge University Press,
2010.
7. Genz A. A Package for Testing Multiple Integration Subroutines // Numer. Integr.: Recent Dev., Softw. and Appl. Vol. 203. 1987. P. 337-340.
8. Donald E. Knuth. Seminumerical Algorithms // The Art of Computer Programming. Vol. 2. Addison-Wesley, Reading, Massachusetts, Second, 1981.
9. Kuipers L., Niederreiter H. Uniform Distribution Of Sequences. Wiley. New York. 1974.
10. Lemieux C. Monte Carlo and Quasi-Monte Carlo Sampling. New York: Springer, 2009.
11. Niederreiter H. Error bounds for Quasi-Monte Carlo integration with uniform point sets // J. of Comput. and Appl. Math. 2003. 150. P. 283-292.
12. Niederreiter H. Random Number Generation and Quasi-Monte Carlo Methods. Dover: Dover Publications, 2006.
13. Owen A. B. Scrambled net variance for integrals of smooth functions // Annals of Statistics. 1997. Vol.25, N4. P. 1541-1562.
14. Schiirer R. A Comparison between (Quasi-)Monte Carlo and Cubature Rule Based Methods for Solving High-dimensional Integration Problems // Mathematics and Computers in Simulation. 2003. Vol.62, N3-6. P. 509-517.
Статья поступила в редакцию 23 октября 2014 г.
Сведения об авторах
Антонов Антон Александрович — аспирант; [email protected]
QUASI MONTE CARLO INTEGRATION ALGORITHM WITH A POSTERIORI ERROR ANALYSIS
Anton A. Antonov
St.Petersburg State University, Universitetskaya nab., 7-9, St.Petersburg, 199034, Russian Federation; [email protected]
A possibility of a probabilistic approach to a deterministic error estimation procedure of quasi Monte-Carlo method is presented. Existing estimation methods of the aforementioned estimation are non-constructive.
A previously obtained result, which was deduced from the theory of random cubature formulae, is expanded with new theoretical results. Qint algorithm, which is presented in a convenient for direct
application step-by-step form, is discussed in terms of parametrization. Several key features of Qint are discussed, including monotonicity by a partition parameter and accuracy dependency on a partition rule.
New practical results are presented, which are obtained with the help of Sobol sequences and TESTPACK integrand functions. In all cases a better estimate compared to a traditional Monte Carlo confidence interval is demonstrated. Refs 14. Figs 4.
Keywords: Monte Carlo and quasi Monte Carlo method, confidence interval, random cubature formulae, Haar functions, Sobol sequences.
References
1. Antonov A. A., Ermakov S. M., "Empirically estimating error of integration by Quasi-Monte Carlo method", Vestnik St.Petersburg University. Mathematics 47(1), 1—8 (2014).
2. Sobol' I. M., "Distribution of points in a cube and approximate evaluation of integrals", Zh. Vych. Mat. Mat. Fiz. 7(4), 784-802 (1967).
3. Sobol' I.M., Multidimens. Quadrature Formulas and Haar Funct. (Nauka, Moscow, 1969) [in Russian].
4. Ermakov S. M., Die Monte Carlo Methode und Verwandte Fragen (VEB Deutscher Verlag der Wissenschaften, Berlin, 1975).
5. Cochran W. G., Sampling Techniques, 3rd ed. (John Wiley, 1977).
6. Dick J., Pillichshammer F., Digital Nets and Sequences (Cambridge University Press, New York, 2010).
7. Genz A., "A Package for Testing Multiple Integration Subroutines", Numer. Integr.: Recent Dev., Softw. and Appl. 203, 337-340 (1987).
8. Donald E. Knuth., "Seminumerical Algorithms", The Art of Computer Programming 2 (Addison-Wesley, Reading, Massachusetts, Second, 1981).
9. Kuipers L., Niederreiter H., "Uniform Distribution Of Sequences" (Wiley, New York, 1974).
10. Lemieux C., "Monte Carlo and Quasi-Monte Carlo Sampling" (Springer, New York, 2009).
11. Niederreiter H., "Error bounds for Quasi-Monte Carlo integration with uniform point sets", J. of Comput. and Appl. Math. 150, 283-292 (2003).
12. Niederreiter H., "Random Number Generation and Quasi-Monte Carlo Methods" (Dover Publications, Dover, 2006).
13. Owen A. B., "Scrambled net variance for integrals of smooth functions", Annals of Statistics. 25(4), 1541-1562 (1997).
14. Schürer R., "A Comparison between (Quasi-)Monte Carlo and Cubature Rule Based Methods for Solving High-dimensional Integration Problems", Mathematics and Computers in Simulation. 62(3-6), 509-517 (2003).
