Снижение дисперсии оценок Монте-Карло тема диссертации и автореферата по ВАК РФ 01.01.05, кандидат наук Иосипой Леонид Сергеевич

• Иосипой Леонид Сергеевич
• кандидат наук
• 2021

Оглавление диссертации кандидат наук Иосипой Леонид Сергеевич

Введение

Список литературы

Приложение

A Статья 1. Variance Reduction in Monte Carlo Estimators via Empirical Variance

Minimization

B Статья 2. Variance reduction for Markov chains with application to MCMC

C Статья 3. Об оценке плотности распределения с помощью ряда Фурье

Введение

В данной работе изучается задача вычисления математического ожидания п(/) := Е [/(X)] некоторой функции / : X ^ М от случайного вектора X с плотностью распределения п(ж) и значениями в множестве X С М^ в предположении, что / € Ь2(п). Мы будем предполагать, пока не сказано иное, что функция /(ж) и плотность п(ж) известны.

Сформулированную задачу можно рассмотреть как задачу вычисления интеграла от функции /(ж) • п(ж) по множеству X. Мы будем предполагать, что функция /(ж) и плотность п(ж) достаточно сложные, и данный интеграл явно не вычисляется. Поэтому мы будем оценивать его с помощью алгоритмов численного интегрирования. Известно, что любые детерминированные (даже адаптивные) алгоритмы численного интегрирования «страдают» от так называемого «проклятия размерности», которое означает, что сложность данных алгоритмов экспоненциально растет с ростом размерности см., например, классическую работу Н.С. Бахвалова [В71] и современный обзор Е.Новака [N16]. В свою очередь, рандомизированные алгоритмы численного интегрирования, к которым относятся методы Монте-Карло, оказываются эффективнее других подходов к решению данной задачи, если размерность ! является большой и/или функция / является сложно вычислимой, см., например, обзор Г. Робертса и Д. Розенталя [КИ04] или книгу К. Роберта и Д. Каселлы [И.С99].

В случае, когда можно получить независимую выборку Х1 ..., Хп из распределения с плотностью п(ж), математическое ожидание п(/) можно оценить с помощью оценки Монте-Карло

1 п

пп(/) := П Е /(Хк).

к=1

Данная оценка является несмещенной, то есть Е[пп(/)] = п(/), и имеет дисперсию, равную Уаг(пп(/)) = V(/)/п, где, здесь и далее, V(/) обозначает дисперсию функции / по отношению к п, то есть

V(/) := Уаг(/(X)), X ~ п.

С помощью центральной предельной теоремы можно построить асимптотический доверительный интервал для п(/), который будет иметь вид

Пп(/) ± , (1)

V п

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

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

По методам снижения дисперсии доступна обширная литература, см., например, классические книги И. Димова [D08], П. Глассермана [G13], Р. Рубинштейна и Д. Кроуза [RK16]. Стоит отметить, что большая часть методов снижения дисперсии явно используют структуру функции f (x) либо плотности п(х), поэтому не применимы в самых общих постановках. Одним из немногих общих методов снижения дисперсии является метод контрольных функций. Его идея заключается в следующем. Допустим, нам дан класс функций G С L2(n), который обладают свойством, что для любой g G G выполняется п(д) = 0. Такие функции называются контрольными. Тогда можно рассмотреть новую оценку Монте-Карло

1 n

nn(f - g) := (f (Xk) - )),

k =

где g — функция из G, на которой достигается (или почти достигается) минимум дисперсии V(f - g) = Var(f (X) - g(X)), X - п, то есть

g G argmin V(f - g). (2)

geG

В случае, когда класс контрольных функций G выбран хорошо, такой поход может значительно снизить дисперсию. Здесь возникают два ключевых вопроса. Как конструктивно построить класс контрольных функций g G G и как численно решить задачу (2), так как дисперсия V(f - g) не предполагается известной, а ее вычисление — более сложная задача, чем исходная.

Вопрос построения класса контрольных функций G был изучен ранее в работе А. Миры, Р. Солги, Д. Импарато [MSI13] и работах К. Оутса, М. Джиролами, Н. Шопена [OGC17] и [OGC19], поэтому на нем подробно мы останавливаться не будем. Мы лишь упомянем, что популярным методом построения функций g G G является подстановка в

gф = (Ф, V log п) + &у(Ф) (3)

различных функций Ф : X ^ Rd; здесь (•, •) — стандартное скалярное произведение в Rd, а div^) — дивергенция Ф. При достаточно слабых условиях на п и Ф, с помощью интегрирования по частям можно показать, что п^ф) = 0, см. Предложения 1 и 2 в [MSI13]. Контрольные функции вида (3) называются контрольными функциями Стейна.

Что касается вопроса о численном решении задачи (2), кажется естественным заменить истинную дисперсию V(f - g) на выборочную дисперсию Vn(f - g), построенную по выборке

XI,..., Xn, то есть на

1 п

К(/ - д) := — Е(/(Xk) - - Пп(/ - 5))2. (4)

к =

Кроме того, мы также заменим класс 0 на его дискретную конечную аппроксимацию 0', так как численная оптимизация по бесконечному множеству функций является сложной задачей, которую в общем случае решить невозможно. Возникает естественный вопрос: сколько мы «проиграем» за счет замены истинной дисперсии на эмпирическую и множества 0 на его аппроксимацию 0'? Ответ на этот вопрос является основной целью данной работы.

Настоящая диссертация состоит из введения и приложения. Введение поделено на три пронумерованных раздела. В Разделе 1 будут сформулированы базовые определения, которые будут использоваться в следующих разделах. Описание основных результатов будет дано в Разделе 2 и Разделе 3. В Разделе 2, мы оценим с большой вероятностью разность

V(/ - д) - 1п£ V(/ - д), (5)

9&3

где

д € arginf Vn (f — g), G' является е-сетью для G,

geG>

и дисперсия V(f — д) является условной по выборке Xi,...,Xn, на которой вычислялась оценка д. Подобные оценки часто встречаются в статистической литературе и называются оценками на «решетках» (в английской литературе skeleton estimates или sieve estimates), см., например, работы В. Вонг и К. Шен [WS95], Л. Девроя, Л. Гёрфи, Г. Лугоши [DGL96] или С. ван де Гир [G00]. Далее, в Разделе 3 мы рассмотрим случай, когда независимую выборку Xi..., Xn из распределения с плотностью п(х) явно получить нельзя, но сама п(х) является известной (возможно, с точностью до нормировочной константы). В таких случаях часто можно воспользоваться алгоритмами Монте-Карло на основе марковских цепей (MCMC в английской литературе) и получить марковскую цепь {Xk }£=1, которая сходится по распределению к распределению с плотностью n(x). Для марковских цепей, при некоторых условиях, также верна центральная предельная теорема и, следовательно, с ее помощью можно построить асимптотический доверительный интервал

/V TO(f)

nn(f) ± -f, (6)

n

где VTO(f) — асимптотическая дисперсия, которая определяется как

. (7)

V~(f) := lim n • E

n

1 n \2'

- £ f (Xfc) — n(f)

n

ч fc =

Ввиду возможных попарных корреляций у цепи асимптотическая дисперсия )

не равна в общем случае дисперсии V(/). Следовательно, в данном случае целесообразно подбирать контрольные функции д € 0, оптимизируя уже асимптотическую дисперсию Vто(/—д). Так как асимптотическая дисперсия так же, как и обычная дисперсия, не предполагается известной, мы ее заменим на ее оценку по выборке. Итак, в Разделе 3 мы оценим с большой вероятностью разность

где уже для некоторой оценки асимптотической дисперсии V

д € ащшт (/ — д), 0' является е-сетью для 0,

дед'

а дисперсия Vто(/ — д ) является условной по выборке Х\,..., Хп, на которой вычислялась оценка д. Данное обобщение оказывается нетривиальным по двум причинам. Во-первых, в связи с попарными корреляциями, оценивание асимптотической дисперсии возможно несколькими способами и требует специальных техник, таких как спектральные методы или разделение марковской цепи на слабо коррелирующие блоки, см. обзор методов оценки асимптотической дисперсии в работе Д. Флегала и Д. Джонса [ЕЛО]. Во-вторых, неасиптотический анализ оценки / — дд требует обобщения используемых в доказательстве неравенств концентрации на цепи Маркова, которое также не является тривиальным. Мы проведем данный анализ для достаточно широкого класса марковских цепей, который включают в себя многие МСМС алгоритмы. Наконец, в приложении будут приведены копии статей, которые выносятся на защиту.

Введение диссертации (часть автореферата) на тему «Снижение дисперсии оценок Монте-Карло»

Цель работы

Цель настоящей работы — изучение теоретических свойств описанного выше метода снижения дисперсии оценок Монте-Карло (метода контрольных функций) и, в частности, получение неасимптотических оценок на избыточный риск (5) и (8) в независимом и зависимом случае соотвественно.

1. Обозначения и определения

Перед тем как сформулировать основные результаты, введем еще некоторые обозначения. Обозначим класс функций Н(ж) = /(ж) — д(ж) для д € 0 через Н:

где 0 — класс контрольных функции (то есть таких, что п(д) = 0 для всех д € 0). Заметим, что у всех функций Н € Н будет одинаковое математическое ожидание, равное п(/). Как уже

(8)

Н := {/ — д : д € 0},

упоминалось, вместо оптимизации по всему классу H мы будем искать минимум в его конечной аппроксимации. Зафиксируем е > 0 и положим r = 1 или r = 2. Предполагая, что класс H вполне ограничен, обозначим через He,r произвольную минимальную по мощности е-сеть множества H в Lr (п)-норме, то есть такой минимальный набор функций He,r = {hi,..., hm} С H, что для произвольной h € H существует h* € He,r, что расстояние между h и h* в Lr (п)-норме не превосходит е. Обозначим метрическую энтропию H через HLr(n)(H,e) := log |He,r|, где |He,r | — это мощность He,r. Далее, определим величину, которая будет использована в основных результатах данной работы:

1ьг(n)(H,n) := inf{n > 0 : HLrw(H,r) < -rf).

Обратим внимание, что r > 0, удовлетворяющее неравенству H^r(n)(H,r) < -nr, существует и конечно, так как метрическая энтропия HLr(n)(H, n) является убывающей функцией по n, а отображение r ^ -nr является строго возрастающим. Значение YLr(n)(H, n) будет использовано для того, чтобы оценить мощность класса He,r. Действительно, выбирая е > YLr(n)(H,n) мы получим |He,r | < ene . Легко также видеть из определения, что 7L^(n)(H, n) является убывающей функцией по n.

Приведем теперь оценку на yl(n)(H, n) в достаточно общем случае, когда H является подмножеством взвешенного соболевского пространства. Напомним, что соболевское пространство Ws,p(X) определяется как Ws,p(X) = {h € Lp(A) : Dah € Lp(A), V|a| < s}, где А является мерой Лебега, a = (ai,..., ad) является мульти-индексом с |a| = ai + ... + ad, а Da является дифференциальным оператором вида Da = d |a|/d xai ... дж^. Здесь все производные понимаются в слабом смысле. Взвешенное соболевское пространство Ws,p(X, (x)e) для полиномиальной весовой функции (x)e = (1 + ||х||2)в/2, в € R, определяется следующим образом:

Ws,p(X, (х)в) = {h € Lp(A) : h • (x)e € Ws,p(X)}. Норма в данном взвешенном соболевском пространстве определяется как

Ww- = £||Da(h )|

'...........МАГ

| a| <s

Мы будем говорить, что Н С ШЯ'Р(Х, (х)в) является ограниченным подмножеством пространства ШЯ,Р(Х, (х)в), если Н вполне ограничено и существует константа с > 0, такая что

|w|,p < с для всех h € H. Верно следующее предложение.

Предложение 1. Пусть Н будет (непустым) ограниченным подмножеством Ш(х)в), где 1 < р < то, в € К и в - ¿/р > 0. Пусть также для некоторого к > 0, ||(х)к-в< то.

Тогда для г = 1, 2,

п при к > в — ¿/р,

(п)(Н,п)

<

П г+(к/а+1/р) 1 при К < в — ¿/р,

где символ < означает неравенство с точностью до некоторой константы, не зависящей

от п.

Данный результат получен подстановкой в определение (П)(Н, п) оценки на энтропию ограниченных подмножеств ШЯ'Р(Х, (ж)в) из Следствия 4 работы Р. Никля и Б. Петшера [МР07]. Предложение 1 не является основным результатом диссертационной работы и приведено здесь в качестве примера возможных значений величины (П)(Н, п), которая далее будет использоваться в основных результатах.

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

2. Снижение дисперсии в независимом случае

Перед тем как сформулировать основные результаты для независимого случая, напомним, что для класса Н = {/ — д : д € 0} оценку, свойства которой мы будем изучать, можно записать как

Нег € ащтш ^П(Н),

ненЕ,г

где ^П(Н) — эмпирическая дисперсия, определенная в (4), а Не,г — е-сеть множества Н в норме Ьг(п), г =1, 2. Подробные определения даны в Разделе 1.

Теперь мы можем перейти к основным результатам раздела. Мы начнем со случая, когда функции Н € Н являются ограниченными и класс Н является замкнутым и выпуклым множеством.

Теорема 2. Пусть класс Н является замкнутым и выпуклым множеством. Пусть также вирнеч 8иРхех |Н(ж)| < Ь для некоторого Ь > 0 и пусть п(Н) = с для всех Н € Н и некоторой с € М. Тогда для некоторого е > 0 (значение которого определено в доказательстве) и любого 5 € (0,1), с вероятностью не менее 1 — 5,

V(Нел — чV(н) <ь27ь1(п)(н,п) + ,

V / нен п

где первая дисперсия является условной по выборке Х1,..., Хп, на которой считалась оценка Нед, а символ < означает неравенство с точностью до некоторой абсолютной константы.

Теперь сформулируем подобный результат, когда вместо выпуклости класса Н, мы потребуем, чтобы Н содержал функцию, тождественно равную константе. Так как по построению п(Н) = п(/) для всех Н € Н, эта константа должна быть равна п(/) и, следовательно, в этом случае Мнен V(Н) = 0.

Теорема 3. Пусть класс Н содержит константную функцию, то есть Н*(х) = с для некоторой Н* € Н и некоторого с € М. Пусть также вирнен 8иРхех |Н(ж) | < Ь для некоторого Ь > 0 и пусть п(Н) = с для всех Н € Н и некоторой с € М. Тогда для некоторого е > 0 (значение которого определено в доказательстве) и любого 5 € (0,1), с вероятностью не менее

1 — 5,

V (НеД) < ь27Ь1(п)(н, п) +

и

V(Не,2) < Ь2Ы(П)(Н,п))2 +,

где дисперсии являются условными по выборке Х1,..., Хп, на которой считались оценки Не1 и Не,2, а символ < означает неравенство с точностью до некоторой абсолютной константы.

Предполагая дополнительно к условиям Теоремы 2 или Теоремы 3, что Н является ограниченным подмножеством некоторого взвешенного соболевского пространства ШЯ'Р(Х, (х)в), и используя результат Предложения 1, мы получим окончательную оценку с большой вероятностью вида:

V (Нм) — 1Пн V(Н) < п-а, а € (0,1), (9)

где значение параметра а будет зависеть от размерности й и параметров соболевского пространства в, р и в. Подобные порядки сходимости называются «быстрыми» в статистической литературе, так как позволяют получать скорость сходимости быстрее стандартных порядков 1/л/п при параметрическом оценивании. В качестве следствия стоит отметить, что порядок длины асимптотического доверительного интервала (1) с помощью данной процедуры снижения дисперсии можно уменьшить с п-1/2 до п-1/2-а/2, если класс Н выбран так, что Мнен V(Н) < п-а.

Насколько автору известно, Теорема 2 и Теорема 3 являются первыми доступными в литературе результатами об оценках, полученных минимизацией эмпирической дисперсии ^П(Н). Отличительной особенностью данных результатов является тот факт, что эмпирическая дисперсия ^(Н) не является суммой независимых случайных величин Н(Х1),..., Н(ХП), как часто бывает в теории эмпирических процессов. Вместо этого, ^П(Н) является и -статистикой. «Быстрые» порядки в задачах минимизации и -статистик были ранее получены, насколько известно автору, только в работе С. Клеменсона, Г. Лугоши, Н. Ваятиса [СЬУ08] (и некоторых других последующих работах авторов), где рассматривалась задача ранжирования случайных

величин. В доказательствах Теоремы 2 и Теоремы 3 автором используется другая техника, основанная на дискретизации множества Н и концентрации и-статистик, что позволяет избежать многих технических проблем, которые присутствовали в [СЬУ08].

Приведем теперь здесь два следствия Теоремы 3, в которых мы явно вычислим значение параметра а из оценки (9). Для простоты мы рассмотрим одномерный случай X = К ив качестве контрольных функций мы возьмем контрольные функции Стейна (3), которые в одномерном случае можно записать как

дф(ж) := ф(ж) • п(ж))' + ф'(ж) = -(-у (ф(ж)п(ж))',

где ф : К ^ К — функции из некоторого класса Ф, который мы определим позднее. Обозначим класс в раз непрерывно дифференцируемых функций на К через СЯ(К) и класс функций из СЯ(К), производные которых растут не быстрее некоторого полинома через Сро1у(К), через

Сьро1у (К) := {ф € С8 (К) : 3 т € М, такой что |ф(к)(ж)| < |ж|т при |ж| ^ то, Ук = 0,. ..,в}.

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

Следствие 4. Пусть п(ж) ж е-с|х|' для некоторых с € К и д € N. Пусть также / € С8о1у(К) для некоторого в € N. Зафиксируем произвольное Ф С С®+у(К), такое что соответствующий класс Н = {/-дф : ф € Ф} является ограниченным, содержит константную функцию и 8иРхех 1^(ж)| < Ь для некоторого Ь > 0. Тогда для некоторого е > 0 и любого 5 € (0,1), с вероятностью не менее 1 - 5,

V (Ч < (1) ^ + ^,

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

Рассмотрим теперь случай, когда плотность п убывает полиномиально. Для этого сначала обозначим класс функций из СЯ(К), производные которых растут не быстрее полинома порядка т € М, через

Сро1у<т(К) := {ф € С8(К) : |ф(к)(ж)| < |ж|т при |ж| ^ то, Ук = 0,..., в}.

Следствие 5. Пусть п(ж) ж (1 + ж2)-9 для д € N. Пусть также / € С®о1у<т(К) для некоторых в,т € N. Зафиксируем произвольное Ф С С®+у(К), такое что соответствующий класс Н = {/ - дф : ф € Ф} является ограниченным, содержит константную функцию и вирьеИБир^х |^(ж)| < Ь для некоторого Ь > 0. Тогда для некоторого е > 0 и любого 5 € (0,1),

с вероятностью не менее 1 — S,

iА Т+Т77 iog(i/s) .

— +----при s < 2q — m + 1/p — 3,

Vi <{ \nl г n

1\ Г + (2,-т+Г7Р-3)-Г log(l/S)

— +--при s > 2q — m + 1/p — 3,

n J n

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

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

3. Снижение дисперсии в зависимом случае

В этом разделе мы сформулируем аналоги Теоремы 2 и Теоремы 3 в зависимом случае. Пусть {Хк}£=о будет марковской цепью (в этом разделе мы нумеруем последовательность для удобства с 0) с переходным ядром Р и начальным значением Хо = хо для некоторого хо € X. Мы будем предполагать, что марковское ядро Р имеет единственную инвариантную меру п.

Как упоминалось во Введении, известно несколько оценок для асимптотической дисперсии Vсм. работу Д. Флегала и Д. Джонса [ЕЛ0]. Для простоты изложения мы рассмотрим только спектральную оценку, так как она является наиболее общей. Тем не менее, результаты, аналогичные приведенным ниже, можно получить и для других оценок. Начнем с определения. Пусть Н = {/ — д : д € 0}, где 0 — класс контрольных функций, см. Раздел 1. Обозначим автоковариационную функцию процесса {Л.(Хк)}£=о через

р(5) := Е[ЩХо) — п(^))(Ь(Х*) — п(^))],

где предполагается, что Хо ~ п. Обозначим также выборочную автоковариационную функцию через

1 п — в — 1

Р«(5):= П Е ) — пп(Н))(Н(Хк+в) — Пп(М).

к =о

Легко видеть, что асимптотическую дисперсию Vопределенную в (7), можно переписать как

+^

V~(ь) = Е р(И).

Спектральная оценка асимптотической дисперсии основана на взвешенных выборочных автоковариационных функциях:

Ь„-1

К°°(Н) := £ ™п(*)рп(М), (10)

8 = -(Ь„-1)

где и>„ — скользящее окно, а Ьп — ширина окна. Ширина окна Ьп при любом значении п является натуральным числом, а скользящее окно тп является ядром вида тп(в) = т(в/Ьп), где т — симметричная неотрицательная функция с носителем [-1,1], равная

{2в + 2, -1 < в< -1/2, 1, -1/2 < в < 1/2, -2в + 2, 1/2 <в < 1.

Возможен и другой выбор т(в), см. подробнее в [ЕЛ0]. В данном разделе мы будем изучать свойства оценки

Не 2 := ащштЪ"°(Н),

ьене,2

где Не,2 — е-сеть множества Н в норме Ь2(п), см. подробные определения в Разделе 1.

Напомним, что для двух марковских ядер Р и Q на X х X, где X — борелевская а-алгебра на X, их произведение тоже является марковским ядром и определяется для ж € X и А € X как PQ(ж, А) = / Р(ж, А). Кроме того, ядро Р" определяется рекуррентным соотношением

Р" = РР"-1. Пусть теперь Ш : X ^ [1, то) — измеримая функция. Ш-норма произвольной функции Н : X ^ К определяется как ||Н||^ = 8иРхех{|Н(ж)|/Ш(ж)}. Для двух вероятностных мер р и V на (X, X), удовлетворяющих ^(Ш) < то и ^(Ш) < то, Ш-нормой разности р - V является по определению - V= Бирц^ц^<1 |м(Н) - V(Н)|.

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

(GE) Марковское ядро Р имеет единственную инвариантную меру п. Кроме того, для некоторой измеримой функции Ш : X ^ [1, то), такой что п(Ш) < то, и некоторых чисел ? > 0, р € (0,1) выполняется для всех ж € X и всех п € N

||Р"(ж, •) - п||^ < (ж)р".

Известно, что условие (СЕ) влечет центральную предельную теорему для последовательности {Н(Хк)}£=0, Н € Н, если дополнительно п(Н2+к) < то для некоторого к > 0, см. работу И. А. Ибрагимова и Ю.В. Линника [!Ь71] или Д. Джонса [ДО4]. Чтобы асимптотический

доверительный интервал (6) был корректно определен, далее мы будем предполагать, что Н С Ь2+к для некоторого к > 0.

Напомним, что Ь2-метрикой Канторовича (в зарубежной литературе — Ь2-метрикой Вас-серштейна) между вероятностными мерами м и V называется

1/2

Ух — у||2 ас(х,у)

С \JXxX

где || • || — евклидова норма в М^, а инфимум взят по всем вероятностным мерам £ на декартовом произведении X х X, у которых маргинальными мерами являются м и V соответственно. Расстояние Кульбака-Лейблера между мерами м и V определяется как

II ^ Л 1о®(^если М < V,

то, иначе.

Мы будем говорить, что вероятностная мера м удовлетворяет (информационному) транспортному неравенству Т2(С), если существует С > 0, такая что для любой вероятностной меры V выполняется

"ЪМ < У2СКЦ^Й.

Для доказательства основных результатов нам необходима концентрация величины ^П°(Н). Важно отметить, что ^П—(Н) является квадратичной формой по {Н(Х^)}П=о1. Нами была доказана гауссовская концентрация ^П°(Н) в случае, когда функций Н € Н являются липшицевыми, а марковское ядро Р является сжимающим отображением в Ь2-метрике Канторовича.

(Ь) Функции Н € Н являются Ь-липшицевыми, то есть |Н(х) — Н(у)| < Ь||х — у|| для некоторого Ь > 0 и всех х, у € X.

(CW) Марковское ядро Р(х, •) удовлетворяет Т2(С) для всех х € X и некоторого С > 0. И, кроме того, существует г € (0,1), такое что "^(Р(х,-),Р(у, •)) < г||х — у|| для всех х, у € X.

Теперь мы можем сформулировать следующую теорему.

Теорема 6. Предположим (ОЕ), (Ь) и (СШ). Положим Ьп = 2^(п)/^(1/р). Тогда для некоторого е > 0 (значение которого определено в доказательстве), любого начального значения хо € Х и любого 5 € (0,1), с вероятностью не менее 1 — 5,

У-ГНЛ — ш! V-(Н) < А11оё(п)7ь2(п)(Н,п) + А21ое(п)1;ё(1/5), V / ней л/п

где первая асимптотическая дисперсия является условной по выборке Х1,..., Хп, на которой вычислялась оценка Не,1, символ < означает неравенство с точностью до некоторой абсолютной константы и

УСь2 V?(п(Ш) + Ш(жо)) (УСь2 , ..,.,2

А1 = п-V!—ТТУТ, А2 = л-1—/ 0 1-+ 8ИР ||Н||ж 1

(1 - г)^(1/р) а/1 - Р!о^(1/р) - г ьен

Предполагая дополнительно к условиям Теоремы 6, что Н является ограниченным подмножеством некоторого взвешенного соболевского пространства Ш^^ (ж)в), и подставляя результат Предложения 1, мы получим окончательную оценку с большой вероятностью вида:

V~ (Не,2) - т£ V~(Н) < п-а, а € (0,1/2),

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

Порядок сходимости может быть улучшен до п-а для а € (0,1), если дополнительно предположить, что Н содержит константную функцию. Так как п(Н) = п(/) для всех Н € Н, эта константа должна быть равна п(/) и, следовательно, в этом случае ш^ен VTO(Н) = 0.

Теорема 7. Предположим (ОЕ), (Ь) и (СШ). Предположим также, что Н содержит константную функцию, то есть Н* (ж) = с для некоторой Н* € Н и некоторого с € К. Положим Ьп = 2^(п)/^(1/р). Тогда для некоторого е > 0 (значение которого определено в доказательстве), любого начального значения ж0 € X и любого 5 € (0,1), с вероятностью не менее 1 - 5,

2 . ^(п)^(1/5)

V~ (Не,2) < А1 log(n)(7ь2(п)(Н, п))2 + А2

(п)( . ,

п

где асимптотическая дисперсия является условной по выборке Х1,..., Хп, на которой вычислялась оценка Не1, символ < означает неравенство с точностью до некоторой абсолютной константы и

. = Сь2 = Сь2 + ? (п(Ш)+ Ш (жо))

11 = (1 - г)2 ь^1/р), А2 = (1 - г)2 ^(1/р) + (1 - р)1/2 ^(1/р) 172.

2

Предполагая теперь, что Н является ограниченным подмножеством Ш^^ (ж)в), мы получим оценку с большой вероятностью вида:

V~ (\,2) < п-а, а € (0,1).

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

В доказательстве Теоремы 6 и Теоремы 7 использовалось следующее неравенство концентрации, которое, как кажется автору, представляет самостоятельный интерес.

Предложение 8. Пусть {Хк}—=о является марковской цепью с переходным ядром Р и начальным значением Хо = хо для некоторого хо € X. Предположим, что существует С > 0, такое что Р(х, •) удовлетворяет Т2(С) для всех х € X. Предположим также, что существует число г € (0,1), такое что для всех х, у € X выполняется

"2(Р(х, 0,Р(у, •)) < г||х — у||.

Для произвольной Ь-липшицевой функции Н : X ^ М обозначим ^„(Н) = (Н(Хо),..., Н(ХП-1))Т. Тогда для любой матрицы А размера п х п и любого £ > 0,

а

p(|Z„(h)TAZ„(h) - E[Z„(h)TAZ„(h)] | >t) < 2exp(^--

t2

сК2 (Е||А^„(Н)||2 + ¿||А||) где К2 = СЬ2/(1 — г)2, а с > 0 — некоторая абсолютная константа.

Доказательство данного утверждения состоит из двух этапов. Сначала, используя результат из статьи Х. Джеллоута, А. Гийена и Л. Ву [DGW04], мы показываем, что совместное распределение {Хкпринадлежит классу Т2(С/(1 — г)2), что, в свою очередь, влечет гаус-совскую концентрацию для всех липшицевых функций. Далее, используя технику из работы Р. Адамчака [А15], мы доказываем концентрацию для квадратичных форм.

Алгоритм Ланжевена

Простым примером MCMC алгоритма, для которого выполняются условия (GE) и (CW) является алгоритм Ланжевена (в английской литературе Unadjusted Langevin Algorithm или ULA). Пусть плотность п(ж) известна с точностью до нормировочной константы и выражается в виде п(ж) ж e-Uдля некоторой неотрицательной функции U(ж), которую принято называть потенциалом. В алгоритме Ланжевена марковская цепь }k>o строится с помощью следующего рекуррентного соотношения

Xfc+1 = Xk - 7VU(Xk) + y2Yzfc+1, (11)

где {Zk}k>i — последовательность независимых стандартных нормальных векторов в Rd, а Y > 0 — шаг алгоритма. Данный алгоритм является дискретизацией первого порядка (по схеме Эйлера-Маруямы) стохастического дифференциального уравнения Ланжевена:

dYt = -VU (Yt)dt + V2dBt,

где {Bt}t>o — d-мерное броуновское движение. При некоторых условиях, единственным стационарным распределением динамики Ланжевена {Yi}t>o является п, а марковская цепь {Xk}k>o

сходится к единственному стационарному распределению п7, которое близко к распределению п в некоторых метриках, таких как расстояние по вариации или Ь2-метрике Канторовича. Подробное описание алгоритма и теоремы сходимости доступны в работах Г. Робертса и Р. Твиди [RT96], А. Далаляна [D17], А. Дурмуса и Э. Мулине [DM17]. В случае, когда плотность п(ж) не является известной, для использования алгоритма Ланжевена можно воспользоваться методами из статей Д. В. Беломестного и автора [BI19] и [BI20].

Рассмотрим следующие условия на потенциал U(ж), которые будут гарантировать выполнение условий (GE) и (CW):

(U1) Потенциал U является дважды непрерывно дифференцируемым U G C2(Rd) и имеет липшицевый градиент, то есть существует L > 0, такое что для всех ж, y G

||VU(ж) -VU(y)|| < L||x - y||.

(U2) Потенциал U является строго выпуклым, то есть существует константа m > 0, такая что для всех ж, y G выполняется

U(y) > U(ж) + (VU(x),y - ж) + m||ж - y||2.

Предложение 9. Пусть выполнены условия (U1) и (U2). Тогда при любом y G (0, 2/(m+L)), марковское ядро PY, соответствующее цепи }k>0 из (11), удовлетворяет условиям (GE) и (CW) при некотором р G (0,1),

W(ж) = ||ж||2, C = 2y и r = - 2YmL/(m + L).

Данный факт является непосредственным следствием результатов из двух работ А. Дур-муса и Э. Мулине [DM17] и [DM19].

Байесовская логистическая регрессия. Рассмотрим также численный пример для алгоритма Ланжевена. Мы сравним два описанных метода выбора контрольных функций: минимизацию эмпирической дисперсии (EVM метод), см. (4), и минимизацию спектральной дисперсии (ESVM метод), см. (10). Сравнение будем проводить на двух наборах данных для задачи классификации из базы данных UCI. Первый набор данных называется Pima и содержит N = 768 наблюдений в размерности d = 9. Второй набор данных называется EEG и содержит N = 14 980 наблюдений в размерности d = 15.

В логистической регрессии вероятность i-ого ответа yj G {-1,1} для i = 1,..., N, вычисляется как р(у^х^в) = (1 + где Xj G Rd — это вектор признаков (факторов) и в G — вектор неизвестных регрессионных коэффицииентов. Сначала мы разобьем данные на обучающую выборку Ttrain = {(yj,Xj)}^-^ и тестовую выборку Ttest = {(y-,xj)}K=1,

выбрав случайно К = 100 наблюдений из данных (и удалив их из обучающей выборки). Далее, мы наложим априорное стандартное нормальное распределение на параметр в, обозначим его Ро(в), и используем алгоритм ИЬА для получения выборки из апостериорного распределения р(в|Г1гаш) ГС Ро(в)П(ур(у* 1х,в). С помощью полученной выборки {вк}„1(1 мы оценим предсказание для фиксированной точки из тестовой выборки (у', х!), то есть р(у'|х') = /к<1 р(у'|х', в) р(в|Т4га'")ав, вычислив оценку п-1 /№) для /(в) = р(у'|х',в).

Чтобы избавиться от случайности, мы также оценим среднее предсказаний по всему тестовому набору , вычислив оценку для f (в) = К-1 XI(у х)ет*ев* р(у^ |х^, в). Кроме того, аналогичные оценки мы посчитаем для ЕУМ и ЕЯУМ методов, в которых мы будем использовать контрольные функции Стейна, см. (3), где Ф(х) = Ь для некоторого Ь € ММ.

Графики с разбросом полученных оценок даны на Рис. 1. Обратим внимание, что оба метода ЕУМ и ЕЯУМ дают значительный выигрыш в дисперсии, но метод ЕЯУМ работает лучше, так как учитывает марковскую структуру цепи.

Рис. 1: Оценивание среднего предсказания в байесовской логистической регрессии. Первые два графика относятся к набору данных Pima, а последние два — к набору данных EEG. В первых графиках для каждого набора данных отражены 3 «ящика с усами» для 100 оценок полученных соответсвенно обычным алгоритмом ULA, алгоритмом ULA с EVM и алгоритмом ULA с ESVM, а следующий график — увеличение масштаба для последних двух алгоритмов.

Основные результаты диссертационной работы

В рамках диссертационной работы получены следующие основные результаты:

1) Неасимптотические оценки на избыточный риск (5) в независимом случае для произвольного класса контрольных функций 0. Данные результаты сформулированы в Теореме 2 и Теореме 3. Кроме того, продемонстрировано, какие получаются оценки в случае контрольных функций Стейна, см. Следствие 4 и Следствие 5.

2) Неасимптотические оценки на избыточный риск (8) в зависимом случае для произвольного класса контрольных функций 0. Данные результаты сформулированы в Теореме 6 и Теореме 7. Кроме того, в данной работе впервые было предложено минимизировать вместо эмпирической дисперсии (4) спектральной оценку асимптотической дисперсии (10). Данное решение лучше работает на практике в случае зависимых случайных величин.

3) Новое концентрационное неравенство для квадратичной формы от марковской цепи, когда ее марковское ядро является сжимающим отображением в Ь2-метрике Канторовича, см. Предложение 8.

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

Личный вклад автора

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

Апробация результатов

Результаты диссертационной работы докладывались на следующих конференциях и научно-исследовательских семинарах:

• Конференция «Информационные технологии и системы» в ИППИ РАН, Санкт-Петербург, сентябрь 2016 года. Тема доклада: «Концентрация нормы изотропного логарифмически вогнутого случайного вектора».

• Аспирантский коллоквиум кафедры теории вероятностей в МГУ, Москва, декабрь 2017 года. Тема доклада: «Снижение дисперсии оценки Монте-Карло».

• Конференция «New frontiers in high-dimensional probability and statistics» в НИУ ВШЭ, Москва, февраль 2018 года. Тема доклада: «Снижение дисперсии оценок Монте-Карло с помощью минимизации выборочной дисперсии».

• Конференция «Ломоносов» в МГУ, Москва, апрель 2018 года. Тема доклада: «Снижение дисперсии оценок Монте-Карло с помощью минимизации выборочной дисперсии».

• Городской семинар по теории вероятностей и математической статистике в ПОМИ РАН, Санкт-Петербург, сентябрь 2018 года. Тема доклада: «Снижение дисперсии оценок Монте-Карло с помощью минимизации выборочной дисперсии».

• Зимняя школа «New frontiers in high-dimensional probability and statistics 2» в НИУ ВШЭ, Москва, февраль 2019 года. Тема доклада: «MCMC алгоритмы для распределений с известной характеристической функцией».

• Семинар центра прикладной математики в Политехнической школе (Ecole polytechnique), Париж, июнь 2019 года. Тема доклада: «Снижение дисперсии оценок Монте-Карло с помощью минимизации выборочной дисперсии».

Список литературы диссертационного исследования кандидат наук Иосипой Леонид Сергеевич, 2021 год

Литература

1. ГОЛУБЕВ Г.К. Непараметрическое оценивание гладких плотностей распределения в L2 // Проблемы передачи информации. - 1992. - №28. - С. 52-62.

2. ЕФРОИМОВИЧ С.Ю. Непараметрическое оценивание плотности неизвестной гладкости // Теория вероятностей и ее применения. - 1985. - №30. - С. 524-534.

3. ЗОРИЧ В.А. Математический анализ. - М.: МЦНМО, 2012.

4. ИБРАГИМОВ И.А., ХАСЬМИНСКИЙ Р.З. Об оценке плотности распределения // Записки научных семинаров ЛОМИ. - 1986. - №153. - С. 45-59.

5. ЛАГУТИН М.Б. Наглядная математическая статистика. -М.: Бином. Лаборатория знаний, 2009.

6. BORWEIN J.M., HUANG W.Z. A fast heuristic method for polynomial moment problems with Boltzmann-Shannon entropy // SIAM J.Opt. - 1995. - Vol. 5. - P. 68-99.

7. BRETAGNOLLE J., HUBER C. Estimation des densites: risque minimax // Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete - 1979. - Vol. 47. - P. 119-137.

8. DEVORE R.A., LORENTZ G.G. Constructive Approximation. -Springer, Berlin, 1993.

9. FASINO D. Approximation of nonnegative functions by means of exponential trigonometric polynomials // J. of Computational and Applied Mathematics. - 2002. - Vol. 140. - P. 315-329.

10. JIANG H. Uniform Convergence Rates for Kernel Density Estimation // Proc. of the 34th Int. Conference on Machine Learning, PMLR. - 2017. - Vol. 70 - P. 1694-1703.

11. HOEFFDING W. Probability inequalities for sums of bounded random variables // J. of the American Statistical Association. -1963. - Vol. 58, No. 301. - P. 13-30.

12. SERRA S. On the extreme spectral properties of Toeplitz matrices generated by L1 functions with several minima (maxima) // BIT Numerical Mathematics. - 1996. - Vol. 36. - P. 135-142.

13. SERRA S. On the extreme eigenvalues of Hermitian (block) Toeplitz matrices // Linear Algebra and its Applications. - 1998. -Vol. 270. - P. 109-129.

14. TSYBAKOV A. Introduction to Nonparametric Estimation. Springer Series in Statistics // Springer Series in Statistics. - 2009.

15. TYRTYSHNIKOV E. A Unifying Approach to Some Old and New Theorems on Distribution and Clustering // Linear Algebra and its Applications - 1996. - Vol. 232, No. 1. - P. 2-43.

ON DENSITY ESTIMATION VIA FOURIER SERIES

Denis Belomestny, National Research University Higher School of Economics, Moscow, Cand.Sc., professor ([e-mail]). Leonid Iosipoi, National Research University Higher School of Economics, Moscow, junior research fellow ([e-mail]).

Abstract: In this paper, we consider the classical statistical problem of probability density estimation based on a sample from this distribution. This problem naturally arises in many applications when one aims at investigation of a probability structure in a random process. For instance, it is possible to identify some structure in a complex system using density estimation. In this paper, a new approach to estimate a density function is proposed. This approach is based on approximation of a log-density via Fourier series with coefficients obtained by solving a system of linear equations. Analysis of theoretical properties of such estimate is the main purpose of this work. As the main results, we prove bounds on the difference between target density and its approximation in the supremum norm and the Kullback-Leibler divergence. Obtained rates are parametric and have order 0(1/VN) with high probability, which is a standard rate in parametric estimation problems. The constants in the rates are obtained up to an absolute factor, which means that we investigated the dependence on all parameters. As a numerical example, we consider a problem of Cauchy density estimation.

Keywords: density estimation, Fourier series, Kullback-Leibler divergence.

УДК 519.23 ББК 22.172

DOI: 10.25728/ubs.2019.82.2

Статья представлена к публикации членом редакционной коллегии Э.Ю. Калимулиной.

Поступила в редакцию 03.09.2019. Опубликована 30.11.2019.