Содержание к диссертации
Введение
Глава 1. Задачи, использующие распределение А2 6
1.1 Приближенное вычисление интегралов (задача 1) 6
1.2 Разделение ошибок в регрессионном анализе (задача 2) 14
1.3 Нахождение точных D-оптимальных планов (задача 3) 20
Глава 2. Моделирование распределения А2 в непрерывном случае 25
2.1 Подходы и известные результаты. Метод отбора 26
2.2 Оценка трудоемкости метода отбора 28
2.3 Условные плотности как смеси распределений 31
2.4 Оценка трудоемкости моделирования условных плотностей 35
2.5 Сравнение величины трудоемкости методов отбора и условных плотностей 41
2.6 Примеры вычисления интегралов 43
Глава 3. Моделирование распределения А2 в дискретном случае 53
3.1 Подходы и известные результаты. Метод обращения 53
3.2 Использование метода отбора 57
3.3 Использование условных плотностей 60
3.4 Сравнение трудоемкостей предложенных методов 65
3.5 Примеры разделения ошибок 66
Глава 4. Моделирование обобщенного распределения А2 и вычисление точных D-оптимальных планов 72
4.1 Вид условных плотностей обобщенного А2-распределения 73
4.2 Моделирование обобщенного распределения А2 и оценка его трудоемкости 77
4.3 Качество приближения 81
4.4 Замечание о моделировании в области сложного вида 84
4.5 Примеры нахождения D-оптимальных планов 84
Заключение 93
Список литературы
- Разделение ошибок в регрессионном анализе (задача 2)
- Условные плотности как смеси распределений
- Использование условных плотностей
- Моделирование обобщенного распределения А2 и оценка его трудоемкости
Введение к работе
Распределение А2 играет важную роль по крайней мере в трех задачах: уменьшения дисперсии при вычислении интегралов методом Монте-Карло (условимся называть ее в дальнейшем задачей 1), анализа систематической погрешности в регрессионном анализе (задача 2), а также построения точных jD-оптимальных планов в теории планирования эксперимента (задача 3). Первая задача рассматривается подробно в [4], а вторая в работах [4] и [16].
Использование распределения А2 требует его моделирования. Однако, до недавнего времени были предложены методы (в непрерывном случае -метод отбора, а в дискретном - метод обращения), чьи трудоемкости растут очень быстро с ростом п размерности моделируемого вектора. Это обстоятельство создавало значительные трудности при решении задач, в которых п велико.
Целями данной работы являются, во первых, построение эффективных алгортимов моделировния распределения А2 во всех его модификациях -непрерывной, дискретной и обобщенной. Во вторых, проиллюстрировать на важных примерах преимущество использования распределения А2 в задачах 1 и 2. В третьих, рассмотреть неисследованную в контексте А2 задачу 3, получить результаты о близости полученной аппроксимации точного D-оптимального плана и проиллюстрировать предложенную
процедуру на примерах. В частности, в двумерной области в случае квадратичной регрессии.
В Главе 1 изложены постановки и известные решения задач 1 и 2, а также исследована задача 3 в контексте обобщенного распределения А2. Главы 2, 3 и 4 посвящены моделированию А2 в непрерывном, дискретном и обобщенном случае, соответственно, причем в каждой из них приводятся содержательные иллюстративные примеры.
Разделение ошибок в регрессионном анализе (задача 2)
В регрессионном анализе разделение систематической ошибки от случайной проводится для того, чтобы оценить качество построенной регрессии. Систематическая ошибка приближения указывает на целесообразность составления регрессионной модели с помощью того или иного набора линейно независимых функций. Поэтому процедура ее выделения представляет большой интерес. Ермаковым и Швабе, [16], предложен подход, основанный на рандомизации оценок метода наименьших квадратов (МНК) с помощью дискретного распределения А2. При этом удается выделить систематическую погрешность как в случае независимых, так и в случае зависимых наблюдений. Остановимся на постановке задачи и изложим коротко вышеуказанную схему рассуждений.
Пусть /J, дискретная равномерная мера, сосредоточенная в N точках: ц(х{) = ±, i = l,N Экспериментально можно получить значения ъ ,/v, где & = f(xt) + ЄІ, ЕЄІ = 0; Ее] = а2; СОУ(ЄІ? Sj) = 0, і ф j. п Аппроксимируем f линейной комбинацией 5 еда, где { i}iLi линейно 1=1 независимы в X. Тогда оценки метода наименьших квадратов сі коэффициентов СІ имеют вид ([16], с.З)
Тогда случайная величина, определяемая соотношением является несмещенной оценкой сі, т.е. ЕХІ = сі (проверяется непосредственно). Более того, случайные величины Xi некоррелированы и имеют конечную дисперсию ([16], с.4): cov(x/,Xr) = N N 3=1 п г=1 + а2 , 1 = г О , 1фт N где (/, w) = j? Е f(xj) pi{xj). i=i С другой стороны, можно заметить, что по своему виду, в силу несмещенности, Xi( ij ? «) являются решениями системы линейных алгебраических уравнений г=1
Таким образом, мы обнаруживаем, что распределение (5) является Д2-распределением по отношению к мере цп, где \х - мера, сосредоточенная равномерно в точках х\,..., х .
Исследование вопроса о дисперсии рандомизированных оценок МНК Также как и в пункте 1.1., возможны два способа проведения эксперимента, в которых оценки метода наименьших квадратов имеют разные дисперсии:
1) Активный эксперимент: точки ж ,..., выбираются согласно плотности (5) и в каждой из них наблюдаем значения ,-.. = (ж ), j: = 1, п. Если взять новые точки x{i,..., х{ , то наблюдаемые значения . не будут зависеть от раньше полученных ..
2) Пассивный эксперимент: рандомизация производится после того, как мы уже получили все N значений функции / (bootstrap).
В первом случае, если получить N\ серий (т.н. репликаций) точек то получаем следующий результат для среднего оценок 1-ого регрессионного коэффициента Q ([16], с.5): есть систематическая погрешность аппроксимации / линейной комбинацией функций из набора (pi,..., рп.
Более того, если при каждом получении точек х\ ,..., х\ (к = 1, Ni) мы измеряем в них iV"2 раз, то получим:
Во втором случае у нас имеются наблюдения і, .., &v в точках х\,..., х . Пусть {(fii}f=l - .D-регулярная ортонормированная система функций, а ii,...,in и г1?... ,гп обозначают две независимые реализации дискретного распределения А2. Тогда известен следующий результат ([16], с.9):
Условные плотности как смеси распределений
Метод моделирования распределения Д2, предложенный в этой работе, основан на использовании условных плотностей. Чтобы получить их явный вид, воспользуемся хорошо известным лапласовским разложением определителя (см. [1]), а также вспомогательным утверждением относительно совместной плотности всевозможных наборов точек (xi,...,xr), г =Т/п:
Лемма 2 (Разложение Лапласа) Пусть D - определитель п-ого порядка. Фиксируем любые т его строк i\,... ,im, т п. Тогда Іі - Іт П где Miu„.tim-ju„.jm - минор D, получающийся при пересечении строк с номерами ii,...,im и столбцов с номерами ji,...,jm, a - гь...,гт;я,...,;? т алгебраическое дополнение М{ъ_ т- ъ_ т.
Утверждение леммы верно, если вместо т строк фиксировать т столбцов. Лемма 3 Пусть функции ь 2, п образуют ортонормированную систему в пространстве L2(X,//); где X - измеримое множество, а ц, - мера Лебега. Обозначим 1 / \2 pn(xi,. ..,хп) = — (det \\ Pi(xj)\\lJ=1) , а Pn-k \Xl, 5 ХП—к) = I Pn(a?i,..., xn)dxn-k+i... UiXfi X Тогда Рп-.к{хЪ...,Хп_к)=— 2 \\ is j)\\ns,j=l) 1 Н ... in-k n
Доказательство: Воспользуемся разложением Лапласа для определителя Рп —Рп{х\...хп), фиксируя первых п — к его столбцов: l ii ... jn_fc n где ji,..., jk - дополнение набора индексов {ц ... гп-к} до {1, 2,... п}. Обозначив fdet /2is( )"jii) =: Aiv„in_k: будем иметь nlpn-k(xi,..., хп-к) = п\ ] рп(осъ..., xn)dxn-k+i - dxn = X / An...in_fc (detyjim(a;r) =ljr n_fc+1) ) dxn-k+i dxn J \1 H — n-fc J I к - -2 dxn-k+i -dx E Дь.. -. I E(-i)"+VA(xn)A ,i)..J, l H ... in_fc n V=l x (0 - Л \\т.(гг\ т. (rr\ in. ( r \ in . (n- \ П s=n—k+l где Aiwfe = det I\Vh(xs), , Vjt-iM,4 ji+x{xa), , Pi (se)
Проинтегрируем сначала по xn. Возведем в квадрат подынтегральное выражение и воспользуемся тем, что функции pi,tp2,- ,Фп образуют ортонормированную систему в пространстве L2(X). Благодаря этому, попарные произведения в подынтегральном выражении, содержащие множители (pi(xn)(fj(xn), г ф j, обратятся в ноль. Более того, определитель Aii„.in_fc под знаком внешней суммы не зависит от хп. Поэтому его можно вынести за знак интеграла. Наборов (ji,... :jk) С„ штук. Каждому из них соответствует к слагаемых, СОСТОЯЩИХ ИЗ ЧЛеНОВ С Элементами (Si, . . . , Sfc-l) С (ji,..., Jfc), в подынтегральном выражении. Следовательно, общее количество членов становится равным кС%. Из них ровно С% различных. Таким образом, пользуясь ортогональностью функций /?i, /?2, , Фп в L2(X), получим, что pn-k(xi,..., хп-к) равняется
Сосчитаем теперь трудоемкость моделирования распределения А2, если последовательно моделируются условные плотности (12). Эта схема включает в себя следующие подходы: а) Метод обращения: представляем A2(Q) = РІ(ЖІ)Р2(Ж2ІЖІ)РЗ(ЖЗЖІ,Ж2) - -Рп(ж„жі,Ж2, - - - ,Жп-і) б) Метод композиций: каждая из плотностей Рі(хі),р2{х2\хі),рз{ха\хі,Х2),..., pn(xn\xi,X2...,xn-.i) моделируется как смесь распределений с плотностями А2. . «1 2—«к V Д2 (mi ... mfc_1)c(ii,...,ifc) к = 1,п (13) где Vitixi) ... (fi Xk-i) Рь(хк) д2. . »1»2- к Vi (a?l) ... tfc(»fc-l) Рік(Хк) есть функция только от хк (поскольку х\,..., Xk-i уже промоделированы и, следовательно, нам известны). Это позволяет записать (13) в следующем виде E(-i)j+1cgji2)...iife}\{ii} ,( ) т Д2 (mi ... mfc_i)c(ii,...,ifc) к = 1,п, (14) ГДЄ Коэффициенты С г; г- г }\{г} ОТНОСЯТСЯ К fc-ОЙ УСЛОВНОЙ плотности (верхний индекс) и представляют определитель (к — 1) го порядка, элементами которого являются значения функций й, , у_1} Уік в точках жі,..., Xk-i. в) Метод отбора: каждая из последних плотностей моделируется методом отбора. Поэтому сложность ее моделирования равна ( max 1 ІІ(ЖА;)1 ) Cfc, где Ck суммарная сложность нахождения коэффициентов С/, . wax- Покажем, что эти коэффициенты получаются итеративно и укажем процедуру их нахождения. Заметим, что в некоторых частных случаях более удобным может оказаться моделирование плотностей (14) методом обращения. Однако в общем случае такая схема рассуждений неприменима и, поэтому, нами не рассматривается.
Использование условных плотностей
В пунктах 3.1, 3.2 и 3.3 изложены три схемы моделирования дискретного распределения А2. При обращении нужно сосчитать С определителей п-го порядка, а при отборе моделируется набор точек ж ,... ,жг-п, в котором проверяется условие остановки Неймана с грубо завышенной оценкой максимума определителя. Что касается использования условных плотностей, нужно промоделировать п дискретных распределений, вероятности которых связаны итеративно.
Сравним полученные трудоемкости (21), (22) и (26) методов обращения, отбора и условных плотностей соответственно. Напомним, что это сравнение ведется в зависимости от значений трех параметров модели - N (числа точек, вкоторых сосредоточена дискретная мера /І), п (количества ортонормированных функций) и R (числа повторений). Заметим, что в (21) определяющим является член С 2п+п+6, в (22) - это пп, а в (26) - (n + 4)2re_1. Таким образом, с ростом N (при фиксированном п) трудоемкость метода обращения заметно возрастает, в отличии от трудоемкостей методов отбора и условных плотностей, которые слабо зависят от N. Так, например, при п = 5, R = 1 и Ср = 10 методом обращения выгоднее пользоваться (по сравнению с методом отбора) до N = 19. Если же п = 10, R = 1, а Ср = 10, то обращением стоит пользоваться до N = 54. Эти результаты получены в предположении, что мы не учитываем множитель МАХ2п в (22). В обоих случаях, однако, наименьшее количество операций совершается при использовании условных плотностей.
С ростом N использование метода отбора предпочтительнее методу обращения. Например, при N = 1000, n = 5, R = 1 и Cv = 10: Cinv = 0(1011), a Crej = О(106). Этого и следовало ожидать, поскольку при "больших"N дискретная модель аппроксимирует непрерывную, в которой и используется отбор (в той или иной форме). Стоит отметить, что в этом же случае C d « 302. Это и иллюстрирует огромное преимущество использования условных плотностей при больших N и R.
Пример 5 Активный эксперимент (одномерный случай) Пусть N = 10, X = {xi}, Х{ = і, /Л(ХІ) = 1/Ю, і = 1,10 п = 3. В качестве ортонормированной системы функций на L2(X,fi) возьмем !( ) = 1 рг(х) = -щ (х - у) рзМ = V Ji(:"2 - ПХ + 22) Рассмотрим 3 выборки (т.е. считаем Ni = 3) = 179,3 41} = 180 6 $] = 182 7 d1} = 175,1 = 182,5 Й1} - 178,4 Й1} = 175,8 $ = 178,3 1} = 178,8 $ = 180,3 р = 177,1 2) = 182,9 2) - 183,9 2) = 179, 5 2) = 177,6 U2) = 176,2 2) = 177,8 Й2) = 175,2 = 179,1 $ = 179,2 f) = 177,6 ф = 182,7 3) = 177,3 f3) = 181,3 3) = 180,0 Й3) = 179,1 3) = 181,3 3) = 183,7 Й3) = 180,2 3) = 180,3 из генеральной совокупности с нормальным распределением N(a, а2), где а = 180, а2 = 2,25.
Разделение ошибок в стандартной линейной полиномиальной регрессии на X = {1,2,..., 10}2: N = 100, п = 6, JVi = 10, N2 = 3 (активный эксперимент)
Пример 7 Пассивный эксперимент (двумерный случай) Проведем разделение ошибок пользуясь результатом (9) Ермакова и Швабе. Рассмотрим случай Примера 5, а именно, пусть N = 100, X = {(xi',Xj)}(, ХІ = i, }j,{(xi,Xj)} = 1/100, г, j — 1,10, п = 6. В качестве ортонормированной системы функций на L2(X, ц) возьмем опять
Моделирование обобщенного распределения А2 и оценка его трудоемкости
Условные плотности обобщенного распределения А2, также как и в случае т — п, имеют вид смеси. Однако число распределений в ней равно ьк Y1 Сп_кС1т, что существенно увеличивает сложность вычислений, когда 1=ак тип велики. Отметим, что коэффициенты этой смеси считаются легко, поскольку это биномиальные коэффициенты, разделенные на условную плотность (п — к — 1)-го порядка, которая подсчитана на предыдущем шаге. Поэтому алгоритм моделирования обобщенного распределения А2 будет представлять лишь модификацию Алгоритма 1, а именно:
Алгоритм 2 Моделирование обобщенного распределения А2 Для к = 1,п: 1. Пересчитываем коэффициенты с\ ... ifc_x по формуле (15). 2. Моделируем дискретное распределение, определяющее смесь в (28). 3. Выбрав в 2. однозначно набор г15...,гг, моделируем методом отбора соответствующее ему распределение в смеси (28). Таким образом, мы промоделировали Хк.
Сосчитаем трудоемкость моделирования обобщенного распределения А2. Разница по сравнению с непрерывным А2 (Глава 1) заключается в том, что в смеси, получающейся для (п — &)-ой условной плотности, имеется ьк X) Сп-кРт слагаемых. Поэтому на моделирование этого дискретного i=a,k распределения, Например С ПОМОЩЬЮ ДИХОТОМИИ, ОТВОДИТСЯ log2 Y1 п-к т операций. Следовательно, сложность моделирования смесей всех п п Ьк п Ьк условных плотностей равна X) 1оё2 С1п_кС1т — log2 П 2 С1п_кС1т. п—к—1 /=а п—к=\ 1=а,к Величину под знаком логарифма, хоть и грубо, можно оценить сверху следующим образом: n bfc n n—k n П Еі-Л П С _Л П 2»- C 2 (С!І )П(29) n—A;=l i=afc n—&=1 —0 n—k—1 Следовательно, n—k=l /=0 С другой стороны, число операций для подсчета коэффициентов вышеуказанных смесей не превосходит E CUci, J2 СЙ]" 2К_к= J2 сТ0(п2«) = cf]0(n22") n—/s=l l=ak n—k=l 1=0 n—k=l
Что касается коэффициентов crJ { -,\и.\ они получаются так же, как и в Главе 1. Поэтому на их пересчет отводится 0{п2) операций. Поэтому, если Сф1 = max (рі(ж), общая трудоемкость моделирования обобщенного х±Х, t=l,n А2-распределения не будет превосходить Сдеп = (о(п2) + 0(п22")СІ?] + nlog2 СІ?1) Cn) v »
Определяющим для величины трудоемкости моделирования в этом общем случае является член О(n22п)Cm . Поэтому, как и следовало ожидать, с ростом тип, она заметно возрастает. Например, без учета С\п при п — 10, т = 3 имеем Сдеп & 3 105, при п — 10, т = 5 -Сдеп « 106, а при п = 20, m = 5 величина трудоемкости не превосходит СдЄп 4 109. В ряде частных случаев, когда оценку (29) можно заметно улучшить, сложность моделирования значительно меньше. Так, например, при п = 4, т = 3 имеем 4, 8, 7, и 4 распределения в смеси для условных
Следовательно, на моделирование дискретных распределений в смесях отводится log2(4 -8-7-4) & 10 операций, а на пересчет с і \\ил -не больше 16. Полагая, как и сверху, СІр1 = 1, и учитывая трудоемкость применения метода отбора при моделировании 4( 4 1 2 3), Рз(#зІі#2), Р2(#2І#і) и рі(хі), которая не превосходит 3\/3 « 5, получаем итоговую сложность моделирования А 3, не превосходящую 130. Важно еще заметить, что моделирование обобщенного распределения А2 усложняется с ростом размерности s множества X лишь в плане вычисления значений функций y?i,... 5 фт- Однако, это обстоятельство не имеет большого влияния на величину Сдеп Чтобы проиллюстрировать выигрыш использования Апт по сравнению с простейшим случайным поиском или случайным поиском с помощью сеток, рассмотрим простейший пример полиномиальной регрессии на отрезке X = [0,1], т = гг, в котором максимум определителя Апп известен ([2]).
