Робастная регрессия с применением t-распределения и EM-алгоритма Текст научной статьи по специальности «Математика»

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

и EM-алгоритма

Шведов A.C.

В работе рассматривается линейная регрессионная модель. EM-алгоритм представляет собой распространенный подход к оценке параметров таких моделей на основе общего принципа максимизации правдоподобия. Известно, что этот метод оценки параметров является робастным, если ошибки независимы, одинаково распределены и имеют многомерное i-распределение.

В предыдущих работах такой подход к оценке параметров регрессионных моделей применялся лишь при условии, что ошибки имеют многомерное i-распределение с числовым параметром степеней свободы. В настоящей работе рассматривается более общая ситуация, когда ошибки могут иметь многомерное t-распределение с векторным параметром степеней свободы. Ненаблюдаемые величины в EM-алгоритме при этом оказываются случайными матрицами.

На численных примерах при различных распределениях ошибок исследованы преимущества такого подхода по сравнению с методом наименьших квадратов.

Ключевые слова: робастная регрессия; многомерное {-распределение; EM-алгоритм.

1. Введение

Регрессионные модели являются основным инструментом для выявления зависимостей между различными показателями практически во всех областях экономической науки. Однако классические и наиболее распространенные методы построения регрессионных моделей не обладают свойством робастности, что, в ряде случаев, может приводить к неверным результатам. Робастные методы в эконометрике известны достаточно давно (см., например, [11]). Но все же в настоящее время, скорее, можно говорить об усилении тенденции к применению робастных методов, а не об обязательном требовании использования таких методов хотя бы наряду с классическими для контроля качества результатов.

Идея использовать регрессионные модели, у которых ошибки имеют не нормальные (гауссовские) распределения, а t-распределения, возникает, даже если не связы-

Шведов A.C. - д. физ.-мат. н., профессор кафедры математической экономики и эконометрики НИУ «Высшей школы экономики», e-mail: ashvedov@hse.ru

Статья поступила в Редакцию в декабре 2010 г.

вать этот подход с робастностью процедур оценки параметров. Еще в XIX в. ученые знали «об опасностях, порождаемых длинными хвостами функций распределения ошибок» (см. [2, с. 7]). (В настоящее время чаще используется термин не «длинные хвосты», а «тяжелые хвосты».) Изначальное использование при анализе предположения о нормальности ошибок включает и предположение о легких хвостах функций распределения. Это может не отвечать существу дела и приводить к искажениям в выводах. Одним из самых распространенных подходов к моделированию тяжелых хвостов является использование t-распределения.

Хотя для линейных регрессионных моделей, у которых ошибки имеют t-распре-деление, и не существует такой замкнутой и красивой теории, как для линейных рег-рессионныгх моделей, у который: ошибки имеют нормальное распределение, можно говорить и о преимуществах моделей с t-распределением. Так, фактически, регрессионные модели, у которых ошибки имеют t-распределение, включают в себя в качестве частного случая регрессионные модели, у которых ошибки имеют нормальное распределение, поскольку при стремлении числа степеней свободы к бесконечности t-распре-деления переходят в нормальное распределение. И результаты в предположении, что ошибки имеют t-распределение с достаточно большим числом степеней свободы, и в предположении, что ошибки имеют нормальное распределение, оказываются практически неотличимыми. Наконец, при оценке методом максимального правдоподобия параметров регрессионных моделей, у которых ошибки имеют t-распределения, можно использовать процедуры, обладающие свойством робастности.

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

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

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

В какой мере можно говорить, что эти претензии соответствуют действительности? Существуют различные подходы к построению робастных алгоритмов. Иногда резко выщеляющиеся наблюдения автоматически игнорируются. Для тех методов, которые изучаются в настоящей работе, вклад таких наблюдений только уменьшается. Для каждого класса алгоритмов слова «хорошее соответствие общей структуре и при наличии резко выделяющихся наблюдений» наполняются своим содержанием.

Среди предшествующих работ, в который: изучаются алгоритмы того же класса, что и у нас, назовем [7, 12, 19]. К перечисленным можно было бы добавить и интересную работу [13], однако в [8, с. 165] указывается на неправильные выводы, имеющиеся в этой работе. Подробнее о «подводных камнях», возникающих, если включать в со-

став аргументов функции правдоподобия число степеней свободы {-распределения, см. [8, 14].

Робастные алгоритмы других классов представлены, например, в книгах [15, 18]. Применяется и байесовский подход (см., например, [9]). Из работ прикладной направленности назовем [17, 22].

Содержание настоящей работы следующее. В параграфе 2 на примере множественной регрессии (наблюдения одномерные, объясняющих факторов несколько) обсуждается связь М-оценок и метода наименьших квадратов с итерационно модифицируемыми весами. В параграфе 3 приводится описание ЕМ-алгоритма, специализированного метода нахождения точки максимума именно функции правдоподобия. В параграфе 4 излагаются некоторые результаты, относящиеся к оценке параметров множественной регрессии с ошибками, имеющими одномерное {-распределение. Объясняется, почему применение ЕМ-алгоритма в данном случае дает робастный метод оценки параметров регрессии. Также в этом параграфе приводятся результаты численного исследования по методу Монте-Карло. В параграфе 5 устанавливаются две новые теоремы о матричном гамма-распределении. Затем эти теоремы используются в параграфе 6, где результаты, изложенные в параграфе 4, обобщаются на случай многомерной регрессии (наблюдения многомерные, объясняющих факторов несколько). При этом ошибки имеют {-распределение с векторным параметром степеней свободы (введенное в [3, 4]). Прием, применяемый в параграфе 6, когда в ЕМ-алгоритме в качестве ненаблюдаемых переменных берутся случайные матрицы, видимо, используется впервые. Также в этом параграфе приводятся результаты одного расчета.

2. Робастность М-оценок

Рассмотрим обычную линейную регрессию

(1) Уг =е X «Р«+ е , ' = .

а=1

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

(2) 1фЗ е

ст ист,

где ф( х) - некоторая известная функция плотности; ст> 0 - масштабирующий множитель. (Как обычно, используется одно и то же обозначение ег и для случайной величины, и для аргумента функции плотности.) Задача состоит в нахождении параметров Р1,...,Р9 и о.

Если ввести в рассмотрение 9 -мерные вектора хг = (хг1,...,хгд) и Р = (р1,...,Р ), штрих означает транспонирование, то сумму, входящую в правую часть (1), можно обозначить хг'Р.

Для определения вектора Р может быть использован метод наименьших квадратов с весами, когда оценка вектора Р строится путем минимизации выражения

(3) е^ (у- - x'p)2,

где w1,...,wn - заранее выбранные положительные числа. В частности, при wl = ... = wn = 1 данный метод является обычным методом наименьших квадратов. (Мы сейчас не касаемся теоретических свойств метода наименьших квадратов с весами. Подчеркнем только, что речь не идет об обобщенном методе наименьших квадратов, см., например, [1, гл. 5], хотя там и возникают сходные уравнения.) Приравнивание к нулю частных производных функции (3) по р1,...,рq, что является необходимым условием минимума, дает систему уравнений

n

(4) е X aWi (У- - X'P)= 0 , а = и. q .

-=1

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

(5) - n logs + еlogф| Уг - X;P

г=1

Если ввести в рассмотрение функцию

/ \ 1 ф'(x)

(6) w( x) =--^-,

X ф(X)

то приравнивание к нулю частных производных функции (5) по P1,...,Pq дает систему уравнений

(7) е.Хг^[^-^jy, - xJP)= 0 , а = 1,...,q .

Уравнения (7), хотя и похожи на уравнения (4), отличаются от них тем, что веса зависят от искомого параметра P . Введем в рассмотрении n х q матрицу X, i -я строка которой - это X , и диагональную n х n матрицу W, у которой i -й элемент на главной диагонали - это . Тогда уравнения (4) можно записать в форме X ' W(y - XP) = 0 ,

где y = (y1,..., Уп ) . Отсюда

(8) P = (X'WX )-1 X 'Wy ,

если матрица X' WX невырожденная. Этот же прием может быть использован и для решения уравнений (7).

Предположим, что функция w(х) обладает следующими свойствами. Во-первых, она принимает только неотрицательные значения. Во-вторых, эта функция монотонно не убывает при х < 0 и монотонно не возрастает при х > 0 . В-третьих, w(х)

стремится к нулю и при х ® да, и при х ® -да. Тогда можно говорить о робастности оценок параметра Р, полученных при помощи системы уравнений (7), поскольку резко выделяющимся наблюдениям у, , как правило, соответствуют большие по абсолютной величине разности у, - хг'Р и, соответственно, малые веса м>{^У—.

1 — 2/2

Нетрудно увидеть, что если ф(х) =_е х - функция плотности стандартно-

V 2 л

го нормального распределения, то w( х) = 1. И в этом случае система уравнений (7) не приводит к робастным оценкам параметра р. Отсюда возникает идея рассмотреть так называемые М-оценки, которые включают в себя в качестве частного случая оценки максимального правдоподобия (см., например, [2, 6]). В этом случае функция w(х), используемая в системе уравнений (7), не обязательно связана с функцией ф( х) соотношением (6). Зато можно потребовать, чтобы функция w(х) обладала тремя перечисленными выше свойствами. Или даже более сильными свойствами, например, обращалась в ноль при достаточно больших по абсолютной величине х.

Аргументация в пользу М-оценок может быть и такой. Если функция ф(х) на практике все равно не известна, то почему нужно начинать с выбора функции ф( х), а не с выбора функции w(х) ? Может быть, правильнее начинать с выбора функции w(х), а функцию ф(х) определять из уравнения (6)?

Но мы все же начнем с выбора функции ф( х), а не с выбора функции w( х). При любом а > 0 можно рассмотреть функцию

(9) j( x) =

1 Г (а + 0,5)

л/2ла Г(а)

( x2 ^ 1 + —

и 2а ш

-а-0,5

функцию плотности t -распределения с 2а степенями свободы. Известно, что при а ® да функция ф(х) переходит в функцию плотности стандартного нормального распределения. Нетрудно увидеть, что если функция ф( х) задается формулой (9),

, , ч , ч 2а +1 тт

то определяемая соотношением (6) функция w( х) имеет вид w( х) =-. И в этом

2а + х2

случае можно ожидать, что полученная путем решения системы уравнений (7) оценка параметра Р будет обладать свойством робастности (даже если не переходить от оценок максимального правдоподобия к более общим М-оценкам), поскольку w( х) стремится к нулю и при х ® да, и при х ® —да.

Если записать систему уравнений (7) в виде (8), то W - это диагональная

Ж у, — х'рЦ ^

п х п матрица, у которой , -й элемент на главной диагонали равен wl —-— I. Ре-

шить уравнение (8) можно попытаться методом простой итерации. Выбрав начальное приближение Р(0), например, при помощи метода наименьших квадратов, т.е. решив уравнение (8) с единичной матрицей Ж, затем принимаем

(10) р(г+1) =(х' Ж(г)Х)-1 X' Ж(г) у ,

где Ж(г) - это диагональная п х п матрица, у которой I -й элемент на главной диа-

ж

гонали равен м>

( уг - хг'Р(г) ц и а(г)

Приравнивание к нулю частной производной функции (5) по ст дает выражение

(11) а2=Iее (у, - х'рш^к - х;Р)

„у, -х;р

п ,=1

что равносильно соотношению

а2 = - (у - Хр)Ж (у - ХР). п

Последнее уравнение позволяет находить значения а итерациями:

(12) а(г+1)2 = I(у - ХР(г+1))Ж(г)(у - ХР(г+1)).

п

Для определения а(0) может быть использована единичная матрица Ж. Хорошо известно, что последовательность значений, получаемых при помощи метода простой итерации, может быть как сходящейся, так и не быть сходящейся. На практике этот метод обычно применяют для небольших расчетов, руководствуясь принципом «раз сошлось, значит, решение получено». Хотя и существуют условия, гарантирующие сходимость метода простой итерации. Кроме того, в правой

части (12) стоит не Р(г), а Р(г+1), т.е. в данном случае производится некоторое усложнение метода простой итерации.

В параграфе 4 показано, что для множественной линейной регрессии, когда ошибки имеют t -распределение, применение ЕМ-алгоритма приводит к тому же

итерационному процессу (10), (12) для определения Р(г+1), а(г+1) . А тогда применимы теоремы о сходимости итерационного процесса, построенного на основе ЕМ-алгоритма. В рамках ЕМ-алгоритма также можно провести обобщение на случай, когда наблюдения у1,...,уп не одномерные, а т -мерные (см. параграф 6).

Более подробно связь метода наименьших квадратов с итерационно модифицируемыми весами и робастных процедур освещается, например, в работе [21].

3. EM-алгоритм

EM-алгоритм предназначен для поиска точки, в которой достигает максимума функция правдоподобия, путем построения некоторого итерационного процесса. Каждый шаг итерационного процесса состоит из двух подшагов. E-подшаг заключается в нахождении ожидания (expectation) некоторой функции от случайных величин. При этом ожидание само оказывается функцией интересующего параметра. M-подшаг -это максимизация (maximization), определение того значения параметра, при котором данная функция достигает максимума. Первые буквы приведенных английских слов и дают название алгоритма. Обзор различных задач, для решения которых применяется EM-алгоритм, можно найти, например, в работе [16].

Пусть у = (y1v.., yn) - это набор наблюдений, вообще говоря, многомерных. z = (z1,..., zn) - набор ненаблюдаемых величин также, вообще говоря, многомерных. С одной стороны, предполагается, что ненаблюдаемые величины заметно влияют на наблюдаемые, и привлечение их для анализа отвечает существу дела. С другой стороны, нахождение точек максимума функций в рамках EM-алгоритма может оказаться значительно более простой и надежной с вычислительной точки зрения процедурой, чем непосредственное нахождение точки максимума исходной функции правдоподобия. В каких-то задачах не вызывает сомнений, что именно следует взять в качестве ненаблюдаемых величин. В других задачах ответ на этот вопрос не столь очевиден.

Обозначим через h совместную функцию плотности случайных векторов y и z, через g - условную функцию плотности случайного вектора z при заданном у, через f - маргинальную функцию плотности случайного вектора у. Все эти функции плотности считаются зависящими от некоторого параметра 6, вообще говоря, многомерного. Имеет место соотношение

f (у; 6) = h( у,z; 6)

g(z | у; 6)

Переходя к логарифмам, получаем соотношение для логарифмических функций правдоподобия

(13) /(61 у) = /(6 | у, z) - log g(z | у; 6).

Задача состоит в нахождении точки 6, в которой достигает максимума функция /(6 | у). Пусть 6(r) - значение параметра 6, найденное при r -й итерации. Умножим левую и правую части (13) на g (z | у; 6( r)) и проинтегрируем по z. Введем обозначения

U(6 | у;6(r)) = J/(61 у,z) g(z | у;6(r)) dz, P(61 у; 6(r)) = J log g (z | у; 6) g (z | у;6(г)) dz.

Тогда

(14) /(61 у) = U(6 | у;6(r))- p(6 | у;6(r)).

Заметим, что

(Г )

f logg(zly>q(r)) & J g( z|y; 6)

- это расстояние Кульбака - Лейблера между функциями плотности g(z | y; 6( r)) и g(z | y;6), которое, как известно, всегда неотрицательно. Поэтому при любом 6

(15) р(61 y;6(r)) <р(6(r) | y;6(r)).

E-подшаг состоит в нахождении ожидаемого логарифмического правдоподобия U(6 | y;6(r)).

M-подшаг состоит в нахождении точки

= argmax U(6 | y;6(r)).

6

Из (14), (15) и способа определения точки 6(r+1) следует, что

(16) /(6(r+1)|y)> l(6(r)|y).

Соотношения (16) показывают, что движение идет «в правильном направлении»,

но еще не гарантируют, что последовательность 6(r) сходится. Условия общего характера, из которых следует сходимость этой последовательности к точке максимума функции l(6 | y), даются в работе [20] теоремами 1 и 4.

4. Множественная линейная регрессия

В этом параграфе рассматривается набор одномерных наблюдений y1,...,yn , и будем предполагать, что ненаблюдаемые величины z1,...,zn также одномерные и, кроме того, положительные. Двумерные случайные величины (y1,z1) ,..., (yn,zn) считаем независимыми.

Предположим, что совместная функция плотности случайных величин yi и zi имеет вид

1/2 Ж 2 Ц

g (Zi) ,

z

(17) -!— exp

av 2 л

и 2а2 y

где ег = уг - хг'Р в соответствии с (1). Будем использовать обозначение 6 для пары Р, о. Тогда ^у,z;Р,о) - это произведение функций (17) при г = 1,...,п . Следовательно,

П 1 П 1 П П

log h(y, z;b, а) =--log2p-n log а+ - V log z;--- V zrf + V log g (z;).

2 2 2а2

i=1 i=1

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

1 n

(18) l(0|y, z) = -n log g--- У z-e"2.

2s2 1=1

В соответствии с определением, данным в параграфе 3,

и (Э | у; 0(г) )= Е (?(0|у, г) | у; 0(г)),

где 1(0 | у, г) рассматривается как функция случайного вектора г. Из (18) следует, что эта функция линейна относительно г1,..., 2п.

Распределение вероятностей с совместной функцией плотности (17) называется нормальным-гамма распределением, если

Aa

(19) g (z) =-exp(- Az) za-1

G(a)

где а > 0 , А > 0. Тогда, как нетрудно увидеть, (выкладки для т -мерного случая приводятся в параграфе 6) маргинальная функция плотности случайной величины у, имеет вид (2), где

(20) Ф<» = (2.)-"2 . пА 2 ) *+0,5 -

Г(а) (А + 0,5х2 ) ,

Из (6) следует, что в этом случае

, ч а + 0,5

(21) w( х) = '

A + 0,5x2

Введем обозначение e(r) = yi - x\b(r) и воспользуемся тем, что

(22) 4i-|y;0(r))= w [g^e(r)] .

Доказательство соотношения (22) для m -мерного случая приводится в параграфе 6. Для одномерного случая (22) доказано в работе [7], и в этом доказательстве не требуется, чтобы функция g (z) обязательно имела вид (19). Из (18) и (22) следует, что

U (01 y; 0(r) )= -n log g- ¿w i-^ e(r) Ц (yt - x b)2.

2g -=1 ug Ш

Таким образом, E-подшаг EM-алгоритма выполнен. Выполнение M-подшага аналогично процедуре, описанной в параграфе 2 (ср. (3), (5), (7), (11), (12)).

Пример 1. Сокращения М1 и М2 в этом примере используются для обозначения одного из двух методов, применяемых для определения параметра Р и стандартного отклонения возмущений 5.

М1 - метод максимального правдоподобия в предположении нормальности возмущений.

М2 - ЕМ-алгоритм в предположении, что возмущения имеют t -распределение с тремя степенями свободы.

Рассматриваются три различных вида сгенерированных рядов наблюдений.

Н0 - наблюдения соответствуют модели с нормальными возмущениями.

Н1 - наблюдения соответствуют модели с нормальными е -засоренными возмущениями, е = 0,1.

Н2 - наблюдения соответствуют модели с возмущениями, имеющими t -распределение с тремя степенями свободы.

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

Пусть п = 15 , q = 1, хп = 1 при каждом i, 1 < i < п . Случайная величина ei имеет функцию плотности (2), где ф(х) - либо функция плотности стандартного нормального распределения, либо функция плотности t -распределения при 2а = 3 ; обе эти функции плотности приведены в параграфе 2. Для генерации наблюдений yi используются значения Р = 1, 5 = 0,3 . При использовании нормальных возмущений ст = 5. При использовании возмущений, имеющих t -распределение с V степенями свободы,

V-2

СТ = 5

V

При генерации ряда с е -засоренными наблюдениями считается, что с вероятностью (1 - е) стандартное отклонение возмущения равно 5, и с вероятностью е равно 55.

Для каждого из трех видов генерируется L = 300 рядов наблюдений длины п . Для I -го эксперимента, I = 1,...,L значения параметров Р1 и 51 определяются и методом М1, и методом М2. Затем определяются средние значения

1 L 1 L

р=7 &, 5-=7 е 51

L 1=1 7 г=1

и среднеквадратические отклонения

Г, 7 . .. У" Г, 7

1-7 ее (Рг-Р)2 , 7£(5г - 5)

г=1

7 е 5 ^2

и г=1 Ш

Результаты для средних значений приведены в табл. 1 и 2. В скобках даются среднеквадратические отклонения.

Таблица 1.

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

НО

Н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.

Результаты для параметра ж, истинное значение ж = 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) дает лучшие результаты и в смысле меньшего разброса (т.е. среднеквадратического отклонения), и в смысле близости среднего значения найденных параметров к истинным значениям (за исключением значения Р= 1,012, которое несколько хуже, чем значение р= 1,011).

На первый взгляд, результаты для параметра Р для рядов вида Н2 противоречат теореме Гаусса - Маркова. Метод М1 совпадает с методом наименьших квадратов, и, казалось бы, дисперсия оценки должна быть наименьшей. А среднеквадратиче-ское отклонение 0,067 существенно больше, чем среднеквадратическое отклонение 0,050, полученное при использовании метода М2. На самом деле, этот пример всего лишь показывает, что требование линейности и несмещенности оценки, содержащееся в теореме Гаусса - Маркова, не может быть отброшено. Метод М2 не является линейным.

Для рядов вида Н1 в результатах для параметра Р видна робастность метода М2. Полученное среднеквадратическое отклонение примерно в 1,65 раза меньше, чем для метода М1. Различие в средних, 1,016 и 1,019, носит случайный характер. Так, для другой серии из 300 экспериментов для рядов вида Н1 для параметра Р при

1,023

использовании метода М1 получены результаты , а при использовании мето-

да М2 - результаты 1,019 . Проявляется робастность метода М2 и в результатах для параметра 5 для рядов этого вида.

5. Две теоремы о матричном гамма-распределении

Теория матричных гамма-распределений излагается, например, в [10]. Наиболее известными из этих распределений являются, видимо, распределения Уишарта. Также в работе [10, с. 122-124] рассматриваются матричные гамма-распределения с векторным параметром - некоторое естественное обобщение распределений Уишарта. При том, что легкими и короткими формулировки здесь быть не могут, обозначения, используемые в [10] и других работах для распределений с векторным параметром, с нашей точки зрения, несколько избыточны, возможно, из-за этого остались неисследованными свойства этих распределений. Более прозрачные обозначения для матричных гамма-распределений с векторным параметром используются в работе [3]. Здесь мы повторим только самые необходимые определения.

Для т х т матрицы С = (сг;- при k = 1,...,т рассматриваются подматрицы С^] = {Сг; = и С№] = {Сг; }т].

Рассматривается также вектор а = (а1,...,ат) такой, что а у > 0,5(у -1) при у = 1,..., т. Многомерная гамма-функция определяется следующим образом:

т

С (а) = рт(т-1)/4 Хг(а; - 0,5( у -1)),

У=1

где Г(-) - обычная гамма-функция. Дополнительно считается, что

а0 = 0 , ат+1 = 0,5(т +1).

Параметрами рассматриваемых гамма-распределений являются вектор а указанного вида и положительно определенная т х т матрица А. Функция плотности имеет вид

т

(23) g (Г) = у а, а etr(-Az) Х 2 у ] |ау-ау+1,

у=1

где 2 - положительно определенная т х т матрица; егг(С) = ехр(гтС); | С | - определитель матрицы С . Коэффициент уаА задается формулой

( т-1

(24) У а,А =

гт(а) Хи[т-у] а

У=о

(ср. (23) и (24) с (19)).

Пусть Т - симметричная т х т матрица такая, что матрица А - Т положительно определенная. (Изначально симметричная матрица Т может быть взята произвольно. При некотором е > 0 матрица А -еТ будет положительно определенной. В этом смысле условие, что матрица А - Т положительно определенная, не является ограни-

чительным.) Пусть Z - положительно определенная случайная матрица с функцией плотности (23). Определим производящую функцию моментов

Из формулы

М (Т) = ЕеГ (Т2).

1т(Т1) = £т„2„ +

к1

1=1

к=1 1>к

следует, что

(25)

а при 1 > к

(26)

Теорема 1.

М (Т) =

дМ (Т)

дТ„

дМ (Т)

Т=0

дТ

к1

--е 2„),

: 2 Е(Zkl).

У а,А

т-1

(

т-1

а,А = ХХ| (А - Т)[т-1 ] |а1-а1+! | А[т-1] ^

У а,А-Т 1 =0 и 1=0

Доказательство следует из выражения

т

М(Т) = Уа,А Т е*(-(А - Т)7) 1] I

г>0

1=1

и из формулы (24).

При 1 = 0,..., т -1 определим т х т матрицу Ст-1 такую, что

(Ст-; ) [т-1 ]=(А[т-1 ] )-1

остальные элементы матрицы Ст-1 равны нулю. Теорема 2.

т-1

Е (2) = еа+1 - а} )с

т- 1

1=0

Доказательство. Чтобы найти ожидание каждого элемента 2Ы , к < 1, случайной матрицы 2, воспользуемся формулами (25), (26) и теоремой 1.

Т=0

-1

Зафиксируем j такое, что I < т - j . Минор элемента Ак1 - Тк1 матрицы (А - Т)[т обозначим Мк1. Через Е обозначим определитель этой матрицы. Введем обозначение

•и = Ак1 - Тк1. Тогда

m- j

(27) Е = е("1)И+Ч^.

и=1

Дифференцирование функции (27) по •п не вызывает затруднений, поскольку ни один из миноров от ••ц не зависит. Имеем

(28) — = Мп = |(А - Т)[т-п|((А - Т)[т-п)- .

5 %

Дифференцирование функции (27) по •к1 при к < I несколько труднее. Все миноры кроме Мп могут зависеть от •к1, поскольку •1к = •к1. Имеем

9Е , , т^, -*„+1 5

(29) — = (-Dk+lMkl + е(-1Г'su,—Mul .

9 skl u=1 9 skl

u№l

Обозначим через Mul lv определитель матрицы, получающейся из (A - T)[m-выкидыванием строк с номерами u и / и столбцов с номерами / и v . Тогда при u < l

l-1 m- j

Mul = е (-1)V+l-1 SlvMul ,lv - е (-1)V+l-1 SlvMul ,lv ; v=1 v=l+1

при u > l

l-1 m- j

Mul = е (-1)v+l SlvMul,lv - е (-1)v+l SlvMul, lv . v=1 v=l+1

Ни один из определителей Mul lv, входящих в две последние формулы, от skl

не зависит. Вторыми суммами в правых частях также, очевидно, можно пренебречь, поскольку k < l. Поэтому при u < l

9 -Mul = (-1)k+l-1 Mul, ik;

при u > l

9 skl

9

- Mul = (-1)k+lMul llk .

9 skl

Из (29) при k < l находим

1-1 m-1 ^=(-\f+iMkl+е (-ir'sul (-i)k+ е (~i)u+'sui (-Dk+iMikui.

5 skl u=1 u=l+1

То есть при k < l

(30) — = 2(-1)k+lMkl = 2\(A - T)[m-j]|((A - T)[m-j])-l.

5 skl

Изменения, которые нужно внести в приведенные выкладки при m - j < 2 или при l = да - j, очевидны.

Те же выражения (28) и (30) получаются при дифференцировании |(A - T)[m-j]| по Tll и по Tkl, только в правые части добавляется знак минус.

Воспользовавшись теоремой 1, формулами (28) и (30), получаем

5M (T)

5T„

m-l / \

■ееa - j)(A[m-j])-;:

j=0

T=0

и при k< l

m-l

5M (T)

5Tkl

2ё<?1 -ji) (A[m-1 ])-l.

j=0

T=0

Воспользовавшись (25) и (26), при k < l получаем

m-1

E (Zkl ) = е(й!+1 - a! )(cm- j )kl .

1=0

Теорема 2 доказана.

Теорема 2 для случая а1 = ... = ат известна. Доказательство приводится, например, в [10]. При этом и в [10], и в других работах используется не производящая функция моментов, а характеристическая функция. Использование производящей функции моментов позволяет избежать рассмотрения матриц с комплексными элементами. Теорема 1, хорошо известная для одномерного случая, при т > 1, по-видимому, является новой даже для случая а1 = ... = ат.

6. Многомерная линейная регрессия

Рассмотрим уравнения, аналогичные (1):

(31) уг =Р'хг + ег, г = 1,..., п .

Объясняющие переменные xiа - известные числа, х{ = (хг1,...,х{С1); Р - дхт

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

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

(32) (2л)-

1/2

1/2

1

*

ехРГ ^ е 18 )

где 8^) определяется формулой (23), и ei = yi -Р'х± в соответствии с (31).

Рассмотрим вектор Ь = (Ь1,...,Ьт), где — = а}- + 0,5 при у = 1,..., т . Пусть Ь0 = 0 , Ьт+1 = 0,5(т +1).

При х е Rm рассмотрим функцию

г*(Ь) т-1 Г 1 '

(33) ф(х) = (2л)-т/2-Гт(—^ Ы|-1/2 П1 1 + -х[т-;] (А[т-;])-1 х[т-I

гт(а) и 2 ш

Она является функцией плотности многомерного t -распределения с векторным параметром степеней свободы (см. [3; 4]); ср. (20). Здесь для т -мерного вектора х через

х[4г] обозначается k -мерный вектор, состоящий из первых k компонент вектора х, k = 1,...,т.

Теорема 3. Маргинальная функция плотности случайного вектора yi имеет вид

(34) -!тФр-е;

ст

где функция ф задается формулой (33).

Доказательство. Маргинальная функция плотности случайного вектора yi получается интегрированием по области zi > 0 совместной функции плотности (32). Чтобы несколько сократить формулы, внутри доказательства теоремы будем использовать обозначение * вместо zi и обозначение е вместо ei. Во-первых, заметим, что

е' = ее' * ).

Используя (23) и (24), получаем следующее выражение для маргинальной функции плотности:

Ы-"2т„,атI -На+^4I П1*"1^" =

1 т-1

= (2р)-т/2Г» П

а 1 =0

1=1

А +-ее '

2а2

[т-Л

Ь-ь.

Если воспользоваться тем, что

1 ;

А[т-1] + е[т-1 ]е[т-1]

2а2

= А[т-1]||1 + 1 е[т-1 ] (А[т- 1])-1е[т-1] | 1 I 2а2 Ш'

(см., например, лемму 5 в [3]), то получаем выражение

(2р)-т/2 |А|-1/2 П(1 + -1-е[т-1]'(А[т-1])

Г>) ат 1=|1 *-2 ( )

Ь. -ь.

2а2

1=0 ■

Теорема 3 доказана.

Теорема 4. Условная функция плотности случайного вектора г, при условии

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

1 ;

ром Ь и с матричным параметром А +--е,е, .

2а2

Доказательство. Условная функция плотности случайного вектора г, при условии у, - это отношение совместной функции плотности (32) к маргинальной функции плотности (34). Как и в доказательстве предыдущей теоремы, будем использовать обозначение г вместо г, и обозначение е вместо е,. Искомая условная функция плотности представима в виде дроби с числителем

е1г

-| А +-

1 Ц Ц т

ее'\г ] П1г'

Ь-Ь

[ Л

и со знаменателем

т-1 Ь Ь т-1 ( 1 ;

ГГ(Ь) П|А[т-1]| '- '+1 П| 1 + "Ге[т-1] (А[т-1])-1е[т-1]

1=0 1=0 и 2а 2

Ь -Ь,

Используя, как и в доказательстве предыдущей теоремы, лемму 5 из [3], получаем следующее выражение для знаменателя:

Г*(Ь) П

1=0

1

А[т-1 + е[т-1 ]е[т-1] 2а 2

(

\

-1

У 1

Ь, А+-- ее

4 2а Ш

Теорема 4 доказана.

Ь +1

Будем использовать обозначение 6 для пары Р , ст . Из (32) следует, что и в многомерном случае логарифмическая функция правдоподобия имеет вид, сходный с (18):

(35)

1 П

l(0|y, z) =- n log sm--- У e,

2s2 tf

2iii

Ввиду линейности 1(6 | у,2) по 2Х,...,2п , чтобы построить функцию и(6 | у;6(г)) достаточно знать условное ожидание Е(2{ | у{; 6(г)).

При у = 0,..., т -1 и при х е Вт определим т х т матрицу Ст- .■ (х) такую, что

Cm-J (X))

[m-J] _

A +—xx '

[m- J] I 1

остальные элементы матрицы Cm- j ( x) равны нулю. Положим

w(X) = У(Ь+1 - bJ )Cm-J (X)

j=o

(ср. (21)). Тогда на основании теорем 2 и 4 получаем E(z, | y,;0(r))= w\—e\r)Ц (ср. (22)). Определим m х m матрицы

w(r) = w\— e(r) I, i = 1,...,n .

г U(r) i ш,

В соответствии с определением функции U(0 | y; 0(r)) в параграфе 3 (см. также параграф 4) и с (35) получаем

n

U(Р,s | y;P(r),s(r)) = -nmlogs- —-У[ y' - х, 'p)w(r) (y -p 'х).

-s2 t1U Ш

Дифференцирование функции U по Рак , а = 1,..., q , к = 1,...,m и приравнивание производной к нулю с учетом симметричности матрицы w дает уравнение

n m

(

\

(36)

У X аУ W]k y,j -У X yP

,=1 J=1

yFy J

g=1 ш

= о.

Дифференцирование функции и по ст и приравнивание производной к нулю дает уравнение

(37)

= — Уи у,' - Xi Р^ w(r) (y, -Р' X )

nm ,=1 v J

2

ст

(ср. (11), (12)). Найденные из уравнений (36) и (37) значения Р и ст принимаются в качестве р(г+1) и ст(г+1). Значения Р(0) и ст(0) рассчитываются с единичными матрицами wi.

Другой вариант применения ЕМ-алгоритма в линейной регрессии, когда ошибки имеют многомерное t -распределение, представлен, например, в работе [12].

г

Пример 2. Пусть т = 3 , п = 1000 , q = 2 . При каждом / = 1,...,п вектор х^ = (1,/), а трехмерная случайная величина в1 имеет функцию плотности (33), (34) с параметрами а = (2,5,9),

ж 1 1,3 1,5Ц

А = 1,3 2 2

Л5 2 4 ш

ст = 12 . При генерации трехмерных возмущений с указанным t -распределением применяется алгоритм Метрополиса (см., например, [5]). Для генерации наблюдений используется матрица

_,10 20 - 30 Р = и 0,3 - 0,2 0,4

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

При применении ЕМ-алгоритма, т.е. при использовании итерационного процесса, основанного на формулах (36), (37), сходимость в проведенном эксперименте была достигнута после 20 итераций. Получены значения параметров

,'10,185 20,268 - 29,681

Р = ,

1 0,299 - 0,201 0,400

ст = 11,533 . Результаты являются удовлетворительными. Отметим, что при применении метода наименьших квадратов получено значение параметра

,10,191 20,399 - 29,634

Р = ,

1 0,299 - 0,201 0,400 ^

Оно же использовалось в качестве начального значения Р(0) в ЕМ-алгоритме. В данном расчете метод наименьших квадратов уступает ЕМ-алгоритму.

* * *

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

1. Магнус Я.Р., Катышев П.К., Пересецкий А.А. Эконометрика. Начальный курс. М.: Дело, 2004.

2. Хъюбер П. Робастность в статистике. М.: Мир, 1984.

3. Шведов А.С. Бета-распределение случайной матрицы и его применение в модели состояние-наблюдение: препринт. WP2/2009/01. М.: ГУ ВШЭ, 2009.

4. Шведов А.С. t-распределение случайной матрицы и его применение в регрессионной модели: препринт. WP2/2010/01. М.: ГУ ВШЭ, 2010.

5. Шведов А.С. О методах Монте-Карло с цепями Маркова // Экономический журнал Высшей школы экономики. 2010. Т. 14. № 2. С. 227-243.

6. Andrews D.F. A Robust Method for МиШр1е 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. Р. 213255.

12. Lange K.L., Little R.J.A., Taylor J.M.G. Robust Statistical Modelling Using the t-distri-bution // 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 // Encyc^edia 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.