Диссертация (1150625), страница 4
Текст из файла (страница 4)
Так, дляувеличения точности интегрирования на один порядок необходимо увеличить количество вычислений подынтегральной функции на два порядка, что далеко не всегда приемлемо. Существует целый спектр методов понижения дисперсии (обзор таких приёмов дан у Ермакова С.М. [6], Кохрана В.Г. [7]), каждый из которых ставит своей целью ускорить убывание остатка. Одним из наиболее известных таких методов является расслоение (расслоенная выборка,stratified sampling).Рассмотрим разбиение на непересекающихся областей 1 , . .
. , , покрывающих весь гиперкуб . Для каждой из этих областей проведём стандартное интегрирование методом МонтеКарло с использованием 1 , . . . , точек соответственно, так что общее количество точек попрежнему равно : 1 + . . . + = . Обозначим точки, принадлежащие , за 1 , . . . , , и15тогда общая оценка интеграла будет выглядеть как=∑︁ ( ) ∑︁=1 ( ).(1.9)=1В этом случае дисперсия оценки из 2 / превращается вVar( )=∑︁2 ( )=12,,(1.10)где⎞2⎛2,=1 ( )∫︁ 2 () −∫︁12 ( )⎜⎝⎟ ()⎠ .(1.11)Степень уменьшения получившейся дисперсии (1.10) зависит от разбиения и количества точек1 , .
. . , в каждом из 1 , . . . , .Простейший вариант расслоения мы получим, если объёмы областей 1 , . . . , равны междусобой, и в каждой из них берётся одинаковое количество точек: 1 = . . . = . Обозначим этоколичество за 0 , тогда для всех имеем ( ) = 1/ и(,0)01 ∑︁ ∑︁= ( ).0 =1 =1(1.12)Как мы можем заметить, оценка интеграла в точности та же, что и для наивного МонтеКарло (1.2). Дисперсия также примет простой вид:Var((,)=0)1 ∑︁ 2 ,2 0 =1 ,где⎞2⎛2,=∫︁⎜ 2 () − 2 ⎝(1.13)∫︁⎟ ()⎠ .(1.14)2Как и для наивного Монте-Карло, на практике нужно строить оценки величин ,. Это делаетсяаналогично оценке (1.7):2ˆ,=⎛⎞21 ∑︁ ⎝1 ∑︁ ( ) − ( )⎠ , − 1 =1 =122при этом, естественно, E(ˆ,) = ,.(1.15)16В упрощённом варианте расслоения для практического применения имеем несмещённуюоценку дисперсии (1.13):̂︂ ) =Var((,0 )2 10 (0 − 1) ∑︁∑︁=1 =1⎞2∑︁⎝ ( ) − 1 ( )⎠ .
=1⎛(1.16)В случае, если разбиение 1 , . . . , выбрано удачно, дисперсия расслоенной выборки (1.10)(или (1.13) в упрощённом случае) будет значительно меньше, чем дисперсия (1.5) наивного MC.В противном случае, однако, дисперсия может даже увеличиться. Гарантированное уменьшениедисперсии достигается в некоторых частных случаях (наиболее известен случай, при которомдоля точек в каждой области пропорциональна её объёму, см.
Ермаков С.М. [6], Гельфанд И.М.,Фролов А.С. и Ченцов Н.Н. [25]). Важно также отметить, что уменьшение дисперсии можетбыть двояким: речь может идти как об улучшении асимптотики убывания дисперсии, так и обуменьшении постоянной при сохранении асимптотики наивного Монте-Карло.1.000.750.500.250.000.000.250.500.751.00Рисунок 1.1: Расслоенная выборкаНа рисунке 1.1 представлена расслоенная выборка из 32 точек.
Расслоение проводится пошестнадцати квадратам равного объёма, в каждый из которых попадает ровно по две точки.171.1.4Случайные квадратурные формулыИспользование случайных квадратурных формул является одним из универсальных приёмов уменьшения дисперсии. Эта техника разработана, в первую очередь, Ермаковым С.М. иЗолотухиным В.Г. [8], Ермаковым С.М.
[9, 10, 13], Ермаковым С.М. и Грановским Б.Л. [12] с дополнениями Хэндскомба Д. [11]. Теория случайных квадратур обладает мощным вычислительным потенциалом и гибкостью, особенно в тех случаях, когда известна некоторая информацияо подынтегральной функции (например, какая ортонормированная система может являться еёхорошим приближением). Мы будем использовать аппарат случайных квадратур для полученияформулы с определёнными необходимыми нам свойствами, которая ляжет в основу предлагаемого метода Qint. Приведём основные результаты из книги Ермакова С.М.
[6] в несколькоупрощённой форме для удобства изложения. В частности, мы рассматриваем только частнуюпостановку задачи численного интегрирования.Случайной квадратурной формулой мы будем называть формулу вида=∑︁∫︁ ( ) ≈=1 (),(1.17)где коэффициенты = (1 , . . . , ) суть заданные функции случайных величин .Далее мы будем использовать следующие обозначения.
Пусть { }=1 – произвольное семейство функций, а = (1 , . . . , ) – набор точек в . В таком случае ∆() обозначаетопределитель альтернантной матрицы⃒⃒⃒ 1 (1 ) 2 (1 )⃒⃒⃒ 1 (2 ) 2 (2 )∆() = ⃒⃒ ...⃒ ...⃒⃒ ( ) ( )⃒ 1 2............⃒⃒ (1 ) ⃒⃒⃒ (2 ) ⃒⃒.. ⃒ ,. ⃒⃒ ( )⃒⃒(1.18)а ∆ () – определитель похожей матрицы, полученной заменой функции 1 на :⃒⃒⃒ (1 ) 2 (1 )⃒⃒⃒ (2 ) 2 (2 )∆ () = ⃒⃒ ...⃒ ...⃒⃒ ( ) ( )⃒2............⃒⃒ (1 ) ⃒⃒⃒ (2 ) ⃒⃒.. ⃒ .. ⃒⃒ ( )⃒⃒(1.19)Мы будем называть систему { }=1 регулярной, если ( : ∆() = 0).Следующий результат является центральным для теории случайных квадратур (теорема 4.1,Ермаков С.М., [6]).Теорема 1.
Пусть система { }=1 состоит из линейно независимых для почти всех и ортогональных функций, причём первая из них тождественно равна единице: 1 ≡ 1. Рассмотрим18случайную квадратурную формулу, задаваемую =∆ ().∆()(1.20)Возьмём в качестве плотности совместного распределения узлов 1 , . . . , функцию(1 , . . . , ) =∆2 ().!(1.21)В этом случае является несмещённой оценкой интеграла .Дополнительно отметим, что квадратурная сумма в этом случае точна для всего семейства { }=1 .
В самом деле, для функции 1 сумма постоянна и равна единице, значениеинтеграла также равно единице. Для любой из 2 , . . . , в силу ортогональности значение интеграла равно нулю, в ноль обращается и квадратурная сумма, поскольку определитель ∆ ()содержит два одинаковых столбца.Приведём также следующий важный результат о дисперсии формулы, рассматриваемой втеореме 1 (теорема 4.2, Ермаков С.М., [6]).Теорема 2. Если система { }=1 является ортонормированной, то в предположениях теоремы 1 справедливо неравенство⎛Var( ) 6∫︁⎝ () −∑︁⎞2(, ) ()⎠ .(1.22)=1Если при этом система { }=1 регулярна, то (1.22) обращается в строгое равенство.Теорема 1 даёт возможность строить такие случайные квадратуры, которые точны для наперёд заданного семейства функций, хорошо приближающих подынтегральную функцию . Чемлучше такое приближение, тем меньше будет дисперсия формулы, что гарантируется теоремой 2.1.1.5Метод квази-Монте-КарлоМетод квази-Монте-Карло2 (quasi-Monte Carlo, QMC) получил широкое распространение втрудных с точки зрения вычислительной сложности задачах.
Такого рода задачи возникают, например, в стохастической финансовой математике или при рендеринге трёхмерных сцен со сложными моделями освещённости. Использование метода Монте-Карло даже с учётом понижениядисперсииможет оказаться неприемлемым в силу достаточно медленной сходимости порядка(︁ )︁1 √ . Основное преимущество квази-Монте-Карло заключается в том, что он имеет сходи(︀ )︀мость, близкую к 1 .2В русскоязычной литературе нет устоявшегося написания метода; мы будем использовать вариант, наиболеесоответствующий нормам русского языка (в частности, именно такое написание рекомендовано справочником Розенталя Д.Э.
[26]).19Оценкой искомой величины при использовании квази-Монте-Карло является среднее значение функции для некоторой детерминированной последовательности 1 , . . . , : 1 ∑︁= ( ), =1(1.23)что в точности совпадает с оценкой (1.2) наивного Монте-Карло с той лишь разницей, что (1.23)не содержит случайных величин. В связи с этим набор 1 , . . . , часто называют квазислучайнойпоследовательностью.Одним из самых значимых результатов для теории квази-Монте-Карло является неравенствоКоксмы-Хлавки (Koksma-Hlawka inequality), связывающее ошибку интегрирования с величинойдискрепанса3 (discrepancy).Определение 1.
Дискрепансом набора {1 , . . . , } называется величина⃒⃒⃒ #{ : ∈ }⃒⃒ (1 , . . . , ) = sup ⃒− ()⃒⃒ ,∈(1.24)где #{ : ∈ } обозначает количество точек из набора, попавших в множество , а определяется как объединение всевозможных множеств вида∏︁[ , ) = { ∈ : 6 < }.(1.25)=1Определение 2. Стар-дискрепансом набора {1 , . . . , } называется величина* (1 , . . .
, )⃒⃒⃒ #{ : ∈ }⃒⃒− ()⃒⃒ ,= sup ⃒∈ *(1.26)где – объединение множеств вида∏︁[0, ) = { ∈ : 0 6 < }.(1.27)=1Теорема 3 (Неравенство Коксмы-Хлавки). Пусть {1 , . . . , } – произвольный набор точек вгиперкубе . Тогда для ошибки интегрирования выполнено⃒⃒⃒∫︁⃒⃒⃒1 ∑︁⃒⃒ ( )⃒ 6 ( )* (1 , . . . , ),⃒ () −⃒⃒ =1⃒⃒3(1.28)В работах Соболя И.М. (например, [27]) близкий по сути термин носит название отклонения.
Мы будем использовать прямую транслитерацию для того, чтобы избежать смысловой коллизии с понятием стандартного отклонения.20где ( ) – вариация функции в смысле Харди-Краузе4 . Кроме того, это неравенство являетсяточным в следующем смысле: для любого набора {1 , . . . , } и для любого > 0 существуеттакая функция , что () = 1 и⃒⃒⃒∫︁⃒⃒⃒1 ∑︁⃒⃒( )⃒ > * (1 , . . . , ) − .⃒ () −⃒⃒ =1⃒⃒(1.29)Интерпретация этого фундаментального неравенства состоит в том, что верхняя границадля ошибки интегрирования распадается на произведение двух множителей, первый из которых зависит только от свойств подынтегральной функции, а второй – только от того, насколькоравномерно покрывается элементами последовательности 1 , .