Адаптивный метод существенной выборки при ограниченном числе шагов бисекционного процесса Текст научной статьи по специальности «Математика»

12. Жгутов, В.М. Математическая модель, алгоритм исследования и анализ устойчивости нелинейно-упругих ребристых оболочек при больших перемещениях [Текст] / В.М. Жгутов //

Научно-технические ведомости СПбГПУ. Сер. «Физико-математические науки» - 2009. - № 4. -С. 24 - 30.

УДK 519.6

Д.Г. Арсеньев, В.М. Иванов, Н.А. Берковский

АДАПТИВНЫМ МЕТОД СУЩЕСТВЕННОЙ ВЫБОРКИ ПРИ ОГРАНИЧЕННОМ ЧИСЛЕ ШАГОВ БИСЕКЦИОННОГО ПРОЦЕССА

Адаптивный метод существенной выборки (АСВ) является одним из проявлений идеи последовательных методов Монте-Карло для вычисления интегралов; впервые он предложен в работе [4]. Наиболее полно теория этого метода развита в монографии [1]. Там представлены оценки сходимости, гарантирующие при определенных условиях большую эффективность адаптивных методов по сравнению с классическими методами Монте-Карло, такими как метод существенной выборки и метод выделения главной части. Однако теория предполагает возможность сколь угодно точных кусочных аппроксимаций подынтегральной функции, что не может быть реализовано на практике. Особенно в случае интегралов высокой кратности в связи с ограничениями объема памяти вычислительного устройства вычислителя лимитирует число подобластей, на которые следует разбивать исходную область интегрирования.

В данной статье исследуются особенности применения АСВ при условии, когда имеются ограничения в мелкости дробления области интегрирования. Разработан алгоритм, который позволяет эффективно применять адаптивную схему при малом числе шагов процесса бисекции. Возможности алгоритма показаны на примере методической одномерной задачи теории фильтрации и сравниваются с результатами, полученными методом Монте-Карло с гауссовой плотностью распределения, а также методом существенной выборки (СВ). Детально разбирается численный пример, в котором адаптивная схема оказывается эффективнее, чем метод СВ. Заметим, что согласно литературным данным адаптация в стохастическом интегрировании может быть основана и на других идеях [6].

Схема АСВ в одномерном случае

Общая схема АСВ для интегралов любой кратности приведена в работе [1]. Здесь для простоты рассмотрим ее одномерный вариант, предполагая интеграл несобственным с бесконечными пределами. Допустим, что необходимо приближенно вычислить I = J/(х)сЬс, где о функции Дх) известно

следующее:

она интегрируема на [—оо 5оо];

Дх) Ф 0 почти везде на всей числовой оси;

значенияДх) вне фиксированного промежутка [а, Ь] пренебрежимо малы по сравнению с ее значениями внутри него, так что вклад этих значений в интеграл несуществен.

Последнему требованию можно придать более точную форму, но в этом нет необходимости. Кроме того, будем считать, что Дх) такова, что дисперсия оценки Дх)/р(х) конечна при всех плотностях р(х), которые фигурируют в статье. Рассмотрим следующую схему вычислений.

1. Выберем произвольную плотность распределения с1(х) такую, что с1(х) > 0 (строго) на всей числовой прямой.

2. Разделим [а, Ь] на две равных части [а, с] и [с, Ь], возьмем малое число 5 > 0 и найдем значения

+ 8, /2 =

Л =

/

/

+ 8.

Обозначим промежутки [а, с] и [с, Ь] как А! и Д2 соответственно и рассмотрим функцию

' 2

к=1

ш=

где Ха1 (*) обозначает характеристическую функцию промежутка Ак. Далее, введем функцию

а(0

_ Ш

оо

¡т

= /(%) Ры{хм)

N

+1)

'ж

=Хк)2

ЛГ

¿=1

промежутки Д1, и Д^ вместо А4 и получаем новое разбиение

Обозначим за !~к середины промежутков Аь положим = [/(^ I и составим функцию

которая есть первое приближение оптимальной плотности распределения. Этот пункт служит инициализацией следующего итеративного процесса (см. далее п.п. 3, 4)

3. На шаге с номером N генерируем случайную точку хм с плотностью распределения рм (х). Вычисляем первичную оценку

/лм(*)=

N+2 к=1

<1{х), х<£ [а,6],

затем введем

По вычисленным ранее первичным оценкам 51!, Б2,..., находим вторичную оценку

к=1

где а^ - весовые коэффициенты, которые в случае кусочно-постоянной аппроксимации выражаются как

Вычисляется асимптотически несмещенная оценка дисперсии

При больших значениях N вторичная оценка 1Н ~ N{1,0^ в силу центральной предельной теоремы и утверждений, доказанных в работе [1]; при этом - 0(1/А^) при дифференцируемой подынтегральной функции [1]. Если а) N>N0 и б) Одт < е/2 (или стж < е/3), где е - точность вычислений, то процесс останавливается; при этом считаем с требуемой точностью. Если же

конъюнкция а) и б) не имеет места, то переходим к п. 4.

4. Делим промежуток А^, в который попала точка хн, на две равных части Д^ и Д^ (бисекция). Вставляем в уже существующее разбиение

хы ={Д,, Д^.-^Адад}

и перейдем к п. 3 на шаге N + 1.

Поясним роль 5, Ы0, и С(х). Параметр 5 следует вводить для подынтегральных функций, которые могут быть очень близки к нулю на промежутке [а, Ь], иначе метод не будет адекватно работать.

- достаточно большое число (например 500) -вводится для того, чтобы исключить неадекватные оценки на первых шагах метода. Функция С(х) нужна, чтобы определить плотность р1^{х) вне промежутка[а, Ь], на котором ведется адаптивный бисекционный процесс.

Схема АСВ при ограниченном числе шагов бисекционного процесса

Как уже отмечалось, вычислитель всегда ограничен в возможности увеличивать мощность разбиения т^ ={Д15 А2, Алг+1}. Будем предполагать, что мы можем разбить [а, Ь] не более чем на N промежутков. Как показывает вычислительная практика, в некоторых задачах, если действовать согласно схеме, представленной в предыдущем разделе (см. пп. 1-4), и провести N итераций, то иногда удается добиться точного результата, а иногда нет. Заметим, что если дисперсия оценки 1Ы =1н(рх,х2,..., хм) мала, то можно перейти к оценке

1 м

ш к=1

к

N,

где - к-ая реализация /ЛГ=/ЛГ(х1 ,х2,..., %), а М - число реализаций

Но в случае большой дисперсии

этот способ не будет эффективен. Обозначим за ■Рлт(0 - случайную плотность, полученную при к-ой реализации 1^. Рассмотрим для каждой плотности случайные величины

1 1 ^42

Ж)

которые являются стандартными оценками

=

и требуется меньше обращений к подынтегральной функции, чем при стандартном методе Монте-Карло и методе СВ. Алгоритм применения АСВ при ограниченном числе N шагов бисекционного процесса может быть следующим.

1. Генерируется небольшое число (3-5) случайных плотностей Рх (лс), и по малой выборке

(200-300 точек) оцениваются величины СУ-! I.

Запоминается функция плотности р™е й(л:), которая соответствует медиане

(

Меё - тес!

\ж\ м4

методом Монте-Карло на плотности (дс). В зависимости от вида Рх{х) их дисперсия будет варьироваться от реализации к реализации. Кроме того, по-разному себя будет вести выборочное

I /(01 и -

стандартное отклонение СГ-{ —тгт^г. Чрезвычайно

неблагоприятно для вычислений, когда оно совершает относительно большие скачки на итерациях с большими номерами, препятствуя монотонному убыванию стандартных ошибок О^/*}. Это может сильно увеличивать время работы программы. Требования к оптимальной плотности р^'(х), выбираемой из случайных плотностей р^ (¡с), можно сформулировать как следующие:

случайная величина с) должна

иметь минимальную дисперсию среди прочих плотностей р^ (х);

оценка о{/*} должна монотонно убывать с ростом Ь.

Практически это означает, что если с какой-то (не очень малой) вероятностью при определенном N генерируется эффективная для вычислений плотность (х), то можно попытаться «поймать» эту плотность и провести вычисления на ней стандартным методом Монте-Карло. Метод отслеживания плотностей может быть разным, и возможно, что это отдельная тема исследований.

В данной работе предлагается производить отбор статистическими методами, которые на первый взгляд занимают много времени. Однако, как это показано в разделе численных экспериментов, в отдельных ситуациях статистический подход оправдывает себя; именно при его использовании значительно быстрее вычисляется интеграл

Плотность

/(О

минимальной

величиной

[К (01

брать опасно, так как эта малость мо-

жет быть неадекватной оценкой вследствие малости выборки. Вычисляется верхний порог тах для

7(0'

последующих значений С<

равный этой

[р№)\

медиане (или 120-150 % от нее для страховки).

2. Начинаются вычисления стандартным методом Монте-Карло с плотностью /?™ей(0, т. е. вычисляется оценка

" /(о ] 1/1а " /(0 1

1/ Г"

и ее выборочное стандартное отклонение <т{у™е(1}. На каждом шаге, начиная с некоторого Ь0 (300-500) вычисляется отношение

Если Я < 0$, то переходим к п. 3; Я0 - контрольное число ((< ^ < 1,5). Кроме того, каждые 100 итераций вычисляется

если < 0, переходим к п. 3.

Этот двойной контроль отменяется, если требуемая точность уже почти достигнута, например, если

о{/Т1}<»яе/2,

где 1 <т < 1,5.

с

Если на Ь-й итерации

а{/Гй}<е/2,

т-шсс! т

то вычисления заканчиваются и и £ ~ 1.

3. Генерируются новые плотности р^(х) до первой плотности р^™ (дс), для которой

Ах) ]

а , Л<тах.

Величины а-

оцениваются по неболь-

шой выборке (200-300) точек. Далее для плотности Рх™(х\ прошедшей отбор, повторяются действия по п. 2 с плотностью р^(х) вместо ^™е<1(л;) и т. д.

Характерные черты и возможности схемы АСВ в условиях фиксированного числа шагов процесса бисекции

Следует заметить, что процедура, предложенная в предыдущем разделе, не обязана заканчиваться. Программа может перебирать плотность за плотностью. Эта проблема решается с помощью установки ограничений по времени, вывода на монитор различных сообщений о стадиях работы программы и т. п.

Предложенная схема легко комбинируется с методом существенной выборки, так как плотность, вычисляемая в соответствии с последним, может пройти такую же селекцию, как и плотности р^(х) (см. предыдущий раздел).

От стандартной схемы АСВ при неограниченном числе N (см. раздел «Схема АСВ в одномерном случае») остается, по существу, только метод подгонки плотности, главная черта которого - та, что промежутки {Дх, Д2,..., Дж} автоматически сгущаются в областях наибольших значений модуля подынтегральной функции; а это дает преимущества при вычислении интегралов от функций, имеющих резкие изолированные максимумы N

аппроксимациями с использованием, например, формулы Тейлора.

Алгоритм выгодно применять в случае функций с резкими, почти сингулярными максимумами, расположение которых в области определения известно с точностью до достаточно крупной подобласти, содержащей экстремум. Примерами таких функций служат подынтегральные функции задачи оптимального оценивания, рассмотренной далее.

Параметры 5, Л^0, Ь0, ^ и вспомогательную плотность с1(х) вычислитель выбирает по своему усмотрению. Для неопытного пользователя или в том случае, если программа вызывается автоматически как подпрограмма, их можно установить по умолчанию.

Постановка задачи оптимального оценивания

Рассмотрим следующую одномерную задачу оптимального оценивания.

Пусть случайная величина X распределена нормально с плотностью распределения р{х)-Ы{х,х,о), а результат измерения У при X=х распределен с плотностью р(у | х)= Ы(у,х,г). Допустим, что измерение положения X дало результат у. При этом условии требуется найти оптимальную оценку для значения X. Эта оценка находится по формуле

I.

^хр(у | х)р(х)с1х |/[ (х)с1х

\р(у I х)р(х)сЬс |/2 {х)сЬс

-. (1)

Известно точное решение задачи (1), а именно

2 2

(2)

Допустим, что интегралы в формуле (1) вычисляются каким-либо приближенным методом. Пусть - относительная погрешность вычисления 1-, а 82 - вычисления /2. Исходя из общей формулы абсолютной погрешности частного [5] для

[2]. Оценки 1Ы = не вычисляются, хотя, \°»(у)= ^ЬУ ^Щ имеем оценку:

к=1

вероятно, можно составить алгоритм, аналогичный предложенному в предыдущем разделе, где они будут эффективно использоваться.

Кусочно-постоянные аппроксимации в процессе бисекции можно заменить более точными

д^^ООв+бД

(3)

Пусть 1Ъ /2, хор*{у)= 1Х/12 - это приближенные значения величин 1Ь /2, хор>(уУ Учитываем, что

8,=

5,=

Ь ±2

и заменяем в оценке (3) точные значения приближенными; тогда имеем следующую оценку для теоретической погрешности:

К1-Л

. (4)

Подразумевается, что абсолютные погрешности —и |/2 — /2| можно эффективно оценить в процессе вычислений. Заметим, что при использовании стохастических методов оценка (4) носит вероятностный характер. Допустим, что

- А | < е) = Р1 и Р(|/2 -12| < е) = р2.

Тогда при независимых выборках, используемых для числителя и знаменателя, с вероятностью Р1Р2 будет верна оценка

7 ( - . ^

А

1 1

\

■ /

р(|/1-71|<е||/2-72|<е)

■ 1,

функций числителя и знаменателя при различных значениях г представлены на рис. 1. Видно, что с уменьшением значения г обе подынтегральные функции становятся все более похожими на дельта-функцию, и соответственно, становятся проблемой для численного интегрирования. В данной задаче оба интеграла вычисляются аналитически (точные значения см. в табл. 1).

Таблица 1

Зависимость точных значений оценки и интегралов формулы (1) от значения г

г х°»(у)

1 1,05000 0,29546 0,2813904

0,1 1,099010 0,43411 0,3950021

0,01 1,0999900 0,43662 0,3969524

0,001 1,09999990 0,43655 0,3969524

Однако, как показала вычислительная практика, использовать независимые выборки неэффективно. Лучше использовать одну выборку для числителя и знаменателя, при этом оказывается, что погрешности ^ — и |/2 — /2| коррелируют таким образом, что для достаточно малых е условная вероятность

следовательно, формулой (4) можно пользоваться с той же степенью достоверности, что и оценками для числителя и знаменателя по методу Монте-Карло. К сожалению, оценка (4) оказывается сильно завышенной (см. далее), и альтернативы ей, кроме различных эмпирических подходов, пока не разработано.

Обозначения: 1Х = | (х)ск, /2 = | /2 (х)с1х

Вычисления проводились в системе МЛТЬЛБ, на тактовой частоте 3 ГГц. Рассмотрено три метода, обозначенные как МК Гаусс, СВ и АСВ.

МК Гаусс - это стандартный метод Монте-Карло, где плотность распределения вбрасываемых точек суть р(х) = Щ(х, 1, 1).

СВ - это метод существенной выборки, в котором I — J/(х)с1х с плотностью

где

Х(Л+8)Ха». хе 4];

к=1

АГ(х,1,1), хе[-2,4];

Раздел численных экспериментов

Основные результаты. Задача оптимального оценивания, сформулированная в предыдущем разделе, решалась при случайной величине X, распределенной нормально с плотностью р(х) = Щ(х, 1, 1), при условии, что известен результат измерения у = 1,1. Графики подынтегральных

при этом

д* =

„ 6(к — 1) „ 6 к

-2 + —--,-2 + —

n n

а /к - значения модуля функции в серединах Ак. Таким образом, имеем кусочно-постоянную ап-роксимацию на промежутке [-2, 4], где равны

промежутки постоянства Ак. Промежуток [-2, 4] выбран из того соображения, что по правилу «трех сигм» случайная величина X почти достоверно в нем находится. Как и прежде, величина 5 > 0 введена во избежание больших дисперсий при подынтегральных функциях, близких к нулю на промежутке [-2, 4] (в данных примерах это имеет место).

АСВ - адаптивный метод существенной выборки, применяемый согласно схеме, описанной во втором разделе данной статьи, причем ё(х) = Щ(х, 1, 1). И в АСВ, и в СВ число N = 1000 (ограничение на мощность разбиения), 5 = 0,005; это обеспечивало корректное сравнение методов. Представляло интерес выяснить, в каких случаях переход от метода СВ к более алгоритмически

сложному АСВ оправдан в отношении трудоемкости вычислений. В табл. 2 и 3 приведены результаты численных экспериментов. Ввиду того, что на месте простых интеграндов могли оказаться более сложные функции, вычисление которых заняло бы больше времени, чем вычисление функций данной задачи, сравнение методов проведено по времени и числу обращений к подынтегральным функциям. В табл. 3 представлено сравнение параметров выигрыша. При этом в столбце «Параметр выигрыша по /» записаны отношения числа обращений к интеграндам методом МК Гаусс к числу обращений рассматриваемым методом (СВ или АСВ). Аналогичные показатели для времени содержатся в столбце «Параметр выигрыша по 1».

Рис. 1. Графики подынтегральных функций ^ (х) (1) и /2 (х) (2) (интеграндов числителя и знаменателя соответственно) при различных значениях среднеквадратичного отклонения г: 1 (а); 0,1 (б); 0,01 (в); 0,001 (г)

(Следует обратить внимание на разные масштабы по обеим осям)

Таблица 2

Результаты численных экспериментов, проведенных тремя методами

Значение параметра для различных значений г

Параметр 1 0,1 0,01 0,001

сравнения МК Гаусс СВ АСВ МК Гаусс СВ АСВ МК Гаусс СВ АСВ МК Гаусс СВ АСВ

Абсолютная 1,66 10-3 1,0410-3 1,03 10-3 0,000 1,6-Ю-4 1-10-4 110-5 110-5 3-10-5 0,000 0,000 0,000

погрешность

Число обращений к функциям, тыс. 378 86 90 9337 116 123 106957 168 171 1081095 110342 305

Время, с 9 48 45 208 35 35 2738 47 48 24252 30281 92

Теоретическая 7,3-10-3 7,3-10-3 7,3-10-3 5,3 10-3 4,9-10-3 5,3 10-3 5,3 10-3 5,3 10-3 5,3 10-3 5,3 10-3 5,3 10-3 5,310-3

погрешность

Обозначения: МК Гаусс - метод Монте-Карло с гауссовой плотностью распределения, СВ и АСВ - традиционный и адаптивный методы существенной выборки; г - среднеквадратичное отклонение измерений

Таблица 3

Сравнение параметров выигрыша в зависимости от среднеквадратичного отклонения г для трех вычислительных методов

г Значение параметра выигрыша

по / по 1

МК Гаусс СВ АСВ МК Гаусс СВ АСВ

1 1,0 4,4 4,2 1,0 0,2 0,2

0,1 80,3 75,5 5,9 5,9

0,01 633,3 623,4 50,6 49,5

0,001 9,8 3533,6 0,8 263,6

Примечание. Параметр выигрыша - это число обращений к функции или время счета методом МК Гаусс, деленное на соответствующий показатель рассматриваемого метода

Счет по всем методам (см. табл. 2 и 3) прекращался по достижении точности 0,001 для обоих интегралов. Отметим, что в методах СВ и АСВ при г = 1 аппроксимировались плотности и в числителе, и в знаменателе. При остальных значениях г аппроксимировалась только оптимальная плотность знаменателя ввиду подобия интеграндов. Это позволило сэкономить время, хотя и незначительно. Для метода АСВ взяты средние значения по 100 просчетам.

Из данных табл. 3 видно, что только при г = 1 стандартный метод МК Гаусс имеет преимущество перед методом АСВ во времени за счет простоты программной реализации. При уменьшении значения г это преимущество теряется; кроме того, при этом резко увеличивается преимущест-

во метода АСВ перед методом МК Гаусс как по времени, так и по малости числа обращений к функции. Более сложная ситуация возникает при сравнении методов АСВ и СВ: когда г > 0,001, результаты примерно одинаковы, а при г = 0,001 АСВ оказывается значительно эффективнее. Объяснение этого результата дано в следующем подразделе. Наконец, отметим, что значения теоретической погрешности (см. табл. 2), полученные по формуле (4), оказываются сильно завышенными, так что оценка (4) в данном случае неэффективна. Поиск формулы, которая бы без избытка оценивала погрешность частного в данной задаче, остается актуальной проблемой, очевидного решения которой пока не найдено. Данные табл. 3 убедительно показывают преимущества метода АСВ в случае малых значений величины г.

Анализ результатов численных экспериментов. По итогам вычислительной практики была получена следующая зависимость среднего количества плотностей распределения, перебираемых методом АСВ в процессе решения задачи, от значения среднеквадратичного отклонения г:

г Количество плотностей

>0,01.....................................1

0,01.......................................3

0,001....................................30

Другими словами, при ограничении мощности разбиения значением N = 1000 при г > 0,01 достоверно генерируется «удачная» плотность распределения, при г = 0,01 она генерируется с вероятностью 1/3, а при г = 0,001 - такая вероятность мала и составляет всего 1/30.

На рис. 2. показаны гистограммы двух плотностей, сгенерированных методом АСВ при г = 0,01: одна - «неудачная», которая удаляется согласно алгоритму отбора, другая - «удачная», благоприятная для вычислений. Кроме того, на рис. 2 показана гистограмма плотности, полученной по методу СВ.

Заметим, что на рис. 2 в качестве «неудачной» плотности для наглядности показан очень неблагоприятный случай, который встречается с малой вероятностью, однако критерий отбора, заложенный в АСВ, удаляет плотности даже более похожие на нормальное распределение, т. е. отбор ведется весьма строго. Видно, что плотность, полученная методом СВ (рис. 2, в), менее точно подогнана под интегрируемую функцию, чем «удач -ная» плотность АСВ (рис. 2, б). Однако, за счет

равномерного измельчения промежутка метод СВ достаточно точно повторяет поведение интегран-да, в частности, симметрию, поэтому метод СВ эффективен при значении г = 0,01. Иная картина наблюдается при г = 0,001 (рис. 3). Видим, что из-за недостаточно мелкого разбиения плотность, построенная в соответствии с методом СВ (рис. 3, б), не повторяет важнейшую особенность поведения интегранда, а именно резкий максимум при х = 1,10, из-за чего метод СВ становится неэффективным при г = 0,001. На рис. 3, а показана гистограмма одной из немногочисленных плотностей, прошедших отбор в методе АСВ; видно, что она воспроизводит поведение подынтегральной функции гораздо точнее. Аналогичная ситуация наблюдалась в работе [2] для многомерного случая. На процедуру отбора тратится немалое

а)

а н

2

б)

1,07

1,09

1,10 в)

200

150

1,11

1,12

1,07

1,09

1,10

1,11

1,12

л н

2 =т

100

50

1,07

1,08

1,09

1,10 1,11

1,12

Рис. 2. Гистограммы плотностей распределения, полученных при г = 0,01 методами

АСВ (а, б) и СВ (в);

а, б - «неудачное» и «удачное» распределения плотности

х

х

х

а) 30

25

20

л н

8 15

ö =т

10

■ I.....I

II.....i_I

0

1,090 1,095 1,100 1,105

л н 2 У =т

б) 20

15 10

0

1,090

1,095

1,100

1,105

Рис. 3. Гистограммы плотностей распределения, полученных при г = 0,001 методами АСВ (а) и СВ (б); а, б - «удачное» и «неудачное» распределения плотности

машинное время, однако в среднем это время оказывается значительно меньшим, чем затраченное на счет для одной плотности методами МК Гаусс или СВ.

Итак, нами предложен практический алгоритм численного стохастического интегрирования, основанный на идее адаптации плотности распределения. По сравнению с основными предшествующими источниками акцент был сделан не на скорости сходимости адаптивного метода

в бесконечности, а на возможностях адаптивного подхода в естественных условиях ограниченной мощности разбиения области интегрирования.

Этот алгоритм следует применять в сложных задачах, не поддающихся быстрому решению при использовании стандартной схемы метода существенной выборки; последняя может входить в данный алгоритм как частный случай. На примере задачи оптимального оценивания показаны преимущества предложенной схемы и объяснен механизм возникновения этого преимущества.

5

5

х

x

СПИСОК ЛИТЕРАТУРЫ

1. Арсеньев, Д.Г. Адаптивное управление в стохастических методах вычислительной математики и механики [Текст] / Д.Г. Арсеньев, В.М. Иванов, М.Л. Кореневский. - 2-е изд., испр. и доп. - СПб.: Наука, 2008 - 423 с.

2. Арсеньев, Д.Г. Анализ эффективности адаптивного метода существенной выборки [Текст] / Д.Г. Арсеньев, В.М. Иванов, Н.А. Берковский // Научно-технические ведомости СПбГПУ. - 2009. -№ 4 (88). - С. 43-50.

4. Halton, J.H. Sequential Monte Carlo [Text] / J.H. Halton// Proc. Cambridge Philos. Soc. - 1962. -Vol. 58. - № 1. - P. 57-78.

5. Демидович, Б.П. Основы вычислительной математики [Текст]: Учеб. пос. для студ. / Б.П. Демидович, И. А. Марон. - СПб.: «Лань», 2009. - 672 с.

6. Ohl, T. Vegas revisited: Adaptive Monte Carlo integration beyond factorization [Text] / T. Ohl // Comput. Phys. Commun. - 1999. - Vol.120. - № 1. -P. 13-19.