close

Вход

Забыли?

вход по аккаунту

?

Робастная регрессия с применением t-распределения и EM-алгоритма.

код для вставкиСкачать
68
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
№1
Робастная регрессия с применением t-распределения
и EM-алгоритма
Шведов А.С.
В pаботе pассматpивается линейная pегpессионная модель. EM-алгоpитм пpедставляет собой pаспpостpаненный подход к оценке паpаметpов таких моделей на основе общего пpинципа максимизации пpавдоподобия. Известно, что этот метод оценки паpаметpов является pобастным,
если ошибки независимы, одинаково pаспpеделены и имеют многомеpное
t-pаспpеделение.
В пpедыдущих pаботах такой подход к оценке паpаметpов pегpессионных моделей пpименялся лишь пpи условии, что ошибки имеют многомеpное t-pаспpеделение с числовым паpаметpом степеней свободы. В настоящей pаботе pассматpивается более общая ситуация, когда ошибки
могут иметь многомеpное t-pаспpеделение с вектоpным паpаметpом степеней свободы. Ненаблюдаемые величины в EM-алгоpитме пpи этом оказываются случайными матpицами.
На численных пpимеpах пpи pазличных pаспpеделениях ошибок исследованы пpеимущества такого подхода по сpавнению с методом наименьших квадpатов.
Ключевые слова: робастная регрессия; многомерное t-распределение; EM-алгоритм.
1. Введение
Регpессионные модели являются основным инструментом для выявления зависимостей между pазличными показателями пpактически во всех областях экономической науки. Однако классические и наиболее pаспpостpаненные методы постpоения
pегpессионных моделей не обладают свойством pобастности, что, в pяде случаев, может
пpиводить к невеpным pезультатам. Робастные методы в эконометpике известны достаточно давно (см., напpимеp, [11]). Но все же в настоящее вpемя, скоpее, можно говоpить об усилении тенденции к пpименению pобастных методов, а не об обязательном
тpебовании использования таких методов хотя бы наpяду с классическими для контpоля качества pезультатов.
Идея использовать pегpессионные модели, у котоpых ошибки имеют не ноpмальные (гауссовские) pаспpеделения, а t-pаспpеделения, возникает, даже если не связы____________________
Шведов А.С. ? д. физ.-мат. н., профессор кафедры математической экономики и эконометрики
НИУ «Высшей школы экономики», e-mail: ashvedov@hse.ru
Статья поступила в Редакцию в декабре 2010 г.
2011
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
69
вать этот подход с pобастностью пpоцедуp оценки паpаметpов. Еще в XIX в. ученые
знали «об опасностях, поpождаемых длинными хвостами функций pаспpеделения ошибок» (см. [2, с. 7]). (В настоящее вpемя чаще используется теpмин не «длинные хвосты»,
а «тяжелые хвосты».) Изначальное использование пpи анализе пpедположения о ноpмальности ошибок включает и пpедположение о легких хвостах функций pаспpеделения. Это может не отвечать существу дела и пpиводить к искажениям в выводах.
Одним из самых pаспpостpаненных подходов к моделиpованию тяжелых хвостов является использование t-pаспpеделения.
Хотя для линейных pегpессионных моделей, у котоpых ошибки имеют t-pаспpеделение, и не существует такой замкнутой и кpасивой теоpии, как для линейных pегpессионных моделей, у котоpых ошибки имеют ноpмальное pаспpеделение, можно говоpить и о пpеимуществах моделей с t-pаспpеделением. Так, фактически, pегpессионные модели, у которых ошибки имеют t-pаспpеделение, включают в себя в качестве
частного случая pегpессионные модели, у котоpых ошибки имеют нормальное pаспpеделение, поскольку пpи стpемлении числа степеней свободы к бесконечности t-pаспpеделения пеpеходят в ноpмальное pаспpеделение. И pезультаты в пpедположении, что
ошибки имеют t-pаспpеделение с достаточно большим числом степеней свободы, и в
пpедположении, что ошибки имеют ноpмальное pаспpеделение, оказываются практически неотличимыми. Наконец, пpи оценке методом максимального пpавдоподобия паpаметpов pегpессионных моделей, у котоpых ошибки имеют t-pаспpеделения, можно
использовать пpоцедуpы, обладающие свойством pобастности.
Неpедко статистические данные, по котоpым стpоится pегpессионная модель,
содеpжат pезко выделяющиеся наблюдения (outliers). Эти наблюдения существенно
отделены от основной части и не подчиняются общей стpуктуpе. В каких-то случаях
такие выбpосы являются пpосто следствием ошибок, допущенных пpи сбоpе или обработке инфоpмации, но могут отpажать и pеальные эффекты.
Пpи использовании многих общепpинятых пpоцедуp для оценки паpаметpов даже одно pезко выделяющееся наблюдение может оказать очень сильное и часто искажающее пpавильную каpтину действие. Это легко понять на пpимеpе выбоpочного
сpеднего или выбоpочной диспеpсии. То же относится и к методу наименьших квадратов пpи опpеделении коэффициентов в линейной pегpессионной модели.
Робастные пpоцедуpы оценки паpаметpов пpетендуют на то, чтобы давать хорошее соответствие общей стpуктуpе и пpи наличии pезко выделяющихся наблюдений, как и в случае, когда pезко выделяющиеся наблюдения отсутствуют. Выявленная таким обpазом стpуктуpа, в свою очеpедь, может быть использована для обнаружения pезко выделяющихся наблюдений даже пpи pаботе с многомеpными статистическими данными.
В какой меpе можно говоpить, что эти пpетензии соответствуют действительности? Существуют pазличные подходы к постpоению pобастных алгоpитмов. Иногда
pезко выделяющиеся наблюдения автоматически игноpиpуются. Для тех методов, которые изучаются в настоящей pаботе, вклад таких наблюдений только уменьшается.
Для каждого класса алгоpитмов слова «хоpошее соответствие общей стpуктуpе и пpи
наличии pезко выделяющихся наблюдений» наполняются своим содеpжанием.
Сpеди пpедшествующих pабот, в котоpых изучаются алгоpитмы того же класса,
что и у нас, назовем [7, 12, 19]. К перечисленным можно было бы добавить и интересную pаботу [13], однако в [8, с. 165] указывается на непpавильные выводы, имеющиеся
в этой pаботе. Подpобнее о «подводных камнях», возникающих, если включать в со-
70
№1
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
став аpгументов функции пpавдоподобия число степеней свободы t-pаспpеделения, см.
[8, 14].
Робастные алгоpитмы дpугих классов пpедставлены, напpимеp, в книгах [15, 18].
Применяется и байесовский подход (см., напpимеp, [9]). Из работ прикладной направленности назовем [17, 22].
Содеpжание настоящей pаботы следующее. В паpагpафе 2 на пpимеpе множественной pегpессии (наблюдения одномеpные, объясняющих фактоpов несколько) обсуждается связь M-оценок и метода наименьших квадpатов с итеpационно модифицируемыми весами. В паpагpафе 3 пpиводится описание EM-алгоpитма, специализированного
метода нахождения точки максимума именно функции пpавдоподобия. В паpагpафе 4
излагаются некотоpые pезультаты, относящиеся к оценке паpаметpов множественной
pегpессии с ошибками, имеющими одномеpное t-pаспpеделение. Объясняется, почему
применение EM-алгоpитма в данном случае дает pобастный метод оценки паpаметpов pегpессии. Также в этом паpагpафе пpиводятся pезультаты численного исследования по методу Монте-Каpло. В паpагpафе 5 устанавливаются две новые теоремы о матpичном гамма-pаспpеделении. Затем эти теоpемы используются в паpагpафе 6,
где pезультаты, изложенные в паpагpафе 4, обобщаются на случай многомерной
pегpессии (наблюдения многомеpные, объясняющих фактоpов несколько). Пpи этом
ошибки имеют t-pаспpеделение с вектоpным паpаметpом степеней свободы (введенное в [3, 4]). Пpием, пpименяемый в паpагpафе 6, когда в EM-алгоpитме в качестве ненаблюдаемых пеpеменных беpутся случайные матpицы, видимо, используется впервые. Также в этом паpагpафе пpиводятся pезультаты одного pасчета.
2. Робастность M-оценок
Рассмотрим обычную линейную регрессию
q
(1)
yi =
?x
i ab a
+ ei , i = 1,..., n .
a =1
Объясняющие переменные xi a считаются известными числами. Через y1 ,..., yn
обозначаются и одномерные наблюдения, и случайные величины, представляющие собой вероятностную модель для этих наблюдений. Предполагается, что случайные величины e1,..., en независимы, одинаково распределены, и каждая из них имеет функцию плотности
(2)
1 ? ei ?
j?? ?? ,
s ?s?
где j( x ) ? некоторая известная функция плотности; s > 0 ? масштабирующий множитель. (Как обычно, используется одно и то же обозначение ei и для случайной величины, и для аргумента функции плотности.) Задача состоит в нахождении параметров b1 ,...,b q и s .
(
)
(
)
Если ввести в рассмотрение q -мерные вектора xi = x i1 ,..., xiq ? и b = b1 ,...,b q ?,
штрих означает транспонирование, то сумму, входящую в правую часть (1), можно
обозначить xi?b .
2011
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
71
Для определения вектора b может быть использован метод наименьших квадратов с весами, когда оценка вектора b строится путем минимизации выражения
n
? w (y
(3)
i
i
- xi?b )2 ,
i =1
где w1,..., wn ? заранее выбранные положительные числа. В частности, при w1 = ... = wn = 1
данный метод является обычным методом наименьших квадратов. (Мы сейчас не касаемся теоретических свойств метода наименьших квадратов с весами. Подчеркнем
только, что речь не идет об обобщенном методе наименьших квадратов, см., например, [1, гл. 5], хотя там и возникают сходные уравнения.) Приравнивание к нулю
частных производных функции (3) по b1 ,..., b q , что является необходимым условием минимума, дает систему уравнений
n
?x
(4)
i a wi
( yi - xi?b ) = 0 ,
a = 1,..., q .
i =1
С другой стороны, оценка параметров b и s может быть произведена методом
максимального правдоподобия, когда ищется максимум функции
n
(5)
- n log s +
? y i - xi?b ?
?? .
s ?
? log j???
i =1
Если ввести в рассмотрение функцию
w( x ) = -
(6)
1 j?( x )
,
x j( x )
то приравнивание к нулю частных производных функции (5) по b1 ,...,b q дает систему уравнений
n
?x
(7)
i =1
? yi - xi?b ?
??( y i - xi?b ) = 0 , a = 1,..., q .
? s ?
i a w?
?
Уравнения (7), хотя и похожи на уравнения (4), отличаются от них тем, что веса
зависят от искомого параметра b . Введем в рассмотрении n ? q матрицу X , i -я строка которой ? это x?i , и диагональную n ? n матрицу W , у которой i -й элемент на главной диагонали ? это wi . Тогда уравнения (4) можно записать в форме X ? W ( y - Xb ) = 0 ,
где y = ( y ,..., y )? . Отсюда
1
(8)
n
b = ( X ? WX ) -1 X ? W y ,
если матрица X ? WX невырожденная. Этот же прием может быть использован и для
решения уравнений (7).
72
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
№1
Предположим, что функция w ( x ) обладает следующими свойствами. Во-первых, она принимает только неотрицательные значения. Во-вторых, эта функция монотонно не убывает при x < 0 и монотонно не возрастает при x > 0 . В-третьих, w( x )
стремится к нулю и при x ® ? , и при x ® -? . Тогда можно говорить о робастности
оценок параметра b, полученных при помощи системы уравнений (7), поскольку резко выделяющимся наблюдениям yi , как правило, соответствуют большие по абсо? y - xi?b ?
?? .
лютной величине разности yi - xi?b и, соответственно, малые веса w?? i
? s ?
1 -x2 / 2
Нетрудно увидеть, что если j( x ) =
e
? функция плотности стандартно2p
го нормального распределения, то w( x ) ? 1 . И в этом случае система уравнений (7) не
приводит к робастным оценкам параметра b. Отсюда возникает идея рассмотреть
так называемые М-оценки, которые включают в себя в качестве частного случая оценки максимального правдоподобия (см., например, [2, 6]). В этом случае функция w( x ) ,
используемая в системе уравнений (7), не обязательно связана с функцией j( x ) соотношением (6). Зато можно потребовать, чтобы функция w ( x ) обладала тремя перечисленными выше свойствами. Или даже более сильными свойствами, например, обращалась в ноль при достаточно больших по абсолютной величине x .
Аргументация в пользу М-оценок может быть и такой. Если функция j( x ) на
практике все равно не известна, то почему нужно начинать с выбора функции j( x ), а
не с выбора функции w ( x ) ? Может быть, правильнее начинать с выбора функции
w ( x ), а функцию j( x ) определять из уравнения (6)?
Но мы все же начнем с выбора функции j( x ), а не с выбора функции w ( x ).
При любом a > 0 можно рассмотреть функцию
(9)
j( x ) =
G (a + 0,5) ?
x 2 ??
?1 +
?
?
2pa
G ( a ) ? 2a ?
1
- a - 0 ,5
? функцию плотности t -распределения с 2a степенями свободы. Известно, что
при a ® ? функция j( x ) переходит в функцию плотности стандартного нормального распределения. Нетрудно увидеть, что если функция j( x ) задается формулой (9),
2a + 1
то определяемая соотношением (6) функция w ( x ) имеет вид w( x ) =
. И в этом
2a + x 2
случае можно ожидать, что полученная путем решения системы уравнений (7) оценка
параметра b будет обладать свойством робастности (даже если не переходить от оценок максимального правдоподобия к более общим М-оценкам), поскольку w ( x ) стремится к нулю и при x ® ? , и при x ® -? .
Если записать систему уравнений (7) в виде (8), то W ? это диагональная
? y - xi?b ?
n ? n матрица, у которой i -й элемент на главной диагонали равен w? i
? . Ре? s ?
2011
73
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
шить уравнение (8) можно попытаться методом простой итерации. Выбрав начальное
приближение b ( 0 ) , например, при помощи метода наименьших квадратов, т.е. решив
уравнение (8) с единичной матрицей W , затем принимаем
(
b ( r +1) = X ? W ( r ) X
(10)
)
-1
X ? W (r) y ,
где W (r ) ? это диагональная n ? n матрица, у которой i -й элемент на главной диа? y - xi?b ( r ) ?
?.
гонали равен w?? i
?
(r)
? s
?
Приравнивание к нулю частной производной функции (5) по s дает выражение
(11)
s2 =
n
1
?(y
n
i
i =1
? y - xi?b ?
??( y i - xi?b ) ,
- xi?b ) w?? i
? s ?
что равносильно соотношению
s2 =
1
( y - Xb )?W ( y - Xb ) .
n
Последнее уравнение позволяет находить значения s итерациями:
(12)
2
s( r +1) =
1
n
(y - Xb )?W (y - Xb ) .
( r +1)
(r)
( r +1)
2
Для определения s( 0 ) может быть использована единичная матрица W .
Хорошо известно, что последовательность значений, получаемых при помощи
метода простой итерации, может быть как сходящейся, так и не быть сходящейся.
На практике этот метод обычно применяют для небольших расчетов, руководствуясь принципом «раз сошлось, значит, решение получено». Хотя и существуют условия, гарантирующие сходимость метода простой итерации. Кроме того, в правой
части (12) стоит не b (r ) , а b ( r +1) , т.е. в данном случае производится некоторое усложнение метода простой итерации.
В параграфе 4 показано, что для множественной линейной регрессии, когда
ошибки имеют t -распределение, применение EM-алгоритма приводит к тому же
2
итерационному процессу (10), (12) для определения b ( r +1) , s ( r +1) . А тогда применимы теоремы о сходимости итерационного процесса, построенного на основе EMалгоритма. В рамках EM-алгоритма также можно провести обобщение на случай,
когда наблюдения y1,..., y n не одномерные, а m -мерные (см. параграф 6).
Более подробно связь метода наименьших квадратов с итерационно модифицируемыми весами и робастных процедур освещается, например, в работе [21].
74
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
№1
3. EM-алгоритм
EM-алгоритм предназначен для поиска точки, в которой достигает максимума
функция правдоподобия, путем построения некоторого итерационного процесса. Каждый шаг итерационного процесса состоит из двух подшагов. E-подшаг заключается в
нахождении ожидания (expectation) некоторой функции от случайных величин. При
этом ожидание само оказывается функцией интересующего параметра. M-подшаг ?
это максимизация (maximization), определение того значения параметра, при котором
данная функция достигает максимума. Первые буквы приведенных английских слов и
дают название алгоритма. Обзор различных задач, для решения которых применяется
EM-алгоритм, можно найти, например, в работе [16].
Пусть y = ( y1,..., yn ) ? это набор наблюдений, вообще говоря, многомерных.
z = (z1 ,..., zn ) ? набор ненаблюдаемых величин также, вообще говоря, многомерных.
С одной стороны, предполагается, что ненаблюдаемые величины заметно влияют на
наблюдаемые, и привлечение их для анализа отвечает существу дела. С другой стороны, нахождение точек максимума функций в рамках EM-алгоритма может оказаться
значительно более простой и надежной с вычислительной точки зрения процедурой,
чем непосредственное нахождение точки максимума исходной функции правдоподобия. В каких-то задачах не вызывает сомнений, что именно следует взять в качестве
ненаблюдаемых величин. В других задачах ответ на этот вопрос не столь очевиден.
Обозначим через h совместную функцию плотности случайных векторов y и
z, через g ? условную функцию плотности случайного вектора z при заданном y ,
через f ? маргинальную функцию плотности случайного вектора y. Все эти функции плотности считаются зависящими от некоторого параметра q, вообще говоря,
многомерного. Имеет место соотношение
f ( y ; q) =
h( y , z; q)
.
g ( z | y; q)
Переходя к логарифмам, получаем соотношение для логарифмических функций
правдоподобия
(13)
l (q | y ) = l (q | y , z ) - log g ( z | y; q) .
Задача состоит в нахождении точки q, в которой достигает максимума функция l (q | y ) . Пусть q(r ) ? значение параметра q, найденное при r -й итерации. Умножим левую и правую части (13) на g ( z | y; q( r ) ) и проинтегрируем по z. Введем обозначения
U (q | y; q( r ) ) = l (q | y , z ) g ( z | y; q( r ) ) dz ,
?
r(q | y; q( r ) ) = log g ( z | y; q) g ( z | y; q( r ) ) dz .
?
Тогда
(14)
l (q | y ) = U (q | y; q(r ) ) - r(q | y; q(r ) ) .
2011
75
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
Заметим, что
?
log
g ( z | y ; q( r ) )
g ( z | y; q( r ) ) dz
g ( z | y ; q)
? это расстояние Кульбака ? Лейблера между функциями плотности g ( z | y; q( r ) ) и
g ( z | y; q) , которое, как известно, всегда неотрицательно. Поэтому при любом q
(15)
r(q | y; q(r ) ) ? r(q( r ) | y; q( r ) ) .
E-подшаг состоит в нахождении ожидаемого логарифмического правдоподобия
U (q | y; q(r ) ) .
M-подшаг состоит в нахождении точки
q(r +1) = arg max U (q | y; q(r ) ) .
q
Из (14), (15) и способа определения точки q( r +1) следует, что
(16)
l (q( r +1) | y ) ? l (q( r ) | y ) .
Соотношения (16) показывают, что движение идет «в правильном направлении»,
но еще не гарантируют, что последовательность q(r ) сходится. Условия общего характера, из которых следует сходимость этой последовательности к точке максимума
функции l (q | y ), даются в работе [20] теоремами 1 и 4.
4. Множественная линейная регрессия
В этом параграфе рассматривается набор одномерных наблюдений y1,..., y n , и
будем предполагать, что ненаблюдаемые величины z1 ,..., zn также одномерные и, кроме
того, положительные. Двумерные случайные величины ( y1 , z1 ) ,?, ( yn , zn ) считаем независимыми.
Предположим, что совместная функция плотности случайных величин yi и zi
имеет вид
(17)
? z e2 ?
exp? - i i ? g (zi ) ,
? 2s 2 ?
s 2p
?
?
z1i / 2
где ei = yi - xi?b в соответствии с (1). Будем использовать обозначение q для пары b, s .
Тогда h( y , z; b, s) ? это произведение функций (17) при i = 1,..., n . Следовательно,
n
1
log h ( y, z; b , s) = - log 2p - n log s +
2
2
n
?
i =1
log zi -
1
2s 2
n
?
i =1
n
zi ei2 +
? log g (z ) .
i
i =1
76
№1
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
Пренебрегая слагаемыми, не зависящими ни от b, ни от s, получаем выражение для
логарифмической функции правдоподобия
l (q | y , z ) = - n log s -
(18)
1
2s 2
n
2
i i
?z e
.
i =1
В соответствии с определением, данным в параграфе 3,
(
) (
)
U q | y ; q( r ) = E l ( q | y , z ) | y ; q( r ) ,
где l (q | y , z ) рассматривается как функция случайного вектора z. Из (18) следует,
что эта функция линейна относительно z1 ,..., z n .
Распределение вероятностей с совместной функцией плотности (17) называется
нормальным-гамма распределением, если
g ( z) =
(19)
Aa
G( a )
exp(- Az ) z a -1 ,
где a > 0 , A > 0 . Тогда, как нетрудно увидеть, (выкладки для m -мерного случая приводятся в параграфе 6) маргинальная функция плотности случайной величины yi
имеет вид (2), где
j( x ) = ( 2p) -1 / 2
(20)
Aa
G(a + 0,5)
(A + 0,5x )
2
G( a )
a + 0 ,5
.
Из (6) следует, что в этом случае
w( x ) =
(21)
a + 0,5
A + 0,5x 2
.
Введем обозначение ei( r ) = y i - xi?b ( r ) и воспользуемся тем, что
? 1
?
E zi | yi ; q( r ) = w ? ( r ) ei( r ) ? .
?s
?
(
(22)
)
Доказательство соотношения (22) для m -мерного случая приводится в параграфе 6. Для одномерного случая (22) доказано в работе [7], и в этом доказательстве не
требуется, чтобы функция g (z ) обязательно имела вид (19). Из (18) и (22) следует, что
(
)
U q | y; q( r ) = - n log s -
n
1
2s
2
? 1
? w ?? s
i =1
(r)
?
ei( r ) ? ( y i - xi?b )2 .
?
Таким образом, E-подшаг EM-алгоритма выполнен. Выполнение M-подшага
аналогично процедуре, описанной в параграфе 2 (ср. (3), (5), (7), (11), (12)).
2011
77
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
Пример 1. Сокращения М1 и М2 в этом примере используются для обозначения
одного из двух методов, применяемых для определения параметра b и стандартного
отклонения возмущений s.
М1 ? метод максимального правдоподобия в предположении нормальности возмущений.
М2 ? EM-алгоритм в предположении, что возмущения имеют t -распределение
с тремя степенями свободы.
Рассматриваются три различных вида сгенерированных рядов наблюдений.
Н0 ? наблюдения соответствуют модели с нормальными возмущениями.
Н1 ? наблюдения соответствуют модели с нормальными e -засоренными возмущениями, e = 0,1 .
Н2 ? наблюдения соответствуют модели с возмущениями, имеющими t -распределение с тремя степенями свободы.
Целью является исследовать поведение каждого из методов для «своих» и «чужих» рядов наблюдений. Для простоты мы ограничиваемся лишь возмущениями,
имеющими конечные дисперсии.
Пусть n = 15 , q = 1 , xi 1 = 1 при каждом i , 1 ? i ? n . Случайная величина ei имеет функцию плотности (2), где j( x ) ? либо функция плотности стандартного нормального распределения, либо функция плотности t -распределения при 2a = 3 ; обе эти
функции плотности приведены в параграфе 2. Для генерации наблюдений yi используются значения b = 1 , s = 0,3 . При использовании нормальных возмущений s = s .
При использовании возмущений, имеющих t -распределение с n степенями свободы,
s=s
n-2
.
v
При генерации ряда с e -засоренными наблюдениями считается, что с вероятностью (1 - e ) стандартное отклонение возмущения равно s, и с вероятностью e равно 5s.
Для каждого из трех видов генерируется L = 300 рядов наблюдений длины n .
Для l -го эксперимента, l = 1,..., L значения параметров b l и sl определяются и методом М1, и методом М2. Затем определяются средние значения
b=
1
L
L
?
bl , s =
l =1
1
L
L
?s
l
l =1
и среднеквадратические отклонения
?1
?
?L
?
L
1/ 2
?
2
? (bl - b ) ??
l =1
?
?1
, ?
?L
?
L
1/ 2
?
? (sl - s ) 2 ??
l =1
.
?
Результаты для средних значений приведены в табл. 1 и 2. В скобках даются
среднеквадратические отклонения.
78
№1
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
Таблица 1.
Результаты для параметра ?, истинное значение ? = 1
Н0
Н1
Н2
М1
1,011
(0,071)
1,016
(0,133)
1,011
(0,067)
М2
1,017
(0,074)
1,019
(0,081)
1,012
(0,050)
Таблица 2.
Результаты для параметра s, истинное значение s = 0,3
Н0
Н1
Н2
М1
0,273
(0,052)
0,472
(0,228)
0,244
(0,093)
М2
0,370
(0,079)
0,454
(0,130)
0,276
(0,070)
Сравнивая результаты для рядов вида Н0 и Н2, мы видим, что «свой» метод
(т.е. метод М1 для рядов Н0 и метод М2 для рядов Н2) дает лучшие результаты и
в смысле меньшего разброса (т.е. среднеквадратического отклонения), и в смысле близости среднего значения найденных параметров к истинным значениям (за исключением значения b = 1,012, которое несколько хуже, чем значение b = 1,011 ).
На первый взгляд, результаты для параметра b для рядов вида Н2 противоречат теореме Гаусса ? Маркова. Метод М1 совпадает с методом наименьших квадратов, и, казалось бы, дисперсия оценки должна быть наименьшей. А среднеквадратическое отклонение 0,067 существенно больше, чем среднеквадратическое отклонение
0,050, полученное при использовании метода М2. На самом деле, этот пример всего
лишь показывает, что требование линейности и несмещенности оценки, содержащееся
в теореме Гаусса ? Маркова, не может быть отброшено. Метод М2 не является линейным.
Для рядов вида Н1 в результатах для параметра b видна робастность метода
М2. Полученное среднеквадратическое отклонение примерно в 1,65 раза меньше,
чем для метода М1. Различие в средних, 1,016 и 1,019, носит случайный характер.
Так, для другой серии из 300 экспериментов для рядов вида Н1 для параметра b при
использовании метода М1 получены результаты
да М2 ? результаты
1,019
(0,085)
1,023
(0,140)
, а при использовании мето-
. Проявляется робастность метода М2 и в результатах
для параметра s для рядов этого вида.
2011
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
79
5. Две теоремы о матричном гамма-распределении
Теория матричных гамма-распределений излагается, например, в [10]. Наиболее
известными из этих распределений являются, видимо, распределения Уишарта. Также в работе [10, с. 122?124] рассматриваются матричные гамма-распределения с векторным параметром ? некоторое естественное обобщение распределений Уишарта.
При том, что легкими и короткими формулировки здесь быть не могут, обозначения,
используемые в [10] и других работах для распределений с векторным параметром, с
нашей точки зрения, несколько избыточны, возможно, из-за этого остались неисследованными свойства этих распределений. Более прозрачные обозначения для матричных гамма-распределений с векторным параметром используются в работе [3].
Здесь мы повторим только самые необходимые определения.
{ }im, j=1
Для m ? m матрицы C = ci j
при k = 1,..., m рассматриваются подматрицы
{ }ik, j =1 и C[k ] = {ci j }im, j =m-k +1 .
C [ k ] = ci j
Рассматривается также вектор a = (a1,..., am ) такой, что a j > 0,5( j - 1) при
j = 1,..., m. Многомерная гамма-функция определяется следующим образом:
m
Gm* ( a ) = p m( m-1) / 4
? G(a
j
)
- 0,5( j - 1) ,
j =1
где G(?) ? обычная гамма-функция. Дополнительно считается, что
a0 = 0 , am+1 = 0,5( m + 1) .
Параметрами рассматриваемых гамма-распределений являются вектор a указанного вида и положительно определенная m ? m матрица A. Функция плотности
имеет вид
m
(23)
g ( z ) = g a , A etr( - Az )
?| z
[ j]
a j - a j +1
|
,
j =1
где z ? положительно определенная m ? m матрица; etr (C ) = exp( tr C ) ; | C | ? определитель матрицы C . Коэффициент g a, A задается формулой
(24)
g a,A
m -1
?
?
a -a
= ? Gm* ( a)
| A[ m- j ] | j j +1 ?
?
?
j =0
?
?
-1
?
(ср. (23) и (24) с (19)).
Пусть T ? симметричная m ? m матрица такая, что матрица A - T положительно определенная. (Изначально симметричная матрица T может быть взята произвольно. При некотором e > 0 матрица A - eT будет положительно определенной. В этом
смысле условие, что матрица A - T положительно определенная, не является ограни-
80
№1
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
чительным.) Пусть Z ? положительно определенная случайная матрица с функцией
плотности (23). Определим производящую функцию моментов
M (T ) = E etr (T Z ) .
Из формулы
m
tr(T Z ) =
m
?
Tl l Z l l + 2
l =1
??T
k l Zk l
k =1 l > k
следует, что
¶M (T )
(25)
¶Tl l
( )
= E Zl l ,
T =0
а при l > k
¶ M (T )
(26)
¶Tk l
( )
= 2 E Zk l .
T =0
Теорема 1.
M (T ) =
g a,A
g a , A-T
m -1
=
?| ( A - T )
[ m- j ] a j - a j +1
|
j =0
?
?
?
?
m -1
?
[ m - j ] a j - a j +1 ?
?| A
|
j =0
?
?
-1
.
Доказательство следует из выражения
m
M (T ) = g a , A
?
etr(-( A - T ) z )
?| z
[ j]
|
a j - a j +1
dz
j =1
z >0
и из формулы (24).
При j = 0,..., m - 1 определим m ? m матрицу Cm- j такую, что
(Cm- j ) [m- j ]= (A[m- j ] )-1 ,
остальные элементы матрицы Cm- j равны нулю.
Теорема 2.
m-1
E(Z ) =
? (a
j +1 - a j
)Cm- j .
j =0
Доказательство. Чтобы найти ожидание каждого элемента Z k l , k ? l , случайной матрицы Z , воспользуемся формулами (25), (26) и теоремой 1.
2011
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
81
Зафиксируем j такое, что l ? m - j . Минор элемента Ak l - Tk l матрицы ( A - T )[ m- j ]
обозначим M k l . Через S обозначим определитель этой матрицы. Введем обозначение
sk l = Ak l - Tk l . Тогда
m- j
(27)
? (-1)
S=
u+ l
su l M u l .
u =1
Дифференцирование функции (27) по sl l не вызывает затруднений, поскольку
ни один из миноров от sl l не зависит. Имеем
(28)
¶S
¶ sl l
(
= M l l = ( A - T )[ m - j ] ( A - T )[ m - j ]
)
-1
ll .
Дифференцирование функции (27) по sk l при k < l несколько труднее. Все миноры кроме M l l могут зависеть от sk l , поскольку sl k = sk l . Имеем
(29)
¶S
¶ sk l
m- j
= (-1)k +l M k l +
? (-1)
u+ l
¶
su l
¶ sk l
u =1
u ?l
M ul .
Обозначим через M u l ,l v определитель матрицы, получающейся из ( A - T )[ m- j ]
выкидыванием строк с номерами u и l и столбцов с номерами l и v . Тогда при u < l
m- j
l -1
M ul =
?
( -1)v + l -1 sl v M u l ,l v -
v =1
? (-1)
v + l -1
sl v M u l ,l v ;
v = l +1
при u > l
m- j
l -1
M ul =
?
( -1)v +l sl v M u l ,l v -
v =1
? (-1)
v +l
sl v M u l ,l v .
v = l +1
Ни один из определителей M u l ,l v , входящих в две последние формулы, от sk l
не зависит. Вторыми суммами в правых частях также, очевидно, можно пренебречь,
поскольку k < l . Поэтому при u < l
¶
¶ sk l
M u l = (-1) k +l -1 M u l ,l k ;
при u > l
¶
¶ sk l
M u l = (-1) k +l M u l ,l k .
82
№1
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
Из (29) при k < l находим
¶S
¶ sk l
m- j
l -1
= (-1)k +l M k l +
?
( -1)u +l su l ( -1)k +l -1 M l k ,u l +
u =1
? (-1)
u+ l
su l (-1) k +l M l k ,u l .
u = l +1
То есть при k < l
(30)
¶S
¶ sk l
(
= 2 ( -1)k +l M k l = 2 ( A - T )[ m- j ] ( A - T )[ m- j ]
)
-1
kl .
Изменения, которые нужно внести в приведенные выкладки при m - j ? 2
или при l = m - j , очевидны.
Те же выражения (28) и (30) получаются при дифференцировании ( A - T )[ m- j ]
по Tl l и по Tk l , только в правые части добавляется знак минус.
Воспользовавшись теоремой 1, формулами (28) и (30), получаем
m-l
¶M (T )
¶Tl l
= -
? (a
- a j +1
j
) (A[m- j ] ) -l l1 ,
j =0
T =0
и при k < l
m- l
¶M (T )
¶Tk l
= -2
? (a
j
- a j +1
) (A[m- j ] ) -k1l .
j =0
T =0
Воспользовавшись (25) и (26), при k ? l получаем
m-1
E ( Zk l ) =
? (a
j +1 - a j
)(Cm- j )k l .
j =0
Теорема 2 доказана.
Теорема 2 для случая a1 = ... = am известна. Доказательство приводится, например, в [10]. При этом и в [10], и в других работах используется не производящая
функция моментов, а характеристическая функция. Использование производящей
функции моментов позволяет избежать рассмотрения матриц с комплексными элементами. Теорема 1, хорошо известная для одномерного случая, при m > 1 , по-видимому, является новой даже для случая a1 = ... = am .
6. Многомерная линейная регрессия
Рассмотрим уравнения, аналогичные (1):
(31)
yi = b? xi + ei , i = 1,..., n .
2011
83
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
Объясняющие переменные xi a ? известные числа, xi = xi1 ,..., xiq ? ; b ? q ? m
(
)
матрица. Через y1 ,..., yn обозначаются и m -мерные наблюдения, и случайные вектора,
представляющие собой вероятностную модель для этих наблюдений. Предполагается,
что m -мерные случайные вектора e1,..., en независимы и одинаково распределены.
Ненаблюдаемые величины z1 ,..., zn являются положительно определенными m ? m
случайными матрицами. Как и в одномерном случае, ( y1 , z1 ) ,?, ( yn , zn ) независимы.
Совместное распределение yi и zi будем считать нормальным-гамма, т.е. совместная
функция плотности имеет вид
(32)
(2p)- m / 2
1
s
2
1/ 2
zi
? 1 ? ?
exp?? ei zi ei ?? g (zi ) ,
? 2s 2
?
где g (zi ) определяется формулой (23), и ei = yi - b? xi в соответствии с (31).
Рассмотрим вектор b = (b1,..., bm ) , где b j = a j + 0,5 при j = 1,..., m . Пусть b0 = 0 ,
bm +1 = 0,5( m + 1) .
При x ? R m рассмотрим функцию
(33)
j( x ) = (2p )-m / 2
Gm* (b)
Gm* ( a)
A
-1 / 2
m -1
? 1 [ m- j ]? [ m- j ] -1 [ m- j ] ?
??1 + x
(A
) x
??
2
?
?
j =0
b j - b j +1
?
Она является функцией плотности многомерного t -распределения с векторным параметром степеней свободы (см. [3; 4]); ср. (20). Здесь для m -мерного вектора x через
x[k ] обозначается k -мерный вектор, состоящий из первых k компонент вектора x ,
k = 1,..., m .
Теорема 3. Маргинальная функция плотности случайного вектора yi имеет вид
(34)
?1 ?
j?? ei ?? ,
s ?s ?
1
m
где функция j задается формулой (33).
Доказательство. Маргинальная функция плотности случайного вектора yi получается интегрированием по области zi > 0 совместной функции плотности (32). Чтобы несколько сократить формулы, внутри доказательства теоремы будем использовать обозначение z вместо zi и обозначение e вместо ei . Во-первых, заметим, что
e?ze = tr (ee?z ) .
Используя (23) и (24), получаем следующее выражение для маргинальной функции
плотности:
84
№1
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
(2p )-m / 2 g a, A
? ?
? ?
1
etr?? - ?? A +
ee? ?? z ??
m
2
s z >0 ? ?
2s
? ?
1
-m / 2
= (2p)
?
ga ,A
m -1
1
s
G * ( b)
m m
?
j =0
m
?z
b j - b j +1
[ j]
dz =
j =1
?
?
1
?? A +
ee? ??
2
2s
?
?
[ m - j ] bi - b j +1
.
Если воспользоваться тем, что
A[ m- j ] +
?
?
1 [ m- j ]? [ m- j ] -1 [ m- j ] ?
?? ,
e[ m- j ]e[ m- j ] = A[ m- j ] ??1 +
e
(A
) e
2s
? 2s 2
?
1
2
(см., например, лемму 5 в [3]), то получаем выражение
*
(2p)- m / 2 Gm* (b)
1
Gm (a ) s m
A
-1 / 2
m -1
?
1 [ m- j ]? [ m- j ] -1 [ m- j ] ?
??
??1 +
e
(A
) e
2s 2
?
j=0 ?
b j - b j +1
?
.
Теорема 3 доказана.
Теорема 4. Условная функция плотности случайного вектора zi при условии
yi ? это функция плотности матричного гамма-распределения с векторным парамет1
ром b и с матричным параметром A +
ei ei ? .
2s 2
Доказательство. Условная функция плотности случайного вектора zi при условии yi ? это отношение совместной функции плотности (32) к маргинальной функции
плотности (34). Как и в доказательстве предыдущей теоремы, будем использовать обозначение z вместо zi и обозначение e вместо ei . Искомая условная функция плотности представима в виде дроби с числителем
? ?
? ?
1
etr?? - ?? A +
ee? ?? z ??
2
2s
? ?
? ?
m
?z
b j - b j +1
[ j]
j =1
и со знаменателем
m -1
Gm* (b)
?
A[ m- j ]
bi - b j +1
m -1
?
1 [ m- j ]? [ m- j ] -1 [ m- j ] ?
??
??1 +
e
(A
) e
2s 2
?
j=0 ?
b j - b j +1
?
j =0
.
Используя, как и в доказательстве предыдущей теоремы, лемму 5 из [3], получаем
следующее выражение для знаменателя:
m -1
Gm* (b)
[ m- j ]
?A
j =0
Теорема 4 доказана.
+
1
2s 2
e
[ m - j ] [ m - j ]?
e
b j - b j +1
-1
?
?
? .
= ?g
1
? b , A + 2 ee? ?
2s
?
?
2011
85
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
Будем использовать обозначение q для пары b , s . Из (32) следует, что и в многомерном случае логарифмическая функция правдоподобия имеет вид, сходный с (18):
(35)
n
1
l (q | y , z ) = - n log s m -
2s
2
? e ?z e
.
i i i
i =1
Ввиду линейности l (q | y , z ) по z1 ,..., zn , чтобы построить функцию U (q | y; q(r ) ) ,
достаточно знать условное ожидание E ( zi | yi ; q(r ) ) .
При j = 0,..., m - 1 и при x ? R m определим m ? m матрицу Cm- j (x ) такую, что
-1
(Cm- j ( x ))
[ m- j ]
[ m- j ] ?
??
1 ?
? ,
= ? ?? A + xx ? ??
?
??
2
?
?
?
остальные элементы матрицы Cm- j (x ) равны нулю. Положим
m -1
? (b
w( x ) =
j +1 - b j
)Cm- j ( x )
j =0
? 1 ( r) ?
(ср. (21)). Тогда на основании теорем 2 и 4 получаем E zi | yi ; q( r ) = w??
ei ?? (ср. (22)).
? s( r )
?
Определим m ? m матрицы
(
)
? 1 ( r) ?
wi( r ) = w??
ei ?? , i = 1,..., n .
? s( r )
?
В соответствии с определением функции U (q | y; q(r ) ) в параграфе 3 (см. также
параграф 4) и с (35) получаем
U (b, s | y; b( r ) , s( r ) ) = -nm log s -
n
1
2s 2
? ??? y ? - x ?b ??? w
i
i
(r)
i
( yi - b? xi ) .
i =1
Дифференцирование функции U по ba k , a = 1,..., q , k = 1,..., m и приравнивание
производной к нулю с учетом симметричности матрицы w дает уравнение
n
(36)
m
?x ?
ia
i =1
j =1
?
wi(,rj)k ? yi j
?
q
-
?
?x
?
?=0.
i gb g j ?
g =1
?
Дифференцирование функции U по s и приравнивание производной к нулю
дает уравнение
(37)
s2 =
1
n
?? y ? - x ?b ?? w
?
?
?
nm
i
i =1
i
(r)
i
( yi - b? xi )
86
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
№1
(ср. (11), (12)). Найденные из уравнений (36) и (37) значения b и s принимаются в качестве b( r +1) и s( r +1) . Значения b( 0 ) и s(0 ) рассчитываются с единичными матрицами wi .
Другой вариант применения EM-алгоритма в линейной регрессии, когда ошибки
имеют многомерное t -распределение, представлен, например, в работе [12].
Пример 2. Пусть m = 3 , n = 1000 , q = 2 . При каждом i = 1,..., n вектор xi ? = (1, i ) ,
а трехмерная случайная величина ei имеет функцию плотности (33), (34) с параметрами a = (2,5,9) ,
? 1 1,3 1,5 ?
?
?
A = ?1,3 2 2 ? ,
?1,5 2 4 ?
?
?
s = 12 . При генерации трехмерных возмущений с указанным t -распределением применяется алгоритм Метрополиса (см., например, [5]). Для генерации наблюдений используется матрица
?10 20 - 30 ?
?? .
b = ??
? 0,3 - 0,2 0,4 ?
Проводя аналогию с временными рядами, можно сказать, что рассматривается модель со свободным членом и с трендом, а возмущения представляют собой трехмерный белый шум, имеющий t -распределение с векторным параметром степеней
свободы.
При применении EM-алгоритма, т.е. при использовании итерационного процесса,
основанного на формулах (36), (37), сходимость в проведенном эксперименте была достигнута после 20 итераций. Получены значения параметров
?10,185
b = ??
? 0,299
20,268
- 29,681?
?,
- 0,201
0,400 ??
s = 11,533 . Результаты являются удовлетворительными. Отметим, что при применении метода наименьших квадратов получено значение параметра
?10,191
b = ??
? 0,299
20,399
- 29,634 ?
?.
- 0,201
0,400 ??
Оно же использовалось в качестве начального значения b ( 0 ) в EM-алгоритме.
В данном расчете метод наименьших квадратов уступает EM-алгоритму.
2011
ЭКОНОМИЧЕСКИЙ ЖУРНАЛ ВШЭ
87
* *
*
С ПИ С ОК Л И Т Е Р А Т У Р Ы
1. Магнус Я.Р., Катышев П.К., Пересецкий А.А. Эконометрика. Начальный курс. М.:
Дело, 2004.
2. Хьюбер П. Робастность в статистике. М.: Мир, 1984.
3. Шведов А.С. Бета-pаспpеделение случайной матpицы и его пpименение в модели
состояние-наблюдение: пpепpинт. WP2/2009/01. М.: ГУ ВШЭ, 2009.
4. Шведов А.С. t-pаспpеделение случайной матpицы и его пpименение в pегpессионной
модели: пpепpинт. WP2/2010/01. М.: ГУ ВШЭ, 2010.
5. Шведов А.С. О методах Монте-Каpло с цепями Маpкова // Экономический журнал Высшей школы экономики. 2010. Т. 14. № 2. С. 227?243.
6. Andrews D.F. A Robust Method for Multiple Linear Regression // Technometrics. 1974.
16. Р. 523?531.
7. Dempster A.P., Laird N.M., Rubin D.B. Iteratively Reweighted Least Squares for Linear
Regression When Errors are Normal/Independent Distributed // Multivariate Analysis ? V / ed.
by P.R. Krishnaiah. Amsterdam: North-Holland, 1980. Р. 35?57.
8. Fernandez C., Steel M.F.J. Multivariate Student-t Regression Models: Pitfalls and
Inference // Biometrika. 1999. 86 (1). Р. 153?167.
9. Fonseca T.C.O., Ferreira M.A.R., Migon H.S. Objective Bayesian Analysis for the
Student-t Regression Model // Biometrika. 2008. 95 (2). Р. 325?333.
10. Gupta A.K., Nagar D.K. Matrix Variate Distributions. N.Y.: Chapman & Hall, 1999.
11. Koenker R. Robust Methods in Econometrics // Econometric Reviews. 1982. 1. Р. 213?
255.
12. Lange K.L., Little R.J.A., Taylor J.M.G. Robust Statistical Modelling Using the t-distribution // Journal of the American Statistical Association. 1989. 84. Р. 881?896.
13. Liu C.H., Rubin D.B. ML Estimation of the t-distribution Using EM and its Extensions, ECM and ECME // Statistica Sinica. 1995. 5. Р. 19?39.
14. Lucas A. Robustness of the Student-t Based M-estimator // Communications in
Statistics ? Theory and Methods. 1997. 26 (5). Р. 1165?1182.
15. Maronna R.A., Martin R.D., Yohai V.J. Robust Regression ? Theory and Methods.
N.Y.: Wiley, 2006.
16. Meng X.L., van Dyk D.A. The EM Algorithm ? An Old Folk-song Sung to a Fast
New Tune (with discussion) // Journal of the Royal Statistical Society. 1997. B. 59. Р. 511?567.
17. Preminger A., Franck R. Forecasting Exchange Rates ? A Robust Regression
Approach // International Journal of Forecasting. 2007. 23(1). Р. 71?84.
18. Rousseeuw P.J., Leroy A.M. Robust Regression and Outlier Detection. N.Y.: Wiley,
1987.
19. Rubin D.B. Iteratively Reweighted Least Squares // Encyclopedia of Statistical
Sciences. N.Y.: Wiley, 1983. Vol. 4. Р. 272?275.
20. Wu C.F.J. On the Convergence Properties of the EM Algorithm // Annals of
Statistics. 1983. 11. Р. 95?103.
21. Yuan K.-H., Bentler P.M. Robust Mean and Covariance Structure Analysis through
Iteratively Reweighted Least Squares // Psychometrika. 2000. 65(1). Р. 43?58.
22. Zellner A., Ando T. Bayesian and Non-Bayesian Analysis of the Seemingly Unrelated
Regression Model with Student-t Errors, and its Application for Forecasting // International
Journal of Forecasting. 2010. 26. Р. 413?434.
Документ
Категория
Без категории
Просмотров
8
Размер файла
749 Кб
Теги
алгоритм, регрессии, робастная, применению, распределение
1/--страниц
Пожаловаться на содержимое документа